Zagadnienia

4. Równania z małym parametrem

Mały parametr w równaniu różniczkowym może pojawiać się zasadniczo na dwa sposoby: albo z prawej strony albo z lewej strony. W pierwszym przypadku mamy do czynienia z małym zaburzeniem układu, o którym na ogół sporo wiemy a w drugim przypadku pojawiają się tzw. drgania relaksacyjne. Oba przypadki są omawiane kolejno w następnych rozdziałach.

4.1. Uśrednianie

Przykładem układu z pierwszej grupy jest znany układ van der Pola

\acute{x}=y,\text{ \ }\dot{y}=-x+\varepsilon(1-x^{{2}})y,

gdzie \varepsilon jest naszym małym parametrem. Jest to szczególny przypadek zaburzenia układu hamiltonowskiego o jednym stopniu swobody

\dot{x}=H_{{y}}^{{\prime}}+\varepsilon P(x,y),\text{ \ \ }\dot{y}=-H_{{x}}^{{\prime}}+\varepsilon Q(x,y). (4.1)

W zastosowaniach często pojawiają się układy hamiltonowskie z wieloma stopniami swobody postaci

\dot{q}_{{i}}=H_{{p_{{i}}}}^{{\prime}},\text{ \ \ }\dot{p}_{{i}}=-H_{{q_{{i}}}}^{{\prime}},\text{ \ \ \ }i=1,\ldots n, (4.2)

gdzie H(q_{{1}},\ldots,q_{{n}},p_{{1}},\ldots,p_{{n}}) jest funkcją Hamiltona, lub hamiltonianem (Zadania 4.11 i 4.12). Na ogół układ (4.2) nie daje się rozwiązać. Jednak istnieje klasa układów hamiltonowkich w pełni rozwiązywalnych.

Definicja 4.1. Układ (4.2) nazywa się zupełnie całkowalnym jeśli istnieje układ funkcjonalnie niezależnych całek pierwszych  F_{{1}}=H, F_{{2}},\ldots, F_{{n}} taki, że każda funkcja F_{{j}} jest całką pierwszą dla innych układów hamiltonowskich generowanych przez inne funkcje F_{{i}}. Mówi się też, że funkcje F_{{j}} są w inwolucji.

Przykładami układów zupełnie całkowalnych jest zagadnienie Keplera i potok geodezyjny na powierzchni elipsoidy (patrz [4]); oba mają dwa stopnie swobody.

Dla układów spełniających warunek z Definicji 4.1 zachodzi następujące twierdzenie, które przytaczamy bez dowodu (patrz [4]).

Twierdzenie 4.2 (Liouville–Arnold). Jeśli wspólne poziomice \{ F_{{1}}=c_{{1}},\ldots, F_{{n}}=c_{{n}}\} zupełnie całkowalnego układu hamiltonowskiego są zwarte i gładkie, to są one torusami \mathbb{T}^{{n}}.

Ponadto  w otoczeniu danego takiego torusa istnieje układ współrzędnych \left(I_{{1}},\ldots,I_{{n}},\varphi _{{1}},\ldots,\varphi _{{n}}\right), tzw. zmienne działanie–kąt, w których układ (4.2) przyjmuje następującą postać hamiltonowską

\dot{I}_{{j}}=0,\text{ \ \ }\dot{\varphi}_{{j}}=\omega _{{j}}(I)=\partial H_{{0}}/\partial I_{{j}},\text{ \ \ \ \ }j=1,\ldots,n, (4.3)

gdzie H(q,p)=H_{{0}}(I_{{1}},\ldots,I_{{n}}) jest hamiltonianem po zamianie. W szczgólności ruch na torusach \left\{ I_{{1}}=d_{{1}},\ldots,I_{{n}}=d_{{n}}\right\}, które są parametryzowane przez kąty \varphi _{{j}} mod 2\pi, jest okresowy lub prawie-okresowy (patrz Rysunek 4.1):

\varphi _{{j}}(t)=\varphi _{{j}}(0)+\omega _{{j}}(I)t.

Przykład 4.3. Dla układu van der Pola z \varepsilon=0 i H=\frac{1}{2}(x^{{2}}+y^{{2}}) zmienne działanie–kąt są następujące: I=H i \varphi=\arg(x+iy).

Rozważmy teraz następujące zaburzenie układu (4.3)

\dot{I}=\varepsilon g(I,\varphi),\text{ \ \ }\dot{\varphi}=\omega(I)+\varepsilon f(I,\varphi), (4.4)

gdzie I=(I_{{1}},\ldots,I_{{n}}), \varphi=(\varphi _{{1}},\ldots,\varphi _{{n}}) and \omega=\left(\omega _{{1}},\ldots,\omega _{{n}}\right). Naturalne jest spodziewać się, że rozwiązanie układu (4.4) po czasie rzędu O(1) różni sę od rozwiązania układu (4.3) z tymi samymi warunkami początkowymi o wielkość rzędu O(\varepsilon). Tymczasem poniższe Twierdzenie 4.4 mówi, że taką samą wielkość O(\varepsilon) można uzyskać po czasie dążącym do nieskończoności przy \varepsilon\rightarrow 0. Tego rodzaju zjawisko ma miejsce dzięki tzw. uśrednieniu.

Rys. 4.1. Dynamika prawie okresowa na torusie.

Idea uśrednienia wiąże się z faktem, że na większości torusów \mathbb{T}^{{n}}=\left\{ I=d\right\} trajektori układu niezaburzonego jest gęsta (jak na Rysunku 4.1). Zatem średnie odchylenie działania I(t) można wyliczyć (w przybliżeniu) poprzez uśrednienie po torusie.

Definijemy układ uśredniony

\dot{J}=\varepsilon G(J), (4.5)

gdzie

G(J)=\left(\frac{1}{2\pi}\right)^{{n}}\int _{{0}}^{{2\pi}}\ldots\int _{{0}}^{{2\pi}}g(J,\varphi)d\varphi _{{1}}\ldots d\varphi _{{n}}

jest uśrednioną po \mathbb{T}^{{n}} wielkością prędkości zmian działania.

Twierdzenie 4.4 (O uśrednianiu). Niech n=1 i funkcje \omega,f,g będą klasy C^{{1}} i \omega(I)>0 na otwartym podzbiorze \mathbb{R}^{{1}}\times\mathbb{T}^{{1}}. Jeśli \left(I(t),\varphi(t)\right) i \left(J(t),\psi(t)\right) są rozwiązaniami układów (4.4) i (4.5) takimi, że I(0)=J(0), to dla

0<t<1/\varepsilon

amy

\left|I(t)-J(t)\right|<C\cdot\varepsilon,

gdzie stała C zależy tylko od\omega,f,g.

Dowód. Dokonajmy zamiany

K=I+\varepsilon k(I,\varphi) (4.6)

tak, aby zachodziło \dot{K}=O(\varepsilon^{{2}}). Wyliczenie k(J,\varphi) przebiega następująco

\begin{array}[]{lll}\dot{K}&=&\dot{I}+\varepsilon\frac{\partial k}{\partial I}\dot{I}+\varepsilon\frac{\partial k}{\partial\varphi}\dot{\varphi}=\varepsilon\left\{ g+\frac{\partial k}{\partial\varphi}\omega\right\}+O(\varepsilon^{{2}})\\
&=&\varepsilon\left\{ g(K,\varphi)+\frac{\partial k}{\partial\varphi}(K,\varphi)\omega(K)\right\}+O(\varepsilon^{{2}}).\end{array}

Zatem chcemy rozwiązać równanie

\frac{\partial k}{\partial\varphi}(K,\varphi)\omega(K)=-g(K,\varphi),

z oczywistym rozwiązaniem g(K,\varphi)=\frac{-1}{\omega(K)}\int _{{0}}^{{\varphi}}g(K,\psi)d\psi. Niestety, na ogół to rozwiązanie nie jest jednoznaczną (czyli okresową) funkcją od \varphi. Przeszkodą jest wielkość \int _{{0}}^{{2\pi}}g(h,\psi)d\psi, która może być niezerowa.

Ale, zapisując

g(K,\varphi)=G(K)+\tilde{g}(K,\varphi)

tak, że \int _{{0}}^{{2\pi}}\tilde{g}(h,\psi)d\psi=0, możemy zdefiniować jednoznaczną funkcję

g(K,\varphi)=\frac{-1}{\omega(K)}\int _{{0}}^{{\varphi}}\tilde{g}(K,\psi)d\psi.

Dostajemy równanie

\dot{K}=\varepsilon G(K)+O(\varepsilon^{{2}}).

Widać, że po czasie O(1/\varepsilon) różnica pomiędzy J(t) i K(t) jest rzędu O(\varepsilon). Z drugiej strony, różnica pomiędzy K(t) i I(t) jest rzędu O(\varepsilon), dzięki zamianie (4.6). ∎

Dla zaburzeń typu (4.4) zupełnie całkowalnych układów hamiltonoskich z wieloma stopniami swobody oszacowania są słabsze niż w tezie Twierdzenia 4.4. Okazuje się, że po czasie czędu O(1/\varepsilon^{{a}}), dla warunków początkowych spoza zbioru o mierze Lebesque's O(\varepsilon^{{b}}), odchylenie J(t) od I(t) nie przekracza O(\varepsilon^{{c}}), gdzie a,b,c>0 są wykładnikami zależnymi od \omega,f,g. Po więcej informacji odsyłam czytelnika do [5].

Przykład 4.5 (Całki abelowe). Rozważmy następujące zaburzenie dwuwymiarowego układu hamiltonowskiego

\dot{x}=H_{{x}}^{{\prime}}+\varepsilon P(x,y),\text{ \ \ \ }\dot{y}=-H_{{x}}^{{\prime}}+\varepsilon Q(x,y).

Dla \varepsilon=0 krzywe fazowe leżą w poziomicach funkcji Hamiktona H(x,y). W pewnym obszarze przestrzeni fazowej te krzywe są zamknięte. Jak już robiliśmy to kilkakrotnie, badanie cykli granicznych układu zaburzonego polega na analizie przekształcenia powrotu Poincarégo z\  S (cięcie transweralne do krzywych fazowych) do S. Parametryzując S za pomocą H|_{{S}} warunek cyklu granicznego to \Delta H=H(B)-H(A)=0 (patrz Rysunek 4.2). Mamy

Rys. 4.2. Przekształcenie powrotu dla zaburzenia układu hamiltonowskiego.
\begin{array}[]{lll}\Delta H&=&\int _{{0}}^{{T}}\frac{dH}{dt}dt=\int _{{0}}^{{T}}\left\{ H_{{x}}^{{\prime}}\left(H_{{y}}^{{\prime}}+\varepsilon P\right)+H_{{y}}^{{\prime}}(-H_{{x}}^{{\prime}}+\varepsilon Q\right\} dt\\
&=&\varepsilon\int\left(PH^{{\prime}}x+QH_{{y}}^{{\prime}}\right)dt=\varepsilon\int\left\{ P(H_{{x}}^{{\prime}}-\varepsilon Q)+Q\left(H_{{y}}^{{\prime}}+\varepsilon P\right)\right\}\\
&=&\varepsilon\int _{{\Gamma(h)}}Qdx-Pdy=\varepsilon\oint _{{H=h}}\left(Qdx-Pdy\right)+O(\varepsilon^{{2}}),\end{array}

gdzie T jest czasem powrotu do S a \Gamma(h) jest krzywą fazową układu zaburzonego startującą z A\in S takiego, że H(A)=h.

Wyrażenie

I(h)=\oint _{{H=h}}Qdx-Pdy (4.7)

jest tzw. całką abelową.15Pojęcie całki abelowej wywodzi się z zespolonej geometrii algebraicznej. Są to całki z 1-form meromorficznych wzdłuż pewnych zamkniętych krzywych na zespolonych krzywych algebraicznych (powierzchniach Riemanna). Gdy H(x,y), P(x,y) i Q(x,y są wielomianami, to powierzchnia Riemanna jest zespoloną krzywą \left\{ H(x,y)=h\right\}\subset\mathbb{C}^{{2}} a 1-forma to \omega=\left(Qdx-Pdy\right)|_{{H=h}}. Z Twierdzenia o funkcjach uwikłanych wynika, że jeśli I(h_{{0}})=0 i I^{{\prime}}(h_{{0}})\not=0, to dla \varepsilon\not=0 i małego istnieje cykl graniczny \gamma _{{\varepsilon}}, który dąży do krzywej H=h_{{0}} przy \varepsilon\rightarrow 0. To podejście do problemu cykli granicznych jest szeroko stosowane w Jakościowej Teorii.

Nietrudno zauważyć, że funkcja I(h) jest odpowiednikiem całki uśrednienia G(J), która występuje we wzorze (4.5).

4.2. Teoria KAM

Rozważmy układ hamiltonowski

\dot{q}=H_{{p}}^{{\prime}},\text{ \ \ \ }\dot{p}=-H_{{q}}^{{\prime}},

p=\left(p_{{1}},\ldots,p_{{n}}\right),q=\left(q_{{1}},\ldots,q_{{n}}\right), z hamiltonianem postaci

H(p,q)=H_{{0}}(p,q)+\varepsilon H_{{1}}(p,q),

gdzie H_{{0}} jest hamiltonianem układu zupełnie całkowalnego, czyli układu typu (4.3) w zmiennych działanie–kąt. Z tą sytuacją wiąże się jedno z najważniejszych twierdzeń matematycznych drugiej połowy XX wieku. Przed jego sformułowaniem musimy wprowadzić jeszcze dwa założenia dotyczące niezdegenerowania niezaburzonego hamiltonianu H_{{0}}:

\det\left(\frac{\partial\omega _{{i}}}{\partial I_{{j}}}\right)=\det\left(\frac{\partial^{{2}}H_{{0}}}{\partial I_{{i}}\partial I_{{j}}}\right)\not=0, (4.8)
\det\left[\begin{array}[]{ll}\frac{\partial^{{2}}H_{{0}}}{\partial I_{{i}}\partial I_{{j}}}&\frac{\partial H_{{0}}}{\partial I_{{i}}}\\
\frac{\partial H_{{0}}}{\partial I_{{j}}}&0\end{array}\right]. (4.9)

Warunek (4.8) oznacza, że częstości \omega _{{i}}(I) znieniają się niezależnie i dosyć szybko wraz ze zmianą działań I_{{j}}, natomiast warunek (4.9) oznacza, że te częstości zmianiają się szybko i w miarę niezależnie po ograniczeniu do poziomic \left\{ H_{{0}}=const\right\}.

Twierdzenie 4.6 (Kołmogorow–Arnold–Moser). Jeśli są spełnione warunki niezdegenerowania (4.8) i (4.9) dla H_{{0}}, to dla małego zaburzenia H=H_{{0}}+\varepsilon H_{{1}} większość torusów niezmienniczych \left\{ I=const\right\} nie znika, ale tylko lekko deformuje się i ruch na nich jest dalej prawie okresowy.

To twierdzenie zostało sformułowane w 1954 roku na Międzynarodowym Kongresie Matematyków w Amsterdamie, ale na ścisły dowód musiało czekać do początku lat sześdziesiątych. Podali go niezależnie V. Arnold (w przypadku analitycznym) i J. Moser (w przypadku gładkim klasy C^{{333}}). Później klasa gladkości została obniżona do C^{{3}}. Oczywiście nie jestem w stanie przedstawić tego dowodu tutaj.

Przykład 4.7. (Płaski ograniczony problem trzech ciał) Płaskie ograniczone zagadnienie trzech ciał jest to układ w którym dwa ciała (oddziałujące na siebie siłą grawitacji) obracają się za stała prędkością kątową wokół ich środka masy (w początku układu współrzędnych) a trzecie ciało porusza się w płaszczyźnie obrotu dwu ciał i ma masę tak małą, że nie zakłóca ich ruchu. Na Rysunku 4.3 mamy taki układ, w którym S oznacza Słońce, J -Jowisz, A zaś jest Asteroidem. Jednostki czasu, długości i masy można dobrać tak, aby prędkość kątowa, suma mas S i J oraz stała grawitacyjna były równe 1. Wtedy też odległość między S i J też równa się 1. Jedynym parametrem charakteryzującym układ jest masa Jowisza \mu.

Równania ruchu Asteroidu są hamiltonowskie z hamiltonianem

\frac{1}{2}(p_{{1}}^{{2}}+p_{{2}}^{{2}})-\frac{1-\mu}{\rho _{{1}}}-\frac{\mu}{\rho _{{2}}},

gdzie \rho _{{1}} i \rho _{{2}} są odległościami A od S i J odpowiednio. Zauważmy, że położenia S i J zmieniają się z czasem: J=\left(1-\mu\right)(\cos t,\sin t), S=(-\mu)\left(\cos t,\sin t\right); zatem hamiltoniam zależy bezpośredno od czasu.

Rys. 4.3. Problem trzech ciał i niezmiennicze torusy.

Aby pozbyć się tej zależności od czasu, dokonujemy następującej zamiany (jednoczesny obrót współrzędnych i pędów)

q^{{\prime}}=M(t)q,\text{ \ \ }p^{{\prime}}=M(t)p,\text{ \ \ }M(t)=\left(\begin{array}[]{ll}\cos t&\sin t\\
-\sin t&\cos t\end{array}\right).

Okazuje się, że w nowych zmiennych układ nadal jest hamiltonowski z nowym hamiltonianem

H=\frac{1}{2}\left(p_{{1}}^{{\prime}}+q_{{2}}^{{\prime}}\right)^{{2}}+\frac{1}{2}\left(p_{{2}}^{{\prime}}-q_{{1}}^{{\prime}}\right)^{{2}}-V(q_{{1}}^{{\prime}},q_{{2}}^{{\prime}}), (4.10)
V=\frac{q_{{1}}^{{\prime 2}}+q_{{2}}^{{\prime 2}}}{2}+\frac{1-\mu}{\rho _{{1}}}+\frac{\mu}{\rho _{{2}}},

\rho _{{1}}^{{2}}=\left(q_{{1}}^{{\prime}}+\mu\right)^{{2}}+q_{{2}}^{{\prime 2}}, \rho _{{2}}^{{2}}=\left(q_{{1}}^{{\prime}}+\mu-1\right)^{{2}}+q_{{2}}^{{\prime 2}} (Zadanie 4.13). W nowych zmiennych q_{{1}}^{{\prime}},q_{{2}}^{{\prime}} ciała\  S i J spoczywają.

Punkty równowagi układu hamiltonowsiego to punkty krytyczne funkcji hamiltona (Zadanie 4.14). W przypadku hamiltonianu (4.10) te punkty, które nazywamy względnymi położeniami równowagi, zadane są przez

p_{{1}}^{{\prime}}=-q_{{2}}^{{\prime}},\text{ \ \ }p_{{2}}^{{\prime}}=q_{{1}}^{{\prime}},\text{ \ \ }\partial V/\partial q_{{1}}^{{\prime}}=\partial V/\partial q_{{2}}^{{\prime}}=0.

Mamy

\displaystyle\frac{\partial V}{\partial q_{{2}}^{{\prime}}} \displaystyle= \displaystyle q_{{2}}^{{\prime}}\left(1-\frac{1-\mu}{\rho _{{1}}^{{3}}}-\frac{\mu}{\rho _{{2}}^{{3}}}\right)=q_{{2}}^{{\prime}}f,
\displaystyle\frac{\partial V}{\partial q_{{2}}^{{\prime}}} \displaystyle= \displaystyle q_{{1}}^{{\prime}}f-\mu(1-\mu)\left(\frac{1}{\rho _{{1}}^{{3}}}-\frac{1}{\rho _{{2}}^{{3}}}\right).

Mamy dwie możliwości:

1. q_{{2}}^{{\prime}}=0; tutaj znajdujemy trzy punkty tzw. współliniowe punkty libracji L_{{1}}, L_{{2}}, L_{{3}} (Zadanie 4.15), które okazują się niestabilne.

2. f=0 i \rho _{{1}}=\rho _{{2}}=1; tutaj mamy dwa tzw. trójkątne punkty libracji L_{{4}} i L_{{5}}, które leżą w wierzchołkach dwóch trójkątów równobocznych o podstawie \overline{SJ}.

Wyliczenia, których nie przeprowadzamy, pokazują, że dla 27\mu(1-\mu)>1 punkty L_{{4,5}} są niestabilne natomiast w przeciwnym przypadku, tj. dla \mu<\mu _{{1}}=\frac{1}{2}(1-\sqrt{23/27})\approx 0.03852, wartości własne części liniowej układu hamiltonowskiego są postaci \pm i\omega _{{1}}, \pm i\omega _{{2}}, gdzie \omega _{{1}}<0<\omega _{{2}}\not=\omega _{{1}}. Jesteśmy na granicy obszaru stabilności.

Ponadto część kwadratowa H w punkcie L_{{4}} przyjmuje postać

H_{{0}}=\frac{1}{2}\omega _{{1}}(\tilde{p}_{{1}}^{{2}}+\tilde{q}_{{1}}^{{2}})+\frac{1}{2}\omega _{{2}}\left(\tilde{p}_{{2}}^{{2}}+\tilde{q}_{{2}}^{{2}}\right)

w odpowiednim układzie wspólrzędnych w otoczeniu L_{{4}} (patrz [19]). Jest to Hamiltonian układu zupełnie całkowalnego ze zmiannymi działanie–kąt I_{{1}}=\frac{1}{2}(\tilde{p}_{{1}}^{{2}}+\tilde{q}_{{1}}^{{2}}), I_{{2}}=\frac{1}{2}\left(\tilde{p}_{{2}}^{{2}}+\tilde{q}_{{2}}^{{2}}\right), \varphi _{{1}}=\arg\left(\tilde{q}_{{1}}+i\tilde{p}_{{1}}\right), \varphi _{{2}}=\arg\left(\tilde{q}_{{2}}+i\tilde{p}{}_{{2}}\right) i z H_{{0}}=\omega _{{1}}I_{{1}}+\omega _{{2}}I_{{2}} (Zadanie 4.16).

Mamy sytuację jak w Twierdzeniu KAM: H=H_{{0}}+H_{{1}}, gdzie H_{{0}} jest zupełnie całkowalny a H_{{1}} zawiera wyrazy rzędu >2 ze względu na I_{{j}} (które są małe). Niestety, to nie wystarcza, ponieważ częstości \omega _{{j}}=\partial H_{{0}}/\partial I_{{j}} są stałe, a z warunku niezdegenerowania (4.8) powinny się zmianiać wraz z I_{{j}}. Należy więc uwzględnić jeszcze dalsze wyrazy rozwinięcia H w otoczeniu L_{{4}}.

Dokładniej, dokonujemy uproszczenia wyrazów rzędu trzeciego i czwartego w hamiltonianie H. To uproszczenie jest analogiem formy normalnej Poincarégo–Dulaca i zostało udowodnione przez G. Birkhoffa w Twierdzeniu 4.9 poniżej. Ta forma normalna  Birkhoffa w naszym przypadku ma następującą postać

H=H_{{0}}+H_{{1}},\text{ \ }H_{{0}}=\omega _{{1}}I_{{1}}+\omega _{{2}}I_{{2}}+\sum\omega _{{ij}}I_{{i}}I_{{j}},\text{ \ }I_{{j}}=\frac{1}{2}(P_{{j}}^{{2}}+Q_{{j}}^{{2}}), (4.11)

gdzie P_{{j}}=\tilde{p}_{{j}}+\ldots, Q_{{j}}=\tilde{q}_{{j}}+\ldots są nowymi zmiennymi a H_{{1}} zawiera wyrazy rzędu \geq 5 (oraz H_{{0}} i H_{{1}} są inne niż powyżej). W założeniu twierdzenia Birkhoffa pojawia się warunek braku relacji rezonansowych rzędu 4 i 3. Okazuje się, że takie relacje zachodzą dla wartości \mu _{{2}}=\frac{1}{2}\left(1-\sqrt{1833}/45\right)\approx 0.02429 i \mu _{{3}}=\frac{1}{2}\left(1-\sqrt{213}/15\right)\approx 0.01352; zatem te wartości parametru \mu należy wykluczyć.

Hamiltonian H_{{0}}=H_{{0}}(I_{{1}},I_{{2}}) jest hamiltonianem zupełnie całkowalnym i ma szansę na spełnienie warunków niezdegenerowania (4.8) i (4.9). Okazuje się, że tylko warunek (4.9) jest istotny. A. Leontowicz pokazał, że może on zostać naruszony tylko dla dyskretnego zbioru wartości parametru \mu.\footnote{W \@@cite{[\@@bibref{Number}{Zol1}{}{}]} można dowiedzieć się, że warunek (4.9) zostaje
naruszony dla dokłanie jednej konkretnej wartości parametru $\mu.$ Ta
wartość wynikała ze wzoru na wyznacznik w równaniu (4.9)
podanego przez francuskich astronomów A. Deprit i A. Deprit-Bartholom\'{e} (i cytowanego w bardzo poważnych monografiach). Ostatnio z moją
magistrantką W. Barwicz odkryliśmy, że ten wzór jest
nieprawdziwy, a nawet sprzeczny z wyliczeniami Leontowicza. W istocie, tenże wyznacznik jest bardzo skomplikowaną funkcją algebraiczną
od $\mu,$ króra nie jest tożsamościowo równa zeru.} Załóżmy zatem, że \mu spełnia wszystkie warunki wypisane powyżej, czyli jest prawdziwa teza twierdzenia KAM.

Jak z twierdzenia KAM wynika stabilność? Otóż znajdujemy się w przestrzeni 4-wymiarowej w otoczeniu punktu równowagi. Ponieważ układ jest hamiltonowski z hamiltonianem niezależnym od czasu, więc ruch odbywa się po powierzchniach H=\textrm{const}. Są one trójwymiarowe. Z twierdzenia KAM wynika, że każda taka powierzchnia jest prawie zapełniona torusami niezmienniczymi \mathbb{T}^{{2}}, których jest tym więcej im bliżej jesteśmy torusa I_{{1}}=I_{{2}}=0. Każdy torus niezmienniczy rozbija powierzchnię H=const na dwie części, swoje wnętrze i zewnętrze. Zaden punkt z wnętrza nie wychodzi zeń w trakcie ewolucji. Ponieważ w przestrzeni zmiennych P,Q torusy mogą być dowolnie małe, to wynika stąd stabilność w sensie Lapunowa.

Uzupełnimy powyższy przykład. Załóżmy, że mamy hamiltonian w postaci

H=\sum\omega _{{j}}\cdot\frac{1}{2}(p_{{j}}^{{2}}+q_{{j}}^{{2}})+\ldots.

Definicja 4.8. Mówimy, że `częstości' \omega _{{j}} spełniają relację resonansową rzędu d, jeśli isnieją liczby całkowite k_{{1}},\ldots,k_{{n}} z \sum\left|k_{{j}}\right|=d takie, że

\sum k_{{j}}\omega _{{j}}=0.

Twierdzenie 4.9 (Birkhoff). Jeśli częstości \omega _{{j}} nie spełniają żadnej relacji rezonansowej rzędu \leq 2m, to istnieje kanoniczna zamiana zmiennych \left(p,q\right)\longmapsto\left(P,Q\right)=\left(p+\ldots,q+\ldots\right) prowadząca do hamiltonianu

H=\sum _{{\left|l\right|\leq m}}a_{{l}}I^{{l}}+O\left(\left|\left(p,q\right)\right|^{{2m+1}}\right),

gdzie I_{{j}}=\frac{1}{2}(P_{{j}}^{{2}}+Q_{{j}}^{{2}}) i sumowanie przebiega po wielowskaźmikach \left(l_{{1}},\ldots,l_{{n}}\right) z \left|l\right|=l_{{1}}+\ldots+l_{{n}} iI^{{l}}=I_{{1}}^{{l_{{1}}}}\ldots I_{{n}}^{{l_{{n}}}}.

Uwaga 4.10. Zamiana \left(p,q\right)\longmapsto\left(P,Q\right), występująca w powyższym twierdzeniu jest kanoniczna jeśli

\sum dp_{{j}}\wedge dq_{{j}}=\sum dP_{{j}}\wedge dQ_{{j}}.

Okazuje się, że po kanonicznej zamianie zmiennych układ hamiltonowski przechodzi w układ hamiltonowski (patrz [4]).

ZADANIA

Zadanie 4.11. Pokazać, że jeśli funkcja hamiltona H nie zależy bezpośrednio od czasu, to jest całką pierwszą dla układu (4.2).

Zadanie 4.12. Pokazać, że pole wektorowe zadane wzorem (4.2) ma zerową dywergencję. Wywnioskować stąd, że odpowiedni potok fazowy zachowuje objętość.

Zadanie 4.13. Udowodnić wzór (4.10).

Zadanie 4.14. Pokazać, że jeśli H nie zależy bezpośrednio od czasu, to punkty równowagi układu (4.2) są dokładnie punktami krytycznymi funkcji H.

Zadanie 4.15. Pokazać, że istnieją dokładnie trzy współliniowe punkty libracji.

Zadanie 4.16. Pokazać, że hamiltonian postaci H_{{0}}=\omega _{{1}}I_{{1}}+\omega _{{2}}I_{{2}} (lub jak we wzorze (4.11)) jest hamiltonianem układu zupełnie całkowalnego.

Zadanie 4.17. Zastosować metodę całek abelowych (Przykład 4.5) do pokazania, że układ van der Pola \dot{x}=y, \dot{y}=-x-a(x^{{2}}-1)y dla małego parametru a>0 posiada dokładnie jeden cykl graniczny.

4.3. Drgania relaksacyjne

Zacznijmy od znanego przykładu.

Przykład 4.18 (Układ van der Pola).

\dot{x}=y-x^{{3}}+x,\text{ \ \ }\dot{y}=-\varepsilon x.

(Gdy \varepsilon=1 i położyć y_{{1}}=y-x^{{3}}+x, to dostaje się \dot{x}=y_{{1}}, \dot{y}_{{1}}=-x-(3x^{{2}}-1)y_{{1}}; z dokładnością do przeskalowania jest to układ z Przykładu 2.34.)

Rys. 4.4. Układ van der Pola typu `wolny–szybki'.

Widać, że x zmienia się szybko w porównaniu z y; mówimy, że x jest szybką zmienną a y wolną. Dla \varepsilon=0 mamy y=const i w istocie mamy równanie na x zależne od parametru y (teoria bifurkacji się kłania, patrz Rysunek 4.4). Gdy \varepsilon\not=0 (ale małe), to fizycy powiedzieliby, że parametr y `płynie'. Oczekuje się istnienia cyklu granicznego \gamma _{{\varepsilon}} (w istocie \gamma _{{\varepsilon}} jest stabilny) dążącego do kawałkmi gładkiej krzywej \gamma _{{0}} przedstawionej na Rysunku 4.5. Cykl \gamma _{{0}} składa się z:

— kawałków ruchu powolnego wzdłuż krzywej y=x^{{3}}-x (gdzie \dot{x}=0),

— odcinków skoku wzdłuż prostych y=const.

Taki ruch jest przykładem dragań relaksacyjnych (jak bicie serca).

Rys. 4.5. Drgania relaksacyjne.

Rozważmy teraz ogólną sytuację. Mamy układ niezaburzony

\dot{x}=f(x,y),\text{ \ \ \ }\dot{y}=0,

(x\in\mathbb{R}^{{k}}, y\in\mathbb{R}^{{l}}); tutaj x to szybkie współrzędne a y to wolne współrzędne. Mamy też układ zaburzony

\dot{x}=F(x,y;\varepsilon),\text{ \ \ }\dot{y}=\varepsilon G(x,y;\varepsilon),\text{ \ \ \ }F(x,y;0)=f(x,y).

Definicja 4.19. Powierzchnia S=\left\{ f(x,y)=0\right\} nazywa się powolną powierzchnią.

Powolna powierzchnia dzieli się na obszary stabilności i niestabilności układu niezaburzonego; odpowiadają one sytuacjom gdy \textrm{Re}\lambda _{{j}}(A)<0j=1,\ldots,k, A=\frac{\partial f}{\partial x}, i gdy istnieje \textrm{Re}\lambda _{{j}}(A)>0.

Na powolnej powierzchni mamy pole wektorowe definiowane następująco. Bierzemy pole

\frac{\partial}{\partial\varepsilon}\left(F\partial _{{x}}+\varepsilon G\partial _{{y}}\right)|_{{\varepsilon=0}}=f_{{1}}(x,y)\partial _{{x}}+g(x,y)\partial _{{y}}

w punkcie \left(x,y\right)\in S i rzutujemy je na T_{{(x,y)}}S wzdłuż zmiennych y. Jest to pole ruchu powolnego.

Przypomnę, że na początku tego rozdziału mówiliśmy, że drgania relaksacyjne charakteryzują się własnością, ze mały parametr występuje po lewej stronie. Aby się o tym przekonać wprowadzamy czas powolny \tau=\varepsilon t. Wtedy dostajemy układ

\varepsilon\frac{dx}{d\tau}=f(x,y)+O(\varepsilon),\text{ \ \ }\frac{dy}{d\tau}=g(x,y)+O(\varepsilon).

Teraz równanie ruchu powolnego na S (lokalnie parametryzowanej przez y) jest postaci

\frac{dy}{d\tau}=h(y)+O(\varepsilon)

(z odpowiednią funkcją h).

Przeanalizujmy ruch typowego punktu \left(x_{{0}},y_{{0}}\right). Składa się on z kawałków trzech rodzajów: dochodzenie do powierzchni powolnej, ruch wzdłuż powierzchnii powolnej i ruch w obszarze przejściowym.

Rys. 4.6. Dochodzenie do powierzchni powolnej.

4.20. Dochodzenie do powierzchni powolnej. Niech punkt (x_{{0}},y_{{0}}) spoza S rzutuje się (wzdłuż współrzędnych y) na punkt \left(x_{{\ast}},y_{{0}}\right), x_{{\ast}}=x_{{\ast}}(y_{{0}}) , na S w obszarze stabilności (patrz Rysunek 4.6). To znaczy, że punkt x_{{0}} leży w basenie przyciągania punktu x_{{\ast}} dla równania \dot{x}=f(x,y_{{0}}) (y_{{0}} stałe). Rozważmy obszar U=\left\{\left|x-x_{{\ast}}(y_{{0}})\right|<\delta,\text{ }y_{{0}}\in V\right\}, gdzie V jest pewnym obszarem odpowiadającum podzbiorowi obszaru stabilności w S. Okazuje się, że powolny czas dochodzenia rozwiązania z warunkiem początkowym \left(x_{{0}},y_{{0}}\right) do U jest rzędu \tau _{{1}}\sim C_{{1}}\varepsilon\left|\ln\varepsilon\right|, co odpowiada rzeczywistemu czasowi

t_{{1}}\sim C_{{1}}\left|\ln\varepsilon\right|

(stała C_{{1}} zależy od U i od F,G).

4.21. Ruch powolny. W obszarze U mamy ruch powolny, opisywany równaniem dy/d\tau=h(y)+O(\varepsilon). Trwa on do momentu \tau _{{2}}=T=O(1), co odpowiada długiemu czasowi rzeczywistemu t_{{2}}=T/\varepsilon.

4.22. Ruch w obszarze przejściowym. Obszar przejściowy leży blisko granicy pomiędzy obszarami stabilności i niestabilności w S. Mamy dwie typowe możliwości (jak w teorii bifurkacji):

A.\lambda _{{1}}(A)=0 (gdzie A=\frac{\partial f}{\partial x}|_{{f=0}});

B.\textrm{Re}\lambda _{{1,2}}=0.

A. Zryw. Ten przypadek (który odpowiada bifurkacji siodło–węzeł) zanalizujemy dla sytuacji gdy x\in\mathbb{R} i y\in\mathbb{R} (można do tego wszystko zredukować). Po odpowiednich przeskalowaniach mamy następujący układ

Rys. 4.7. Zjawisko `zrywu'.
\dot{x}=x^{{2}}-y+\ldots,\text{ \ \ }\dot{y}=-\varepsilon+\ldots

Dokonujemy normalizacji

\varepsilon=\mu^{{3}},\text{ \ \ }x=\mu X,\text{ \ \ }y=\mu^{{2}}Y.

łatwo sprawdzić, że prowadzi to do pola

\dot{X}=\mu\left\{ X^{{2}}-Y+O(\mu)\right\},\text{ \ \ }\dot{Y}=\mu\left\{-1+O(\mu)\right\}

orbitalnie równoważnemu polu \left(X^{{2}}-Y\right)\partial _{{X}}-\partial _{{Y}}. Jego portret fazowy jest zadany równaniem Riccatiego

dX/dY=X^{{2}}-Y (4.12)

i jest przedstawiony na Rysunku 4.717Równanie (4.12) jest chyba najprostszym przykładem równania różniczkowego, którego nie można rozwiązać w tzw. kwadraturach.. Zjawisko, które tutaj obserwujemy nosi nazwę zrywu.

B. Opóźnienie utraty stabilności. W tym przypadku, który odpowiada bifurkacji Andronowa–Hopfa, problem redukuje się do następującego modelowego układu

\dot{z}=\left(y+i\omega\right)z+cz\left|z\right|^{{2}},\text{ \ \ }\dot{y}=\varepsilon, (4.13)

zx_{{1}}+ix_{{2}}\in\mathbb{C}\simeq\mathbb{R}^{{2}},y\in\mathbb{R}. Oczywiście y=\varepsilon t jest `płynącym' parametrem. Załóżmy jeszcze, że

c=-1;

przypadek c>0 jest mniej ciekawy. Dla amplitudy r=\left|z\right| dostajemy równanie Bernoulliego

\dot{r}=r\left(\varepsilon t-r^{{2}}\right).

Połóżmy warunek początkowy

y(t_{{0}})=-\mu,\text{ \ \ \ }r(t_{{0}})=r_{{0}},\text{ \ \ }t_{{0}}=-\mu/\varepsilon,

gdzie \mu>0 jest ustaloną (nie za dużą i nie za małą) stałą. To zagadnienie początkowe ma następujące rozwiązanie

r(t)=r_{{0}}\left\{ e^{{\varepsilon\left(t_{{0}}^{{2}}-t^{{2}}\right)}}+2r_{{0}}^{{2}}\int _{{t_{{0}}}}^{{t}}e^{{\varepsilon\left(s^{{2}}-t^{{2}}\right)}}ds\right\}^{{-1/2}} (4.14)

(Zadanie 4.24). Zbadamy asymptotyczne zachowanie się tego rozwiązania przy \varepsilon\rightarrow 0 dzieląc zakres czasu t na cztery obszary:

(a) 0<t-t_{{0}}<O(1), czyli 0<y+\mu<O(\varepsilon).

Niech u=t-t_{{0}}. Wtedy \varepsilon\left(t_{{0}}^{{2}}-t^{{2}}\right)=\varepsilon(t_{{0}}+t)u\approx 2\mu u i \varepsilon\left(s^{{2}}-t^{{2}}\right)\approx 2\mu(u-v), gdzie v=s-t_{{0}}. Zatem

\int _{{t_{{0}}}}^{{t}}e^{{\varepsilon\left(s^{{2}}-t^{{2}}\right)}}ds\approx\int _{{0}}^{{u}}e^{{2\mu(u-v)}}dv=\frac{1}{2\mu}(e^{{2\mu u}}-1)

oraz

r(t)\approx r_{{0}}\left\{ e^{{2\mu u}}+r_{{0}}^{{2}}(e^{{2\mu u}}-1)/\mu\right\}{}^{{-1/2}}

jest malejącą funkcją od u.

(b) y=\varepsilon t jest ustalone tak, że -\mu<y<\mu.

Tutaj e^{{\varepsilon\left(t_{{0}}^{{2}}-t^{{2}}\right)}}\approx e^{{\left(\mu^{{2}}-y^{{2}}\right)/\varepsilon}}\rightarrow\infty. Zatem

r(t)<C_{{1}}e^{{-C_{{2}}/\varepsilon}}\rightarrow 0,

przy czym jest to bardzo szybkie dążenie do zera.

(c) 0<\left|t_{{0}}\right|-t<O(1), czyli 0<\mu-y<O(\varepsilon).

Wprowadźmy zmienną w=\left|t_{{0}}\right|-t. Jak w punkcie (a) mamy e^{{\varepsilon\left(t_{{0}}^{{2}}-t^{{2}}\right)}}\approx e^{{2\mu w}}.

Obszar całkowania dla całki we wzorze (4.14) podzielimy na trzy odcinki: od t_{{0}} do t_{{0}}/2<0, od t_{{0}}/2 do \left|t_{{0}}\right|/2 i od \left|t_{{0}}\right|/2 do t. Przez I_{{1}}, I_{{2}} i I_{{3}} oznaczymy odpowiednie całki. Podobnie jak w punkcie (a) pokazuje się, że I_{{1}}=O(1) i I_{{3}}=O(1). Z rachunków w punkcie (b) wynika, że I_{{2}}\rightarrow 0 bardzo szybko. Zatem

r(t)=O(1).

(d) \left|t_{{0}}\right|<t, czyli y>\mu i jest ustalone. Teraz \exp\left\{\varepsilon\left(t_{{0}}^{{2}}-t^{{2}}\right)\right\}\approx\exp\left\{-(y^{{2}}-\mu^{{2}})/\varepsilon\right\}\rightarrow 0. Następnie \varepsilon\left(s^{{2}}-t^{{2}}\right)\approx(s-t)\cdot 2y dla s bliskich t, tj. dla tych s, dla których wkład do całki jest dominujący. Dostajemy \int^{{t}}e^{{2y(s-t)}}ds\approx\frac{1}{2y}. Stąd

r(t)\approx r_{{0}}\left\{ r_{{0}}^{{2}}/2y\right\}^{{-1/2}}=\sqrt{y}.

Możemy podsumować powyższe obliczenia.

Rys. 4.8. Zjawisko opóźnienia utraty stabilności.

Twierdzenie 4.23.W przypadku B opisywanego układem (4.13) z c<0 zachodzi zjawisko opóźnienia utraty stabilności. Polega ono na tym, przy zmianie zmiennej y (która jest współczynnikiem stabilności ruchu niezaburzonego) od wartości ujemnej y(t_{{0}})=-\mu do wartości dodatniej \mu układ (względem z) jest cały czas stabilny, a zmiana stabilności rozwiązania następuje dla parametru y=\mu, przy czym dalej amblituda oscylacji rośnie jak w zwykłej bifurkacji Andronowa–Hopfa.

Zjawisko opóźnienia utraty stabilności można objaśnić fizycznie.Zmienna y jest ujemna przez bardzo długi czas, rzędu 1/\varepsilon. Wtedy układ fizyczny zdąży podejść bardzo blisko położenia równowagi; na tyle blisko, że potrzeba potem tyle samo czasu, aby od położenia równowago odejść (patrz Rysunek 4.8).

ZADANIA

Zadanie 4.24. Udowodnić wzór (4.14).

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.