Symulowanie i analiza układów dynamicznych

Spis treści

1 Układy dynamiczne

Układ dynamiczny będziemy modelować za pomocą układu autonomicznych równań różniczkowych zwyczajnych. To są takie równania, które w sposób jawny nie uwzględniają czasu i zawierają tylko jedną zmienną niezależną. Układ autonomiczny nazywamy jednowymiarowym, jeśli składa się z jednej zmiennej zależnej od czasu np. x(t); do opisania takiego układu stosuje się równania różniczkowe pierwszego rzędu; dwu-wymiarowy składa się z dwóch zmiennych np. x(t) i y(t) i tutaj stosuje się równania różniczkowe drugiego rzędu. Analogicznie można opisać układy wyższych rzędów.

Równanie różniczkowe zwyczajne jest równaniem zawierającym jedną lub więcej funkcji z tylko jedną zmienną niezależną oraz pochodne tych funkcji, np. równanie: \(\ddot x + \epsilon \dot x (x^2-1) + x=0\), w którym zmienną niezależną jest czas, a zmienną zależną jest x. Jednym z głównych celów analizy jest poszukiwanie nieznanej krzywej całkowej określającej zależność x od t. W przypadku układów dynamicznych chcemy wiedzieć, jak badane zjawisko zmienia się w każdej chwili czasu.

Równanie różniczkowe zwyczajne drugiego rzędu postaci: \[\ddot x = f(t,x,\dot x)\] np. równanie oscylatora harmonicznego:
\[\text{ postać nieautonomiczna: }m\ddot x + b\dot x + kx = f(t)\] \[\text{ postać autonomiczna: }m\ddot x + b\dot x + kx = 0\]

można sprowadzić do postaci układu dwuwymiarowego przez podstawienie: \[\dot x = y\] \[\dot y = -by - kx\]

Układ jednowymiarowy autonomiczny: \[\dot x = f(x)\]

Układ dwuwymiarowy autonomiczny: \[\dot x = f(x,y)\] \[\dot y = g(x,y)\]

Omawianie układów dynamicznych rozpoczniemy od jednowymiarowych tzn., takich które poruszają się po linii bądź okręgu. Złożoność ruchu jest bardzo ograniczona - prędkość w takich układach monotoniczne rośnie bądź maleje. W układach dwuwymiarowych linowych będziemy mogli zauważyć dodatkowe typy punktów stałych, takich jak: środek (orbita okresowa) i wiry, ale nie zauważymy cykli granicznych, które są charakterystyczne dla nieliniowych układów. Dopiero w trójwymiarowych układach nieliniowych pojawia się chaos.

2 Układy jednowymiarowe

2.1 Liniowe i nieliniowe równania różniczkowe pierwszego rzędu

Równanie różniczkowe zwyczajne pierwszego rzędu jest równaniem różniczkowym zawierającym jedną lub więcej funkcji z tylko jedną zmienną niezależną oraz pochodne pierwszego rzędu tych funkcji, np. równanie: \(\dot x - x^2 + k=0\), w którym zmienną niezależną jest np. czas t, a zmienną zależną jest x. To równanie moglibyśmy zapisać takiej formie, aby pokazać zależność od czasu: \(\frac{dx}{dt} = x(t)^2 - k\). W dalszej części będę korzystał z wersji uproszczonej: \(\dot x = x^2 - k\). Postać ogólna takiego równania: \(\dot x = f(x)\)

f jest funkcją rzeczywistą lub inaczej - funkcją o wartościach rzeczywistych (ang. real-valued function). Zakładamy, że jej pochodna każdego rzędu istnieje i jest ciągła - zapisujemy to w taki sposób: \(f \in C^{\infty}\)

Funkcję f nazywamy funkcją regularną rzędu k na zbiorze U, jeżeli wszystkie pochodne cząstkowe funkcji f do rzędu k włącznie istnieją i są ciągłe w całej dziedzinie U. Mówimy wtedy, że funkcja jest klasy \(C^k(U)\) i piszemy \(f \in C^k(U)\)

Przykłady:

  • \(f \in C^0\) - funkcja ciągła
  • \(f \in C^1\) - pochodne cząstkowe pierwszego rzędu istnieją i są ciągłe
  • \(f \in C^\infty\) - funkcja gładka

Dla układu jednowymiarowego autonomicznego równanie \(\dot x = f(x)\) definiuje nam pole nachyleń. Dla dowolnego miejsca na linii fazowej (ang. phase-line) możemy obliczyć pochodną, która jest styczną do szukanej krzywej całkowej. Linia fazowa jest to linia wyznaczona przez pochodne równania dla konkretnego czasu t, czyli na rysunku poniżej jest to linia wzdłuż osi rzędnych.

pole_nachylen.png

Figure 1: Pole nachyleń dla \(\dot x = -x\)

par a = -1
x' = a*x
s'=1
init x=-2

@ total=5, bounds=10000, xp=S,yp=X,xlo=0,ylo=-2,xhi=5,yhi=2
done

GUI
# Initialconds –> Range —> dobranie odpowiednich parametrów –> OK
# Dir.field/flow —> Direct Field —> dobranie siatki —> ENTER

Rozwiązaniem równania w sensie geometrycznym jest krzywa całkowa x(t), którą w przestrzeni fazowej układu dla wymiarów większych od 1 możemy przedstawić jako krzywą parametryczną zależną od czasu t: \(x=f(t),\;y=g(t)\). Tej zależności od czasu właśnie szukamy, a której nie widać w równaniu ruchu. Nasza krzywa musi spełniać równanie: \(\dot x = f(x)\), \(\dot x = -x\). Rozwiązaniem ogólnym tego równania jest \(x=x(t)=x_{0}e^{-t}\). To rozwiązanie jest poprawne, jeśli jego pochodna będzie równa -x. \(\dot x = - (x_{0}e^{-t}) = -(x)\). Wszystko się zgadza.

Punkt \(x_*\), dla którego zachodzi \(f(x_*)=0\) nazywamy punktem równowagi (ang. equalibrium) lub inaczej punktem stałym (ang. fixed point). Punkt ten odpowiada rozwiązaniu stałemu \(x(t) = x_*\). W tym miejscu nasz układ się nie rusza. Wiedza na temat istnienia punktów stałych w układzie pozwala również na określenie obszarów w przestrzeni fazowej, gdzie pochodna zmienia znak! Gdy pochodna zmienia znak to kierunek ruchu się zmienia!

Szukanie rozwiązania równania różniczkowego pierwszego rzędu, którego celem jest znalezienie krzywej całkowej x(t), możemy podzielić na rozwiązania w punkcie równowagi, rozwiązanie szczególne, rozwiązanie ogóle. W przypadku dwuwymiarowych układów będziemy również szukać rozwiązań wzdłuż wektorów własnych.

Podstawowe twierdzenie algebry: każde równanie algebraiczne stopnia n ma n pierwiastków. k-krotne pierwiastki liczymy k razy np. dla równania kwadratowego, w którym Δ=0, mamy jedno rozwiązanie rzeczywiste (2 pierwiastki równe)

Zgodnie z tym twierdzeniem to równanie algebraiczne pierwszego rzędu (liniowe), posiada tylko jedno rozwiązanie:
Postać uporządkowana: \(ax+b=0\)
Rozwiązanie: \(x_{*}=-\frac{b}{a}\)

Równanie algebraiczne rzędu większego niż 1 jest równaniem nieliniowym. Ma rozwiązań w punkcie równowagi minimum jedno. Już na tym etapie widzimy, że w przypadku równań nieliniowych zachowanie układu, opisane tymi równaniami, rośnie.

Rozwiązanie szczególne jest rozwiązaniem dla konkretnych warunków początkowych np. \(x(t)=e^{-t}\) - jest rozwiązaniem dla \(x_{0} = 1\)

Rozwiązanie ogólne jest rozwiązaniem dla dowolnych warunków początkowych np. \(x(t)=x_{0}e^{-t}\), rozwiązanie to jest odpowiedzią na problem wartości początkowej t.j. \(\dot x = f(x)\), \(x(0)=x_{0}\) Rozwiązanie istnieje dla dowolnej wartości początkowej i jest jednoznaczne.

\(\dot x - ax^2 = 0\) - jest to nieliniowe równanie różniczkowe pierwszego rzędu, gdyż równanie algebraiczne jest kwadratowe. Jest to układ autonomiczny, gdyż w sposób jawny nie zależy od czasu t. Dodatkowo jest to równanie jednorodne, gdyż nie ma wolnych współczynników i przez to posiada minimum jedno rozwiązanie x(t)=0, które nazywamy trywialnym.

Przykład równania niejednorodnego: \(\dot x = ax^2 - k\), które posiada wolny współczynnik k. Przykład równania liniowego jednorodnego \(\dot x = ax\).

Liniowe równanie niejednorodne nieautonomiczne:
\(f(x,t) = t,\; \dot x = t,\; x=0.5t^2 + x_{0},\; x(t) \ne 0\)

Liniowe równanie niejednorodne autonomiczne:
\(f(x,t) = x-1,\; \dot x = x - 1,\; \ln{\bigl|\frac{x-1}{x_{0}-1}\bigr|}=t,\;x=(x_{0}-1)e^t+1,\;x(t) \ne 0\)

Liniowe równanie jednorodne autonomiczne:
\(f(x,t) = x,\; \dot x = x,\; ln{|\frac{x}{x_{0}}|}=t,\; x=x_{0}e^t,\;\text{dla}\;x_{0}=0 \rightarrow x(t) = 0\)

Liniowe równanie jednorodne nieautonomiczne:
\(f(x,t) = xt,\; \dot x = xt,\; ln{|\frac{x}{x_{0}}|}=0.5t^2,\;x=x_{0}e^{0.5t^2},\;\text{dla}\;x_{0}=0 \rightarrow x(t) = 0\)

2.1.1 Przykład 1: \(\dot x = ax-k, \quad k=0\)

Rozwiązanie analityczne liniowego równania różniczkowego:

$$\frac{dx}{dt} = ax$$ $$\frac{1}{ax}dx = dt$$ $$\int \frac{1}{ax}dx = \int dt$$ $$\frac{1}{a}\ln{|x|} = t + C$$ $$x(0) = x_0 \rightarrow C = \frac{1}{a}\ln{|x_{0}|} - warunek\quad początkowy$$ $$\frac{1}{a}\ln{|x|} - \frac{1}{a}\ln{|x_{0}|} = t$$ $$t = \frac{1}{a}\ln{\biggl|\frac{x}{x_{0}}\biggr|}$$ $$\biggr|\frac{x}{x_{0}}\biggl| = e^{at}$$ $$dla \;x<0,\;\biggr(\frac{x}{x_{0}}\biggl) = e^{at}$$ $$dla \;x>0,\;\biggr(\frac{x}{x_{0}}\biggl) = e^{at}$$ $$x(t)=x_{0}e^{at}$$

x(t) jest rozwiązaniem ogólnym równania i określa nam trajektorię - wiemy jak się zmienia x w zadanym czasie w zależności od warunków początkowych. Trajektoria punktu fazowego znajduje się natomiast w przestrzeni fazowej. Dla układów pierwszego rzędu jest to przestrzeń jednowymiarowa, czyli linia fazowa.

linear1d.png

Figure 2: Rozwiązanie numeryczne równania dla a = 1

par a=1
x'=a*x 
init x=-2
@ total=5, bounds=10000, xlo=0,ylo=-2,xhi=5,yhi=2
done
# W GUI:
# Initialconds –> Range —> dobranie odpowiednich parametrów –> OK

Zachowanie układu zależy od parametru a oraz warunków początkowych x0. Jeśli \(a > 0\) to układ jest niestabilny i punktem stałym jest \(x_{*}=0\) 2.

2.1.2 Przykład 2a: \(\dot x = x^2-k,\quad k=1\)

Rozwiązanie analityczne:

$$\frac{dx}{dt} = x^2 - 1$$ $$\frac{1}{x^2 - 1}dx = dt$$ $$\int \frac{1}{(x-1)(x+1)}dx = \int dt$$ $$\frac{1}{(x-1)(x+1)} \equiv \frac{A}{x-1} + \frac{B}{x+1}$$ $$1 \equiv A(x+1) + B(x-1)$$ $$1 \equiv x(A+B) + (A-B)$$ $$A+B=0, A-B=1, A=1/2, B=-1/2$$ $$\frac{1}{2}\int\frac{1}{x-1} - \frac{1}{2}\int\frac{1}{x+1}$$ $$\frac{1}{2}\ln{|x-1|} - \frac{1}{2}\ln{|x+1|} = t + C$$ $$x(0) = x_0 \rightarrow C = \frac{1}{2}\ln{\biggl|\frac{x_{0}-1}{x_{0}+1}\biggr|} - warunek\quad początkowy$$ $$\frac{1}{2}\ln{\biggl|\frac{x-1}{x+1}\biggr|} - \frac{1}{2}\ln{\biggl|\frac{x_{0}-1}{x_{0}+1}\biggr|} = t$$ $$\biggl|\frac{x-1}{x+1}\biggr| \cdot \biggl|\frac{x_{0}+1}{x_{0}-1}\biggr| = e^{2t}$$ $$dla\; x<-1,\; \biggl(\frac{x-1}{x+1}\biggr) \cdot \biggl(\frac{x_{0}+1}{x_{0}-1}\biggr) = e^{2t}$$ $$dla\; x \in (-1,1),\; \biggl(-\frac{x-1}{x+1}\biggr) \cdot \biggl(-\frac{x_{0}+1}{x_{0}-1}\biggr) = e^{2t}$$ $$dla\; x>1,\; \biggl(\frac{x-1}{x+1}\biggr) \cdot \biggl(\frac{x_{0}+1}{x_{0}-1}\biggr) = e^{2t}$$ $$x(t)=\frac{e^{2t}(x_{0}-1)+x_{0}+1}{e^{2t}(1-x_{0})+x_{0}+1}$$

example221_a.png

Figure 3: Rozwiązanie numeryczne równania

x'=x^2-1 
init x=.8
@ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
done
# GUI
# Initialconds –> Range —> dobranie odpowiednich parametrów –> OK

example221.mp4; XPP: x'=x^2-1; @total=5

Zachowanie układu zależy od warunku początkowego x0. Z rysunku 3 możemy wywnioskować istnienie stabilnego punktu stałego x=-1, gdyż trajektoria zmierza właśnie do osiągnięcia tego punktu; x=1 jest niestabilnym punktem stałym. Punkt stały jest to stan układu, w którym nie zachodzą żadne zmiany: \(\dot x = 0\). Dla \(x^2-1=0\) punktami stałymi są \(x_*=1\) i \(x_{*}=-1\).

example221_anim.gif

Figure 4: Animacja ruchu dla x0 = 0.999

Dynamika układu zmienia się w zależności od odległości od punktów stałych. Gdy układ posiada punkt stabilny to dąży do osiągnięcia tego stanu. Punkt niestabilny natomiast to taki punkt, od którego trajektorie układu się oddalają. Z rysunku 3 możemy wyczytać, że jeśli układ zaczyna działać od punktu \(x_{0} = -1\) to pozostaje cały czas w tym punkcie. Gdy \(x_{0} < 1\) to zawsze układ zmierza do punktu stabilnego, czyli do punktu \(x_* = -1\). Dzieje to się jednak z różną dynamiką, która zależy od warunków początkowych. Jeśli x0 będzie bardzo blisko drugiego punktu stałego np. \(x_* = 0.999\) to układ będzie powoli od niego się oddalał (ten przykład widzimy na animacji). W miarę oddalania się o niego jego prędkość będzie się zwiększała, by znowu zmaleć, gdy układ będzie zbliżał się do punktu stabilnego \(x_* = -1\). Inaczej układ się zachowa, gdy wystartuje z punktu x0 = -3. Na początku będzie poruszał się z dużą prędkością, by w miarę zbliżania się do punktu stabilnego, zwalniać.

Przykład obliczeń symbolicznych za pomocą języka programowania Julia:

using Plots, SymPy
@symfuns x;
@vars t positive=true;
eqn = x'(t) - x(t)^2 + 1;
out = dsolve(eqn);
eq = out.rhs;    # prawa strona równania
C1 = first(setdiff(free_symbols(eq), (t))); 
# warunki początkowe x(0)=0 
c1 = solve(eq(t=>0), C1);
eq(C1 => c1[1]);

# warunki początkowe można podać podczas obliczania całki
t0, x0  = 0, 0;
s1 = dsolve(eqn, x(t), ics=(x, t0, x0));
s1[1]

         -1   
x(t) = ───────
       coth(t)

Ta postać jest równoważna naszemu rozwiązaniu dla x0=0

\(\frac{-1}{coth(t)} = -\frac{e^{2t}-1}{e^{2t}+1}\)

\[x(t)=\frac{e^{2t}(x_{0}-1)+x_{0}+1}{e^{2t}(1-x_{0})+x_{0}+1}\]

2.1.3 Przykład 2b: geometryczne spojrzenie na równanie różniczkowe

Nieliniowe, w przeciwieństwie do linowych równań, często nie mają rozwiązania analitycznego. Podstawowym narzędziem do radzenia sobie z takim problemem w analizie układów dynamicznych jest interpretacja równań różniczkowych jako pola wektorowego: każdemu punktowi P w przestrzeni fazowej układu przyporządkowujemy wektor V:

Dla układów jednowymiarowych np. \(V_{x}=\dot x = x^2-1\), \(V_{x} = V_{x}(x)\), \(\vec V = V_{x}\vec i\). Każdy punkt fazowy na linii fazowej jest reprezentowany przez taki wektor5

Dla układów dwuwymiarowych np. \(V_{x}=\dot x = x -y\), \(V_{y}=\dot y = x^2 - 4\). \(V_{x} = V_{x}(x,y),\quad V_{y} = V_{y}(x,y)\), \(\vec V = V_{x}\vec i + V_{y}\vec j\)

W programie XPP: każdemu punktowi P=(x,y) jest przyporządkowany wektor: \((x + sV_{x}(x, y), y + sV_{y}(x, y))\), gdzie s jest parametrem skalującym.

vector_field_line.png

Figure 5: Pole wektorowe na linii fazowej

Wykres na samym dole rysunku 5 nazywa się portretem fazowym i przedstawia jakościowo różne zachowania badanego układu. Równanie różniczkowe interpretujemy jako pole wektorowe w taki sposób, że rzutujemy krzywą równania ruchu \(\frac{dx}{dt}\) na oś X i tworzymy wektory wokół punktów stałych. Jeśli wektor \(V_{x}=\frac{dx}{dt}\) ma wartość dodatnią to jest skierowany zgodnie z osią X. W przeciwnym razie skierowany jest w drugą stronę. Z wykresu dowiadujemy się, że układ ma dwa punkty stałe tzn. takie dla których: \(\dot x = 0\). Te punkty to \(x_* = -1\) i \(x_{**}= 1\). Punkt \(x_*\) jest stabilnym punktem stałym, gdyż nachylenie krzywej równania ruchu jest ujemne. W punkcie \(x_{**}\) nachylenie krzywej jest dodatnie i dlatego jest to punkt niestabilny. Punkty pomiędzy punktami stałymi są przyciągane do punktu stabilnego lub odpychane od punktu niestabilnego.

Można wyróżnić trzy typy punkty stałych: ściek (ang. sink), źródło (ang. source) i połowicznie stabilny (ang. half-stable); dla układów wyższy wymiarów mamy jeszcze punkt siodło. Punkt jest stabilny, jeśli po obu stronach punktu stałego wektory się skierowane do tego punktu. Takie punkty nazywane są także atraktorami, ściekiem. Punkt stały jest niestabilny, jeśli wektory są zwrócone od punktu - nazywany jest repellerem, źródłem. Punkt jest połowicznie-stabilny, jeśli z jednej strony wektor skierowany jest do punktu, a z drugiej od punktu. Punkt siodło znajduje się na abstrakcyjnym siodle tzn. wzdłuż jednej krzywej jest on stabilny, a w wzdłuż drugiej jest niestabilny. Punkty stałe nazywa się także stanami równowagi. Stany równowagi mogą być stabilne i niestabilne. W stabilnych stanach równwagi układ mimo drobnych odchyleń pozostaje w tym stanie. W przypadku niestabilnych stanów drobna odchyłka od tego punktu prowadzi do zachowań gwałtownych - na przykład takich, jakie są widoczne na rysunku powyżej dla /x/>1.

Układ jest stabilny globalne, jeśli posiada tylko jeden punkt stały i jest on punktem stabilnym. Układ jest stabilny lokalnie, jeśli wśród punktów stabilnych posiada również niestabilne. Nasz układ x2 - 1 jest stabilny lokalnie tzn. dla \(x \leq 1\)

2.1.4 Przykład 3a: model logistyczny

\[\dot N = r \cdot N(1 - \frac{N}{K})\] Jest to model logistyczny liczebności populacji, gdzie:

  • N: liczebność populacji
  • r: współczynnik tempa wzrostu populacji
  • K: graniczna pojemność otoczenia dla danej populacji

Pytanie: dlaczego ten model jest nazywany logistycznym?
Funkcja logistyczna: \[f(x) = \frac{K}{1+e^{-r(x-x_{0})}}\]

  • x0 - wartość zmiennej niezależnej, dla której f(x0) = K/2
  • K - maksymalna wartość funkcji
  • r - wskaźnik wzrostu

Krzywa logistyczna jest rozwiązaniem równania postaci:
\[\dot x = f(x)(1-f(x))\]

dyn_logistic_function.png

Figure 6: Funkcja logistyczna

Rozwiązanie analityczne:

$$\frac{dN}{dt} = rN(1-\frac{N}{K})$$ $$\frac{1}{rN\biggl(1-\frac{N}{K}\biggr)}dN = dt$$ $$\frac{K}{r}\int\frac{1}{N(K-N)}dN = dt$$ $$\frac{1}{N(K-N)} \equiv \frac{A}{N} + \frac{B}{K-N}$$ $$1 \equiv A(K-N) + B(N)$$ $$1 \equiv N(B-A) + AK$$ $$B-A=0, AK=1, A=1/K, B=1/K$$ $$\frac{1}{r}\int\frac{1}{N}dN + \frac{1}{r}\int\frac{1}{K-N}dN = \int dt$$ $$\frac{1}{r}\ln{|N|} - \frac{1}{r}\ln{|K-N|} = t + C$$ $$N(0) = N_{0} \rightarrow C = \frac{1}{r}\ln{\biggl|\frac{N_{0}}{K-N_{0}}\biggr|} - warunek\quad początkowy$$ $$\frac{1}{r}\ln{\biggl|\frac{N}{K-N}\biggr|} - \frac{1}{r}\ln{\biggl|\frac{N_{0}}{K-N_{0}}\biggr|} = t$$ $$\biggl|\frac{N}{K-N}\biggr| \cdot \biggl|\frac{K-N_{0}}{N_{0}}\biggr| = e^{rt}$$ $$dla\; N<0,\; \biggl(-\frac{N}{K-N}\biggr) \cdot \biggl(-\frac{K-N_{0}}{N_{0}}\biggr) = e^{rt}$$ $$dla\; N \in (0,K) ,\; \biggl(\frac{N}{K-N}\biggr) \cdot \biggl(\frac{K-N_{0}}{N_{0}}\biggr) = e^{rt}$$ $$dla\; N>K,\; \biggl(-\frac{N}{K-N}\biggr) \cdot \biggl(-\frac{K-N_{0}}{N_{0}}\biggr) = e^{rt}$$ $$N(t)=\frac{KN_{0}e^{rt}}{N_{0}(e^{rt}-1)+K}$$

Punkty stałe stanowią rozwiązanie równań: \(r \cdot N(1 - \frac{N}{K}) = 0\), \(1 - \frac{N}{K} = 0\). Punkty stałe, więc to: N=0 i N=K

logistic_model_a.png

Figure 7: Rozwiązanie numeryczne równania

par r=1
par K=2
f(x) = r*x*(1 - x/K)
x' = f(x)
s'=1
@ total=10, bounds=10000, xp=S,yp=X,xlo=0,ylo=-4,xhi=10,yhi=4
done
# W GUI:
# Initialconds –> Range —> dobranie odpowiednich parametrów –> OK
# Dir.field/flow —> Direct Field —> dobranie siatki —> ENTER

Na powyższym rysunku widzimy dynamikę modelu logistycznego dla r = 1, K = 2 oraz warunkami początkowymi N0 od -2 do 4. Możemy się przekonać, że nasze obliczenia punktów stałych są poprawne.

Pytanie: dlaczego dla x0=0 wektor nie jest równy 0?

2.1.5 Przykład 3b: model logistyczny - analiza geometryczna

Wizualizacja: https://www.geogebra.org/classic/pjk4pjs5

logistic_model_phase.png

Figure 8: Pole wektorowe na linii fazowej

Punkt N1 = 0 jest punktem niestabilnym, gdyż wektor $Vx=\frac{dx}{dt} jest skierowany od tego punktu; punkt N2 = 2 jest punktem stabilnym, gdyż wektory, z jednej i drugiej strony punktu, są skierowane do tego punktu. To zachowanie układu widzimy na rysunku 7 - trajektorie zmierzają do N2 = 2.

logistic_model_r1.gif

Figure 9: Animacja ruchu dla x0 = 0.26 i r = 1

logistic_model_r2.gif

Figure 10: Animacja ruchu dla x0 = 0.26 i r = 2

2.1.6 Liniowa analiza stabilności

Analizę stabilność możemy przeprowadzić bez tworzenia rysunku. W tym celu potraktujemy zmiany w okolicach punktów stałych jako zmiany liniowe - wykonamy linearyzację wokół punktów stałych. Liniowa analiza stabilności polega na obliczeniu pochodnej funkcji f(x) i sprawdzeniu jej wartości w punktach stałych. Jeśli jest dodatnia to punkt jest niestabilny, jeśli ujemna to punkt jest stabilny. Gdy pochodna jest równa 0, to wtedy konieczna jest interpretacja równania różniczkowego w postaci pola wektorowego. Na przykład mamy takie układy dynamiczne:

Układ dynamiczny:1
$$\dot x = x^2 - k, \quad k=1$$ $$f'(x) = 2x$$ $$f'(1) = 2>0,\text{ punkt niestabilny}$$ $$f'(-1) = -2<0,\text{ punkt stabilny}$$
Układ dynamiczny:2
$$\dot N = rN(1-\frac{N}{K})$$ $$f'(N) = r - \frac{2r}{K}N$$ $$f'(0) = r > 0,\text{ punkt niestabilny}$$ $$f'(K) = r - 2r < 0,\text{ punkt stabilny}$$
Układ dynamiczny:3
$$\dot x = x^2 -k, \quad k=0$$ $$f'(x) = 2x$$ $$f'(0) = 0?$$
Układ dynamiczny:4
$$\dot x = x^3$$ $$f'(x) = 3x^2$$ $$f'(0) = 0?$$

dyn_linear_analysis.png

2.1.7 Potencjał

Zachowanie jakościowe układu możemy jeszcze w inny sposób przedstawić niż tylko za pomocą pola wektorowego. To przedstawienie będzie analogią do przestrzeni i nawiązuje do koncepcji energii potencjalnej. W tym celu opiszemy nasz układ jako potencjał, czyli funkcję zależną od położenia. \(\dot x = f(x) = -\frac{dV}{dx}\). Nasze układy omawiane wcześniej przyjmują postać:

Układ dynamiczny:1
$$\dot x = x^2 - 1$$ $$-\frac{dV_{1}}{dx}=x^2-1$$ $$\int dV_{1} = -\int(x^2-1)dx$$ $$V_{1}=x-\frac{x^{3}}{3}$$ $$V_{1}=x(1-\frac{1}{3}x^{2})$$
Układ dynamiczny:2
$$\dot N = rN(1-\frac{N}{K})$$ $$-\frac{dV_{2}}{dN}=rN(1-\frac{N}{K})$$ $$\int dV_{2}=-\int(rN-\frac{rN^2}{K})dN$$ $$V_{2}=-\frac{rN^2}{2} + \frac{rN^3}{3K}$$ $$V_{2}=rN^2(\frac{N}{3K} - \frac{1}{2})$$

dyn_potential.png

Figure 12: Rozwiązania

2.1.8 Metody numeryczne

Metoda Eulera:
$$x_{n+1} = x_{n} + f(x_{n})\Delta t$$
Ulepszona metoda Eulera:
$$\tilde{x}_{n+1} = x_{n} + f(x_{n})\Delta t$$ $$x_{n+1} = x_{n} + 0.5(f(x_{n})+f(\tilde{x}_{n+1}))\Delta t$$
Metoda Runge–Kutta czwartego rzędu:
$$k_{1} = f(x_{n})\Delta t$$ $$k_{2} = f(x_{n}+0.5k_{1})\Delta t$$ $$k_{3} = f(x_{n}+0.5k_{2})\Delta t$$ $$k_{4} = f(x_{n}+0.5k_{3})\Delta t$$ $$x_{n+1} = x_{n}+\frac{1}{6}(k_{1} + 2k_{2} + 2k_{3} + k_{4})$$

Implementacja algorytmu Eulera dla \(\dot x = x^2-1\):

using Plots
y = zeros(100,21)
for j=1:21 
    y[1,j] = j/10 - 0.1    
    for i=2:100
        y[i,j] = y[i-1,j] + (y[i-1,j]^2-1)*0.1
    end
end
x = 0:0.1:9.9
plot(y,linestyle = :dot)
savefig("dyn_slope_field_euler.png")

Implementacja algorytmu Runge-Kutta dla \(\dot x = x-x^3\):

using Plots
y = zeros(100,10)
delta = 0.1
for j=1:10 
    y[1,j] = j/10 - 0.5    
    for i=2:100
        k1 = (y[i-1,j] - y[i-1,j]^3)*delta
        k2 = (y[i-1,j] + 0.5*k1 - (y[i-1,j] + 0.5*k1)^3)*delta
        k3 = (y[i-1,j] + 0.5*k2 - (y[i-1,j] + 0.5*k2)^3)*delta
        k4 = (y[i-1,j] + 0.5*k3 - (y[i-1,j] + 0.5*k3)^3)*delta
        y[i,j] = y[i-1,j] + 1/6*(k1+2*k2+2*k3+k4)
    end
end
plot(y,linestyle = :dot)
savefig("dyn_slope_field.png")

dyn_slope_field.png

Figure 13: Metoda Runge-Kutta dla \(\dot x = x - x^3\)

runge_kutta.png

Figure 14: Rozwiązanie w XPP dla \(\dot x = x - x^3\)

2.2 Bifurkacje

Najbardziej interesującą cechą, w omawianych do tej pory układach, jest ich zależność od parametrów. Dobierając odpowiednie wartości tych parametrów, dynamika układu jakościowo się zmienia tzn. pojawiają i znikają punkty stałe oraz zmianie może ulec ich stabilność. Te jakościowe zmiany w zależności od parametru nazywane są bifurkacjami a wartości parametrów, przy których zachodzą, nazywają się punktami bifurkacji.

2.2.1 Bifurkacje siodło-węzeł (ang. saddle-node bifurcation)

\[\dot x = x^2 - k, \quad k=1\] \[\dot x = x^2 - k, \quad k=0\]

dyn_bifurcation.png

Figure 15: Punkty bifurkacji

Wyraźnie widzimy, że nasz układ zmienia się w zależności od parametru k. Jeśli będziemy zmieniali parametr k od wartości ujemnych do dodatnich, to dla dla k < 0 - nie ma punktów stałych. Jednak gdy dojdziemy do wartości k = 0, to zachodzi jakościowa zmiana działania układu. Pojawia się punkt stały połowicznie-stabilny - zachodzi bifurkacja siodło-węzeł (inne nazwy: fold bifurcation, turning-point bifurcation, blue sky bifurcation). Gdy przesuwany dalej w kierunku dodatnich wartości k, to widzimy jak zachowanie układu nam się "rozgałęziło" - pojawiły się dwa punkty stałe. Z jednego punktu stałego dostaliśmy dwa - zaszła tutaj bifurkacja. Do graficznego przedstawienie zmian zachowania się układu w zależności od parametru stosujemy diagram bifurkacji:

bifurcation_saddle_node2.png

Figure 16: Diagram bifurkacji

par k = 1
x' = x^2-k

init x=1
@ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
done
# GUI
# Initialconds->Go - obliczamy trajektorie z wartością początkową będącą punktem stałym x=1
# File->Auto - uruchamiamy program Auto
# W programie Auto:
# Axes->hi-lo (ustaw zakres)
# Numerics-> Par Min i Par Max (ustaw zakres)
# Run
# Grab -> za pomocą tabulacji ustaw się na najdalej w lewo wysuniętym punkcie końcowym gałęzi; Enter 
# Numerics -> zmień Ds z 0.02 na -0.02 - zmieniamy kierunek obliczeń
# Run

Na czerwono zaznaczona jest część diagramu, którego układ jest stabilny. Poniżej diagramu możemy to zidentyfikować na podstawie wartości Pt, która jest ujemna. Dla układu niestabilnego Pt jest dodatnie. Wartość Ty określa typ punktu na diagramie. EP jest punkt końcowy gałęzi diagramu (ang. end point of a branch). Program „Auto“ oblicza diagram bifurkacji w przedziałach na gałęziach BP (ang. branch) określonych w sekcji: Numerics –> Nmax - maksymalna liczba punktów obliczeniowych na gałęzi diagramu. Numery oznaczone na diagramie odpowiadają numerom znajdującym się pod Lab. Bifurkacja węzeł-siodło zachodzi w punkcie 5, dla którego Ty jest równe LP (ang. limit point of a branch) - punkt graniczny gałęzi diagramu. To jest nasz szukany punkt bifurkacji siodło-węzeł.

Układ dynamiczny, którego ruch jest określony funkcją postaci \(\dot X = R - X^2\) (prototypem, postać normalna) ma zawsze bifurkacje siodło-węzeł. Na przykład szukając bifurkacji dla takiego układu: \(\dot x = r - x - exp(-x)\). Po rozwinięciu funkcji w szereg Taylora wokół punktu 0 (i dla r=1) dostaniemy \(r - x - exp(-x) = r-1 - \frac{1}{2}x^2 + \frac{1}{6}x^3\dots\). Bifurkacja siodło-węzeł występuje tam, gdzie układ porusza się po paraboli.

dyn_normal_forms.png

T2 jest graficzną reprezentacją rozwinięcia funkcji \(f(x)=r-x-exp(-x)\) w szereg Taylora, aż do elementu kwadratowego; T3 aż do elementu sześciennego. Widać wyraźnie paraboliczny kształt krzywych.

2.2.2 Bifurkacje transkrytyczne (ang. transcritical bifurcation)

Są to bifurkacje, które zachodzą wtedy, gdy punkt stały zmienia się np. ze stabilnego w niestabilny lub odwrotnie, z połowicznie-stabilnego w stabilny lub odwrotnie.

Przykład: \(\dot x = r\ln{x} + x -1\)
Widzimy, że układ posiada punkt stały dla x=1. Czy zachodzi tu bifurkacja transkrytyczna? \[f'(x) = \frac{r}{x} + 1, \quad f'(1) = r + 1\] Z analizy liniowej wynika, że dla r = -1 zachodzi bifurkacja, gdyż pochodna zmienia znak - jest dodatnia dla r > -1; ujemna dla r < -1 a dla r = 0 po zbadaniu portretu fazowego wynika, że jest to punkt połowicznie-stabilny.

dyn_normal_forms3.png

Po rozwinięciu funkcji w szereg Taylora możemy się przekonać, że postać normalna tej funkcji jest \(\dot X = RX - X^2\)

Dla układu \(\dot x = rx - x^2\) diagram bifurkacji jest następujący:

bifurcation_transcritical2.png

Figure 19: Bifurkacja transkrytyczna

par r = 1
x' = r*x - x^2

init x=1
@ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
done
# GUI
# Initialconds->Go - obliczamy trajektorie z wartością początkową będącą punktem stałym x=1
# File->Auto - uruchamiamy program Auto
# W programie Auto:
# Axes->hi-lo (ustaw zakres)
# Numerics-> Par Min i Par Max (ustaw zakres)
# Run
# Grab -> za pomocą tabulacji ustaw się na najdalej w lewo wysuniętym punkcie końcowym gałęzi; Enter 
# Numerics -> zmień Ds z 0.02 na -0.02 - zmieniamy kierunek obliczeń
# Run

Punkt nr 3 jest punktem bifurkacji transkrytycznej. Pod wykresem mamy zidentyfikowany ten punkt jako BP (ang. bifurcation point; branch point) - punkt bifurkacji. Układ dynamiczny z bifurkacją transkrytyczną w przeciwieństwie do bifurkacji siodło-węzeł zawsze posiada punkt stały. Gdy r jest mniejsze od wartości punktu bifurkacji to w układzie istnieją dwa punkty stałe: stabilny i niestabilny. Gdy zwiększymy r na tyle, że stanie się punktem bifurkacji transkrytycznej r = 0 to w układzie pozostaje tylko jeden punkt stały - połowicznie stabilny. Zachowanie układu jest określone przez parabole. Po przekroczeniu tego punktu w układzie znowu mamy dwa punkty stałe: stabilny i niestabilny.

2.2.3 Bifurkacje widłowe (ang. pitchfork bifurcation)

Bifurkacje tego typu są popularne w układach charakteryzujących się symetrią - równanie różniczkowe jest wtedy symetryczne tzn. jest niezmienne przy zamianie zmiennych z x na -x. Są ich dwa typy: nadkrytyczne (ang. supercritical; nazywana także bifurkacją przednią lub w mechanice statystycznej przejściem fazowym drugiego rzędu) i podkrytyczne (nazywana także bifurkacją wsteczną lub w mechanice statystycznej przejściem fazowym pierwszego rzędu; ang. subcritical). Cechą charakterystyczną jest jednoczesne pojawianie się lub znikanie symetrycznie rozstawionych par punktów stałych. W przypadku bifurkacji nadkrytycznych są to punkty stabilne a dla podkrytycznych - niestabilne.

  1. Bifurkacja widłowa nadkrytyczna

    Przykład: \(\dot x = rx - x^3\), \(\dot x = x(r-x^2)\)
    Punkty stałe są dla x = 0 i x = \(\pm \sqrt{r}\) \[f'(x) = r - 3x^2\] \[f'(0) = r\]

    • r > 0, x* jest niestabilny
    • r < 0, x* jest stabilny
    • r = 0, x* jest stabilny (wynika z analizy portretu fazowego; ale jest to słaba stabilność)

    \[f'(\pm \sqrt{r}) = -2r\]

    • możliwe tylko dla r > 0, x* jest stabilny,

    dyn_pitchfork.png

    Gdy r = -1 układ ma jeden punkt stały w \(x_* = 0\) i jest on stabilny. Gdy zwiększymy wartość r do 0 to nadal mamy układ z jednym punktem stałym, ale zachodzi charakterystyczne dla takiego układu z elementem sześciennym -x3 zjawisko spowolnienia krytycznego (ang. critical slowing down). Działanie takiego układu dąży do osiągnięcia punktu stałego \(x_*=0\) najwolniej w porównaniu działaniem z parametrem r < 0 (dla r > 0 układ dąży do innego punktu stałego). Dzieję się tak, gdyż rozwiązaniem układu dla r = 0 jest funkcja:

    \[x(t) = x_{0}\sqrt{\frac{1}{2x_{0}t+1}}\]

    Zmiany określone taką funkcją po przekroczeniu pewnego punktu zachodzą wolniej niż eksponencjalne. Funkcja eksponencjalna z kolei jest rozwiązaniem układu dla \(r \ne 0\)

    \[x =x_{0}\sqrt{\frac{r}{(\frac{r}{x_{0}^2}-1)e^{-2rt} + 1}}\]

    dyn_ex249_critical_slowing_down.png

    Figure 21: Spowolnienie krytyczne

    par r=-0.01
    x'=r*x - x^3
    y'=-y
    z'= -x^3
    init x=5, y=5, z=5
    
    @ total=40, ylo=0,yhi=5
    @ nplot=3,xp2=t,yp2=y,xp3=t,yp3=z
    @ runnow=1
    

    Na rysunku 21 widzimy rozwiązania dla trzech różnych układów. Najszybciej zmieniający się wartość x jest zaznaczona na czerwono i dotyczy układu liniowego \(\dot x=-x\) - rozwiązaniem tego układu jest funkcja eksponencjalna. Wolniej będzie poruszał się punkt fazowy po czarnej krzywej, która jest rozwiązaniem układu \(\dot x = rx-x^3,\;r=0.01\). Pomarańczowa krzywa stanowi rozwiązanie układu \(\dot x=-x^3\)

    dyn_ex249_critical_slowing_down.gif

    Figure 22: Spowolnienie krytyczne

    Gdy r > 0 to nasz punkt x0 staje się niestabilny i pojawiają się dwa symetrycznie rozstawione punkty stałe. Dla r = 1 są to x1 = 1 i x2 = -1. Oba punkty są stabilne.

    gdy \(t \rightarrow \infty\) to \((\frac{r}{x_{0}^2}-1)e^{-2rt}\rightarrow0\)
    a więc \(x = \pm\sqrt{r}\)

    Diagram bifurkacji:

    bifurcation_pitchfork.png

    Figure 23: Bifurkacja widłowa nadkrytyczna

    par r = -1
    x' = r*x - x^3
    
    init x=0
    @ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
    done
    # GUI
    # Initialconds->Go - obliczamy trajektorie z wartością początkową będącą punktem stałym x=0
    # File->Auto - uruchamiamy program Auto
    # W programie Auto:
    # Axes->hi-lo (ustaw zakres)
    # Numerics-> Par Min i Par Max (ustaw zakres)
    # Run
    # Grab -> za pomocą tabulacji ustaw się na najdalej w prawo wysuniętym punkcie końcowym gałęzi; Enter 
    # Numerics -> zmień Ds z 0.02 na -0.02 - zmieniamy kierunek obliczeń
    # Run
    

    Punkt nr 2 jest punktem bifurkacji widłowej nadkrytycznej. Na diagramie bifurkacji pojawiają się charakterystyczne widły.

  2. Bifurkacja widłowa podkrytyczna

    Wystąpi w układzie, w którym część sześcienna już nie jest czynnikiem stabilizującym:
    Przykład: \(\dot x = rx + x^3\), \(\dot x = x(r+x^2)\)
    Punkty stałe są dla x = 0 i \(x = \pm \sqrt{-r}\)
    \[f'(x) = r + 3x^2\] \[f'(0) = r\]

    • r > 0, x* jest niestabilny
    • r < 0, x* jest stabilny
    • r = 0, x* jest stabilny (wynika z analizy portretu fazowego)

    \[f'(\pm \sqrt{-r}) = -2r\]

    • możliwe tylko dla r < 0, x1 i x2 są niestabilne

    dyn_pitchfork2.png

    Diagram bifurkacji:

    bifurcation_pitchfork2.png

    Figure 25: Bifurkacja widłowa podkrytyczna

    par r = -1
    x' = r*x + x^3
    
    init x=0
    @ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
    done
    # GUI
    # Initialconds->Go - obliczamy trajektorie z wartością początkową będącą punktem stałym x=0
    # File->Auto - uruchamiamy program Auto
    # W programie Auto:
    # Axes->hi-lo (ustaw zakres)
    # Numerics-> Par Min i Par Max (ustaw zakres)
    # Run
    # Grab -> za pomocą tabulacji ustaw się na najdalej w prawo wysuniętym punkcie końcowym gałęzi; Enter 
    # Numerics -> zmień Ds z 0.02 na -0.02 - zmieniamy kierunek obliczeń
    # Run
    

    Dla przypadku niestabilnego r=1 i x0 = 2

    dyn_pitchfork5.png

    Figure 26: Zjawisko wybuchu

    Charakterystyczne dla tego układu jest zjawisko wybuchu (ang. blow-up), czyli układ osiąga wartość nieskończoną w skończonym czasie - to widać na rysunku jako wyraźny skok przebiegu funkcji. Rozwiązaniem układu jest:

    \[x =x_{0}\sqrt{\frac{r}{(\frac{r}{x_{0}^2}+1)e^{-2rt} - 1}}\]
    \[t =-\frac{1}{2r}\log{|\frac{r}{x^2}+1|}+\frac{1}{2r}\log{|\frac{r}{x_{0}^2}+1|}\]
    gdy \(x \rightarrow \infty\) to \(-\frac{1}{2r}\log{|\frac{r}{x^2}+1|}\rightarrow0\)
    a więc dla \(t = \frac{1}{2r}\log{|\frac{r}{x_{0}^2}+1|}\) układ osiąga nieskończoność.

    Aby zabezpieczyć układ przed wybuchem możemy dodać do naszego układu element stabilizujący np. \[\dot x = rx + x^3 - x^5\]

    Dla tego równania diagram bifurkacji wygląda następująco:

    bifurcation_pitchfork3.png

    Figure 27: Bifurkacja widłowa podkrytyczna

    par r = -1
    x' = r*x + x^3 - x^5
    
    init x=0
    @ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
    done
    # GUI - podobnie jak w poprzednich przykładach
    

    Dzięki dodatkowemu elementowi na diagramie bifurkacji pojawiły się dodatkowe stabilne duże amplitudy odpowiadające punktom stabilnym i zabezpieczające nas przed wybuchem. Istnienie różnych stanów stabilnych stwarza możliwość pojawiania się skoków i histerezy. Na przykład działanie układu gdy r < 0 i x* = 0 jest stabilne. Pozostaje takie dopóki r < 0. Gdy zwiększymy r powyżej 0, punkt x* staje się niestabilny i następuje skok do jednej z dwóch dużych amplitud. Układ pozostaje w stanie stabilnym, gdy nadal będziemy zwiększać r. Natomiast gdy zaczniemy zmniejszać ten parametr, to po przekroczeniu wartości 0 nadal będziemy w stanie stabilnym, ale innym niż z tego co wyszliśmy na początku. Skok do początkowego stanu stabilnego nastąpi, gdy r < rc. Taki brak odwracalności działania układu podczas zmiany wartości parametru nazywa się histerezą.

    W punktach 4 i 6 zachodzi bifurkacja siodło-węzeł a w punkcie 2 bifurkacja widłowa. Gdy w układzie zmieniamy r od wartości ujemnych do dodatnich, to zachodzi bifurkacja widłowa nadkrytyczna a gdy odwrotnie to podkrytyczna.

2.2.4 Bifurkacje niedoskonałe i katastrofy

Jeśli do postaci normalnej układu z bifurkacją widłową dodamy parametr zaburzający symetrię np: \(\dot x = h + rx - x^3\), to h nazywamy parametrem niedoskonałości.

Pole wektorowe jest określone przez równanie algebraiczne trzeciego stopnia. Zatem maksymalnie może mieć trzy rozwiązania w stanach równowagi.

dyn_catastrophe.png

Figure 28: Portrety fazowe

Rysunki przedstawiają portrety fazowe z dodatkowymi horyzontalnymi liniami odpowiadającym wartościom parametrowi h: zielona linia h = 0.3849, czerwona h = 0. Dodatkowy parametr powoduje, że analiza dynamiki naszego układu staje się bardziej złożona. Gdy \(r \le 0\) to w układzie jest jeden punkt stały i zmiana parametru h wpływa jedynie na miejsce pojawienia się tego punktu. Ta część jest jakościowa taka sama jak dla układu z bifurkacją widłową nadkrytyczną.

Gdy r > 0 i h = 0 to układ redukuje się do układu z bifurkacją widłową, która charakteryzuje się symetrycznie rozłożonymi punktami stałymi. Gdy \(h \ne 0\) symetria zostaje zniszczona. Na rysunku po lewej stronie na dole linia h tworzy dwa punkty stałe: przesuwa x* = -1 w lewą stronę oraz łączy dwa pozostałe punkty w jeden w miejscu gdzie zielona linia jest styczna do krzywej. Dla układu po prawej stronie parametr h przesunął tylko punkty w inne miejsca. To przesuwanie punktów niszczy symetrię. Gdybyśmy zwiększyli parametr h o dostatecznie dużą wartość to ten tylko po lewej prawej pozostałby, dwa pozostałe najpierw złączyłby się w jeden, tak jak to jest pokazane na rysunku po lewej stronie, a później by zniknął i ten. Ten punkt, gdzie zachodzi łączenie i anihilacja punktów stałych zachodzi bifurkacji siodło-węzeł. Bifurkacje powstające poprzez modyfikowanie w równaniu układu dwóch parametrów są nazywane bifurkacjami z kowymiarem-2 (ang. codimension-2 bifurcation). Do tej pory zajmowaliśmy się bifurkacjami o kowymiarze 1, gdyż modyfikowaliśmy tylko jeden parametr.

dyn_catastrophe2.png

Figure 29: Portrety fazowe

Jeszcze jedno przedstawienie portretów fazowych z uwzględnieniem wpływu parametru h na dynamikę układu.

Pytanie: czy punkty stałe dla fioletowej linii są rozłożone symetrycznie względem osi rzędnych?

cusp_bifurcation_h.png

Figure 30: Diagram bifurkacji dla parametru h prz stałym parametrze r = 1

par r = 1
par h = 0.1
x' = h + r*x - x^3

init x=1
@ total=10, xlo=0,ylo=-2,xhi=5,yhi=2
done
# GUI
# Obliczamy trajektorie z wartością początkową będącą blisko punktu stałego x=1
# Initialconds->Go –> Initialconds->Last
# File->Auto - uruchamiamy program Auto
# W programie Auto:
# Axes->hi-lo (ustaw zakres)
# Numerics-> Par Min i Par Max (ustaw zakres)
# Run
# Grab -> za pomocą tabulacji ustaw się na najdalej w lewo wysuniętym punkcie końcowym gałęzi; Enter 
# Numerics -> zmień Ds z 0.02 na -0.02 - zmieniamy kierunek obliczeń
# Run

Punkty 3 i 4 są punktami bifurkacji siodło-węzeł. Gdy przeprowadzimy linie pionową zaczynając od wartości ujemnych h w kierunku dodatnich to do wartości \(h<-0.3849\) mamy jeden punkt stały, który jest stabilny. Dla \(h=-0.3849\) mamy dwa punkty stałe. Dla \(h>-0.3849\) i \(h<0.3849\) mamy trzy punkty stałe (ta część jest podobna jakościowa do układy z bifurkacją widłową nadkrytyczną). Dla \(h=0.3848\) mamy dwa punkty i dla \(h>0.3849\) mamy znowu jeden punkt stały.

Równania parametrów h i r: \[h + rx - x^3 = 0\] \[h=x^3-rx\] \[r= x^2 - \frac{h}{x}\]

dyn_catastrophe3.png

Figure 31: Diagramy bifurkacji

Na powyższych rysunkach widzimy zmiany parametrów dla kilku przypadków. To są diagramy bifurkacji tylko odwrócone. Na postawie rysunku po prawej stronie na dole, obliczymy diagram bifurkacji dla r w XPP.

cusp_bif_r.png

Figure 32: Diagram bifurkacji

par r = -1
par h = 0.3849
x' = h + r*x - x^3

init x=0.5

cusp_bif_r2.png

Figure 33: Diagram bifurkacji

par r = 2
par h = 0.3849
x' = h + r*x - x^3

init x=-1.5

cusp_bifurcation_hr.png

Figure 34: Diagram stabilności

Zaczynając analizę rysunku od najdalej wysuniętych na lewo wartości h mamy zbiór punktów (h,/r/) dla których w układzie istnieje tylko jeden punkt stały. Przesuwając się w stronę dodatnich wartości h wchodzimy na czerwoną krzywą, która informuje nas, że dla tych punktów (h,/r/) mamy dwa punkty stałe. Pomiędzy krzywymi mamy trzy punkty stałe. Punkt (h/=0,/r/=0) jest punktem wierzchołkowym (ang. /cusp point).

Termin katastrofa ma swoje uzasadnienie w takim sensie, że zmiana parametru r poniżej punktu wierzchołkowego powoduje przeskok do innego stanu stałego. Ten przeskok może być tak duży, że w swoich skutkach katastrofalny.

2.2.5 Przykład: Bifurkacja widłowa nadkrytyczna

Przykład: \(\dot x = -x + \beta \tanh{x}\) - to równianie jest symetryczne i często pojawia się sieciach neuronowych
Dla \(\beta < 1\) mamy tylko jednen punkt stały x*=0 Dla \(\beta \ge 1\) będzie więcej niż jeden punkt stały
\(-x + \beta \tanh{x} = 0\) - jest to równanie z funkcją przestępną i nie jest one trywialne w rozwiązaniu. Możemy jednak potraktować punkty stałe x* jako zmienną niezależną i obliczyć je z równania \(\beta = \frac{x_*}{\tanh{x_*}}\)

W rezultacie naszych obliczeń dostaniemy diagram bifurkacji. Bifurkacja widłowa nadkrytyczna zachodzi dla \(\beta = 1\)

dyn_ex341_bif.png

Figure 35: Diagram bifurkacji

2.3 Dynamika na okręgu

Pole wektorowe na okręgu jest najprostszym układem, który oscyluje. To pole wektorowe jest funkcją, która jednoznacznie określa wektor prędkości w każdym punkcie na okręgu. Każdy punkt na okręgu jest kątem lub fazą. Faza \(\theta = 0\) i \(\theta = 2\pi\) jest tym samym, gdyż dla funkcji okresowej np. \(\sin\) mamy \(f(\theta + k2\pi) = f(\theta)\)

Przykład: \(\dot \theta = sin(\theta)\)

example411.png

Figure 36: Potok na okręgu

theta' = sin(theta)
s'=1

@ total=10, bounds=10000, xp=S,yp=theta,xlo=0,ylo=0,xhi=10,yhi=6.26
# GUI
# Dir.field/flow —> Direct Field —> dobranie siatki —> ENTER
# Dir.field/flow —> Flow —> dobranie siatki —> ENTER

Punkty stałe:
\(\theta_{1} = 0\) - niestabilny \(\theta_{2} = \pi\) - stabilny

2.3.1 Jednorodny oscylator

\[\dot \theta = \omega\]

Jednorodność wynika z faktu, że prędkość kątowa jest stała. Dwa nie oddziałujące na siebie jednorodne oscylatory działające z różną prędkością kątową będą okresowo wchodzić w tę samą fazę i to zjawisko nazywa się ang. beat phenomenon.

Przykład (Zadanie 4.2.1 - Strogatz, 2014): dzwony w dwóch różnych kościołach dzwonią z różną częstotliwością: pierwszy co 3 sekundy, drugi co 4 sekundy. Zakładamy, że zaczęły dzwonić w tym samym czasie. Po jakim czasie dzwony wybrzmią razem?

Wybrzmią razem po upływie 12 sekund, gdyż jest to pierwsza liczba, której wspólnym dzielnikiem są liczby 3 i 4.

\(\dot \theta_{1} = \omega_{1}\), \(\omega_{1}=\frac{2\pi}{T_{1}}\) \(\dot \theta_{2} = \omega_{2}\), \(\omega_{2}=\frac{2\pi}{T_{2}}\)

Rozwiązaniem naszych równań ruchu są: \(\theta_{1} = \omega_{1}t\) \(\theta_{2} = \omega_{2}t\)

Dzwony wybrzmią razem jeśli: \(\omega_{1}t - \omega_{2}t = 2\pi\)

Rozwiązanie tego równania względem t: \(t = \frac{2\pi}{\omega_{1} - \omega_{2}}\) \(t = \frac{T_{1}T_{2}}{T_{2} - T_{1}}\) \(t=(\frac{1}{T_{1}} - \frac{1}{T_{2}})^{-1}\) \(t = 12\)

Ponieważ okres dzwonienia T1 i T2 jest liczbą całkowitą to czas t kiedy oba dzwony wybrzmiewają razem też musi być liczbą całkowitą. W obliczeniach poniżej dostaniemy w wyniku liczbę wymierną, która jest czasem, kiedy jeden dzwon zdubluje w dzwonieniu drugi. Wtedy wielokrotność tej liczby jest czasem, kiedy dzwony wybrzmią razem np.

  • T1 = 3 i T2 = 5 \(t = (\frac{1}{3} - \frac{1}{5})^{-1}=\frac{15}{2}\). Oznacza, że pierwszy dzwoń zdubluje drugiego w t = 7.5 s (\(\omega_{1}t = \omega_{2}t\)) a dla t = 15 s dzwony wybrzmią razem.
  • T1 = 3 i T2 = 6 \(t = (\frac{1}{3} - \frac{1}{6})^{-1}=6\)

ex421.png

Figure 37: Wizualizacja rozwiązania w XPP

x'=2*pi/3
y'=2*pi/4

@ total=12, bounds=10000, xp=X,yp=Y,xlo=0,ylo=0,xhi=6.28,yhi=6.28
done
GUI
# phAsespace —> All
# Initialconds –> Go

Na rysunku widzimy zachowanie się układu dla warunków początkowych \(X=\theta_{1} = 0\) i \(Y=\theta_{2} = 0\). Po upływie t = 12 s układ powraca do punktu początkowego, czyli \(X=\theta_{1} = 2\pi\) i \(Y=\theta_{2} = 2\pi\). Czas potrzeby do tego, aby układ powrócił do tego punktu można wyliczyć jako iloczyn przecięć trajektorii układu z osią X i Y. Przecięcia z osią X oznaczają, że 3 razy oscylator \(\theta_{1}\) powrócił do stanu początkowego a przecięcia z osią Y oznaczają, że 4 razy oscylator \(\theta_{2}\) powrócił do stanu początkowego.

Liczba wymierna:

  • T1 = 1/2 i T2 = 1/4; \(t = (4-2)^{-1}=\frac{1}{2}\) - wybrzmią razem w t = 0.5 s
  • T1 = 1/2 i T2 = 1/5; \(t = (5-2)^{-1}=\frac{1}{3}\) - wybrzmią razem w t = 1 s, gdy T2 dwa razy zdubluje T1

ex421_a.png

Figure 38: Oscylator z prędkościami kątowymi o wartościach wymiernych

x'=2*2*pi
y'=5*2*pi

@ total=1, bounds=10000, xp=X,yp=Y,xlo=0,ylo=0,xhi=6.28,yhi=6.28
done
GUI
# phAsaspace —> All
# Initialconds –> Go

ex421_b2.png

Figure 39: Oscylator drugi ma prędkość kątową o wartości niewymiernej

x'=2*pi
y'=sqrt(2*pi)

@ total=200, bounds=10000, xp=X,yp=Y,xlo=0,ylo=0,xhi=6.28,yhi=6.28
done
GUI
# phAsaspace —> All
# Initialconds –> Go

Pytanie: czy zachowanie układu widoczne na rysunku 39 jest chaotyczne?

2.3.2 Niejednorodny oscylator

\[\dot \theta = \omega - a\sin{\theta}\]

Zastosowanie ma np w modelowaniu oscylujących neuronów, zachowania świetlików.

ex432.png

Figure 40: Portret fazowy

Parametr a wpływa na jakościowe zachowanie się układu. Na powyższym rysunku 40 widzimy zachowanie układu z prędkością kątową \(\omega=2\) i trzema różnymi wartościami parametru a. Gdy \(a < \omega\) układ oscyluje. W miarę jak wartość parametru a będzie się zbliżała do ω to prędkość cały czas maleje. W punkcie \(\theta_{1}=-\frac{\pi}{2}\) jest maksymalna a w punkcie \(\theta_{2}=\frac{\pi}{2}\) jest minimalna. Punkt \(\theta_{2}\) jest punktem bifurkacji siodło-węzeł. Gdy \(a=\omega\) to pojawia się punkt połowicznie stabilny - układ z lewej strony będzie dążył do tego punktu a z prawej odchodził. Gdy \(a>\omega\) to pojawiają się w układzie dwa punktu stałe - stabilny i niestabilny.

ex432_a.png

Figure 41: Cztery trajektorie układu dla x0=0,1,2,3 i a=1.99

ex432.mp4; XPP: x'=2-a*sin(x); @total=100

Stabilność punktów stałych:
\(\omega-a \cdot sin(\theta)=0\)

  1. \(\omega = a,\, sin(\theta)=1\)
    • punt stały: \(\theta=\pi/2\)
    • stabilność: \(f'= a \cdot cos(\theta)\), \(f'(\pi/2)=a \cdot cos(\pi/2)= 0\)
  2. \(\omega < a\), \(sin(\theta)=\frac{\omega}{a}\), \(\theta=\arcsin(\frac{\omega}{a})\), - dwa punkty stałe;

\(sin(\theta)^2=\frac{\omega^2}{a^2},\, 1-cos(\theta)^2=\frac{\omega^2}{a^2},\,cos(\theta)=\pm\sqrt{1-\frac{\omega^2}{a^2}}\),
Stabilność: \(f' = \mp a \sqrt{1-\frac{\omega^2}{a^2}}\) - dwa punkty: stabilny i niestabilny

Pytanie: jaki jest okres oscylacji?

\(\int\limits_{-\pi}^{\pi} \frac{1}{\omega-a \cdot sin(\theta)}\;\mathrm{d}\theta=\int\mathrm{d}t\)

W celu obliczenia tej całki wykonany podstawienie tangensa kąta połowicznego. \(tan(\frac{\theta}{2}) = u\), \(\theta = 2\arctan(u)\), \(\mathrm{d}\theta = \frac{2}{1+u^2}\mathrm{d}u\) Funkcję sinus możemy wyrazić jako iloraz funkcji u. Wykorzystując formułę podwójnego kąta mamy:
\[\sin(\theta)=2\sin\Biggl(\frac{\theta}{2}\Biggr)\cos\Biggl(\frac{\theta}{2}\Biggr)\] \[\sin(\theta)=2\frac{\sin\Bigl(\frac{\theta}{2}\Bigr)}{\cos\Bigl(\frac{\theta}{2}\Bigr)}\cos^2\Biggl(\frac{\theta}{2}\Biggr)\] \[\sin(\theta)=2\tan\Biggl(\frac{\theta}{2}\Biggr)\cos^2\Biggl(\frac{\theta}{2}\Biggr)\] \[\sin(\theta)=\frac{2u}{\sec^2(\frac{\theta}{2})}=\frac{2u}{1+u^2}\] Ponieważ zrobiliśmy podstawienie \(tan(\frac{\theta}{2}) = u\) to musimy zmienić odpowiednio zakres całkowania zgodnie ze wzorem:
\(\int\limits_a^b f(g(x))g'(x)\;\mathrm{d}x=\int\limits_{g(a)}^{g(b)} f(u)\;\mathrm{d}u\)

\(\tan(-\frac{\pi}{2})=-\infty\), \(\tan(\frac{\pi}{2})=\infty\)

\(\int\limits_{-\infty}^{\infty} \frac{\frac{u}{1+u^2}}{\omega-\frac{2ua}{1+u^2}}\;\mathrm{d}u\) \(\int\limits_{-\infty}^{\infty} \frac{2}{\omega u^2 -2ua+\omega}\;\mathrm{d}u\)

Pamiętamy, że liczymy okres oscylacji, czyli dla \(\omega>a\). Ponieważ wyróżnik trójmianu kwadratowego jest mniejszy od zera to rozwiązywanie takiej całki polega to sprowadzaniu równania w mianowniku do postaci kanonicznej.
\[ax^2+bx+c = a\Biggl[\Biggl(x+\frac{b}{2a}\Biggr)^2 + \frac{-\Delta}{4a^2}\Biggr]\] \[\omega u^2 -2ua+\omega = \omega\Biggl[\Biggl(u-\frac{2a}{2\omega}\Biggr)^2+\frac{-\Delta}{4\omega^2}\Biggr]\]

\[\Delta=4a^2-4\omega^2<0\]

Rozwiązywanie całki funkcji wymiernej, czyli ilorazu wielomianów, polega na rozkładaniu funkcji podcałkowej na ułamki proste jeśli wyróżnik trójmianu jest równy bądź większy od zera. Jeśli jest mniejszy od zera to sprowadzamy do postaci kanonicznej i rozwiązaniem będzie funkcja \(arctan\). Gdy licznik jest pochodną mianownika to rozwiązaniem jest funkcja logarytmiczna.

\[\frac{2}{\omega}\int\limits_{-\infty}^{\infty}\frac{1}{(u-\frac{a}{\omega})^2+\frac{\omega^2-a^2}{\omega^2}}\;\mathrm{d}u\] \[u-\frac{a}{\omega} = \sqrt{\frac{\omega^2-a^2}{\omega^2}}k,\; \mathrm{d}u=\sqrt{\frac{\omega^2-a^2}{\omega^2}}\mathrm{d}k\] \[\frac{2}{\omega}\int\limits_{-\infty}^{\infty} \frac{\sqrt{\frac{\omega^2-a^2}{\omega^2}}}{\frac{\omega^2-a^2}{\omega^2}k^2+\frac{\omega^2-a^2}{\omega^2}}\;\mathrm{d}k\]

\[\frac{2\sqrt{\omega^2-a^2}}{\omega^2-a^2}\int\limits_{-\infty}^{\infty} \frac{1}{k^2+1}\] \[2\sqrt{\frac{1}{\omega^2-a^2}}\int\limits_{-\infty}^{\infty} \frac{1}{k^2+1}\]

\[t=2\sqrt{\frac{1}{\omega^2-a^2}}\arctan\left(k\right)\bigg|_{-\infty}^{\infty}\],

\[k=\frac{uw-a}{\sqrt{\omega^2-a^2}}\]

\[t=2\sqrt{\frac{1}{\omega^2-a^2}}\arctan\left(\frac{uw-a}{\sqrt{\omega^2-a^2}}\right)\bigg|_{-\infty}^{\infty}\],

gdy \(u \rightarrow \infty\) \(T_{1}=\frac{\pi}{\sqrt{\omega^2-a^2}}\),
gdy \(u \rightarrow -\infty\) \(T_{2}=-\frac{\pi}{\sqrt{\omega^2-a^2}}\),
Szukany okres oscylacji wynosi:
\(T=T_{1}-T_{2} = \frac{2\pi}{\sqrt{\omega^2-a^2}}\)

ex432_t.png

Figure 42: Wpływ parametru a na okres oscylacji T dla ω = 2

Pytanie: jak szybko zachodzą zmiany w pobliżu punktu bifurkacji?

\(\frac{2\pi}{\sqrt{\omega^2-a^2}}=\frac{2\pi}{\sqrt{(\omega+a)(\omega-a)}}\),

Gdy \(a\rightarrow \omega^-\) - w pobliżu punktu bifurkacji parametr a jest o bardzo niewielką część mniejszy od ω, więc \(\omega \approx a\),

\(\sqrt{(\omega+a)(\omega-a)}\approx \sqrt{2\omega}\sqrt{(\omega-a)}\)

Drugi człon wyrażenia podpierwiastkowego pozostał niezmieniony, gdyż niewielka różnica pomiędzy ω i a ma bardzo duży wpływ na wartość całego wyrażenia.

To czas potrzebny na przejście przez zwężenie szyjkowe (ang. bottleneck) wynosi:
\(T = \sqrt{\frac{2}{\omega}}\frac{\pi}{\sqrt{\omega-a}}\)

Równanie to jest prawem skalowania pierwiastkiem kwadratowym (ang. square-root scaling law). Jest ono powszechne w układach przechodzących bifurkacje siodło-węzeł. Analizując układ dynamiczny możemy skupić się na kluczowym miejscu jego dynamiki. Dla naszego układu oscylatora niejednorodnego takim miejscem jest \(\theta=\pi/2\). W tym miejscu równanie ruchu jest parabolą - taką dynamiką charakteryzują się właśnie układy z bifurkacją siodło-węzeł. Układ zbliżając się do punktu bifurkacji od prawej strony to dwa punkty stabilne łączą się w jeden, a następnie znika i ten punkt stały. Zaraz po anihilacji punktu stałego mamy zwężenie szyjkowe, które wygląda tak, jakby istniał nadal punkt stały i to zjawisko nazywam się duchem.

ex432_t2.png

Figure 43: Portret fazowy układu bliksko punktu bifurkacji

Odległość pomiędzy parabolą a osią X jest równa r. Postać normalna równania ruchu z bifurkacją siodło-węzeł to: \(\dot x = r - x^2\) Gdy układ jest bardzo blisko punkt bifurkacji to dominuje efekt zwężenia szyjowego, który możemy oszacować na podstawie prawa skalowania pierwiastka kwadratowego:
\(T= \int\limits_{-\infty}^{\infty} \frac{1}{r-x^2}\;\mathrm{d}x=\frac{\pi}{\sqrt{r}}\)

Aby przekonać się, że nasz oscylator niejednorodny można opisać za pomocą równania ruchu o postaci normalnej pokazanej powyżej, rozwiniemy go w szereg Taylora:
\(\dot \theta = \omega - a\cdot sin(\theta)\), dla \(\theta=\pi/2\),
\(a\cdot\sin(\theta) = a-\frac{a}{2}(\theta-\pi/2)^2+...+O(3)\),
\(\dot \theta = \omega - a + \frac{a}{2}(\theta^2-\theta\pi + (\pi/2)^2)\),
\(\phi=\theta-\pi/2\),
\(\dot \phi = \omega - a + \frac{a}{2}(\phi)^2\)

ex432_taylor.png

Figure 44: Portret fazowy dla ω = a = 1

Postać normalną uzyskamy po podstawieniach: \(\omega - a = r\), \(x = \sqrt{\frac{a}{2}}\phi\),
\(\dot x=\sqrt{\frac{a}{2}}\dot\phi\), \(\sqrt{\frac{2}{a}}\dot x = r + x^2\)
\(T=\sqrt{\frac{2}{a}}\int\limits_{-\infty}^{\infty}\frac{1}{r+x^2}\;\mathrm{d}{x}\), \(T=\sqrt{\frac{2}{a}}\frac{1}{\sqrt{r}}\arctan{\frac{x}{\sqrt{r}}}\Bigr|_{-\infty}^{\infty}\), \(T=\sqrt{\frac{2}{a}} (\frac{\pi}{2\sqrt{r}} + \frac{\pi}{2\sqrt{r}}) = \sqrt{\frac{2}{a}}\frac{\pi}{\sqrt{r}}\), \(T=\sqrt{\frac{2}{a}}\frac{\pi}{\sqrt{\omega-a}}\), \(a\rightarrow \omega^-\), \(T=\sqrt{\frac{2}{\omega}}\frac{\pi}{\sqrt{\omega-a}}\)
To równanie jest takie samo jakie uzyskaliśmy wcześniej.

square_rot_law_circle.png

Figure 45: Trajektoria oscylatora niejednorodnego \(\dot \theta=1-0.99sin(\theta)\)

square_rot_law_circle.gif

Figure 46: Animacja ruchu oscylatora niejednorodnego \(\dot \theta=1-0.99sin(\theta)\); czarna linia pozioma jest okręgiem 0-2π

Układy z bifurkacją siodło-węzeł poruszające się po linii fazowej np. \(\dot x = x^2-k\) przechodzą bifurkacje dla k = 0. Na podstawie prawa skalowania pierwiastkowego możemy oszacować różnice czasu pomiędzy zachowaniem się dwóch układów o różnych k w pobliżu punktu bifurkacji np. k1 = -0.1 i k2 = -0.01

Pytanie: ile wyniesie różnica?

\(T= \int\limits_{0}^{\infty} \frac{1}{r-x^2}\;\mathrm{d}x\), \(T=\frac{1}{\sqrt{r}}\arctan{\frac{x}{\sqrt{r}}}\Bigr|_{0}^{\infty}\), \(T=\frac{\pi}{2\sqrt{r}}\),

\(T_{diff}\approx \frac{\pi}{2}(\frac{1}{\sqrt{k_{2}}}-\frac{1}{\sqrt{k_{1}}})\approx 10.74\)

square_rot_law_phase_line.png

Figure 47: Cztery trajektorie układu dla x0=0,1,2,3 i a=1.99

Trajektorie układu dla różnych wartości parametru k: pierwsza od lewej k = -0.1 a ostatnia k = -0.01 T1 = 5, T2 = 16, Tdiff = 11

Trajektorie układu dla różnych wartości parametru k: pierwsza od lewej k = -0.001 a ostatnia k = -0.0001 T1 = 50, T2 = 157, Tdiff = 107

\(T_{diff}\approx \frac{\pi}{2}(\frac{1}{\sqrt{0.0001}}-\frac{1}{\sqrt{0.001}})\approx 107.4\)

square_rot_law_phase_line.gif

Figure 48: Animacja ruchu dla k=-0.001

3 Dwuwymiarowe układy liniowe

3.1 Dwuwymiarowe układy liniowe

3.1.1 Układ równań różniczkowych:

\(\dot x_{1} = f_{1}(t,x_{1},x_{2},...,x_{n})\)
\(\dot x_{2} = f_{2}(t,x_{1},x_{2},...,x_{n})\)
\(\ldots\)
\(\dot x_{n} = f_{n}(t,x_{1},x_{2},...,x_{n})\)

fi jest funkcją rzeczywistą - funkcją o wartościach rzeczywistych (ang. real-valued function). Zakładamy, że jej pochodne cząstkowe każdego rzędu n istnieją i są ciągłe - zapisujemy to w taki sposób: fi są \(C^{n}\)

Postać macierzowa: \(\dot X = F(t,X)\)

\(F(t,X) = \begin{pmatrix}{f_{1}(t,x_{1},x_{2},...,x_{n}) \\f_{1}(t,x_{1},x_{2},...,x_{n}) \\ ... \\f_{n}(t,x_{1},x_{2},...,x_{n})}\end{pmatrix}\)

F(t,X)- jest funkcją wektorową - funkcją o wartościach wektorowych (ang. vector-valued function)

Dla układu dwuwymiarowego autonomicznego funkcja wektorowa jest postaci:
\(F(x_{1},x_{2}) = (f(x_{1},x_{2}),g(x_{1},x_{2}))\) - definiuje ona nam pole wektorowe.

Rozwiązaniem równania jest funkcja postaci \(X(t) = (x_{1}(t),x_{2}(t),...,x_{n}(t))\), która spełnia równanie: \(\dot X(t) = F(t,X(t))\)

\(\dot X(t) = (x'_{1}(t),x'_{2}(t),...,x'_{n}(t))\)

Jeśli równania ruchu nie zależą bezpośrednio od czasu to taki układ nazywamy autonomicznym: \(\dot X = F(X)\)

Wektor X0, dla którego zachodzi \(F(X_{0})=0\) nazywamy punktem równowagi (ang. equalibrium) lub inaczej punktem stałym (ang. fixed point). Punkt ten odpowiada rozwiązaniu stałemu \(X(t) = X_{0}\)

Równanie różniczkowe drugiego rzędu postaci: \(\ddot x = f(t,x,\dot x)\)

np. równanie oscylatora harmonicznego: \(m\ddot x + b\dot x + kx = f(t)\)
\(m\ddot x + b\dot x + kx = 0\)

można sprowadzić do postaci układu dwuwymiarowego przez podstawienie:
\(\dot x = y\)
\(\dot y = -\frac{b}{m}y - \frac{k}{m}x\)

3.1.2 Dwuwymiarowy układ liniowych równań różniczkowych

\(\dot x = f(x,y)\)
\(\dot y = g(x,y)\)

Funkcja wektorowa \(F(x,y) = \dbinom{f(x,y)}{g(x,y)}\) określa regułę pola wektorowego. Dla wartości x(t), y(t) mamy wektor \(F(x(t),y(t))\).

Wizualizacja pola wektorowego może być nieczytelna, gdyż wektory mogą nakładać się siebie, dlatego skaluje się wektory tak aby wizualnie obraz był bardziej czytelny.

vector_field2d.png

Figure 49: Pole wektorowe 2D

Każdemu punktowi na rysunku x(t) i y(t) można przyporządkować wektor \(\dbinom{f(x,y)}{g(x,y)}\): na rysunku mamy pole wektorowe: \((sf(x,y),sg(x,y))\), gdzie s jest parametrem skalującym. W XPP dla punktu \((x,y)\) dostajemy wektor: \((x + sf(x, y), y + sg(x, y))\)

Układ 2D:
\(\dot x = y\)
\(\dot y = -x\)

Możemy zapisać jako równanie drugiego rzędu:
\(\ddot x + x = 0\)

Jest to równanie oscylatora harmonicznego bez tłumienia: To równanie jest typu: \(a_{0}(t)x + a_{1}(t)x'+a_{2}(t)x'' + ... + a_{n}(t)x^{n} + b(t) = 0\)

Szczególnym przypadkiem są równania o stałych współczynnikach: \(a_{0}x + a_{1}x'+a_{2}x'' + ... + a_{n}x^{n} + b = 0\)

Możemy jeszcze bardziej zawęzić równania do jednorodnych: \(a_{0}x + a_{1}x'+a_{2}x'' + ... + a_{n}x^{n} = 0\)

Równanie oscylatora harmonicznego można zapisać jako:
\(m\ddot x + b\dot x + kx = 0\)
b - jest współczynnikiem tłumienia
k - jest współczynnikiem sprężystości

Postać macierzowa:
\(\dot X = AX\)

\(A=\begin{pmatrix}a & b\\c & d\end{pmatrix}\), \(X = \begin{pmatrix}x\\y\end{pmatrix}\)

Szukanie miejsc w przestrzeni fazowej, gdzie pochodne są równe zero, sprowadza się do rozwiązania układu liniowego i jednorodnego:
\(ax + by = 0\), \(cx + dy = 0\)

Gdy układ równań jest postaci jednorodnej to zawsze mamy minimum jedno rozwiązanie (dla liniowego tylko jedno) - rozwiązanie zerowe \((x,y)=(0,0)\), wtedy gdy \(det(A) \neq 0\). Gdy \(det(A)=0\) to mamy nieskończenie wiele rozwiązań.

W celu znalezienia rozwiązania układu liniowych równań różniczkowych szukamy takich wektorów i wartości, dla których spełnione jest równanie: \(AV_{0} = \lambda V_{0}\), \(V_{0}\) - jest nazywany wektorem własnym (ang. eigenvector); \(\lambda\) - jest nazywana wartością własną (ang. eigenvalue);

W układach liniowych zawsze znajdziemy rozwiązanie wzdłuż wektorów własnych, gdyż to jest specyfika układów liniowych

Wtedy rozwiązaniem układu jest:
\(X(t) = e^{\lambda t}V_{0}\), To rozwiązanie jest na linii prostej - wzdłuż wektora V0. To jest analogiczne do rozwiązania jednowymiarowego układu dynamicznego. Pierwsza współrzędna vx odpowiada osi X a vy osi Y.
\(V_{0} = \begin{pmatrix}v_{x}\\v_{y}\end{pmatrix}\) Dowód:
\(\dot X(t) = \lambda e^{\lambda t} V_{0}\)
\(\dot X(t) = e^{\lambda t} \lambda V_{0}\)
\(\dot X(t) = e^{\lambda t} A V_{0}\)
\(\dot X(t) = AX(t)\)

Rozwiązaniem jest również:
\(X(t) = \alpha e^{\lambda t}V_{0}\), \(\alpha\) jest dowolną stałą różną od zera.
\(A \alpha V_{0} = \alpha A V_{0} = \lambda \alpha V_{0}\) - można przeprowadzić dodatkowo dowód, co wyżej.

Podejście do rozwiązania równania zatem polega na znalezieniu wektora i wartości własnych układu równań. Aby znaleźć wektor własny musimy znaleźć niezerowe rozwiązanie równania. Mysi być niezerowe, aby znaleźć niezerowy wektor:

\(A\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = \lambda \begin{pmatrix}v_{x}\\v_{y}\end{pmatrix}\)

To równanie możemy również zapisać jako: \((A-\lambda I)V_{0}=0\). Równanie ma niezerowe rozwiązanie wtedy i tylko wtedy, gdy \(det(A-\lambda I)=0\) - jest to równanie charakterystyczne:
\(det\begin{pmatrix}a-\lambda & b\\c & d-\lambda\end{pmatrix}=0\)

\((a-\lambda)(d-\lambda) - bc = ad-d\lambda-a\lambda+\lambda^2 - bc = \lambda^2-(a+d)\lambda+ad-bc\) - jest to wielomian charakterystyczny.

Rozwiązaniami tego wielomianu charakterystycznego są wartości własne \(\lambda_{1}\) i \(\lambda_{2}\)

W celu obliczenia wektorów własnych podstawiamy wartości własne do równania: \((A-\lambda_{1} I)V_{0}=0\) \((A-\lambda_{2} I)W_{0}=0\)

W ten sposób uzyskujemy dwa wektory własne. Na podstawie tych obliczeń jesteśmy w stanie określić rozwiązania układu, które są liniami prostymi (ang. strigth-line solution):

\(X_{1}\):
\(x(t) = \alpha e^{\lambda_{1} t}v_{x}\) \(y(t) = \alpha e^{\lambda_{1} t}v_{y}\)

\(X_{2}\):
\(x(t) = \beta e^{\lambda_{2} t}w_{x}\) \(y(t) = \beta e^{\lambda_{2} t}w_{y}\)

Oprócz tych dwóch mamy jeszcze rozwiązanie w punkcie równowagi, który jest dla układów jednorodnych zawsze (0,0) - dla tego punktu rozwiązaniem jest: \(x(t) = 0\) i \(y(t) = 0\).

Te trzy rozwiązania są oczywiście niepełnym obrazem dynamiki zawartej w równaniach. Nie wiemy jak się zachowa układ dla dowolnego punktu nie należącego do tych trzech rozwiązań. Nie rozwiązaliśmy jeszcze problemu wartości początkowej tzn. nie mamy ogólnego rozwiązania.

Jeśli wektory własne są liniowo niezależne: \(det\begin{pmatrix}v_{x} & w_{x}\\v_{y} & w_{y}\end{pmatrix}\ne 0\)

to stanowią bazę \(\mathbb R^2\) i możemy znaleźć taki wektor: \(\begin{pmatrix}z_{x}\\z_{y}\end{pmatrix} = \alpha \begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} + \beta \begin{pmatrix}w_{x}\\w_{y}\end{pmatrix}\) \(Z=\alpha V_{0} + \beta W_{0}\)

który będzie unikatową kombinacją \(\alpha\) i \(\beta\). To znaczy, że suma naszych rozwiązań wzdłuż wektorów własnych jest także rozwiązaniem naszego układu równań - wynika to z zasady liniowości:
\(Z(t) = \alpha X_{1} + \beta X_{2}\) - jest unikalnym rozwiązaniem

Mamy więc rozwiązania ogólne postaci:
\(Z(t) = \alpha e^{\lambda_{1} t}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} + \beta e^{\lambda_{2} t}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix}\)

Wartości \(\alpha\) i \(\beta\) określamy na podstawie warunków początkowych, czyli np. dla t = 0 \(\begin{pmatrix}x(0)\\y(0)\end{pmatrix} = \alpha\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} + \beta\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix}\)

Pytanie: czy istnieją takie warunku, że dwie trajektorie tego samego układu się przetną?

Układ dwuwymiarowy będzie graficznie reprezentowany na płaszczyźnie fazowej

3.1.3 Liniowa niezależność wektorów własnych

\(V_{0} = (1,0)\), \(W_{0} = (1,1)\) . Wyznacznik jest różny od zera, więc są liniowo niezależne. To od razu widać jeśli zwizualizujemy sobie te dwa wektory we współrzędnych kartezjańskich - zobaczymy, że się nie pokrywają.

Najczęściej układ współrzędnych kartezjańskich jest oparty o bazę standardową, czyli \(E_{x}=(1,0)\), \(E_{y}=(0,1)\) Zawsze mogę znaleźć taki wektor, który będzie unikatową kombinacją \(\alpha\) i \(\beta\) \(\begin{pmatrix}z_{x}\\z_{y}\end{pmatrix} = \begin{pmatrix}\alpha\\\beta\end{pmatrix} = \alpha \begin{pmatrix}1\\0\end{pmatrix} + \beta \begin{pmatrix}0\\1\end{pmatrix}\)

base_normal.png

Figure 50: Baza standardowa

Tak samo jest dla każdych dwóch wektorów liniowo niezależnych.

Gdy wektory są liniowo zależne np. \(\begin{pmatrix}z_{x}\\z_{y}\end{pmatrix} = \alpha \begin{pmatrix}1\\1\end{pmatrix} + \beta \begin{pmatrix}-1\\-1\end{pmatrix}\) \(\begin{pmatrix}z_{x}\\z_{y}\end{pmatrix} = \begin{pmatrix}\alpha-\beta\\\alpha-\beta\end{pmatrix}\)

Wektor Z będzie leżał na prostej jak pozostałe wektory.

base_new.png

Figure 51: Nowa baza

3.1.4 Przykład: źródło

\(\dot x = x+2y\),\(\dot y = 3y\)
\(A = \begin{pmatrix}1 & 2\\0 & 3\end{pmatrix}\), \(\dot X = AX\)
\(det(A)=\lambda_{1}\lambda_{2}=3>0\) \(\tau = \lambda_{1}+\lambda_{2}=4\)

  1. Rozwiązanie w punkcie równowagi: mamy jedno unikatowe rozwiązanie (0,0):
    \(det\begin{pmatrix}1 & 2\\0 & 3\end{pmatrix}> 0\)
  2. Szukanie rozwiązania poza punktem równowagi
    a. Obliczanie wartości własnych
    \(det\begin{pmatrix}1-\lambda & 2\\0 & 3-\lambda\end{pmatrix} = 0\)
    \(\lambda^2 - 4\lambda+3=0\)
    \(\Delta=4>0\),\(\lambda_{1}=1\),\(\lambda_{2}=3\), \(\lambda_{1}\lambda_{2}=det(A)\)
    \(\lambda_{1,2}=\frac{1}{2}(\tau\pm\sqrt{\tau^2-4det(A)})\),
    \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb R^+\) - punkt stały w takim układzie nazywamy źródłem
    b. Obliczanie wektorów własnych
    \(\begin{pmatrix}1-\lambda_{1} & 2\\0 & 3-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\)
    \(0v_{x} + 2v_{y} = 0\) wektor spełniający to równanie: (1,0),(-1,0)…
    \(V_{0} = \begin{pmatrix}1\\0\end{pmatrix}\)
    \(\begin{pmatrix}1-\lambda_{2} & 2\\0 & 3-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\)
    \(-2v_{x} + 2v_{y} = 0\) wektor spełniający to równanie: (1,1),(-1,-1)
    \(W_{0} = \begin{pmatrix}1\\1\end{pmatrix}\)
    Rozwiązanie układu wzdłuż dwóch wektorów jest następujące:
    \(X_{1}(t) = e^{\lambda_{1} t}V_{0} = e^{t}\begin{pmatrix}1\\0\end{pmatrix}\)
    \(X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{3t}\begin{pmatrix}1\\1\end{pmatrix}\)

\(X_{1}(t)\) i \(X_{1}(t)\) są rozwiązaniami po linii prostej.

Wektory własne są liniowo niezależne więc:
\(Z(t) = \alpha e^{t}\begin{pmatrix}1\\0\end{pmatrix} + \beta e^{3t}\begin{pmatrix}1\\1\end{pmatrix}\)
Jest to rozwiązanie ogólne. Dla dowolnych warunków początkowych rozwiązanie istnieje i jest jednoznaczne. To zapewnia nam zasada liniowości. Trajektorie układów dynamicznych określone tymi samymi równaniami nie przecinają się. Jeśli dla dwóch takich układów istniałby taki wspólny punkt początkowy to nie mogłyby być określone takimi samymi równaniami.

  1. Symulacja w XPP:
    na zdjęciu mamy zaznaczone rozwiązania wzdłuż wektorów własnych oraz rozwiązanie dla wartości początkowej (-1,1).

hirsh1203.png

Figure 52: Źródło

W XPP obliczenia punktów stałych wykonuje się za pomocą poleceń Sing pts —> Go. Pojawią się okna dialogowe informujące o możliwych dodatkowych czynnościach. Jeśli wybierzemy opcję Draw Strong Sets to w przestrzeni fazowej naszego układu pojawi się rozwiązanie wzdłuż prostej określone przez wektor własny dla największej wartości własnej tzn. w omawianym przypadku będzie to \(\lambda_{2}\). Prosta będzie koloru żółtego, co jest informacją o niestabilności tej trajektorii układu. Dodatkowo punkt stały, który znajduje się w początku układu współrzędnych, jest zaznaczony kwadratem, co informuje nas, że jest niestabilny.

hirsh1203_b.png

Figure 53: Źródło - oznaczenie w XPP

par a=1
par b=2
par c=0
par d=3
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field
# Sing Pts —> Go

Zadanie: Na podstawie analizy wartości własnych wyjaśnij dlaczego pole wektorowe rozkłada się w ten a nie inny sposób.

3.1.5 Przykład: nieskończenie wiele źródeł 1

\[a=b=c=d=1\] \[\dot x = ax + by \text{, } \dot y = cx + dy\], \[A = \begin{pmatrix}1 & 1\\1 & 1\end{pmatrix}\], \[\dot X = AX\]

\(\det(A)=\lambda_{1}\lambda_{2}=0\), \(\tau = \lambda_{1}+\lambda_{2} = 2\)

  1. Rozwiązanie w punkcie równowagi: mamy nieskończenie wiele rozwiązań:

\[det\begin{pmatrix}1 & 1\\1 & 1\end{pmatrix} = 0\]

  1. Szukanie rozwiązania poza punktami równowagi
    a. Obliczanie wartości własnych \[det\begin{pmatrix}1-\lambda & 1\\1 & 1-\lambda\end{pmatrix} = 0\] \[\lambda(\lambda - 2)=0\], \[\lambda_{1}=0 \text{, } \lambda_{2}=2 \text{, } \lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb R\] Pytanie: jakie konsekwencje dla dynamiki układu ma wartość \(\lambda_{1}=0\)?

    \[\lambda_{1,2}=\frac{1}{2}(\tau\pm\sqrt{\tau^2-4det(A)})\],

    b. Obliczanie wektorów własnych \[\begin{pmatrix}1-\lambda_{1} & 1\\1 & 1-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\]

    \(v_{x} + v_{y} = 0\) wektor spełniający to równanie: (1,-1),(-1,1)…
    \[V_{0} = \begin{pmatrix}1\\-1\end{pmatrix}\]

    \[\begin{pmatrix}1-\lambda_{2} & 1\\1 & 1-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(-v_{x} + v_{y} = 0\) wektor spełniający to równanie: (1,1),(-1,-1)
    \[W_{0} = \begin{pmatrix}1\\1\end{pmatrix}\]

    Rozwiązanie układu wzdłuż dwóch wektorów jest następujące: \[X_{1}(t) = e^{\lambda_{1} t}V_{0} = \begin{pmatrix}1\\-1\end{pmatrix}\] - układ nie zmienia swojego położenia. \[X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{2t}\begin{pmatrix}1\\1\end{pmatrix}\]

Wektory własne są liniowo niezależne, więc rozwiązanie ogólne jest postaci: \[Z(t) = \alpha \begin{pmatrix}1\\-1\end{pmatrix} + \beta e^{2t}\begin{pmatrix}1\\1\end{pmatrix}\]

  1. Symulacja w XPP:

inf_fixed_points.png

Figure 54: Nieskończenie wiele punktów stałych

par a=1
par b=1
par c=1
par d=1
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field (Grid: 16)
# Dir.filed/flow —> Flow (Grid: 8)

Rozwiązaniami są tylko linie proste, gdyż \(\dot x\) i \(\dot y\) są określone takimi samymi równaniami. \(\frac{\dot y}{\dot x} = \frac{x + y}{x + y} = 1\),

3.1.6 Przykład: nieskończenie wiele źródeł 2

\[\dot x = x+2y,\, \dot y = 3x+6y\] \[A = \begin{pmatrix}1 & 2\\3 & 6\end{pmatrix}\] \[\dot X = AX\]

  1. Rozwiązanie w punkcie równowagi: mamy nieskończenie wiele rozwiązań:

\[det\begin{pmatrix}1 & 2\\3 & 6\end{pmatrix}= \lambda_{1}\lambda_{2} = 0,\,\tau=\lambda_{1}+\lambda_{2}=7\]

  1. Szukanie rozwiązania poza punktami równowagi
    a. Obliczanie wartości własnych
    \[det\begin{pmatrix}1-\lambda & 2\\3 & 6-\lambda\end{pmatrix} = 0\] \[\lambda(\lambda-7)=0\] \(\lambda_{1}=0\),\(\lambda_{2}=7\), \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb R\)

    b. Obliczanie wektorów własnych
    \[\begin{pmatrix}1-\lambda_{1} & 2\\3 & 6-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] \(v_{x} + 2v_{y} = 0\) - wektor spełniający to równanie: (1,-0.5),(-1,0.5)…
    \[V_{0} = \begin{pmatrix}1\\-0.5\end{pmatrix}\] \[\begin{pmatrix}1-\lambda_{2} & 2\\3 & 6-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(-6v_{x} + 2v_{y} = 0\) - wektor spełniający to równanie: (1,3),(-1,-3)
    \[W_{0} = \begin{pmatrix}1\\3\end{pmatrix}\] Rozwiązanie układu wzdłuż dwóch wektorów jest następujące:
    \(X_{1}(t) = e^{\lambda_{1} t}V_{0} = \begin{pmatrix}1\\-0.5\end{pmatrix}\)
    \(X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{7t}\begin{pmatrix}1\\3\end{pmatrix}\)
    Wektor własne są liniowo niezależne więc: \(Z(t) = \alpha \begin{pmatrix}1\\-0.5\end{pmatrix} + \beta e^{7t}\begin{pmatrix}1\\3\end{pmatrix}\)

    1. Symulacja w XPP: na zdjęciu mamy zaznaczone na żółto rozwiązania wzdłuż wektora własnego związanego z \(\lambda_{2}\) oraz rozwiązanie dla wartości początkowej (-1,1). Układ ma nieskończenie wiele punktów stałych, która rozkładają się wzdłuż wektora własnego związanego z \(\lambda_{1}\) i kilka tych punktów jest zaznaczonych za pomocą trójkątów i prostokątów - powinny to być same prostokąty, gdyż punkty są niestabilne.

      hirsh1236.png

      Figure 55: Nieskończenie wiele źródeł

par a=1
par b=2
par c=3
par d=6
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Do plotowania prostych należy ustawić warunku począkowe na punkty stałe np 
# Initialconds —> New —> -0.5 ENTER 1 ENTER 
# Sing pts —> Go —> Print eigenvalues –> True (dostaniemy w konsoli informacje) 
# —> Draw Strong Sets (informacja w postaci graficznej)

Rozwiązaniami są tylko linie proste, gdyż: \(\frac{\dot y}{\dot x} = \frac{x + 2y}{3x + 6y} = \frac{x + 2y}{3(x + 2y)}=\frac{1}{3}\)

Pytanie: jakie są możliwe kształty trajektorii układu i dlaczego?

3.1.7 Przykład: ściek

\[\dot x = -2x+y,\, \dot y = 2x-3y\] \[A = \begin{pmatrix}-2 & 1\\2 & -3\end{pmatrix}\] \[\dot x = AX\]

  1. Rozwiązanie w punkcie równowagi: mamy jedno unikatowe rozwiązanie:

\[det\begin{pmatrix}-2 & 1\\2 & -3\end{pmatrix}=\lambda_{1}\lambda_{2} > 0,\,\tau=\lambda_{1}+\lambda_{2}=-5\]

  1. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \[det\begin{pmatrix}-2-\lambda & 1\\2 & -3-\lambda\end{pmatrix} = 0\] \[\lambda^2+5\lambda+4=0\] \(\Delta>0\), \(\lambda_{1}=-1\),\(\lambda_{2}=-4\), \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb{R}^-\) - punkt stały w takim układzie nazywamy ściekiem

    b. Obliczanie wektorów własnych \[\begin{pmatrix}-2-\lambda_{1} & 1\\2 & -3-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] \(-v_{x} + v_{y} = 0\) wektor spełniający to równanie: (1,1),(-1,-1)…
    \[V_{0} = \begin{pmatrix}1\\1\end{pmatrix}\] \[\begin{pmatrix}-2-\lambda_{2} & 1\\2 & -3-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(2v_{x} + v_{y} = 0\) wektor spełniający to równanie: (1,-2),(-1,2)
    \[W_{0} = \begin{pmatrix}1\\-2\end{pmatrix}\] Rozwiązanie układu wzdłuż dwóch wektorów jest następujące: \(X_{1}(t) = e^{\lambda_{1} t}V_{0} = e^{-t}\begin{pmatrix}1\\1\end{pmatrix}\)
    \(X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{-4t}\begin{pmatrix}1\\-2\end{pmatrix}\)
    Wektor własne są liniowo niezależne więc: \(Z(t) = \alpha e^{-t}\begin{pmatrix}1\\1\end{pmatrix} + \beta e^{-4t}\begin{pmatrix}1\\-2\end{pmatrix}\)

    1. Symulacja w XPP:\\ na zdjęciu mamy zaznaczone na niebiesko (oznacza linię stabilną) rozwiązanie wzdłuż wektora własnego z \(\lambda_{2}\); rozwiązanie dla wektora z \(\lambda_{1}\) jest kolejną widoczną prostą; trzecie rozwiązanie jest dla wartości początkowej (4,1). Punkt stały w początku układu współrzędnych jest zaznaczony za pomocą okręgu (oznacza stabilny punkt).

      hirsh-212-3.png

      Figure 56: Ściek

par a=-2
par b=1
par c=2
par d=-3
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Sing pts —> Go —> Print eigenvalues (True) —> Draw Strong Sets
# Do plotowania rozwiązania wzdłuż drugiego wektora własnego należy 
# scałkować dla warunków począkowych np 
# Initialconds —> New —> 10 ENTER 10 ENTER 

Pytanie: czy układ osiąga punkt stały?

3.1.8 Przykład: siodło

\(\dot x = x+2y\),\(\dot y = x\)
\[A = \begin{pmatrix}1 & 2\\1 & 0\end{pmatrix}\] \[\dot X = AX\]

  1. Rozwiązanie w punkcie równowagi: jedno unikatowe rozwiązanie: \[det\begin{pmatrix}1 & 2\\1 & 0\end{pmatrix} = \lambda_{1}\lambda_{2} < 0,\,\tau=\lambda_{1}+\lambda_{2}=1\]
  2. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \[det\begin{pmatrix}1-\lambda & 2\\1 & 0-\lambda\end{pmatrix} = 0\] \[\lambda^2-\lambda-2=0\] \(\Delta>0\), \(\lambda_{1}=-1\),\(\lambda_{2}=2\), \(\lambda_{1}\in\mathbb R^-\), i, \(\lambda_{2}\in\mathbb R^+\) - punkt stały w takim układzie nazywamy siodłem
    Pytanie: jakie konsekwencje dla dynamiki układu mają wartości \(\lambda_{1}\) i \(\lambda_{2}\)?
    b. Obliczanie wektorów własnych \[\begin{pmatrix}1-\lambda_{1} & 2\\1 & 0-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] \(2v_{x} + 2v_{y} = 0\) wektor spełniający to równanie: (1,-1),(-1,1)…
    \[V_{0} = \begin{pmatrix}1\\-1\end{pmatrix}\] \[\begin{pmatrix}1-\lambda_{2} & 2\\1 & 0-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(-v_{x} + 2v_{y} = 0\) wektor spełniający to równanie: (1,0.5),(-1,-0.5)
    \[W_{0} = \begin{pmatrix}1\\0.5\end{pmatrix}\] Rozwiązanie układu wzdłuż dwóch wektorów jest następujące: \(X_{1}(t) = e^{\lambda_{1} t}V_{0} = e^{-t}\begin{pmatrix}1\\-1\end{pmatrix}\)
    \(X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{2t}\begin{pmatrix}1\\0.5\end{pmatrix}\)
    Wektor własne są liniowo niezależne więc: \(Z(t) = \alpha e^{-t}\begin{pmatrix}1\\-1\end{pmatrix} + \beta e^{2t}\begin{pmatrix}1\\0.5\end{pmatrix}\)
    1. Symulacja w XPP:\\ na zdjęciu mamy zaznaczone na żółto rozwiązanie wzdłuż wektora własnego związanego z \(\lambda_{2}\) (rozmaitość niestabilna) oraz na niebiesko z wektorem związanym z \(\lambda_{1}\) (rozmaitość stabilna); dodatkowo jedno rozwiązanie dla wartości początkowej (1,1). Punkt stały znajdujący się w początku układu współrzędnych jest punktem siodłowym i jest zaznaczony trójkątem.

      hirsh1210.png

      Figure 57: Siodło

par a=1
par b=2
par c=1
par d=0
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Sing pts —> Go —> Print eigenvalues (True) —> Draw Invariant Sets
# Initialconds —> New —> 1 ENTER 1 ENTER 

Pytanie: Jak różnice w wartościach własnych wpływają na zachowanie się układu?
Pytanie: Dlaczego zachowanie układu bardziej determinuje rozwiązanie wzdłuż wektora \(W_{0}\)?

3.1.9 Przykład: trzy podstawowe topologiczne typy punktów stałych

\(\dot x = ax+by\),\(\dot y = cx + ay\)
\(A = \begin{pmatrix}a & b\\c & a\end{pmatrix}\), dla \(b \cdot c > 0\)
\[\dot X = AX\]

  1. Rozwiązanie w punkcie równowagi:
    • jeśli \(a^2 = bc\) to \(det(A) = 0\) —> nieskończenie wiele rozwiązań
    • jeśli \(a^2 \ne bc\) to \(det(A) > 0\) —> jedno rozwiązanie
  2. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \[det\begin{pmatrix}a-\lambda & b\\c & a-\lambda\end{pmatrix} = 0\] \[\lambda^2-2a\lambda+a^2-bc=0\] \(\Delta=4bc>0\), \(\lambda_{1}=\frac{2a-2\sqrt{bc}}{2}=a-\sqrt{bc}\), \(\lambda_{2}=a+\sqrt{bc}\) , \(\lambda_{1}, \lambda_{2}\in\mathbb R\)

    • dla \(\sqrt{bc} > a\)
      • dla \(a > -\sqrt{bc}\) —> siodło
      • dla \(a < -\sqrt{bc}\) —> ściek
    • dla \(\sqrt{bc} < a\) —> źródło

    b. Obliczanie wektorów własnych \[\begin{pmatrix}a-\lambda_{1} & b\\c & a-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] \(\sqrt{bc}v_x + bv_{y} = 0\) wektor spełniający to równanie:
    \[V_{0} = \begin{pmatrix}1\\-\frac{\sqrt{bc}}{b}\end{pmatrix}\] \[\begin{pmatrix}a-\lambda_{2} & b\\c & a-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(-\sqrt{bc}v_x + bv_{y} = 0\) wektor spełniający to równanie:
    \[W_{0} = \begin{pmatrix}1\\\frac{\sqrt{bc}}{b}\end{pmatrix}\] Rozwiązanie układu wzdłuż dwóch wektorów jest następujące: \(X_{1}(t) = e^{\lambda_{1} t}V_{0} = e^{(a-\sqrt{cb})t}\begin{pmatrix}1\\-\frac{\sqrt{bc}}{b}\end{pmatrix}\)
    \(X_{2}(t) = e^{\lambda_{2} t}W_{0} = e^{(a+\sqrt{cb})t}\begin{pmatrix}1\\\frac{\sqrt{bc}}{b}\end{pmatrix}\)
    \(Z(t) = \alpha e^{(a-\sqrt{cb})t}\begin{pmatrix}1\\-\frac{\sqrt{bc}}{b}\end{pmatrix} + \beta e^{(a+\sqrt{cb})t}\begin{pmatrix}1\\\frac{\sqrt{bc}}{b}\end{pmatrix}\)

    1. Symulacja w XPP:

      topo_points_type.gif

      Figure 58: Trzy podstawowe topologiczne typy punktów stałych

par a=-3
par b=2
par c=2
par d=-3
x'=a*x+b*y
y'=c*x+a*y

@ total=10,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000, bounds=100000 
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 2 ENTER 0 ENTER 
# Initialconds —> Range —> Range Over: a —> Steps: 249 —> Start: -3 —> End: 3
# Kinescope —> Make AniGif

3.1.10 Nietłumiony oscylator harmoniczny - zasada zachowania energii (zad 5.1.1)

\[m\ddot x + kx =0\] \[\dot x = v\] \[\dot v = -\frac{k}{m}x = -\omega^2 x\] \[\omega=\sqrt{\frac{k}{m}}\], prędkość kątowa oscylatora; ile oscylacji jest na jeden obrót \[T = \frac{2\pi}{\omega}\], czas pojedynczej oscylacji.

\[A = \begin{pmatrix}0 & 1\\-\omega^2 & 0\end{pmatrix}\]

Pytanie: czy w naszym układzie zachodzi zasada zachowania energii?

\[\frac{\dot v}{\dot x}=-\frac{\omega^2 x}{v}\] \[\int vdv=-\omega^2 \int xdx\] \[0.5v^2 = -0.5\omega^2 x^2 +C\] \[v^2 + \omega^2 x^2 = 2C\]

\(E_k=0.5mv^2\),\(E_p=0.5kx^2\) \[\frac{2E_k}{m}+\frac{k}{m}\frac{2E_p}{k}=2C\] \[E_k+E_p=Cm=const\]

  1. Rozwiązanie w punkcie równowagi: jedno unikatowe rozwiązanie: \[det\begin{pmatrix}0 & 1\\-\omega^2 & 0\end{pmatrix}> 0\]

\[det(A)=\lambda_{1}\lambda_{2} = \omega^2>0 \text{, } \tau = \lambda_{1}+\lambda_{2}=0\]

  1. Szukanie rozwiązania poza punktami równowagi
    a. Obliczanie wartości własnych \(\lambda_{1,2}=\pm\frac{1}{2}(\sqrt{-4\omega^2})=\pm i\omega\), \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb C^+\) - gdy \(\tau = 0\) punkt stały jest środkiem (ang. center)

    b. Obliczanie wektorów własnych \[\begin{pmatrix}-\lambda_{1} & 1\\-\omega^2 & -\lambda_{1}\end{pmatrix}\begin{pmatrix}m_{x}\\m_{y}\end{pmatrix} = 0\] \(-i\omega m_{x} + m_{y} = 0\) wektor spełniający to równanie:
    \[M_{0} = \begin{pmatrix}1\\i\omega\end{pmatrix}\]

\[\begin{pmatrix}-\lambda_{2} & 1\\-\omega^2 & -\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{y}\end{pmatrix} = 0\] \(i\omega w_{x} + w_{y} = 0\) wektor spełniający to równanie:
\[W_{0} = \begin{pmatrix}1\\-i\omega\end{pmatrix}\] Wartości własne są sprzężone ze sobą, więc i wektory własne są sprzężone. Wystarczy, że weźmiemy pod uwagę tylko jeden i skorzystamy ze wzoru Eulera: \(e^{ix} = cos(x) + isin(x)\)
\[X(t) = e^{\lambda_{1} t}M_{0} = e^{i\omega t}\begin{pmatrix}1\\i\omega\end{pmatrix}\] \[X(t) = [cos(\omega t)+isin(\omega t)]\begin{pmatrix}1\\i\omega\end{pmatrix}\]
\[X(t) = \dbinom{cos(\omega t) + i sin(\omega t)}{i\omega cos(\omega t) -\omega sin(\omega t)} = \dbinom{cos(\omega t) + i sin(\omega t)}{ -\omega sin(\omega t) + i\omega cos(\omega t)}\] \[X(t) = X_{re}(t) + iX_{im}(t)\] \[X_{re} = \dbinom{cos(\omega t)}{ -\omega sin(\omega t)} \text{, }X_{im} = \dbinom{sin(\omega t)}{\omega cos(\omega t)} \text{ - liczby rzeczywiste!}\]

\(X'_{re}(t) + iX'_{im}(t)=X'(t) = AX(t)=A(X_{re}(t) + iX_{im}(t)) = AX_{re}(t) + iAX_{im}(t)\)
\(X'_{re}(t) = AX_{re}(t)\), \(X'_{im}(t) = AX_{im}(t)\) - są rozwiązaniami równania. Ponadto są liniowo niezależne, gdyż: \(X_{re}(0) = \dbinom{1}{0}\) oraz \(X_{im}(0) = \dbinom{0}{\omega}\) Rozwiązaniem ogólnym równania jest: \(X(t) = c_{1}X_{re}(t) + c_{2}X_{im}(t)\)

Pytanie: od czego zależy kształt orbity? dla jakich parametrów trajektoria układu będzie miała kształt okręgu?

Pytanie: od czego zależy amplituda?

  1. Symulacja w xpp:

    undamped_oscillator.gif

    Figure 59: Przestrzeń fazowa oscylatora

par a=0
par b=1
par c=-4
x'=a*x+b*y
y'=c*x+a*y

@ total=100,xp=x,yp=y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 2 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: c —> Steps: 249 —> Start: -4 —> End: -0.5
# Kinescope —> Make AniGif

3.1.11 Tłumiony oscylator harmoniczny

\[m\ddot x + c\dot x + kx =0\] \[\dot x = v\] \[\dot v = -\frac{c}{m}v - \frac{k}{m} = -2\zeta\omega v -\omega^2 x\] \[\omega_{0}=\sqrt{\frac{k}{m}}\], prędkość kątowa oscylatora nietłumionego; \(\zeta = \frac{c}{2\sqrt{mk}}\) - współczynnik tłumienia;

\[A = \begin{pmatrix}0 & 1\\-\omega_{0}^2 & -2\zeta\omega_{0} \end{pmatrix}\]

  1. Rozwiązanie w punkcie równowagi: jedno unikatowe rozwiązanie: \[det\begin{pmatrix}0 & 1\\-\omega_{0}^2 & -2\zeta\omega_{0} \end{pmatrix}> 0\]

\[det(A)=\lambda_{1}\lambda_{2}=\omega_{0}^2>0 \text{, } \tau = \lambda_{1}+\lambda_{2}=-2\zeta\omega_{0}\]

  1. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \(\lambda_{1}=\frac{-2\zeta\omega_{0}+(\sqrt{4\zeta^2\omega_{0}^2-4\omega_{0}^2})}{2}=\omega_{0}(-\zeta + \sqrt{\zeta^2-1})\), \(\lambda_{2}=\frac{-2\zeta\omega_{0}-(\sqrt{4\zeta^2\omega_{0}^2-4\omega_{0}^2})}{2}=-\omega_{0}(\zeta + \sqrt{\zeta^2-1})\),
    • Dla \(\zeta \in (-1,1)\):
      • \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb C^+\), \(\tau \ne 0\) —> rozwiązaniem jest wir
      • dla \(\zeta \in (-1,0)\): —> wir niestabilny
      • dla \(\zeta \in (0,1)\): —> wir stabilny
    • Dla \(\zeta \geq 1\):
      • \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb R\), \(\tau \ne 0\) —> rozwiązaniem jest ściek
    • Dla \(\zeta \leq -1\):
      • \(\lambda_{1}\,\text{i}\,\lambda_{2}\in\mathbb R\), \(\tau \ne 0\) —> rozwiązaniem jest źródło

Pytanie: dla jakich parametrów oscylator zachowuje się jak nietłumiony?

Pytanie: dla jakich parametrów oscylator jest tłumiony najszybciej? dla jakich parametrów najszybciej osiąga punkt stały?

Pytanie: jak szybko jest tłumiony oscylator?
Pytanie: ile wynosi prędkość kątowa?

b. Obliczanie wektorów własnych \[\begin{pmatrix}-\lambda_{1} & 1\\-\omega_{0}^2 & -2\zeta\omega_{0}-\lambda_{1}\end{pmatrix}\begin{pmatrix}m_{x}\\m_{v}\end{pmatrix} = 0\] Wektor spełniający to równanie:
\[M_{0} = \begin{pmatrix}1\\-\omega_{0}(\zeta-\sqrt{\zeta^2-1})\end{pmatrix}\]

\[\begin{pmatrix}-\lambda_{2} & 1\\-\omega_{0}^2 & -2\zeta\omega_{0}-\lambda_{2}\end{pmatrix}\begin{pmatrix}w_{x}\\w_{v}\end{pmatrix} = 0\] Wektor spełniający to równanie:
\[W_{0} = \begin{pmatrix}1\\-\omega_{0}(\zeta+\sqrt{\zeta^2-1})\end{pmatrix}\]

Dla \(\zeta \in (-1,1)\): \(\lambda_{1}=-\omega_{0}\zeta + i\omega_{0}\sqrt{1-\zeta^2}\) \(X(t) = e^{\lambda_{1} t}M_{0} = e^{(-\omega_{0}\zeta + i\omega_{0}\sqrt{1-\zeta^2})t}\begin{pmatrix}1\\-\omega_{0}(\zeta-i\sqrt{1-\zeta^2})\end{pmatrix}\)
\(X(t) = e^{-\omega_{0}\zeta t}e^{i\omega_{0}\sqrt{1-\zeta^2}t}\begin{pmatrix}1\\-\omega_{0}(\zeta-i\sqrt{1-\zeta^2})\end{pmatrix}\) Wzór Eulera: \(e^{ix} = cos(x) + isin(x)\)
\(X(t) = e^{-\omega_{0}\zeta t}Z(t)\)
1\(Z(t) = [cos(\omega_{0}\sqrt{1-\zeta^2}t)+isin(\omega_{0}\sqrt{1-\zeta^2}t)]\begin{pmatrix}1\\-\omega_{0}(\zeta-i\sqrt{1-\zeta^2})\end{pmatrix}\)
\(Z(t) = \dbinom{cos(\omega_{0}\sqrt{1-\zeta^2}t)+isin(\omega_{0}\sqrt{1-\zeta^2}t)}{-\omega_{0}(\zeta-i\sqrt{1-\zeta^2}) cos(\omega_{0}\sqrt{1-\zeta^2}t) -i\omega_{0}(\zeta-i\sqrt{1-\zeta^2}) sin(\omega_{0}\sqrt{1-\zeta^2}t)}\) \(Z(t) = \dbinom{cos(\omega_{0}\sqrt{1-\zeta^2}t)+isin(\omega_{0}\sqrt{1-\zeta^2}t)}{-\omega_{0}\zeta cos(\omega_{0}\sqrt{1-\zeta^2}t) +i\omega_{0}\sqrt{1-\zeta^2} cos(\omega_{0}\sqrt{1-\zeta^2}t) -i\omega_{0}\zeta sin(\omega_{0}\sqrt{1-\zeta^2}t) -\omega_{0}\sqrt{1-\zeta^2} sin(\omega_{0}\sqrt{1-\zeta^2}t)}\) \(Z(t) = \dbinom{cos(\omega_{0}\sqrt{1-\zeta^2}t)+isin(\omega_{0}\sqrt{1-\zeta^2}t)}{-\omega_{0}\zeta cos(\omega_{0}\sqrt{1-\zeta^2}t)-\omega_{0}\sqrt{1-\zeta^2} sin(\omega_{0}\sqrt{1-\zeta^2}t) +i\omega_{0}\sqrt{1-\zeta^2} cos(\omega_{0}\sqrt{1-\zeta^2}t) -i\omega_{0}\zeta sin(\omega_{0}\sqrt{1-\zeta^2}t) }\) \(Z(t) = \dbinom{cos(\omega_{0}\sqrt{1-\zeta^2}t)+isin(\omega_{0}\sqrt{1-\zeta^2}t)}{-\omega_{0}[\zeta cos(\omega_{0}\sqrt{1-\zeta^2}t)+\sqrt{1-\zeta^2} sin(\omega_{0}\sqrt{1-\zeta^2}t)] +i\omega_{0}[\sqrt{1-\zeta^2} cos(\omega_{0}\sqrt{1-\zeta^2}t) -\zeta sin(\omega_{0}\sqrt{1-\zeta^2}t)] }\)

\(Z(t) = Z_{re}(t) + iZ_{im}(t)\) \(Z_{re}(t) = \dbinom{cos(\omega_{0}\sqrt{1-\zeta^2}t)}{-\omega_{0}[\zeta cos(\omega_{0}\sqrt{1-\zeta^2}t)+\sqrt{1-\zeta^2} sin(\omega_{0}\sqrt{1-\zeta^2}t)]}\) \(Z_{im}(t) = \dbinom{sin(\omega_{0}\sqrt{1-\zeta^2}t)}{\omega_{0}[\sqrt{1-\zeta^2} cos(\omega_{0}\sqrt{1-\zeta^2}t)-\zeta sin(\omega_{0}\sqrt{1-\zeta^2}t)]}\)

\(Z_{re}(0) = \dbinom{1}{-\omega_{0}\zeta}\) oraz \(Z_{im}(0) = \dbinom{0}{\omega_{0}\sqrt{1-\zeta^2}}\) Rozwiązaniem ogólnym równania jest: \(X(t) = e^{-\omega_{0}\zeta t}[c_{1}Z_{re}(t) + c_{2}Z_{im}(t)]\)

Rozwiązanie szczególne:
Dla t = 0, x(0) = 3 i v(0) = 0, \(\omega=2\), \(\zeta = 1/8\) \(3 = c_{1}\cdot 1 + c_{2}\cdot 0\) \(c_{1} = 3\) \(0 = 3\cdot -\omega_{0}\zeta + c_{2} \cdot \omega_{0}\sqrt{1-\zeta^2}\) \(0 = 3 \cdot -0.25 + c_{2} \cdot 2\sqrt{\frac{63}{64}}\) \(0 = 3 \cdot -0.25 + c_{2} \cdot \frac{\sqrt{63}}{4}\) \(c_{2} = \frac{3}{\sqrt{63}} = \frac{\sqrt{63}}{21}\)

\(\dbinom{x(t)}{v(t)} = e^{-0.25 t}[3 \dbinom{cos(\frac{\sqrt{63}}{4}t)}{-0.25 cos(\frac{\sqrt{63}}{4}t)-\frac{\sqrt{63}}{4} sin(\frac{\sqrt{63}}{4}t)} + \frac{\sqrt{63}}{21} \dbinom{sin(\frac{\sqrt{63}}{4}t)}{\frac{\sqrt{63}}{4} cos(\frac{\sqrt{63}}{4}t)-0.25 sin(\frac{\sqrt{63}}{4}t)]}]\)

  1. Symulacja w xpp:

    damped_oscillator.gif

    Figure 60: Przestrzeń fazowa oscylatora

par a=0
par b=1
par c=-4
par d=-4
x'=a*x+b*y
y'=c*x+d*y

@ total=100,xp=x,yp=y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 3 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: d —> Steps: 249 —> Start: -4 —> End: 4
# Kinescope —> Make AniGif

dyn_damped_oscillator.png

Figure 61: Wykres zmiany położenia oscylatora

Pytanie: Czy dla \(\zeta=1\) układ najszybciej osiąga punkt stały?

\(\lambda_{1}=\omega_{0}(-\zeta +\sqrt{\zeta^2-1})\), \(\lambda_{2}=-\omega_{0}(\zeta + \sqrt{\zeta^2-1})\),

\(M_{0} = \begin{pmatrix}1\\-\omega_{0}(\zeta-\sqrt{\zeta^2-1})\end{pmatrix}\)

\(W_{0} = \begin{pmatrix}1\\-\omega_{0}(\zeta+\sqrt{\zeta^2-1})\end{pmatrix}\)

Dla \(\zeta = 1\):
\(X(t) = \alpha e^{-\omega_{0} t} \dbinom{1}{-\omega_{0}} + \beta e^{-\omega_{0} t} \dbinom{1}{-\omega_{0}}\) Wektory własne nie są liniowo niezależne. Możemy się przekonać, że nasza dotychczasowa metoda w takim przypadku zawodzi:

\(\dbinom{x(t)}{v(t)} = (\alpha +\beta)e^{-\omega_{0} t} \dbinom{1}{-\omega_{0}}\)

Rozwiązanie szczególne:
dla t=0, x(0)=3 i v(0)=0, ω=2, ζ = 1 \(\dbinom{3}{0} = (\alpha +\beta) \dbinom{1}{-2}\)

\(0 = -2\alpha - 2\beta\), \(\alpha=-\beta\)
\(3 = \alpha + \beta\), \(3 = 0\)!!! Znajdziemy rozwiązanie tylko dla warunków początkowych leżących na wektorze własnym: \(\dbinom{2}{-4} = (\alpha +\beta) \dbinom{1}{-2}\) \(2 = \alpha + \beta\), \(\alpha = 2 - \beta\)
\(-4 = -2\alpha - 2\beta\), \(-4 = -4 + \beta - 2\beta\),
\(\beta = 0\), \(\alpha=2\)
\(\dbinom{x(t)}{v(t)} = 2e^{-2 t} \dbinom{1}{-2}\)

Dla \(\zeta > 1\):
\(\lambda_{1}=\omega_{0}(-\zeta + \sqrt{\zeta^2-1})\) , \(\lambda_{2}=-\omega_{0}(\zeta + \sqrt{\zeta^2-1})\) ,

\(\dbinom{x(t)}{v(t)} = \alpha e^{-\omega_{0}(\zeta - \sqrt{\zeta^2-1}) t} \dbinom{1}{-\omega_{0}(\zeta-\sqrt{\zeta^2-1})} + \beta e^{-\omega_{0}(\zeta + \sqrt{\zeta^2-1}) t} \dbinom{1}{-\omega_{0}(\zeta+\sqrt{\zeta^2-1})}\)

Rozwiązanie szczególne:
dla \(t=0, x(0)=3\; \text{i}\; v(0)=0, \omega=2, \zeta = \sqrt{2}\) \(\dbinom{3}{0} = \alpha \dbinom{1}{-2(\sqrt{2}-1)} + \beta \dbinom{1}{-2(\sqrt{2}+1)}\) \(0 = 2(\sqrt{2}-1)\alpha - 2(\sqrt{2}+1)\beta\), \(\alpha=\frac{\sqrt{2}+1}{\sqrt{2}-1}\beta\)
\(3 = \alpha + \beta\), \(3 = \frac{\sqrt{2}+1}{\sqrt{2}-1}\beta + \beta\), \(3 = \frac{2(\sqrt{2}+1)}{\sqrt{2}-1}\beta\),
\(\beta = \frac{3(\sqrt{2}-1)}{2(\sqrt{2}+1)}\), \(\alpha=1.5\)

\(\dbinom{x(t)}{v(t)} = 1.5 e^{-2(\sqrt{2}-1) t} \dbinom{1}{-2(\sqrt{2}-1)} + \frac{3(\sqrt{2}-1)}{2(\sqrt{2}+1)} e^{-2(\sqrt{2}+1) t} \dbinom{1}{-2(\sqrt{2}+1)}\)

dyn_damped_oscillator2.png

Figure 62: Wykres zmiany położenia oscylatora

Dlaczego ten jest wolniejszy?

3.1.12 Wir stabilny i niestabilny

\[A = \begin{pmatrix}a & b\\-b & a \end{pmatrix}\]

  1. Rozwiązanie w punkcie równowagi: jedno unikatowe rozwiązanie: \[det\begin{pmatrix}a & b\\-b & a \end{pmatrix}> 0\]

\[det(A)=\lambda_{1}\lambda_{2}=a^2+b^2>0 \text{, } \tau = \lambda_{1}+\lambda_{2}=2a\]

  1. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \(\lambda_{1}=\frac{2a+(\sqrt{4a^2-4a^2-4b^2})}{2} = a+ib\) \(\lambda_{2}=a-ib\),

    • Dla \(a > 0\) —> wir niestabilny
    • Dla \(a < 0\) —> wir stabilny

    b. Obliczanie wektorów własnych \[\begin{pmatrix}a-\lambda_{1} & b\\-b & a-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] Wektor spełniający to równanie:
    \[V_{0} = \begin{pmatrix}1\\i\end{pmatrix}\] \[\begin{pmatrix}a-\lambda_{2} & b\\-b & a-\lambda_{2}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] Wektor spełniający to równanie:

    \[W_{0} = \begin{pmatrix}1\\-i\end{pmatrix}\]

    \[X(t) = e^{\lambda_{1} t}V_{0} = e^{(a+ib)t}\dbinom{1}{i}=e^{at}e^{ib t}\dbinom{1}{i}\] \[X(t) = e^{at}Z(t)\] \[Z(t) = [cos(bt) + isin(bt)]\dbinom{1}{i}\] \[Z(t) = \dbinom{cos(bt) + isin(bt)}{-sin(bt) + icos(bt)}\] \[Z_{re}(t) = \dbinom{cos(bt)}{-sin(bt)}\] \[Z_{im}(t) = \dbinom{sin(bt)}{cos(bt)}\] Pytanie: czy wektory Zre i Zim są liniowo niezależne? \[Z_{re}(0) = \dbinom{1}{0}\] \[Z_{im}(0) = \dbinom{0}{1}\]

    \[X(t) = e^{at}[c_{1}Z_{re}(t) + c_{2}Z_{im}(t)]\]

  2. Symulacja w XPP

spirals.gif

Figure 63: Przestrzeń fazowa

par a=1
par b=1
x'=a*x+b*y
y'=-b*x+a*y

@ total=100,xp=x,yp=y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 2 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: d —> Steps: 249 —> Start: -4 —> End: 4
# Kinescope —> Make AniGif

3.1.13 Wielokrotne wartości własne 1

\[A = \begin{pmatrix}a & 0\\0 & a \end{pmatrix}\]

  1. Rozwiązanie w punkcie równowagi: jedno unikatowe rozwiązanie:

\[det\begin{pmatrix}a & 0\\0 & a \end{pmatrix}> 0\] \[det(A)=\lambda_{1}\lambda_{2}=a^2>0 \text{, } \tau = \lambda_{1}+\lambda_{2}=2a\]

  1. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \(\lambda_{1}=\frac{2a+(\sqrt{4a^2-4a^2})}{2} = a\) \(\lambda_{2}=a\)

    b. Obliczanie wektorów własnych \[\begin{pmatrix}a-\lambda_{1} & 0\\0 & a-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] Dowolny wektor spełnia te równania:
    Pytanie: jakie będzie wyglądała przestrzeń fazowa?

\[X(t) = e^{a t}\dbinom{1}{1}\] \[X(t) = e^{a t}\dbinom{1}{0}\] \[X(t) = e^{a t}\dbinom{3}{1}\] \[X(t) = e^{a t}\dbinom{-3}{1}\]

  1. Symulacja w XPP

star_phase_space.png

Figure 64: Przestrzeń fazowa

par a=1
par b=0
par c=0
x'=a*x+b*y
y'=c*x+a*y

@ total=100,xp=x,yp=y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 2 ENTER 2 ENTER

3.1.14 Wielokrotne wartości własne 2

\(\dot x=ax+ y\), \(\dot y=ay\)

\[A = \begin{pmatrix}a & 1\\0 & a \end{pmatrix}\]

  1. Rozwiązanie w punkcie równowagi:

\[det\begin{pmatrix}a & 1\\0 & a \end{pmatrix}\] \[det(A)=\lambda_{1}\lambda_{2}=a^2>0 \text{, } \tau = \lambda_{1}+\lambda_{2}=2a\]

  1. Szukanie rozwiązania poza punktami równowagi a. Obliczanie wartości własnych \(\lambda_{1}=\frac{2a+(\sqrt{4a^2-4a^2})}{2} = a\) \(\lambda_{2}=a\)

    Ponieważ wartości własne są takie same, więc wektory własne też takie są. Obliczając wektor własny znajdziemy rozwiązanie wzdłuż tego jednego wektora: b. Obliczanie wektorów własnych \[\begin{pmatrix}a-\lambda_{1} & 1\\0 & a-\lambda_{1}\end{pmatrix}\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix} = 0\] \[V = \dbinom{1}{0}\]

Rozwiązanie wzdłuż wektora własnego:
\(X(t) = \alpha e^{at}\dbinom{1}{0}\)

Pytanie: jak obliczyć pozostałe rozwiązania? Pytanie: jak będzie wyglądała przestrzeń fazowa dla a = 0?

\(\dot y = ay\)
\(y(t) = \beta e^{at}\)
\(\dot x = ax + \beta e^{at}\), nieautonomiczne równanie pierwszego rzędu.

Zgadujemy rozwiązanie:
\(x(t) = \alpha e^{at} + \mu t e^{at}\)

Celem teraz jest pokazać, że pochodna naszego rozwiązania jest równa \(ax + \beta e^{at}\)

\(\dot x = a \alpha e^{at} + \mu e^{at} + \mu t ae^{at}\)

\(\dot x = a (\alpha e^{at} + \mu t e^{at}) + \mu e^{at}\)
\(\dot x = ax + \mu e^{at}\), \(\mu = \beta\)

\[X(t) = \alpha e^{at} \dbinom{1}{0} + \beta e^{at}\dbinom{t}{1}\]

  1. Symulacja w XPP

repeated_eigen.gif

Figure 65: Przestrzeń fazowa

par a=2
par b=1
par c=0
x'=a*x+b*y
y'=c*x+a*y

@ total=100,xp=x,yp=y,xlo=-5,ylo=-5,xhi=5,yhi=5 maxstor=20000
done
# GUI 
# Dir.filed/flow —> Direct Field 
# Initialconds —> New —> 2 ENTER 2 ENTER
# Initialconds —> Range —> Range Over: d —> Steps: 249 —> Start: -2 —> End: 2
# Kinescope —> Make AniGif

3.2 Klasyfikacja dwuwymiarowych układów liniowych

3.2.1 Trzy klasy układów - postać kanoniczna

Każdy układ o rzeczywistych wartościach własnych można sprowadzić do układu o zmiennych niesprzężonych, gdy \(det(A)\neq 0\) \[A^k = \begin{pmatrix}\lambda_{1} & 0\\0 & \lambda_{2} \end{pmatrix}\]

Każdy układ o zespolonych wartościach własnych można sprowadzić do postaci, gdy \(det(A)\neq 0\) \[B^k = \begin{pmatrix}\alpha & \beta\\ -\beta & \alpha \end{pmatrix}\]

Każdy układ o wielokrotnych wartościach własnych i niesprzężonych zmiennych można sprowadzić do postaci, gdy \(det(A)\neq 0\): \[C^k = \begin{pmatrix}\lambda & 1\\0 & \lambda \end{pmatrix}\]

  1. Źródło — postać kanoniczna Ak

    \[A = \begin{pmatrix}1 & 2\\0 & 3\end{pmatrix}\] \[\begin{pmatrix}1-\lambda_{1} & 2\\0 & 3-\lambda_{2}\end{pmatrix}=(1-\lambda_{1})(3-\lambda_{3})=0\] \[\lambda_{1}=1 \text{, } \lambda_{2}=3\]

    Transformacja złożona z wektorów własnych \[T = \begin{pmatrix}1 & 1\\0 & 1\end{pmatrix}\] \(T^{-1}AT=\begin{pmatrix}\lambda_{1} & 0\\0 & \lambda_{2} \end{pmatrix}\), \(T^{-1}\) - macierz odwrotna

    \(T^{-1}= \frac{1}{det(T)}T^D\), TD to macierz dołączona (ang. adjugate matrix)
    \(T_{ij} = (-1)^{i+j}det_{ij}T\) - dopełnienie algebraiczne macierzy T; \(det_{ij}T\) - wyznacznik podmacierzy
    \([T_{ij}]\) - macierz dopełnień algebraicznych
    \(T^D = [T_{ij}]^T = [T_{ji}]\), macierz dołączona powstaje poprzez transpozycje macierzy dopełnień

    \[[T_{ij}] = \begin{pmatrix}1 & 0\\-1 & 1\end{pmatrix}\] \[T^{-1} = \begin{pmatrix}1 & -1\\0 & 1\end{pmatrix}\]

    \[A^k= \begin{pmatrix}1 & -1\\0 & 1\end{pmatrix}\begin{pmatrix}1 & 2\\0 & 3\end{pmatrix}\begin{pmatrix}1 & 1\\0 & 1\end{pmatrix}=\begin{pmatrix}1 & 0\\0 & 3\end{pmatrix}=\begin{pmatrix}\lambda_{1} & 0\\0 & \lambda_{2}\end{pmatrix}\]

    using LinearAlgebra
    T = [1 1;0 1]; A=[1 2;0 3];
    Aᵏ = inv(T)*A*T
    

    transormation_real.png

    Figure 66: Przestrzeń fazowa przed i po transformacji

  2. Ściek - postać kanoniczna Ak

    \[A = \begin{pmatrix}-2 & 1\\2 & -3\end{pmatrix}\]

    \[\lambda_{1}=-1 \text{, } \lambda_{2}=-4\]

    Transformacja złożona z wektorów własnych:
    \[T = \begin{pmatrix}1 & 1\\1 & -2\end{pmatrix}\] \[[T_{ij}] = \begin{pmatrix}-2 & -1\\-1 & 1\end{pmatrix}\] \[T^{-1} = -\frac{1}{3}\begin{pmatrix}-2 & -1\\-1 & 1\end{pmatrix}\]

    \[A^k = -\frac{1}{3}\begin{pmatrix}-2 & -1\\-1 & 1\end{pmatrix}\begin{pmatrix}-2 & 1\\2 & -3\end{pmatrix}\begin{pmatrix}1 & 1\\1 & -2\end{pmatrix}=\begin{pmatrix}-1 & 0\\0 & -4\end{pmatrix}=\begin{pmatrix}\lambda_{1} & 0\\0 & \lambda_{2}\end{pmatrix}\]

    using LinearAlgebra
    T = [1 1;1 -2]; A=[-2 1;2 -3];
    Aᵏ = inv(T)*A*T
    

    transormation_real_sink.png

    Figure 67: Przestrzeń fazowa przed i po transformacji

  3. Siodło - postać kanoniczna Ak

    \[A = \begin{pmatrix}1 & 2\\1 & 0\end{pmatrix}\]

    \[\lambda_{1}=-1 \text{, } \lambda_{2}=2\]

    \[T = \begin{pmatrix}1 & 1\\-1 & 0.5\end{pmatrix}\] \[[T_{ij}] = \begin{pmatrix}0.5 & 1\\-1 & 1\end{pmatrix}\] \[T^{-1} = \frac{2}{3}\begin{pmatrix}0.5 & -1\\1 & 1\end{pmatrix}\]

    \[A^k = \frac{2}{3}\begin{pmatrix}0.5 & -1\\1 & 1\end{pmatrix}\begin{pmatrix}1 & 2\\1 & 0\end{pmatrix}\begin{pmatrix}1 & 1\\-1 & 0.5\end{pmatrix}=\begin{pmatrix}-1 & 0\\0 & 2\end{pmatrix}=\begin{pmatrix}\lambda_{1} & 0\\0 & \lambda_{2}\end{pmatrix}\]

    using LinearAlgebra
    T = [1 1;-1 0.5]; A=[1 2;1 0];
    Aᵏ = inv(T)*A*T
    

    transormation_real_saddle.png

    Figure 68: Przestrzeń fazowa przed i po transformacji

  4. Wielokrotne wartości własne - postać kanoniczna Ak

    \[A = \begin{pmatrix}2 & 0\\0 & 2 \end{pmatrix}\]

    \[T = \begin{pmatrix}2 & 1\\1 & 1\end{pmatrix}\] \[[T_{ij}] = \begin{pmatrix}1 & -1\\-1 & 2\end{pmatrix}\] \[T^{-1} = \begin{pmatrix}\frac{2}{3} & -\frac{1}{3}\\-\frac{1}{3} & \frac{2}{3}\end{pmatrix}\]

    \[A^k = \begin{pmatrix}2 & 0\\0 & 2\end{pmatrix}\]

    using LinearAlgebra
    T = [4 1;1 2]; A=[2 0;0 2];
    Aᵏ = inv(T)*A*T
    
  5. Tłumiony oscylator harmoniczny - postać kanoniczna Bk

    \[A = \begin{pmatrix}0 & 1\\-\omega_{0}^2 & -2\zeta\omega_{0} \end{pmatrix}\]

    \[A = \begin{pmatrix}0 & 1\\-4 & -2 \end{pmatrix}\]

    \[\lambda_{1}=-\omega_{0}\zeta + i\omega_{0}\sqrt{1-\zeta^2} = -1 + i\sqrt{3} = \alpha + i\beta\]

    \[\omega_{0}=2, \zeta=0.5,\; \text{wir ściekowy}\]

    \[Z_{re}(0) = \dbinom{1}{-\omega_{0}\zeta} \text{, } Z_{im}(0) = \dbinom{0}{\omega_{0}\sqrt{1-\zeta^2}}\]

    \[T = \begin{pmatrix}1 & 0\\-1 & \sqrt{3} \end{pmatrix}\] \[[T_{ij}] = \begin{pmatrix}\sqrt{3} & 1\\0 & 1\end{pmatrix}\]

    \[T^{-1} = \frac{\sqrt{3}}{3}\begin{pmatrix}\sqrt{3} & 0\\1 & 1\end{pmatrix}\]

    \[B^k = \frac{\sqrt{3}}{3}\begin{pmatrix}\sqrt{3} & 0\\1 & 1\end{pmatrix}\begin{pmatrix}0 & 1\\-4 & -2\end{pmatrix}\begin{pmatrix}1 & 0\\-1 & \sqrt{3}\end{pmatrix}=\begin{pmatrix}-1 & \sqrt{3}\\-\sqrt{3} & -1\end{pmatrix}=\begin{pmatrix}\alpha & \beta\\-\beta & \alpha\end{pmatrix}\]

    using LinearAlgebra
    T = [1 0;-1 sqrt(3)]; A=[0 1;-4 -2];
    Bᵏ = inv(T)*A*T
    

    transormation_spiral.png

    Figure 69: Przestrzeń fazowa przed i po transformacji

    \(\omega_{0}=2, \zeta=0\) - punkt środkowy

    \[A = \begin{pmatrix}0 & 1\\-\omega_{0}^2 & -2\zeta\omega_{0} \end{pmatrix}\] \[A = \begin{pmatrix}0 & 1\\-4 & 0 \end{pmatrix}\]

    Postać ogólna wartości własnej:\(\lambda_{1,2}=-\omega_{0}\zeta \pm i\omega_{0}\sqrt{1-\zeta^2} = \alpha \pm i\beta\)
    Dla naszego przypadku:\(\lambda_{1}= i2 = i\beta\)

    \[Z_{re}(0) = \dbinom{1}{0} \text{, } Z_{im}(0) = \dbinom{0}{2}\]

    \[T = \begin{pmatrix}1 & 0\\0 & 2 \end{pmatrix}\]

    \[[T_{ij}] = \begin{pmatrix}2 & 0\\0 & 1\end{pmatrix}\]

    \[T^{-1} = \begin{pmatrix}1 & 0\\0 & 0.5\end{pmatrix}\]

    \[B^k = \begin{pmatrix}1 & 0\\0 & 0.5\end{pmatrix}\begin{pmatrix}0 & 1\\-4 & 0\end{pmatrix}\begin{pmatrix}1 & 0\\0 & 2\end{pmatrix}=\begin{pmatrix}0 & 2\\-2 & 0\end{pmatrix}=\begin{pmatrix}\alpha & \beta\\-\beta & \alpha\end{pmatrix}\]

    using LinearAlgebra
    T = [1 0;0 2]; A=[0 1;-4 0];
    Bᵏ = inv(T)*A*T
    

    transormation_center.png

    Figure 70: Przestrzeń fazowa przed i po transformacji

  6. Tłumiony oscylator harmoniczny - postać kanoniczna Ck

    \[A = \begin{pmatrix}0 & 1\\-\omega_{0}^2 & -2\zeta\omega_{0} \end{pmatrix}\]

    \[\omega_{0}=2,\; \zeta=1,\; \text{wielokrotne wartości własne}\]

    \[A = \begin{pmatrix}0 & 1\\-4 & -4 \end{pmatrix}\]

    \[\lambda_{1}=\omega_{0}(-\zeta +\sqrt{\zeta^2-1}) = -2\] \[\lambda_{2}=-\omega_{0}(\zeta + \sqrt{\zeta^2-1}) = -2\]

    \[M_{0} = \begin{pmatrix}1\\-\omega_{0}(\zeta-\sqrt{\zeta^2-1})\end{pmatrix} = \dbinom{1}{-2}\]

    \[C^k = \begin{pmatrix}-2 & 1\\0 & -2\end{pmatrix}=\begin{pmatrix}\lambda & 1\\0 & \lambda\end{pmatrix}\]

    transormation_repeated.png

    Figure 71: Przestrzeń fazowa przed i po transformacji

3.2.2 Klasyfikacja geometryczna

linear_classification.png

Figure 72: Klasyfikacja układów liniowych

\[A = \begin{pmatrix}a & b\\c & d\end{pmatrix}\] \[det(A) = ad - cb\] \[det\begin{pmatrix}a-\lambda & b\\c & d-\lambda\end{pmatrix} = 0\] \[ad - a\lambda - d\lambda + \lambda^2 - cb = 0\] \[\lambda^2 - \lambda(a + d) + ad - cb = 0\] \(\tau = a + d\), \(\tau\) - jest sumą elementów na głównej przekątnej macierzy i nazywamy śladem macierzy (ang. trace); jest niezależny od układu współrzędnych: sprawdź przykłady w poprzednim rozdziale. \[\lambda^2 - \lambda\tau + det(A)\]

\[\Delta = \tau^2 - 4det(A)\]

\[\lambda_{1,2} = \frac{1}{2}(\tau \pm \sqrt{\Delta})\]

Wzory Viete'a wiążą współczynniki wielomianu z jego pierwiastkami:
\[ax^2 + bx + c\] \[x_{1}x_{1} = \frac{c}{a}\], \[x_{1} + x_{1} = -\frac{b}{a}\]

Odnosząc je do naszego wielomianu charakterystycznego dostaniemy:
\(\lambda_{1}\lambda_{2}=det(A)\), \(\lambda_{1} + \lambda_{2} = \tau\)

Trzy podstawowe topologiczne typy punktów stałych to: ścieki, źródła i siodła. Każdy ściek jest asymptotycznie stabilny, a źródła i siodła są niestabilne. Wedle rozróżnienia algebraicznego ściek nazywamy węzłem stabilnym, a źródło węzłem niestabilnym, gdyż wszystkie ich wartości własne są rzeczywiste. W układzie z wartościami zespolonymi istnieją wiry stabilne (ściek spiralny) , niestabilne (źródło spiralne) oraz punkty środkowe (orbity okresowe).

Nieizolowane punkty stałe i środki są neutralnie stabilne, gdyż nie są ściekami, ale są stabilne w sensie Lapunowa. Stabilność Lapunowa oznacza, że trajektoria układu wystarczająco blisko punktu stałego pozostaje blisko tego punktu dla każdego t.

https://experiences.math.cnrs.fr/Pendule-simple.html

4 Dwuwymiarowe układy nieliniowe

4.1 Twierdzenie o istnieniu i jednoznaczności

Rozwiązanie analityczne układu równań nieliniowych zazwyczaj nie jest możliwe, a te rozwiązania, które można osiągnąć, często nie mówią na pierwszy rzut oka zbyt wiele. Nasze podejście będzie skupiało się na określeniu jakościowych cech dynamiki układu.

Istnienie i jednoznaczność rozwiązania jest spełnione jeśli funkcja definiująca zachowanie układu jest różniczkowalnie ciągła t.j. odpowiednio gładka. Jest to warunek istnienia i jednoznaczności nieliniowych równań różniczkowych. Jednoznaczność oznacza, że punkty trajektorii się nie przecinają. Gdyby tak było to rozwiązanie problemu warunku początkowego nie byłoby możliwe, gdyż dawałoby dwie możliwości ruchu punktu fazowego.

4.2 Stany równowagi w układach nieliniowych

Szukanie rozwiązania w punktach równowagi w układach nieliniowych polega na potraktowaniu układu nieliniowego jako liniowy - dokonujemy linearyzacji wokół punktów stałych. W przypadku liniowych mieliśmy tylko jeden albo nieskończenie wiele punktów równowagi. W układach nieliniowych jest więcej możliwości.

\[\dot x = f(x,y)\] \[\dot y = g(x,y)\]

Punkty stałe: \[x_{*},y_{*}\]

Odchylenie od punktów stałych: \[u=x-x_{*}\] \[v=y-y_{*}\]

\[\dot u = \dot x = f(x_{*}+u,y_{*}+v)=\] \[= f(x_{*},y_{*}) + (x_{*}+u-x_{*}) \frac{\partial f}{\partial x} + (y_{*}+v-y_{*}) \frac{\partial f}{\partial y} + \mathcal{O}(u^2,v^2,uv)=\] \[= u\frac{\partial f}{\partial x} + v\frac{\partial f}{\partial y} + \mathcal{O}(u^2,v^2,uv)\]

\[\dot v = \dot y = g(x_{*}+u,y_{*}+v)=\] \[= g(x_{*},y_{*}) + (x_{*}+u-x_{*}) \frac{\partial g}{\partial x} + (y_{*}+v-y_{*}) \frac{\partial g}{\partial y} + \mathcal{O}(u^2,v^2,uv)=\] \[= u\frac{\partial g}{\partial x} + v\frac{\partial g}{\partial y} + \mathcal{O}(u^2,v^2,uv)\]

Ponieważ badamy bardzo bliskie otoczenie punktów x* i y* to u i v są bardzo małe, więc kwadrat ich wartości będzie o wiele mniejszy. Dlatego opuszczamy dalsze rozwiniecie szeregu Taylora i wszystkie nieznaczące (nieliniowe) czynniki umieszczamy w:
\(\mathcal{O}(u^2,v^2,uv)\)

\[\dbinom{\dot u}{\dot v}=\begin{pmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y}\\\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}\end{pmatrix}\dbinom{u}{v} + \mathcal{O}(u^2,v^2,uv)\]

Gdy pominiemy nieznaczący element, otrzymujemy znany nam z poprzedniego rozdziału układ równań liniowych.:

\[\dbinom{\dot u}{\dot v}=\begin{pmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y}\\\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}\end{pmatrix}\dbinom{u}{v}\]

\[A = \begin{pmatrix}\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y}\\\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}\end{pmatrix}_{(x_{*},y_{*})}\]

Macierz A jest nazywaną Jakobianem.

Pytanie: czy pominięcie nieliniowego czynnika może mieć wpływ na wyniki analizy?

Tylko w przypadku, gdy mamy przypadki graniczne t.j. środki oraz układ z minimum jedną zerową wartości własną, wyniki analizy linearyzacji mogą się różnić.

Gdy z analizy liniowej będzie wynikało źródło, ściek, siodło, wir to na pewno taki sam charakter będą miały te punkty w analizowanym układzie nieliniowym (Twierdzenie Hartmana-Grobmana).

4.2.1 Przykład: linearyzacja wokół punktów stałych

\[\dot x = -x + x^3, \dot y = -2y\] \(-x + x^3 =0, x(-1+x^2)=0, x_{0}=0, x_{0} = \pm 1\) \[y_{0}=0\] \[A = \begin{pmatrix}-1 +3x^2 & 0\\ 0 & -2\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}-1 & 0\\ 0 & -2\end{pmatrix}, det(A)=2>0, \tau=-3<0, \text{ściek}\] \[A(-1,0) = \begin{pmatrix}2 & 0\\ 0 & -2\end{pmatrix}, det(A)=-4<0, \tau=0, \text{siodło}\] \[A(1,0) = \begin{pmatrix}2 & 0\\ 0 & -2\end{pmatrix}, det(A)=-4<0, \tau=0, \text{siodło}\]

Wartości własne wynikają bezpośrednio z macierzy, gdyż jest to układ o zmiennych niesprzężonych. To samo dotyczy wektorów własnych. Są one zgodne z bazą standardową tylko siodła są przesunięte o 1 w lewo i prawo.

  1. Symulacja w XPP

example631_points.png

Figure 73: Przestrzeń fazowa

x'=-x+x^3 
y'=-2*y 

@ total=10,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5, maxstor=20000, bounds=100000, nmesh=800 
done
# GUI 
# Dir.filed/flow —> Direct Field —> 16
# Dir.filed/flow —> Flow —> 8
# Sing pts —> monteCar —> Domyślne ustawienia

4.2.2 Przykład: wrażliwość linearyzacji na nieliniowość

\[\dot x = -y + ax(x^2+y^2), \dot y = x+ay(x^2+y^2)\] \(-y + ax(x^2+y^2) = 0, x+ay(x^2+y^2) = 0\) \[a= \frac{y}{x(x^2+y^2)}\] \[x+ \frac{y^2}{x(x^2+y^2)}(x^2+y^2) = 0\] \[x+ \frac{y^2}{x} = 0, y^2 = -x^2, x_{0}=0, y_{0}=0\] \[A = \begin{pmatrix}3ax^2+ay^2 & -1+2axy0\\ 1+2ayx & ax^2+3ay^2\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & -1\\ 1 & 0\end{pmatrix}, det(A)=1>0, \tau=0, \text{punkt środkowy?}\]

Z analizy wynika, że punkt (0,0) jest środkiem. Powiedzieliśmy sobie jednak wcześniej, że w przypadkach granicznych, a takim jest punkt środkowy, możemy dostać błędny wynik. Przy niewielkich wahaniach ślad macierzy może być zarówno większy od zera, bądź mniejszy. W pierwszym przypadku nasz punkt byłby wirem niestabilnym, a w drugim wirem stabilnym.

Zamieńmy układ współrzędnych na polarny i przeprowadźmy jeszcze raz analizę.

\[x=rcos\theta, y=rsin\theta\] \[x^2+y^2=r^2 |\frac{d}{dt}\] \[\frac{d}{dt}r^2 = \frac{d}{dt}x^2 + \frac{d}{dt}y^2\] \[\frac{dr}{dt}r = \frac{dx}{dt}x + \frac{dy}{dt}y\] \[\dot r r = \dot x x + \dot y y\] \[\dot x = -y + ax(x^2+y^2), \dot y = x+ay(x^2+y^2)\] \[\dot r r = [-y + ax(x^2+y^2)] x + [x+ay(x^2+y^2)] y\] \[\dot r r = -yx + ax^2(x^2+y^2) + xy+ ay^2(x^2+y^2)\] \[\dot r r = ar^4cos\theta^2 + ar^4sin\theta^2\] \[\dot r r = ar^4\] \[\dot r = ar^3\] \[r_{0}=0\] \[f'(r) = 3ar^2\] \[f'(0) = 0, \text{analogiczny rezultat jak dla układu 2D}\] \[cos\theta = \frac{x}{r}, sin\theta = \frac{y}{r}\] \[\frac{y}{x} = tan\theta\] \[\frac{\dot y x - \dot x y}{x^2} = \frac{1}{cos\theta^2} \dot \theta\] \[\frac{\dot y x - \dot x y}{x^2} = \frac{r^2}{x^2} \dot \theta\] \[\frac{\dot y x - \dot x y}{r^2} = \dot \theta\] \[\frac{\dot y x - \dot x y}{x^2+y^2} = \dot \theta\] \[\dot x = -y + ax(x^2+y^2), \dot y = x+ay(x^2+y^2)\] \[\frac{[x+ay(x^2+y^2)] x - [-y + ax(x^2+y^2)] y}{x^2+y^2} = \dot \theta\] \[\dot \theta = 1\] \[\dot r = ar^3, \text{układ o zmiennych niesprzężonych}\]

Z analizy przestrzeni fazowej układu we współrzędnych polarnych wynika, że dla a > 0 układ jest niestabilny; dla a < 0 jest stabilny.

  1. Symulacja w XPP

example632.png

Figure 74: Przestrzeń fazowa

par a=-1
x'=-y+a*x*(x^2+y^2)
y'=x+a*y*(x^2+y^2)

@ total=10,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5, maxstor=20000, bounds=100000, nmesh=800 
done
# GUI 
# Dir.filed/flow —> Direct Field —> 16
# Dir.filed/flow —> Flow —> 8
# Sing pts —> Go —> Print eigenvalues (True) 
# W konsoli dostaniemy obliczone wartości właasne. 

example632.gif

Figure 75: Przestrzeń fazowa

par a=-1
x'=-y+a*x*(x^2+y^2)
y'=x+a*y*(x^2+y^2)

@ total=10,xp=X,yp=Y,xlo=-5,ylo=-5,xhi=5,yhi=5, maxstor=20000, bounds=100000, nmesh=800 
done
# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 1 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: a —> Steps: 249 —> Start: -1 —> End: 1
# Kinescope —> Make AniGif

Analiza numeryczna potwierdziła nasze wcześniejsze wyniki, że jednak punkt (0,0) nie jest środkiem a wirem. To ma ogromne znacznie na zachowanie się układu, gdyż wir może być niestabilny! W konsoli mamy obliczenia wartości własnych. Ich część rzeczywista jest tak mała, że widzimy same zera. Nie mogą być równe zero, gdyż w takim przypadku punkt (0,0) byłby środkiem.

Pytanie: dlaczego mamy pustą przestrzeń fazową w środku wiru?

4.2.3 Hiperboliczne punkty stałe, równoważność topologiczna, stabilność strukturalna

Hiperboliczne punkty stałe to takie punkty, których część rzeczywista wszystkich wartości własnych jest różna od zera: \(Re(\lambda)\ne 0\). Taki układ jest odporny na wpływ nielinowego czynnika na linearyzacje. Warunek ten odnosi się do stabilności i zapewnia nam, że stabilność punktów obliczona metodą linearyzacji pozostanie taka sama dla układu nieliniowego. Podczas klasyfikowania układów liniowych rozróżniliśmy klasyfikację topologiczną i algebraiczną. W sensie topologicznym węzeł stabilny, zdegenerowany węzeł stabilny, gwiazda stabilna oraz wir stabilny są takie same. Mówimy, że zachodzi równoważność topologiczna, którą nazywamy homeomorfizmem. To znaczy istnieje funkcja ciągła, która w sposób jednoznaczny mapuje jeden układ na drugi. Mapowanie jednoznaczne zawiera w sobie założenie, że istnieje ciągła funkcja odwrotna. Węzeł stabilny i niestabilny nie są topologicznie równoważne, gdyż nie jest zachowany kierunek pola wektorowego. Punkty, dla których nie jest spełniony warunek, to środek i nieizolowane punkty stałe (ang. non-isolated fixed points). Punkty, które są odporne na drobne perturbacje są określane jako strukturalnie stabilne

4.2.4 Równania konkurencji Lotki-Volterry - rozwinięcie modelu logistycznego liczebności populacji

\[\dot N = r \cdot N\biggl(1 - \frac{N}{K}\biggr)\] Jest to model logistyczny liczebności populacji, gdzie:

  • N: liczebność populacji
  • r: współczynnik tempa wzrostu populacji
  • K: graniczna pojemność otoczenia dla danej populacji

Ten model omawialiśmy podczas dyskusji na temat jednowymiarowych układów dynamicznych. Modele konkurencji Lotki-Volterry (ang. Competitive Lotka–Volterra equations) dotyczą modelowania dynamiki populacji gatunków konkurujących o jakieś wspólne zasoby.

\[\dot x = r_{1} \cdot x\biggl(1 - \frac{x+\alpha_{12}y}{K_{1}}\biggr)\] \[\dot y = r_{2} \cdot y\biggl(1 - \frac{y + \alpha_{21}x}{K_{2}}\biggr)\] Założenia do modelu wrażone w wartościach parametrów: \[r_{1}=K_{1}=3, \alpha_{12}=2, r_{2}=K_{2}=2, \alpha_{21}=1\] Polucja gatunku x rośnie 1.5 razy szybciej i jej pojemność jest 1.5 większa niż populacja Y. Wpływ populacji Y na X jest dwa razy większy niż X na Y.

\[\dot x = x(3 - x-2y)\] \[\dot y = y(2 - y-x)\]

Punkt stałe:
\[x(3 - x-2y)=0,y(2 - y-x)=0\] \[x_{0}=0,y_{0}=0\] \[x_{0}=0,y_{0}=2\] \[x_{0}=3,y_{0}=0\] \[3 - x-2y=0,2 - y-x=0\] \[3 - x-2y=0 \rightarrow x=3-2y \rightarrow 2-y-3+2y=0\] \[y_{0}=1, x_{0}=1\]

\[A = \begin{pmatrix}3 - 2x - 2y & -2x\\ -y & 2-2y-x\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}3 & 0\\ 0 & 2\end{pmatrix}, det(A)=6>0, \tau=5>0, \text{źródło}\] \[A(0,2) = \begin{pmatrix}-1 & 0\\ -2 & -2\end{pmatrix}, det(A)=2>0, \tau=-3<0, \text{ściek}\] \[A(3,0) = \begin{pmatrix}-3 & -6\\ 0 & -1\end{pmatrix}, det(A)=3>0, \tau=-4<0, \text{ściek}\] \[A(1,1) = \begin{pmatrix}-1 & -2\\ -1 & -1\end{pmatrix}, det(A)=-1<0, \tau=-2<0, \text{siodło}\]

lotka_volterra.png

Figure 76: Przestrzeń fazowa

par r1=3
par K1=3
par a12=2
par r2=2
par K2=2
par a21=1
x' = r1*x*(1 - (x+a12*y)/K1)
y' = r2*y*(1 - (y+a21*x)/K2)

@ total=10,xp=X,yp=Y,xlo=0,ylo=0,xhi=3,yhi=3, maxstor=2000000, bounds=100000, nmesh=800 
done
# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Dir.filed/flow —> Flow –> 8
# Sing pts —> monteCar —> Shoot 1; pozostałe domyślne ustawienia 

Pytanie: dla jakich warunków początkowych wygrywa gatunek Y?

Basen przyciągania jest to obszar w przestrzeni fazowej, w której każdy punkt zmierza do tego samego punktu stałego. Na rysunku 76 rozmaitość stabilna zaznaczona na niebiesko dzieli na dwa takie obszary i nazywana jest przez to granicą basenu (ang. basin boundary). Dla punktów powyżej tej krzywej wszystkie one dążą do punktu stałego (0,2); poniżej tej krzywej dążą do punktu (3,0): \(X(t) \rightarrow X_{0}, t \rightarrow \infty\)

4.2.5 Badanie wpływa parametru α na zachowanie się układu

lotka_volterra_a12.gif

Figure 77: Przestrzeń fazowa

par r1=2
par K1=2
par a12=2
par r2=2
par K2=2
par a21=2
x' = r1*x*(1 - (x+a12*y)/K1)
y' = r2*y*(1 - (y+a21*x)/K2)

@ total=10,xp=X,yp=Y,xlo=0,ylo=0,xhi=3,yhi=3, maxstor=2000000, bounds=100000, nmesh=800 
done
# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 2 ENTER 2.5 ENTER
# Initialconds —> Range —> Range Over: a12 —> Steps: 249 —> Start: 1.5 —> End: 2.2
# Kinescope —> Make AniGif

4.2.6 Układy zachowawcze (ang. conservative systems)

W układzie istnieje pewna wielkość np. energia całkowita E(t), która się nie zmienia wzdłuż trajektorii układu:

\[\frac{dE}{dt}= 0 \]

jeśli istnieje funkcja odpowiednio gładka (klasy \(C^k,\; k \ge 1\)) \(U: R^n \rightarrow R\) taka że: \(F(X) = -\Biggl(\frac{\partial U}{\partial x_{1}}(X),\frac{\partial U}{\partial x_{2}}(X),...,\frac{\partial U}{\partial x_{n}}(X)\Biggr)=-grad\,U(X)=-\nabla U(X)\). Taki układ jest nazywany zachowawczym. Funkcja gładka U(X) jest nazywana energią potencjalną układu. Jest ona zależna od położenia w przestrzeni, czyli w naszym przypadku od zmiennej x.

Dla układów mechanicznych ma postać: \[m\ddot x = F(X)\] \[\dot x = v\] \[\dot v = -\frac{1}{m}\nabla U(X)\]

Twierdzenie o nieliniowych środkach w układach zachowawczych:

Mamy dany układ \(\dot X = F(X), X \in R^2, F \in C^1\), w którym istnieje zachowana wielkość E i punkt stały X0 jest izolowany. Jeśli X0 jest lokalnym minimum U to wszystkie trajektorie wystarczająco blisko X0 są zamknięte.

Wniosek z twierdzenia: nieliniowe środki w układach zachowawczych występują w lokalnych minimach U

  1. Przykład: potencjał podwójnej studni (ang. double-well potential)

    Energia potencjalna układu mechanicznego:
    \[U(x) = -\frac{1}{2}x^2 + \frac{1}{4}x^4\]

    example_652_potential.png

    Obliczanie pola zachowawczego:
    \[F(x) = -\frac{dU}{dx}=x - x^3\]

    Równanie ruchu:
    \[m\ddot x = -\frac{dU}{dx}, m=1\] \[\ddot x - x + x^3 = 0\]

    \[\dot x = v, \dot v = x-x^3\] \[v=0,x=0\] \[v=0,x=1\] \[v=0,x=-1\]

    \[A = \begin{pmatrix}0 & 1\\ 1-3x^2 & 0\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}, det(A)=-1<0, \tau=0, \text{siodło}\] \[A(1,0) = \begin{pmatrix}0 & 1\\ -2 & 0\end{pmatrix}, det(A)=2>0, \tau=0, \text{środek?}\] \[A(-1,0) = \begin{pmatrix}0 & 1\\ -2 & 0\end{pmatrix}, det(A)=2>0, \tau=0, \text{środek?}\]

    Jest to układ zachowawczy, więc wzdłuż każdej trajektorii energia jest stała 79. Jeśli nie byłby to środek to musiałby to być wir stabilny bądź niestabilny. Aby spełniona była zasada zachowania energii to dla każdej trajektorii energia musi być równa tej w punkcie stałym E(x*). Co oznacza, że funkcja U(X) w układzie z wirami byłaby funkcją stałą. Zasada zachowania zachodzi w układach, w których U(X) nie jest funkcją stała. W przeciwnym razie każdy układ byłby zachowawczy.

    example652_phase.png

    Figure 79: Przestrzeń fazowa

    Pytanie: po jakiej trajektorii koszt energetyczny jest najmniejszy?

    Na rysunku79 widzimy trzy punkty stałe: siodło i dwa nieliniowe środki. Orbita, która rozpoczyna się i kończy w tym samym punkcie stałym nazywa się orbitą homokliniczną. Poruszając się po tej orbicie E = 0, ale ruch trwa w nieskończoność!. Orbity homokliniczne często występują w układach zachowawczych; w pozostałych są rzadkością. Ujemny koszt energetyczny jest w punktach środkowych E=-0.25.

    example652_homoclinic.png

    Figure 80: Orbita homokliniczna

    Pytanie: po jakiej trajektorii koszt czasowy jest najmniejszy?

    Pytanie: ile wynosi energia całkowita?:
    \[\ddot x - x + x^3 = 0 \;/ \cdot \dot x\] \[\dot x \ddot x + \dot x \frac{dU(x)}{dx} = 0 \] \[\frac{d}{dt}\biggl(\frac{1}{2}\dot x^2\biggr) + \frac{dx}{dt}\frac{dU}{dx}= 0 \]

    \[\frac{d}{dt}\biggl(\frac{1}{2}\dot x^2 +U(x)\biggr)= 0 \] \[\frac{d}{dt}(E_{k} + E_{p})= 0 \]

    \[E = \frac{1}{2}\dot x^2 -\frac{1}{2}x^2 + \frac{1}{4}x^4\] dla x = 1 i v = 1 \[E = \frac{1}{2} -\frac{1}{2} + \frac{1}{4} = \frac{1}{4}\]

    t=5; x=-1.532816; v=0.298982;
    E = 0.5*v^2-0.5*x^2+0.25*x^4
    
    
    0.24999947524721233
    

    example652_total_energy.png

    Figure 81: Orbita dla x0=1 i v0=1

    example652_energy3d.gif

    Figure 82: Powierzchnia energii

    x'=y
    y'=x-x^3
    aux e=0.5*y^2 -0.5*x^2 + 0.25*x^4
    
    init x=0.3, y=0
    
    @ total=100,axes=3,xp=x,yp=y,zp=e,xlo=-2,ylo=-2,xhi=2,yhi=2,
    @ xmax=2,xmin=-2, ymax=2,ymin=-2, zmax=0.5,zmin=-0.5
    @ maxstor=2000000, bounds=100000, nmesh=800 
    done
    
    # GUI 
    # Initialconds —> Range —> Range Over: Y —> Steps: 249 —> Start: 0 —> End: 0.5
    # Kinescope —> Make AniGif
    
  2. Twierdzenie o nieliniowych środkach i nieizolowane punkty stałe

    \[\dot x = xy, \; \dot y = -x^2\]

    \[\frac{\dot x}{\dot y} = -\frac{y}{x}\] \[xdx = -ydy\] \[\frac{1}{2}x^2 = -\frac{1}{2}y^2 + C\] \[E = x^2 + y^2 = 2C,\;\text{wielkość zachowywana}\]

    Rozwiązanie w punkcie równowagi: \[xy=0 \text{, } -x^2=0\] \[y=0,\;x\ne0\] \[x=0,\;y\ne0\]

    \[A = \begin{pmatrix}y & x\\ -2x & 0\end{pmatrix}\] \[A(0,a) = \begin{pmatrix}a & 0\\ 0 & 0\end{pmatrix}=0, \tau=a,\;\text{nieizolowane punkty stałe}\]

    \[\lambda_{1} = \frac{1}{2}(a-a)=0\] \[\lambda_{2} = \frac{1}{2}(a+a)=a\] Ponieważ w układzie są nieizolowane punkty stałe to nasze twierdzenie nie obowiązuje. Nie możemy powiedzieć, że dla x0 = (0,0), które jest lokalny minimum \(E=x^2+y^2\), wszystkie trajektorie wystarczająco blisko x0 są zamknięte.

    ex6512_conservative.png

    Figure 83: Przestrzeń fazowa

    x'=xy
    y'=-x^2
    
    @ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2
    @ yhi=2, xmax=2,xmin=-2, ymax=2,ymin=-2
    @ maxstor=2000000, bounds=100000,  nmesh=800 
    done
    
    # GUI 
    # Dir.filed/flow —> Direct Field –> 16
    # Initialconds —> New —> 1 ENTER —> 0 ENTER
    # Initialconds —> New —> -1.5 ENTER —> 0 ENTER
    

    Układ, w którym są izolowane ścieki nie jest zachowawczy. W tym układzie mamy nieizolowane ścieki, więc może to być układ zachowawczy.

4.2.7 Układy odwracalne (ang. reversible systems)

Układ \(\dot X = f(X)\) jest odwracalnym, jeśli istnieje takie przekształcenie T, że \(T^2(X)=X\). To znaczy, że jeśli dwukrotnie zastosujemy przekształcenie T to powrócimy do punktu wyjściowego. Układ odwracalny jest niezmienny przy zmianie zmiennych \[t \rightarrow -t,\;X \rightarrow T(X)\] np.

\[t \rightarrow -t\] \[t \rightarrow -t,\;y \rightarrow -y\] \[t \rightarrow -t,\;x \rightarrow -x\] \[t \rightarrow -t,\;x \rightarrow -x,\;y \rightarrow -y\]

Układy mechaniczne są często odwracalne typu \(t \rightarrow -t,\;y \rightarrow -y\). W takich układach dynamika wygląda tak samo niezależnie, czy czas biegnie do przodu, czy do tyłu.

Nietłumiony układ mechaniczny: \[m\frac{d^2x}{dt^2}=F(x)=m\frac{d^2x}{d(-t)^2}\] \[m\ddot x=F(x)\]

\[\frac{dx}{dt}=y\] \[\frac{dy}{dt}=\frac{1}{m}F(x)\]

\[t \rightarrow -t,\;\text{układ będzie poruszał się w kierunku -t}\] \[m\frac{d}{dt}\frac{dx}{dt} = F(x)\] \[m\frac{d}{d(-t)}\frac{dx}{d(-t)} = F(x)\] \[\frac{dx}{d(-t)}=y,\;-\frac{dx}{dt}=y, \dot x = -y\] \[\frac{d}{d(-t)}y=\frac{dy}{d(-t)}=\frac{1}{m}F(x),\;\dot y=-\frac{1}{m}F(x)\]

Jeśli wiem, że układ po zastosowaniu przekształcenia: \(t \rightarrow -t\) zmienia kierunek prędkości, to po dodatkowym podstawieniu: \(y \rightarrow -y\) układ musi pozostać niezmienny

\[t \rightarrow -t,\;y \rightarrow -y,\;\text{układ jest niezmienny}\] \[\frac{dx}{d(-t)}=-y,\; \dot x = y\] \[\frac{d}{d(-t)}(-y)=\frac{d(-y)}{d(-t)}=\frac{1}{m}F(x),\;\dot y=\frac{1}{m}F(x)\]

Jeśli rozwiązaniem są punkty stałe \(X_{0}=(x_{0},y_{0})\), to w takim układzie punktem stałym będzie także \(X_{0}'=(x_{0},-y_{0})\)

Nietłumiony oscylator harmoniczny
\[m\ddot x + kx =0\]

\[\dot x = v\] \[\dot v = -\frac{k}{m}x = -\omega^2 x\]

\(t \rightarrow -t,\;\text{układ będzie poruszał się w kierunku -t}\), \[\frac{dx}{d(-t)} = v,\;\frac{dx}{dt} = -v\]\[ \]

\[\frac{d}{d(-t)} \frac{dx}{d(-t)} = \frac{dv}{d(-t)},\; \dot v= \frac{k}{m}x = \omega^2 x\]

\(t \rightarrow -t,\; v \rightarrow -v,\;\text{układ jest niezmienny}\) \[\frac{dx}{d(-t)} = -v,\;\dot x = v\]

\[\frac{d}{d(-t)} \frac{dx}{d(-t)} = \frac{d(-v)}{d(-t)},\; \dot v= -\frac{k}{m}x = -\omega^2 x\]

W takich układach przestrzeń fazowa jest symetryczna względem osi X z odwróconym kierunkiem pola wektorowego.

reversible_ty.png

Figure 84: Mechaniczne układy odwracalne

par a=-1
par b=1

x'=a*y
y'=b*x

@ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2
@ yhi=2, xmax=2,xmin=-2, ymax=2,ymin=-2
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Makewindow —> Create
# W lewym górny rogu czarny kwadrat oznacza aktywne okno
# Przejdź na drugie okno i zmień paramtery 
# Param —> a=1 i b=-1
# Dir.filed/flow —> Direct Field –> 16

Jeśli układ drugiego rzędu \(\dot x = f(x,y),\;\dot y = g(x,y)\) jest niezmienny po zastosowaniu zmiany zmiennych \(t \rightarrow -t,\;y \rightarrow -y\), to jest odwracalny. Na podstawie parzystości i nieparzystości funkcji f i g możemy to sprawdzić:

  • f jest nieparzysta w y: \(f(x,-y)=-f(x,y)\)
  • g jest parzysta w y: \(g(x,-y) = g(x,y)\)

Na rysunku widzimy, że przestrzeń fazowa nie tylko jest symetryczna względem osi X, ale również względem osi Y. Musi zatem układ być niezmienny po zastosowaniu przekształcenia \(t \rightarrow -t,\;x \rightarrow -x\):

\(t \rightarrow -t,\; x \rightarrow -x\) \[\frac{d(-x)}{d(-t)} = v,\;\dot x = v\] \[\frac{d}{d(-t)} \frac{d(-x)}{d(-t)} = -\omega^2 (-x)\] \[\dot v= -\omega^2 x\]

Twierdzenie o nieliniowych środkach w układach odwracalnych:

Niech \(X_{0}=0\) będzie nieliniowym środkiem w układzie odwracalnym \(\dot X = F(X), X \in R^2, f \in C^1\). Wtedy wszystkie trajektorie wystarczająco blisko \(X_{0}=0\) są zamknięte.

  1. Przykład zastosowania twierdzenia o nieliniowych środkach

    \[\dot x = f(x,y),\; \dot y = g(x,y)\] \[\dot x = y-y^3,\; \dot y = -x-y^2\] \[y(1-y^2)=0, x=-y^2\] \[x=0,y=0\] \[x=-1,y=1\] \[x=-1,y=-1\] Na podstawie obliczonych punktów stałych możemy już przewidzieć, czy może być to układ odwracalny. W naszym przypadku widzimy symetrie wzdłuż osi X.

    \[A = \begin{pmatrix}0 & 1-3y^2\\ -1 & -2y\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & 1\\ -1 & 0\end{pmatrix}, det(A)=1>0, \tau=0, \text{środek?}\] \[A(-1,1) = \begin{pmatrix}0 & -2\\ -1 & -2\end{pmatrix}, det(A)=-2<0, \tau=0, \text{siodło}\] \[A(-1,-1) = \begin{pmatrix}0 & -2\\ -1 & 2\end{pmatrix}, det(A)=-2<0, \tau=0, \text{siodło}\]

    \[\dot x = y-y^3, \dot y = -x-y^2\]

    \[t \rightarrow -t\]: \[\frac{dx}{d(-t)} = y-y^3, \frac{dy}{d(-t)} = -x-y^2\] \[\dot x = -y+y^3, \dot y = x+y^2\] \[y(-1+y^2)=0, x=-y^2\] \[x=0,y=0\] \[x=-1,y=1\] \[x=-1,y=-1\]

    To jest układ odwracalny, gdyż po zastosowaniu przekształcenia \(t \rightarrow -t\), \(y \rightarrow -y\): układ pozostał taki sam:

    \[\frac{dx}{d(-t)} = -y+y^3, \frac{d(-y)}{d(-t)} = -x-y^2\] \[\dot x = y-y^3, \dot y = -x-y^2\]

    To samo możemy sprawdzić przez podstawienia: \(f(x,-y) = -f(x,y),\;g(x,-y) = g(x,y)\)

    Na podstawie tego faktu można wnioskować, że nasz liniowy środek jest nieliniowym środkiem. Istnienie siodła w punkcie X0 = (-1,1) pociąga w takim układzie istnienie drugiego siodła w punkcie X0'= (-1,-1).

    example661_reversible.png

    Figure 85: Przestrzeń fazowa

    x'=y-y^3
    y'=-x-y^2
    
    @ total=10,xp=x,yp=y,xlo=-2,ylo=-2,xhi=2,yhi=2
    @ xmax=2,xmin=-2, ymax=2,ymin=-2
    @ maxstor=2000000, bounds=100000, nmesh=800 
    done
    
    # GUI 
    # Dir.filed/flow —> Direct Field –> 16
    # Dir.filed/flow —> Flow –> 8
    # Sing pts —> monteCar —> Shoot 1; pozostałe domyślne ustawienia 
    

    Trajektorie łączące punkty siodło są nazwane trajektoriami heteroklinicznymi (ang. heteroclinic trajectories) lub połączeniami siodłowymi (ang. saddle connections). Występują często w układach odwracalnych; w pozostałych są rzadkością.

  2. Przykład na wykorzystanie argumentu odwracalności układu w celu wydedukowania istnienia orbity homoklinicznej

    \[\dot x = y\] \[\dot y = x-x^2\]

    To jest układ odwracalny, gdyż:
    \(f(x,-y) = -f(x,y),\;g(x,-y) = g(x,y)\)

    \[y=0, x(1-x)=0\] \[x=0,y=0\] \[x=1,y=0\]

    \[A = \begin{pmatrix}0 & 1\\1-2x & 0\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}, det(A)=-1<0, \tau=0, \text{siodło}\] \[A(1,0) = \begin{pmatrix}0 & 1\\ -1 & 0\end{pmatrix}, det(A)=1>0, \tau=0, \text{środek}\]

    Ponieważ jest to układ odwracalny to liniowy środek jest nieliniowym środkiem. Układ jest symetryczny względem osi x. W punkcie (0,0) mamy siodło a w punkcie (1,0) środek.

    Pytanie: czy jest tu orbita homokliniczna?
    Na podstawie obliczeń pola wektorowego dla rozmaitości stabilnej, bądź niestabilnej naszkicuj jej trajektorie. Z jej analizy zorientujesz się, że dla x<1 \(\dot y > 0\); dla y>0 \(\dot x > 0\), Gdy x=1 to \(\dot y = 0\) i dla x>0 będzie ujemne, aż y osiągnie zero. Powstały łuk krzywej ma swoje odbicie lustrzane względem osi x, gdyż jest to układ odwracalny. Ponieważ pokazaliśmy, że w układzie istnieje taka orbita, która zaczyna się i kończy w tym samym punkcie, to tym samym udowodniliśmy istnienie orbity homoklinicznej.

    example662_homoclinic.png

    Figure 86: Przestrzeń fazowa

    x'=y
    y'=x-x^2
    
    @ total=10,xp=x,yp=y,xlo=-2,ylo=-2,xhi=2,yhi=2
    @ xmax=2,xmin=-2, ymax=2,ymin=-2
    @ maxstor=2000000, bounds=100000, nmesh=800 
    done
    
    # GUI 
    # Dir.filed/flow —> Direct Field –> 16
    # Dir.filed/flow —> Flow –> 8
    
  3. Różnica pomiędzy odwracalnym a zachowawczym

    \[\dot x = -2cos(x)-cos(y)\] \[\dot y = -2cos(y)-cos(x)\]

    Pytanie: czy jest to układ zachowawczy?
    Jeśli istnieje ściek to możemy to wykluczyć:

    \[y_{0}=\pm\frac{\pi}{2}, x_{0}=\pm\frac{\pi}{2}\]

    \[A = \begin{pmatrix}2sin(x) & sin(y)\\sin(x) & 2sin(y)\end{pmatrix}\] \[A\biggl(\frac{\pi}{2},\frac{\pi}{2}\biggr) = \begin{pmatrix}2 & 1\\ 1 & 2\end{pmatrix}, det(A)=3>0, \tau=4, \text{niestabilny wir}\]

    \[A\biggl(-\frac{\pi}{2},-\frac{\pi}{2}\biggr) = \begin{pmatrix}-2 & -1\\ -1 & -2\end{pmatrix}, det(A)=3>0, \tau=-4, \text{stabilny wir}\] To nie jest układ zachowawczy, gdyż stabilny wir jest ściekiem w sensie topologicznym. Pozostałe punkty stałe to siodła. \[A\biggl(-\frac{\pi}{2},\frac{\pi}{2}\biggr) = \begin{pmatrix}-2 & 1\\ -1 & 2\end{pmatrix}, det(A)=-3<0, \tau=0, \text{siodło}\]

    \[A\biggl(\frac{\pi}{2},-\frac{\pi}{2}\biggr) = \begin{pmatrix}2 & -1\\ 1 & -2\end{pmatrix}, det(A)=-3<0, \tau=0, \text{siodło}\]

    Pytanie: czy jest to układ odwracalny?

    To jest układ odwracalny, po zastosowaniu przekształcenia: \[t \rightarrow -t,\;y \rightarrow -y,\;x \rightarrow -x\], co znaczy, że jest symetryczny względem przekątnej XY, ale nie względem osi X lub Y

    \(t \rightarrow -t,\; y \rightarrow -y\): \[\frac{dx}{d(-t)} = -2cos(x)-cos(-y)\] \[\frac{d(-y)}{d(-t)} = -2cos(-y)-cos(x)\]

    \[\frac{dx}{dt} = 2cos(x)+cos(y)\] \[\frac{dy}{dt} = -2cos(y)-cos(x)\]

    Tak samo będzie z przekształceniem \(t \rightarrow -t,\;x \rightarrow -x\):

    W takim układzie jeśli mamy punkt stały \(X_{0} = (\frac{\pi}{2},\frac{\pi}{2})\) to musi istnieć \(X_{0}' = (-\frac{\pi}{2},-\frac{\pi}{2})\) z odwróconym kierunkiem prędkości. To znaczy jeśli X0 jest niestabilnym wirem to X0' jest stabilnym wirem. To samo dotyczy dwóch pozostałych punktów stałych.

    example663_reversible.png

    Figure 87: Przestrzeń fazowa

    x'=-2*cos(x)-cos(y)
    y'=-2*cos(y)-cos(x)
    
    @ total=10,xp=x,yp=y, xlo=-3,ylo=-3, xhi=3
    @ yhi=3, xmax=3,xmin=-3, ymax=3,ymin=-3
    @ maxstor=2000000, bounds=100000,  nmesh=800 
    done
    
    # GUI 
    # Dir.filed/flow —> Direct Field –> 16
    # Sing Pts —> monteCar —-> można zostawić domyśle ustawienia z jednym wyjątiem: Shoot=1
    
  4. Kolejny przykład

    \[\dot x = y(1-x^2),\; \dot y = 1 - y^2\]

    \[t \rightarrow -t,\;y \rightarrow -y\]

    \[\frac{dx}{d(-t)} = -y(1-x^2),\; \frac{d(-y)}{d(-t)} = 1 - (-y)^2\]

    \[\frac{dx}{dt} = y(1-x^2),\; \frac{dy}{dt} = 1 - y^2\]

    To jest układ odwracalny; jest niezmienny po zastosowaniu przekształcenia:
    \[t \rightarrow -t,\;y \rightarrow -y\] \[f(x,-y) = -f(x,y),\;g(x,-y) = g(x,y)\]

    Nie jest on niezmienny po zastosowaniu przekształcenia: \(t \rightarrow -t,\;x \rightarrow -x\):
    To znaczy nie jest on symetryczny względem osi Y.

    \[\dot x = y(1-x^2),\; \dot y = 1 - y^2\]

    \[t \rightarrow -t,\;x \rightarrow -x\]:

    \[\frac{d(-x)}{d(-t)} = y(1-(-x)^2),\; \frac{dy}{d(-t)} = 1 - y^2\]

    \[\frac{dx}{dt} = y(1-x^2),\; \frac{dy}{dt} = -1 + y^2\]

    \[f(-x,y) = f(x,y),\;g(-x,y) = g(x,y)\]

    ex661.png

    Figure 88: Przestrzeń fazowa

    x'=y*(1-x^2)
    y'=1-y^2
    
    @ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2
    @ yhi=2, xmax=2,xmin=-2, ymax=2,ymin=-2
    @ maxstor=2000000, bounds=100000,  nmesh=800 
    done
    
    # GUI 
    # Dir.filed/flow —> Direct Field –> 16
    # Sing Pts —> monteCar —-> można zostawić domyśle ustawienia z jednym wyjątiem: Shoot=1
    

    ex661_t.png

    Figure 89: Przestrzeń fazowa przed i po przekształceniu t —> -t

    Pytanie: co się stanie z dynamiką układu jeśli usuniemy nieliniowość z funkcji g? \[\dot x = y(1-x^2),\; \dot y = 1 - y\]

4.2.8 Przykład: wahadło

\[\ddot \theta + \frac{g}{L}sin\theta=0\] \[\dot \theta = \omega\] \[\dot \omega = - \frac{g}{L}sin\theta\]

  • g - siła ciężkości
  • L - długość wahadła
  • Rozwiązanie w punktach równowagi \[\omega = 0,\; \theta=0\] \[\omega = 0,\; \theta=\pi\]
  • Linearyzacja wokół punktów równowagi \[A = \begin{pmatrix}0 & 1\\-\frac{g}{L}cos\theta & 0\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & 1\\ -\frac{g}{L} & 0\end{pmatrix}, det(A)=\frac{g}L{}>0, \tau=0, \text{środek?}\] \[A(\pi,0) = \begin{pmatrix}0 & 1\\ \frac{g}{L} & 0\end{pmatrix}, det(A)=-\frac{g}L{}<0, \tau=0, \text{siodło}\] To jest układ zachowawczy, gdyż istnieje funkcja potencjału: \[- \frac{g}{L}sin\theta = -\frac{dU}{d\theta}\] \(U=-\frac{g}{L}cos\theta\) Wielkość zachowywana to:
    \[\dot \theta \ddot \theta + \dot \theta \frac{g}{L}sin\theta=0\] \[\frac{d}{dt}\biggl(\frac{1}{2}\dot \theta^2 - \frac{g}{L}cos\theta\biggr)=0\] \[E=\frac{1}{2}\dot \theta^2 - \frac{g}{L}cos\theta\]

Punkt X0 = (0,0) jest lokalnym minimum U i zgodnie z naszym twierdzeniem wszystkie trajektorie wystarczająco blisko tego punktu są zamknięte. Punkt X0=(0,0) jest nieliniowym środkiem.

Ponieważ jest to nieliniowy środek w X0=(0,0), to układ musi być odwracalny

\[t \rightarrow -t,\;\omega \rightarrow -\omega\] \[\frac{d\theta}{d(-t)} = -\omega,\; \dot \theta = \omega\] \[\frac{d}{d(-t)}\frac{d\theta}{d(-t)} = -\frac{g}{L}sin\theta\] \[\dot \omega = -\frac{g}{L}sin\theta\]

\[t \rightarrow -t,\;\theta \rightarrow -\theta\] \[\frac{d(-\theta)}{d(-t)} = \omega,\; \dot \theta = \omega\] \[\frac{d}{d(-t)}\frac{d(-\theta)}{d(-t)} = -\frac{g}{L}sin(-\theta)\] \[\dot \omega = -\frac{g}{L}sin\theta\]

example67_pendulum.png

Figure 90: Przestrzeń fazowa dla L=1 (lewy) i L=10 (prawy)

par g=9.81
par L=1

theta'=omega
omega'=-g/L*sin(theta)

aux e=0.5*omega^2-g/L*cos(theta)

@ total=100,axes=2,xp=theta,yp=omega,xlo=-6.28,ylo=-6.28,xhi=6.28,yhi=6.28,
@ xmax=6.28,xmin=-6.28, ymax=6.28,ymin=-6.28
@ maxstor=2000000, bounds=100000, nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 32
# Dir.filed/flow —> Flow –> 16
# Makewindow —> Create
# Wybierz jedno okno
# Sing Pts —> monteCar —-> można zostawić domyśle ustawienia z jednym wyjątiem: Shoot=1
# Wybierz drugie okno
# Param —> L=10
# Dir.filed/flow —> Direct Field –> 32
# Dir.filed/flow —> Flow –> 16
# Sing Pts —> monteCar —-> można zostawić domyśle ustawienia z jednym wyjątiem: Shoot=1

example67_energy_surface.png

Figure 91: Powierzchnia energii dla L=1

par g=9.81
par L=1

theta'=omega
omega'=-g/L*sin(theta)

aux e=0.5*omega^2-g/L*cos(theta)

@ total=100,axes=3,xp=theta,yp=omega,zp=e,xlo=-6.28,ylo=-6,xhi=6.28,yhi=6,
@ xmax=6.28,xmin=-6.28, ymax=6,ymin=-6, zmax=15,zmin=-15
@ maxstor=2000000, bounds=100000, nmesh=800 
done

# GUI 
# Initialconds —> Range —> Range Over: THETA —> Steps: 100 —> Start: -4 —> End: 4 —> Movie: No

Z wykresu 3D powierzchni energii wyraźnie widać, że im wyższa energia tym większa orbita.

Pytanie: czy w układzie jest orbita heterokliniczna czy homokliniczna?

4.3 Globalne techniki nieliniowe

4.3.1 Teoria indeksu

Topologiczna metoda dostarczająca globalnych informacji o przestrzeni fazowej układu. Tę metodę możemy wykorzystać do zbadania, czy istnieją izolowane punkty stałe lub orbita zamknięta.

Metoda polega na naszkicowaniu krzywej zamkniętej w przestrzeni fazowej tak, aby nie przechodziła przez punkt stały. Nanosimy krzywą zamkniętą w tym miejscu gdzie chcemy zbadać istnienie punktu stałego, bądź orbity zamkniętej. Po naszkicowaniu takiej krzywej obliczamy wektory w wybranych miejscach na krzywej. Sumujemy kąty nachylenia wektorów do osi X. Kąt jest dostani w stronę przeciwną do ruchu wskazówek zegara. Indeks krzywej liczymy ze wzoru:

\[I_{C} = \frac{1}{2\pi}[\phi_{C}]\]

Jeśli \(I_{C}=0\) to krzywa zamknięta C nie okrążą żadnego punktu stałego. \(I_{C}=-1\) dla siodła, a dla wszystkich pozostałych izolowanych punktów stałych \(I_{C}=1\)

index_theory.png

Figure 92: Zoastosowanie metodu dla dwóch ścieków i siodła

Twierdzenie:

Jeśli krzywa zamknięta C okrąża n izolowanych punktów stałych \(X_{0}^1...X_{0}^n\) to: \(I_{C} = I_{1} + I_{2}+...+I_{n}\)

Twierdzenie:

Każda orbita zamknięta w przestrzeni fazowej musi zawierać w sobie punkty stałe o sumarycznym indeksie równym 1

  1. Przykład: układ z potencjałem podwójnej studni

    \[\dot x = y, \dot y = x-x^3\]

    Są trzy punkty stałe:
    Siodło: \(x=0,y=0\)
    Środek 1: \(x=1,y=0\)
    Środek 2: \(x=-1,y=0\)

    Orbita zawierająca te trzy punkty może istnieć w przestrzeni fazowej układu, gdyż suma indeksów jest równa 1
    Siodło: IC1 = -1
    Środek 1: IC2 = 1
    Środek 2: IC3 = 1

    IC = -1 + 1 + 1 = 1

    W celu weryfikacji spójrz na rysunek 79

  2. Przykład:

    \[x'=y,\;y'=x-x^2\]

    Są dwa punkty stałe:
    Siodło: \(x=0,y=0\)
    Środek: \(x=1,y=0\)

    Orbita zawierająca te dwa punkty nie może istnieć w przestrzeni fazowej układu, gdyż suma indeksów jest równa 0: Siodło: IC1 = -1
    Środek 1: IC2 = 1

    IC = -1 + 1 = 0 (rysunek 86)

  3. Przykład: istnieje orbita zamknięta?

    \[\dot x = x(4-y-x^2),\;\dot y=y(x-1)\]

    Są cztery punkty stałe:
    \[x=0, y=0\] \[x=1, y=3\] \[x=2, y=0\] \[x=-2, y=0\]

    \[A = \begin{pmatrix}4-y-3x^2 & -x\\y & x-1\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}4 & 0\\ 0 & -1\end{pmatrix}, det(A)=-4, \tau=3, \text{siodło}\] \[A(1,3) = \begin{pmatrix}-2 & -1\\ 3 & 0\end{pmatrix}, det(A)=3, \tau=-2, \text{wir stabilny}\] \[A(2,0) = \begin{pmatrix}-8 & -2\\ 0 & 1\end{pmatrix}, det(A)=-8, \tau=-7, \text{siodło}\] \[A(-2,0) = \begin{pmatrix}-8 & 2\\ 0 & -3\end{pmatrix}, det(A)=24, \tau=-11, \text{ściek}\]

    Wniosek: całkowity indeks układu wynosi 0, więc jeśli istnieje orbita zamknięta w układzie to nie zawiera w sobie wszystkich punktów stałych. Orbita zamknięta mogłaby jedynie zaistnieć w okolicach ścieku bądź wiru stabilnego. W pierwszym przypadku w pobliży mamy siodło o rozmaitościach w osiach X i Y. Orbita zamknięta nie może tam istnieć, gdyż wtedy przecinałaby rozmaitość. W drugim przypadku jest podobnie: wir stabilny jest połączony z rozmaitością drugiego siodła. Orbita nie może istnieć z tego samego powodu.

4.3.2 Izokliny zerowego wzrostu

Krzywe w przestrzeni fazowej, które spełniają równania:

\[\dot x =0\] \[\dot y = 0\]

nazywamy izoklinami zerowego wzrostu (ang. nullclines, zero-growth isoclines). Dzielą przestrzeń fazową na jakościowo różne obszary. Punkty przecięcia się krzywych są punktami stałymi.

ex613_nullclines.png

Figure 93: Izokliny zerowego wzrostu

ex613_nullclines_xpp.png

Figure 94: Izokliny zerowego wzrostu - XPPAut

x'=x*(x-y)
y'=y*(2*x-y)

@ total=10,xp=x,yp=y, xlo=-3,ylo=-3, xhi=4
@ yhi=4, maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Nullcline —> New

Izokliny zerowego wzrostu mogą mieć złożone kształty:

example_nullcline.png

Figure 95: Izokliny zerowego wzrostu - XPPAut

x'=-2*cos(x)-cos(y)
y'=-2*cos(y)-cos(x)

@ total=10,xp=x,yp=y, xlo=-3,ylo=-3, xhi=3
@ yhi=3, xmax=3,xmin=-3, ymax=3,ymin=-3
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Nullcline —> New

example_nullcline2.png

Figure 96: Izokliny zerowego wzrostu - XPPAut

x'=y
y'=x-x^2

@ total=10,xp=x,yp=y,xlo=-2,ylo=-2,xhi=2,yhi=2
@ xmax=2,xmin=-2, ymax=2,ymin=-2
@ maxstor=2000000, bounds=100000, nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Nullcline —> New

4.4 Orbity zamknięte i zbiory graniczne

Cykle graniczne są to izolowane orbity zamknięte i występują tylko w układach nieliniowych. W układach liniowych mogą wystąpić punkty środkowe. Poznamy kilka sposobów na zbadanie, czy w przestrzeni fazowej układu znajduje się cykl graniczny - układy gradientowe, twierdzenie Poincare–Bendixsona

Pierwszych kilka przykładów, jakie będziemy omawiać, dotyczy jakościowej analizy tzn. odpowiedzi na pytanie, czy w danym dwu-wymiarowym układzie istnieją rozwiązania periodyczne. Poznamy kryteria, dzięki którym będziemy mogli odrzucić hipotezę istnienia takich orbit oraz metody pozwalające na stwierdzenie, że rzeczywiście takie orbity istnieją w przestrzeni fazowej układu. Kolejne przykłady dotyczą ilościowej analizy tzn. mając w układzie orbitę zamknięta jaki jest jej kształt i okres.

4.4.1 Układ we współrzędnych polarnych

\[\dot r = r(1-r^2)\] \[\dot \theta = 1\]

Układ o zmiennych niesprzężonych dla \(r \ge 0\).

  1. Punkty stałe: \(r_{0}=0,\; r_{0}' = 1\)

exampe_limit_polar.png

Figure 97: Przestrzeń fazowa

r'=r*(1-r^2)
theta'=1

@ total=10,xp=r,yp=theta, xlo=0,ylo=0, xhi=2 yhi=6.28318
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Phasespace —> Choose —> Period: Default ENTER Theta
# Dir.filed/flow —> Direct Field –> 16
# Dir.filed/flow —> Flow –> 8

exampe_limit_polar3.png

Figure 98: Przestrzeń fazowa

r'=0.5*(r-r^3)
theta'=1

aux x=r*cos(theta)
aux y=r*sin(theta)

@ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2 yhi=2
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Initialconds —> Range —-> Steps=100, Start=-2, Stop=2 —> OK

Ten sam układ we współrzędnych kartezjańskich:
\[\dot x=x(1-x^2-y^2) - y\] \[\dot y=y(1-x^2-y^2) + x\]

exampe_limit_polar2.png

Figure 99: Przestrzeń fazowa

x'=x*(1-x^2-y^2) - y
y'=y*(1-x^2-y^2) + x 

@ total=10,xp=x,yp=y, xlo=-3,ylo=-3, xhi=3 yhi=3
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Dir.filed/flow —> Flow –> 8

4.4.2 Oscylator van der Pola

\[\ddot x + \mu(x^2-1)\dot x + x=0\]

\[\dot x = v\] \[\dot v = -\mu(x^2-1)v - x\]

Pytanie: jak wpływa nieliniowy element na tłumienie układu?

  1. Punkty stałe: \(v_{0}=0,\; x_{0}=0\)
  2. Linearyzacja:
    \[A = \begin{pmatrix}0 & 1\\-2\mu vx - 1 & -\mu(x^2-1)\end{pmatrix}\] \[A(0,0) = \begin{pmatrix}0 & 1\\-1 & \mu \end{pmatrix}, \; det(A)=1>0, \tau=\mu\] \[\mu = 0, \; \text{nietłumiony oscylator harmoniczny}\] \[\mu > 0, \; \text{niestabilny wir}\] \[\mu < 0, \; \text{stabilny wir}\]

Pytanie: czy jest możliwe dla \(\mu \ne 0\), aby w przestrzeni fazowej układu pojawiła się orbita zamknięta? Jeśli tak, to czy będzie ona neutralnie stabilna, asymptotycznie stabilna bądź asymptotycznie niestabilna?

4.4.3 Układy Liénarda

\[\ddot x + f(x)\dot x + g(x)=0\]

Twierdzenie Liénarda:

Jeśli:

  • funkcje f(x) i g(x) są C1 w całej dziedzinie
  • g(-x) = -g(x) w całej dziedzinie
  • g(x) > 0 dla x > 0
  • f(-x) = f(x) w całej dziedzinie
  • \(F(x) = \int_0^x f(u)du\)
    • ma dokładnie jedno dodatnie miejsce zerowe x=a,
    • jest ujemna dla 0 < x < a,
    • dodatnia i niemalejąca dla x > a
    • oraz \(F(x) \rightarrow \infty\) gdy \(x \rightarrow \infty\)

wtedy układ Liénarda ma unikatowy i stabilny cykl graniczny otaczający początek układu współrzędnych w przestrzeni fazowej

  1. Przykład: oscylator van der Pola

    \[\dot x = v\] \[\dot v = -\mu(x^2-1)v - x\] \[\ddot x + \mu(x^2-1)\dot x + x=0\]

    \[f(x)=\mu(x^2-1), \; g(x)=x, C^1\] \[g(-x)=-g(x),\;g(x)>0, dla\; x>0,\;f(-x) = f(x)\] \[F(x)=\int_0^x \mu(x^2-1) = \frac{1}{3}\mu x^3 - \mu x\]

    \[\mu x(\frac{1}{3} x^2 - 1)=0\], \[x_{0}=0, x_{0}'= \sqrt{3},x_{0}''= -\sqrt{3}\]

    • dokładnie jedno dodatnie miejsce zerowe \(x=\sqrt{3}\)
    • F(x) < 0 dla \(0 < x < \sqrt{3}\)
    • dodatnie i niemalejące dla x > \sqrt{3}
    • oraz \(F(x) \rightarrow \infty\) gdy \(x \rightarrow \infty\)

    example_van_der_pol_lienard.png

    W przestrzeni fazowej oscylatora van der Pola znajduje się unikatowy i stabilny cykl graniczny otaczający początek układu współrzędnych.

4.4.4 Układy gradientowe

W celu zbadania, czy dany układ ma rozwiązanie, które jest orbitą zamknięta można np. wykazać, że jest układem gradientowym z dowolną funkcją potencjału (na podobieństwa układów mechanicznych z energią potencjalną)

Jeśli mamy funkcję odpowiednio gładką o wartościach skalarnych \(V: R^n \rightarrow R\) taką że:
\(\nabla V = -\Biggl(\frac{\partial V}{\partial x_{1}}(X),\frac{\partial V}{\partial x_{2}}(X),...,\frac{\partial V}{\partial x_{n}}(X)\Biggr)\)

To układ gradientowy na zbiorze Rn jest to układ równań różniczkowych postaci:
\[\dot X = -\nabla V(X)\]

np. dwuwymiarowy układ gradientowy:
\[\dot x = -\frac{\partial V(x,y)}{\partial x}\] \[\dot y = -\frac{\partial V(x,y)}{\partial y}\]

Większość układów dwuwymiarowych to nie są układy gradientowe. Natomiast wszystkie jednowymiarowe są, gdyż można je przedstawić za pomocą potencjału.

  1. Przykład: układ gradientowy

    \[V(x, y) = x^2(x-1)^2 + y^2\]

    \[\dot x = -4x^3+6x^2-2x=-2x(2x^2-3x+1)\] \[\dot x = -2x(x-1)(x-\frac{1}{2})\] \[\dot y = 2y\]

    gradient_saddle_sink_ps.png

    Figure 101: Przestrzeń fazowa

    gradient_saddle_sink.png

    Figure 102: Płaszczyzna potencjału

  2. Przykład: czy w układzie są możliwe oscylacje?

    \[\dot x = siny,\; \dot y = xcosy\]

    Jeśli funkcja potencjału istnieje w tym układzie to musi być funkcją odpowiednio gładką tzn. jej pochodne cząstkowe istnieją i są ciągłe. Pochodne cząstkowe: \[\frac{\partial V(x,y)}{\partial x} = siny\] \[\frac{\partial V(x,y)}{\partial y} = xcosy\]

    Ponadto rozwiązanie w układzie dynamicznym jest unikatowe i jednoznaczne jeśli funkcje opisujące dynamikę są odpowiednio gładkie tzn. pochodna cząstkowe z \(siny\) i \(xcosy\) istnieją i są ciągłe. Dostajemy wtedy pochodne drugiego rzędu. Macierz kwadratowa pochodnych cząstkowych drugiego rzędu jest nazywana Hesjanem (macierzą Hessego).

    \[H = \begin{pmatrix} \frac{\partial^2 V(x,y)}{\partial x^2} & \frac{\partial^2 V(x,y)}{\partial y \partial x}\\ \frac{\partial^2 V(x,y)}{\partial x \partial y} & \frac{\partial^2 V(x,y)}{\partial y^2}\end{pmatrix}\]

    \[H = \begin{pmatrix} 0 & cosy\\ cosy & -xsiny\end{pmatrix}\]

    Twierdzeniem Schwartza:

    Jeśli dla funkcji \(V:\; R^n \rightarrow R\) pochodne cząstkowe drugiego rzędu istnieją i są ciągłe to pochodne cząstkowe mieszane są sobie równe.

    \[\frac{ \partial^2 V(x,y)}{\partial y \partial x} = \frac{ \partial^2 V(x,y)}{\partial x \partial y}\]

    Gdy twierdzenie Schwartza zachodzi to Hesjan jest symetryczny tzn. \(H = H^T\)

    Hesjan w tym przypadku jest wynikiem linearyzacji wokół punktów stałych. Dla układów gradientowych jest on postaci: \[H = \begin{pmatrix} a & b\\ b & c\end{pmatrix}\] \[\lambda_{1,2} = \frac{(a+c) \pm \sqrt{(a+c)^2 - 4(ac-b^2)}}{2}\] \[\lambda_{1,2} = \frac{(a+c) \pm \sqrt{a^2+2ac+c^2 - 4ac+4b^2}}{2}\] \[\lambda_{1,2} = \frac{(a+c) \pm \sqrt{(a -c)^2 + 4b^2}}{2}\]

    Wartości własne \(\lambda_{1,2}\) macierzy symetrycznej są liczbami rzeczywistymi! Pokazaliśmy, że w układzie gradientowym nie mogą wystąpić oscylacje 72.

    Funkcja potencjału: V(x,y) możemy obliczyć całkując poszczególne pochodne cząstkowe:
    \[\dot x = \frac{\partial V(x,y)}{\partial x} = siny\] \[\dot y = \frac{\partial V(x,y)}{\partial y} = xcosy\]

    \[\int \frac{\partial V(x,y)}{\partial x} = \int siny\]

    \[\int \frac{\partial V(x,y)}{\partial y} = \int xcosy\]

    \[\int dV(x,y) = \int siny\, dx = xsiny\]

    \[\int dV(x,y) = \int xcosy\, dy = xsiny\]

    \[V(x,y)= xsiny\]

    example721_gradient2.png

    Figure 103: Przestrzeń fazowa

    gradient_sincos.png

    Figure 104: Płaszczyzna potencjału

    Pytanie: czy istnieje taka przestrzeń fazowa tego układu, w której znajdziemy oscylacje?

  3. Przykład: znajdź funkcję potencjału

    \(\dot x = y^2 + ycosx, \; \dot y = 2xy + sinx\)

    \(\nabla V = -\Biggl(\frac{\partial V}{\partial x_{1}}(X),\frac{\partial V}{\partial x_{2}}(X),...,\frac{\partial V}{\partial x_{n}}(X)\Biggr)\) Przez całkowanie cząstkowe znajdziemy funkcję potencjału:

    \(\int_x \dot x = y^2x + ysinx\) \(\int_y \dot y = xy^2+ysinx\)

    \(V(x,y) = xy^2 + ysinx\)

4.4.5 Układy hamiltonowskie

Jeśli mamy funkcję odpowiednio gładką o wartościach skalarnych \(H: R^n \rightarrow R\) taką że:
\[\dot x = \frac{\partial H}{\partial y}(x,y)\] \[\dot y = -\frac{\partial H}{\partial x}(x,y)\]

To H jest funkcją hamiltonowską, a układ nazywamy hamiltonowskim rzędu 2.

\[ H = \int_y \dot x + \int_x \dot y\]

\[\dot H = \frac{d}{dt}\biggl(\int_y\frac{\partial H}{\partial y}(x,y) - \int_x \frac{\partial H}{\partial x}(x,y)\biggr)\]

\[\dot H = \frac{d}{dt}\biggl(\int\frac{dH}{dy}dy - \int \frac{dH}{dx}dx\biggr)\]

\[\dot H = \frac{d}{dt}\biggl(\int dH - \int dH\biggr)=0\]

Dla układów hamiltonowskich w przestrzeni R2 funkcja H nie zmienia się. Funkcja H jest nazywaną całką pierwszą lub całką ruchu, czyli wielkością fizyczną, która jest stała podczas ruchu.

Jeśli X = (x*,y*) są punktami stałymi w układzie hamiltonowskim to wartości własne zlinearyzowanego układu wynoszą \(\pm \lambda\) lub \(\pm i\lambda\)

  1. Przykład: nietłumiony oscylator harmoniczny

    \[\dot x=y,\; \dot y=-kx\]

    \[\int_y \dot x = \frac{1}{2}y^2\] \[\int_x -\dot y = \frac{1}{2}kx^2\]

    \[H = \frac{1}{2}y^2 +\frac{1}{2}kx^2\]

    \[\dot H = \frac{d}{dt}\biggl(\frac{1}{2}y^2 +\frac{1}{2}kx^2\biggr)\] \[\dot H = \frac{1}{2}\biggl(\frac{d(y^2)}{dy}\frac{dy}{dt} +\frac{d(kx^2)}{dx}\frac{dx}{dt}\biggr)\] \[\dot H = \frac{1}{2}\biggl(2y\frac{dy}{dt} + 2kx\frac{dx}{dt}\biggr)\] \[\dot H = y\dot y + kx \dot x\] \[\dot H = y(-kx) + kx y=0\]

  2. Przykład: wahadło

    \[\dot \theta = \omega\] \[\dot \omega = - \frac{g}{L}sin\theta\]

    \[\int_{\omega} \dot \theta = \frac{1}{2}\omega^2\] \[\int_{\theta} -\dot \omega = -\frac{g}{L}cos\theta\]

    \[H = \frac{1}{2}\omega^2 - \frac{g}{L}cos\theta\]

    \[\dot H = \frac{d}{dt}\biggl(\frac{1}{2}\omega^2 -\frac{g}{L}cos\theta\biggr)\] \[\dot H = \frac{1}{2}\frac{d(\omega^2)}{d\omega}\frac{d\omega}{dt} - \frac{g}{L}\frac{d(cos\theta)}{d\theta}\frac{d\theta}{dt}\] \[\dot H = \omega\frac{d\omega}{dt} + \frac{g}{L}sin\theta\frac{d\theta}{dt}\] \[\dot H = \omega(- \frac{g}{L}sin\theta) + \frac{g}{L}sin\theta\omega=0\]

4.4.6 Funkcje Lapunowa

Funkcja Lapunowa jest C1 o wartościach rzeczywistych posiadająca cechy:

  • \(V(X) > 0\) w całej dziedzinie oprócz punktów stałych (dla V(X*) = 0
  • \(\dot V < 0\) w całej dziedzinie oprócz punktów stałych - zmiana

W takim układzie X* jest globalnie i asymptotycznie stabilny tzn. dla \(t \rightarrow \infty,\; X(t) \rightarrow X_{0}\)

Jeśli układ posiada taką funkcję to możemy wykluczyć istnienie orbit zamkniętych. Ponieważ V(X) > 0 to do zbudowania takiej funkcji może wykorzystać sumę kwadratów.

Pytanie: czy w tym układzie jest orbita zamknięta? \[\dot x = -x +4y,\; \dot y = -x - y^3\]

Sprawdzimy, czy istnieje funkcją Lapunowa, która określimy za pomocą sumy kwadratów:
\(V(x,y) = ax^2 + by^2\), w całej dziedzinie jest dodatnia.

\[\dot V = 2ax\dot x + 2by\dot y = 2ax(-x +4y) +2by(-x - y^3)\] \[\dot V = -2ax^2 + 8axy - 2bxy - 2by^4\]

Ponieważ \(\dot V > 0\) to dobierzmy a=1 i b=4:
\[\dot V = -2x^2 + 8xy - 8xy - 8y^4 = -2x^2 - 8y^4\]

lapunov_function.png

Figure 105: Pochodna funkcji Lapunowa

4.4.7 Twierdzenie Poincarégo-Bendixona

Twierdzenie Poincarégo−Bendixona określa jakie warunki muszą być spełnione, aby w układzie istniała orbita zamknięta.

Twierdzenie Poincarégo-Bendixona:

Jeśli:

  • zbiór R jest zamkniętym i ograniczonym podzbiorem w przestrzeni fazowej
  • pole wektorowe w R jest C1
  • R nie zawiera żadnego punktu stałego
  • w R jest trajektoria C, która zaczyna się i kończy w R

to wtedy C jest orbitą zamkniętą lub spiralą, która dąży do orbity zamkniętej gdy \(t \rightarrow \infty\)

1.png

Figure 106: Twierdzenie Poincare-Bendixsona

Zastosowanie twierdzenia Poincarégo-Bendixona nie jest w praktyce łatwe. Aby znaleźć region R tworzy się region pułapkowania (ang. trapping region) t.j. taki region na krawędziach którego pole wektorowe jest skierowane w stronę tego regionu.

example712_vanderpol.png

Figure 107: Oscylator van der Pola

par mu=1
par r=1

x'=v
v'=-mu*(x^2-1)*v - x

xc' = 1
yc' = 1

aux x1=r*cos(xc)
aux y1=r*sin(yc)

@ total=10,xp=x,yp=v, xlo=-4,ylo=-4, xhi=4 yhi=4
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Graphic Stuff —> Add Curve—> X-axis=x1, Y-axis=y1, Color=1
# Graphic Stuff —> Freeze –> On freeze
# Nulicine —> New
# Initialconds —> New —> X=2,V=0, XC=0, XC=0
# Initialconds —> New —> X=-2,V=3, XC=0, XC=0
# Initialconds —> New —> X=2,V=-3, XC=0, XC=0
# Parameters —> r ENTER 3 ENTER
# Initialconds —> New —> X=2,V=0, XC=0, XC=0
# Dir.filed/flow —> Scaled Dir.Fld –> 16

Przechodzimy do układu współrzędnych polarnych.

example712_vanderpol_polar.png

Figure 108: Oscylator van der Pola

par mu=1
par rc=1
r'=mu*r*(sin(theta))^2*(1-r^2*(cos(theta))^2)
theta'=mu*cos(theta)*sin(theta)*(1-r^2*cos(theta)^2)-1

xc' = 1
yc' = 1

aux x1=rc
aux y1=yc
aux x2=r*cos(theta)
aux y2=r*sin(theta)

@ total=10,xp=r,yp=theta, xlo=0,ylo=0, xhi=4 yhi=6.28
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Phasespace —> Choose —> Period = Default ENTER Theta
# Graphic Stuff —> Add Curve—> X-axis=x2, Y-axis=y2, Color=1
# Initialconds —> New —> R=2,THETA=0, XC=0, XC=0
# Parameters —> r ENTER 3 ENTER
# Initialconds —> New —> X=2,V=0, XC=0, XC=0
# Dir.filed/flow —> Scaled Dir.Fld –> 16

Rezultat jest podobny.

\[\dot r=\mu r sin^2\theta(1-r^2cos^2\theta)\]

\[\mu r sin^2\theta(1-r^2cos^2\theta) > 0\] \[1 - r^2cos^2\theta > 0\] \[r^2cos^2\theta < 1\] \[r < \sqrt{\frac{1}{cos^2\theta}},\;cos\theta=1\] \[cos\theta = \frac{\pi}{2}, \; r < \sqrt{\frac{1}{0}}\] \[r_{min} < 1\]

\[\mu r sin^2\theta(1-r^2cos^2\theta) < 0\] \[1 - r^2cos^2\theta < 0\] \[r > \sqrt{\frac{1}{cos^2\theta}},\;cos\theta=1\] \[cos\theta = \frac{\pi}{2}, \; r > \sqrt{\frac{1}{0}}\] \[r_{max} > \infty\]

Pytanie: co się dzieje w układzie dla \(\theta = \pm \frac{\pi}{2}\)?

example712_vanderpol_polar_check.png

Figure 109: Oscylator van der Pola - weryfikacja obliczeń

par mu=1

r'=mu*r*(sin(theta))^2*(1-r^2*(cos(theta))^2)
theta'=mu*cos(theta)*sin(theta)*(1-r^2*cos(theta)^2)-1

x'=v
v'=-mu*(x^2-1)*v - x

aux x1=r*cos(theta)
aux y1=r*sin(theta)

@ total=10,xp=x,yp=v, xlo=-3,ylo=-3, xhi=3 yhi=3
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

# GUI 
# Initialconds —> New —> R=2,THETA=0, X=2, Y=0
# Graphic Stuff —> Add Curve—> X-axis=x1, Y-axis=y1, Color=1, Line type=-4

4.4.8 Przykład - region pułapkowania dla układu we współrzędnych polarnych

\[\dot r = r(1-r^2) - \mu rcos\theta\] \[\dot \theta = 1\]

Twierdzenie Poincarégo-Bendixona może zostać zastosowane do weryfikacji, czy w danym układzie występuje orbita zamknięta. Sprowadza się do określenie regionu pułapkowania, czyli takiej podprzestrzeni w przestrzeni fazowej naszego układu, na krańcach której wektory są skierowane do tej podprzestrzeni. W układzie współrzędnych polarnych określa się rmin i rmax, dla których odpowiednio \(\dot r >0\) i \(\dot r <0\). W układzie współrzędnych kartezjańskich stosuje się analizę izoklin zerowego wzrostu.

\[r_{min}(1-r_{min}^2) - \mu r_{min}cos\theta>0\]

\[(1-r_{min}^2) > \mu cos\theta\] \[r_{min}^2 < 1 - \mu cos\theta\] \[r_{min} < \sqrt{1 - \mu cos\theta},\; cos\theta = 1\]

dla \(\mu < 1\) istnieje minimalny promień początku regionu pułapkowania: dla \(\mu = 0.5,\; r_{min}<\frac{\sqrt{2}}{2}\)

\[r_{max}(1-r_{max}^2) - \mu r_{max}cos\theta<0\] \[r_{max} > \sqrt{1 - \mu cos\theta},\; cos\theta = -1\] dla \(\mu = 0.5,\; r_{max}>\frac{\sqrt{6}}{2}\)

example731trapping_region.png

Figure 110: Region pułapkowania

Takie poszukiwanie regionu pułapkowania jest wygodne, gdy mamy układ współrzędnych polarnych

4.4.9 Przykład - region pułapkowania dla układu wspołrzędnych kartezjańskich

\[\dot x = -x + ay + x^2y\] \[\dot y = b - ay - x^2y\]

Punkty równowagi:
\[-x + ay + x^2y=0,\;-x + y(a + x^2)=0,\;y=\frac{x}{a + x^2}\] \[b - ay - x^2y=0,\;b - \frac{x}{a+x^2}(a + x^2)=0\] \[x_{0}=b,\;y_{0}=\frac{b}{a+b^2}\]

Linearyzacja:
\[A = \begin{pmatrix} -1+2xy & a+x^2\\ -2xy & -a-x^2\end{pmatrix}\]

\[A\biggl(b,\frac{b}{a+b^2}\biggr) = \begin{pmatrix} -1+\frac{2b^2}{a+b^2} & a+b^2\\ -\frac{2b^2}{a+b^2} & -a-b^2\end{pmatrix}\] \[A\biggl(b,\frac{b}{a+b^2}\biggr) = \begin{pmatrix} \frac{b^2-a}{a+b^2} & a+b^2\\ -\frac{2b^2}{a+b^2} & -a-b^2\end{pmatrix}\]

\[det(A)=\frac{b^2-a}{a+b^2}(-a-b^2)-(a+b^2)(-\frac{2b^2}{a+b^2})=a-b^2+2b^2=a+b^2\]

Pytanie: dla jakich parametrów a i b orbita zamknięta może wystąpić?

Jeśli \(a + b^2 < 0\) to dostaniemy siodło, czyli dla takich warunków nie jest możliwa orbita zamknięta. \(b < \pm \sqrt{-a}\) - ten warunek może być spełniony dla a < 0

Załóżmy, że nasz układ jest zaprojektowany dla a > 0.

\[\tau=\frac{b^2-a}{a+b^2} -(a + b^2)\] Pytanie: jaki znak będzie miało τ?

  • \(-(a + b^2)\) zawsze będzie ujemne
  • \(\frac{(b^2-a)}{a+b^2}\)
    • dla \(b^2-a \le 0,\;b \le \pm \sqrt{a}\)
    • dla \(b^2-a > 0,\;b>\pm \sqrt{a}\)

Z tej analizy wynika, że τ prawie zawsze będzie ujemne, czyli będziemy mieli zazwyczaj układ stabilny. W takim razie jeśli orbita jest w naszym układzie to jest niestabilna. Istnieje jednak region, dla którego τ jest większe od zera: znajdziemy go dla a<1.

\[b^2 = a^2 + ab^2 + ab^2 + b^4 + a\] \[b^2 - b^4= a^2 + 2ab^2 + a\]

dla \(a<1,\; a^2 \ll 1\)

wtedy nasze równanie: \[b^2 - b^4= a(2b^2 + 1)\] \[a \approxeq \frac{b^2(1 - b^2)}{2b^2 + 1}\]

dla niewielkiego zakresu a i b znajdziemy τ dodatnie, czyli układ będzie niestabilny a co za tym idzie orbita zamknięta, jeśli istnieje, będzie stabilna.

example732_tau.png

Orbita zamknięta może wystąpić, gdy \(\tau^2-4\Delta<0\)

\[\biggl(\frac{b^2-a}{a+b^2} -(a + b^2)\biggr)^2 - 4(a^2+b^2)<0\] \[\biggl(\frac{b^2-a}{a+b^2}\biggr)^2 - \frac{b^2-a}{a+b^2}(a + b^2) + (a + b^2)^2 - 4(a^2+b^2)<0\] \[\biggl(\frac{b^2-a}{a+b^2}\biggr)^2 - (b^2-a) + (a + b^2)^2 - 4(a^2+b^2)<0\] \[\biggl(\frac{b^2-a^2}{a+b^2}\biggr)^2 + (a + b^2)^2 - 4a^2 - 5b^2+a<0\] \[\frac{(b^2-a^2)^2+(a+b^2)^4}{(a+b^2)^2} - 4a^2 - 6b^2 + a+b^2<0\] \[\frac{(b^2-a^2)^2+(a+b^2)^4 + (a+b^2)^2}{(a+b^2)^2} < (4a^2 + 6b^2)\]

Dla a=0.05 i b=0.5, τ=0.36666 \(\tau^2-4\Delta = -1.06555 <0\)

Dla a=0.5 i b=0.5, \(\tau=-1.08333\) \(\tau^2-4\Delta = -1.82638<0\)

Dla a=2 i b=3, \(\tau=-10.363636\) \(\tau^2-4\Delta = 63.405>0\)

Pytanie: w przestrzeniach fazowych widzianych na rysunkach poniżej znajdują się orbity zamknięte?

example732_a005b05.png

Figure 112: a = 0.05 i b = 0.5

example732_a05b05.png

Figure 113: a = 0.5 i b = 0.5

Szukamy region pułapkowania dla układu z parametrami a = 0.05 i b = 0.5 \[\dot x = -x + ay + x^2y\] \[\dot y = b - ay - x^2y\]

Izokliny zerowego wzrostu:

  1. \(\dot x = 0\)

\[-x + ay + x^2y = 0,\; -x + y(a + x^2) = 0,\; y=\frac{x}{a+x^2}\] \[\dot y = b - \frac{x}{a+x^2}(a + x^2)=b-x\]

  1. \(\dot y = 0\)

\[b - ay - x^2y=0,\;b - y(a + x^2)=0,\;y=\frac{b}{a+x^2}\] \[\dot x = -x + \frac{b}{a+x^2}(a + x^2)=b-x\]

Dla x,y > 0 i a,b > 0 Region pułapkowania jest ograniczony od dołu - oś X:

  • dla y=0 \(\dot x =-x < 0\) i \(\dot y = b>0\) wektor zawsze będzie skierowany w lewo i do góry. dla x=0 będzie skierowany do góry

Region pułapkowania jest ograniczony od lewej strony - oś Y:

  • dla x = 0 \(\dot x = ay > 0\) wektor zawsze będzie skierowany w prawo.
  • \(\dot y = b-ay\) będzie zmieniało znak w punkcie \(\dot y = b-ay=0,\; y=\frac{b}{a}\)
  • dla \(y > \frac{b}{a}\) wektor będzie skierowany w dół.
  • dla \(y < \frac{b}{a}\) wektor będzie skierowany w górę.

Na tym etapie możemy już odrzucić hipotezę istnienia orbity zamkniętej dla układu z parametrami a = 0.5 i b = 0.5. Wektory na osiach są skierowane do regionu pułapkowania a nasz punkt stały jest stabilny więc na granicach regionu pułapkowania wokół punktu stałego pole wektorowe będzie skierowane w stronę tego punktu.

Aby zamknąć region musimy znaleźć jeszcze jego wszystkie granice. W miejscu \(y=\frac{b}{a}\) izoklina zerowego rzędu \(\dot y\) styka się z osią Y. W tym miejscu możemy ograniczyć nasz zbiór od góry. Pozostaje nam znalezienie ograniczenia od prawej strony.

Szukamy granicy, która będzie skierowana w dół pod kątem. Kierunek pola wektorowego jest jakościowo różnych dla regionów dzielonych przez izokliny zerowego rzędu. Taki kierunek utrzyma się aż do izokliny czerwonej, po przekroczeniu której pole wektorowe zmienia kierunek - w dół i lewo by następnie osiągnąć izoklinę zieloną.

Wiemy zatem, że nasz układ jak rozpocznie działanie w regionie to pole wektorowe zawsze będzie skierowane w dół i prawą stroną. Szukamy prostej wyznaczającej ten region np x=b i \(y=\frac{b}{a}\) . Jakie nachylenie tej krzywej musi być? \[\dot x = -x + ay + x^2y\] \[\dot y = b - ay - x^2y\]

dla \(x \rightarrow \infty\) \[\dot x \approxeq x^2y\] \[\dot y \approxeq -x^2y\]

\(\frac{\dot x}{\dot y} \approxeq -1\)

dla x>b pole wektorowe będzie skierowane w dól i prawą stronę
np. dla x=2b \[\dot x = -2b + a\frac{b}{a} + 4b^2\frac{b}{a}=\frac{4b^3}{a}-b\] to będzie zawsze dodatnie, gdyż a<1 \[\dot y = b - a\frac{b}{a} - 4b^2\frac{b}{a}=-\frac{4b^3}{a}\]

Ten dodatkowy czynnik ‚b‘ w równaniu na \(\dot x\) powoduje, że dla coraz większych x wektor będzie o niewielką część mniejszy od wektora \(\dot y\) co powoduje, że pole wektorowe będzie miało kierunek tuż pod przekątną \(\frac{\dot x}{\dot y}=-1\)

\(x=b\) i \(y=\frac{b}{a}\) \[\dot x = -b + a\frac{b}{a} + b^2\frac{b}{a}=\frac{b^3}{a}\] \[\dot y = b - a\frac{b}{a} - b^2\frac{b}{a}=-\frac{b^3}{a}\]

Szukamy prostej \(y-y_{0} =m(x-x_{0})\)

\[k = \sqrt{\frac{b^6}{a^2} + \frac{b^6}{a^2}}= \frac{b^3}{a}\sqrt{2}\] \[m = - \frac{b^3}{a}\sqrt{2}\] \[y-y_{0} =m(x-x_{0}), y-\frac{b}{a} = -\frac{b^3}{a}\sqrt{2}(x-b)\] \[y = -\frac{b^3}{a}\sqrt{2}(x - b) +\frac{b}{a}\]

example732_trapping.png

Figure 114: Region pułapkowania dla układu z parametrami a = 0.05 i b = 0.5

par a=0.05
par b=0.5

x'=-x+a*y+x^2*y
y'=b-a*y-x^2*y
xc' = 1
yc' = 1

aux x3=b
aux y3=yc
aux y4=-(b^3/a)*sqrt(a)*(xc-b) + b/a
aux x4 = xc

@ total=10,xp=x,yp=y,
@ xlo=0,ylo=0, xhi=20 yhi=12
@ maxstor=2000000, bounds=100000, nmesh=800 
done

# GUI 
# Nulicine —> New
# Graphic Stuff —> Freeze –> On freeze
# Graphic Stuff —> Add Curve—> X-axis=x3, Y-axis=y3, Color=3
# Graphic Stuff —> Add Curve—> X-axis=x4, Y-axis=y4, Color=3
# Dir.filed/flow —> Scaled Dir.Fld –> 32
# Initialconds —> Mouse —> wybierz kilka warunkow początkowcy blisko krawdzie regionu pułakownia 

Z twierdzenia Poincarégo-Bendixona wynika także, że w układach dwuwymiarowych mamy bardzo ograniczony zbiór zachowań tych układów - od punktów stałych do cykli granicznych. Cykl graniczny występuje tylko w układach nieliniowych. Oscylacje w cyklu graniczny są zdeterminowane przez strukturę układu nieliniowego w przeciwieństwie do orbit zamkniętych w układach liniowych, gdzie amplituda oscylacji jest zdeterminowana przez warunki początkowe.

4.4.10 Oscylator relaksacyjny

Przykładem jest oscylator van der Pola

relaxation_oscillator.gif

Figure 115: Oscylacje relaksacyjne

par mu=1

x'=v
v'=-mu*(x^2-1)*v - x

@ total=20,xp=x,yp=v, xlo=-3
@ ylo=-10, xhi=15 yhi=10
@ maxstor=2000000, bounds=100000,  nmesh=800
@ nplot=2, xp2=t,yp2=v

done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 2 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: 1 —> End: 6
# Kinescope —> Make AniGif

Wraz ze wzrostem nieliniowości czynnika wpływającego na tłumienie układu zwiększa się okres pojawiania się relaksacji w postaci nagłych rozładowań. Takie oscylacje z fazą zwężenia szyjkowego i nagłych rozładowań nazywamy relaksacyjnymi. Oscylacje tego typu nazywa się relaksacyjnymi, gdyż nagromadzony "stres" podczas wolniejszej aktywności oscylatora jest uwalniany "zrelaksowany" podczas szybkiej aktywności. Takie zachowanie też wyraźnie pokazuje dwie skale czasowe: wolną i szybką.

relaxation_oscillator2.gif

Figure 116: Oscylacje relaksacyjne - pole wektorowe

par mu=1

x'=v
v'=-mu*(x^2-1)*v - x

@ total=20,xp=x,yp=v, xlo=-3
@ ylo=-10, xhi=15 yhi=10
@ maxstor=2000000, bounds=100000,  nmesh=800
@ dt=0.01
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 2 ENTER 0 ENTER
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: 1 —> End: 6
# Kinescope —> Make AniGif

Porównaj z oscylatorem niejednorodnym: 45, 46

4.4.11 Przykład: oszacowanie okresu oscylacji

\[\ddot x + \mu(x^2-1)\dot x +x=0,\; \mu \gg 1\]

\[\dot x = v,\; \dot v = -\mu(x^2-1)v - x=0\]

Analiza izoklin zerowego wzrostu:

  1. \(\dot x = 0, \; v=0,\; \dot v = -x\)
  2. \(\dot v=0,\; -\mu(x^2-1)v - x=0,\; v=-\frac{x}{\mu(x^2-1)},\;\dot x = -\frac{x}{\mu(x^2-1)}\)

Druga izoklina jest hiperbolą z punktami osobliwymi dla \(x=\pm 1\)

example75_relaxation2.pdf

example751_nullclines.png

Figure 117: Portret fazowy układu z parametrem μ = 6

example75_relaxation_period2.gif

Figure 118: Izokliny zerowego wzrostu

par mu=1

x'=v
v'=-mu*(x^2-1)*v - x

@ total=40,dt=0.01,xp=x,yp=v, xlo=-3
@ ylo=-14, xhi=3 yhi=14
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> -1 ENTER 2 ENTER
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: 1 —> End: 8
# Kinescope —> Make AniGif

Im większa nieliniowość tym izoklina bliżej osi X i układ najwięcej czasu spędza dla x > 1 i x < -1, gdyż wtedy prędkość jest bliska zeru. Ponieważ układ porusza się prawie wzdłuż izokliny to okres możemy oszacować na podstawie scałkowania izokliny w tym zakresie:

\[\dot x = -\frac{x}{\mu(x^2-1)}\] \[-\frac{\mu(x^2-1)}{x} dx = dt\] \[-\mu(x - \frac{1}{x}) dx = dt\] \[T = 2\int_{x_{1}}^{x_{2}}\mu(\frac{1}{x}-x) dx\] \[T = 2\mu(\ln{x} - \frac{1}{2}x^2)\biggr|_{x_{1}}^{x_{2}}\] \[T = 2\mu(\ln{x} - \frac{1}{2}x^2)\biggr|_{1}^{2} dx\]

Ponieważ nasz układ zmienia się od x1 = 2 do x2 = 1 to:
\[T = 2\mu(\ln{x} - \frac{1}{2}x^2)\biggr|_{2}^{1} dx\] \[T = 2\mu(- \frac{1}{2}-\ln{2} + 2)\] \[T = \mu(3 - 2\ln{2})\]

Weryfikacja oszacowań: dla μ=6, T=9.68

\(\delta_{err} = (1 - \frac{9.68}{13.06})\cdot 100 = 26\%\)

dla μ=10, T=16.13

\(\delta_{err} = (1 - \frac{16.13}{19})\cdot 100 = 15\%\)

Aby przekonać się, że istnieje faza szybka i wolna a zakres scałkowania powinien być od x = 1 do x = 2 zastosujmy przekształcenie:

\[\ddot x + \mu(x^2-1)\dot x +x=0 = \frac{d}{dt}[\dot x + \mu(\frac{1}{3}x^3-x)] + x\] \[F(x) = \mu(\frac{1}{3}x^3-x)\] \[w = \dot x + F(x)\] \[\dot w + x=0\]

\[\dot x = w-F(x),\; \text{prędkość będzie bliska zeru jeśli układ będzie poruszał się wzdłuż F(x)}\] \[\dot w = -x\]

Analiza izoklin zerowego wzrostu:

  1. \(\dot x = 0, \; w-F(x)=0,\; w= \mu(\frac{1}{3}x^3-x),\;\dot w = ?\)
  2. \(\dot w=0,\; x=0,\; \dot x = w\)

example75_relaxation3.pdf

Aby uniezależnić izoklinę \(\dot x=0\) od wartości μ zastosujmy dodatkowe przekształcenie:
\(y=\frac{w}{\mu}\)

\[\dot x = y\mu-F(x)\] \[\dot y = -\frac{x}{\mu}\]

Analiza izoklin:

  1. \(\dot x = 0, \; \mu y-F(x)=0,\; y= \frac{1}{3}x^3-x,\;\dot y = ?\)
  2. \(\dot y=0,\; x=0,\; \dot x = \mu y\)

\[\frac{dy}{dt} = \frac{dF(x)}{dx}\frac{dx}{dt}\] \[\frac{dy}{dt} = (x^2-1)\frac{dx}{dt}\] \[-\frac{x}{\mu} = (x^2-1)\frac{dx}{dt}\] \[dt= -\mu \frac{(x^2-1)}{x}{dx}\] \[T = 2\mu\int_{x1}^{x2} (\frac{1}{x}-x)dx=2\mu(ln{|x|}-0.5x^2)\biggl|_{2}^{1}=2\mu(-0.5 - ln{2} + 2)=\mu(3-2ln{2})\]

Obliczenia pola wektorowego:

xf=range(-1,2,length=10);
xs=range(2,1,length=10); 
mu=6; y=2/3;
dotx = mu*y .- mu .* (1/3 .* xf .^3 .- xf);
doty = collect(-xf ./ mu);
dotx2 = xs ./ (mu*(xs .^2 .- 1));
[xf dotx doty xs dotx2]
10×5 Array{Float64,2}:
 -1.0       0.0           0.166667   2.0        0.111111
 -0.666667  0.592593      0.111111   1.88889    0.122596
 -0.333333  2.07407       0.0555556  1.77778    0.137143
  0.0       4.0           0.0        1.66667    0.15625 
  0.333333  5.92593      -0.0555556  1.55556    0.182609
  0.666667  7.40741      -0.111111   1.44444    0.221591
  1.0       8.0          -0.166667   1.33333    0.285714
  1.33333   7.25926      -0.222222   1.22222    0.4125  
  1.66667   4.74074      -0.277778   1.11111    0.789474
  2.0       8.88178e-16  -0.333333   1.0      Inf

Pierwsza kolumna określa zakres zmian x, na którego podstawie zostały obliczone dwie kolejne \(\dot x\) i \(\dot y\). Czwarta kolumna kolejny zakres zmian x, na którego podstawie została obliczona piąta kolumna \(\dot x\). Z tych liczb wyraźnie widać, że układ przechodzi dwie fazy: szybką od x = -1 do x = 2 oraz wolną od x = 2 do x = 1. Podczas fazy szybkiej prędkość \(\dot x\) jest o rząd wielkości większa niż podczas fazy wolnej.

example75_relaxation_calulation.png

Figure 119: Obliczenia prędkości

par mu=6

x'=mu*(y - (1/3*x^3-x))
y'=-x/mu

aux dotx = mu*y - mu*(1/3*x^3-x)
aux dotx2 = x/(mu*(x^2-1))

@ total=40,dt=0.05,xp=x,yp=y, xlo=-3
@ ylo=-3, xhi=3 yhi=3
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 2 ENTER 0 ENTER
# Data

example75_relaxation4.pdf

example75_relaxation_period.gif

Figure 120: Izokliny zerowego wzrostu

par mu=1

x'=mu*(y - (1/3*x^3-x))
y'=-x/mu

@ total=40,dt=0.05,xp=x,yp=y, xlo=-3
@ ylo=-3, xhi=3 yhi=3
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> -1 ENTER 2 ENTER
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: 1 —> End: 8
# Kinescope —> Make AniGif

W oscylatorze relaksacyjnym faza szybka i wolna następują po sobie. W słabo nieliniowych oscylatorach (ang. weakly nonlinear oscillators) te fazy oddziałują dodatkowo na siebie tzn., że w tym samym czasie te dwie skale czasowe są aktywne, tak jak to widać na rysunku poniżej. Na czerwono są zaznaczone oscylacje z fazą szybką i wolną, a na czarno oscylacje z oddziałującymi na siebie skalami czasowymi. Bardzo wolna skala czasowo określa zmiany amplitudy oscylacji, natomiast szybka określa częstotliwość oscylacji. W tłumionym oscylatorze harmonicznym są podobne zmiany amplitudy 61

example_time_scales.png

Figure 121: Skale czasowe

4.4.12 Przykład słabo nieliniowego układu drgającego

\[\ddot x + x + \epsilon h(x,\dot x)=0,\;0 \le \epsilon \ll 1\]

Ze względu na niewielki czynnik nieliniowy te oscylatory wykazują podobne zachowanie do oscylatorów harmonicznych. Jeżeli w równaniu van der Pola μ ustawimy na bardzo małym poziomie to portret fazowy będzie podobny do oscylatora harmonicznego. W oscylatorach liniowych amplituda oscylacji w nietłumionych jest zdeterminowana przez warunki początkowe i pozostaje stała co skutkuje tym, że istnieje tam orbita zamknięta; w tłumionych amplituda zmienia się do zera lub nieskończoności. W słabo nieliniowych oscylatorach istnieje orbita zamknięta, która nie jest zdeterminowana warunkami początkowymi lecz strukturą równania i amplituda oscylacji może się zmieniać o czym świadczą dwie skale czasowe. W słabo nieliniowych oscylatorach mamy minimum dwie skale czasowe.

weakly_nonlinear_oscillator.png

Figure 122: słabo nieliniowy oscylator

par mu=0.01

x'=v
v'=-mu*(x^2-1)*v - x

@ total=1000,dt=0.05,xp=x,yp=v, xlo=-2
@ ylo=-2, xhi=2 yhi=2
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Initialconds —> New —> 0.3 ENTER 0 ENTER
  1. Metoda podwójnej skali czasowej (ang. two-timing)

    Skala szybka: \(\tau = t\)
    Skala wolna: \(T=\epsilon t\)

    \(\tau\; i\; T\) - zmienne niezależne

    Jeśli rozwiązaniem równania układu jest \(x(t,\epsilon)\) to rozwinięcie w szereg potęgowy:
    \(x(t,\epsilon) = x_{0}(t) + \epsilon x_{1}(t) + \epsilon^2 x_{2}(t) + \dotso\)

    W naszych przybliżeniach mamy nadzieję, że rozwiązanie będzie zawarte w dwóch pierwszych elementach szeregu. Dla dwóch skali czasowych:
    \(x(t,\epsilon) = x_{0}(\tau,T) + \epsilon x_{1}(\tau,T) + O(\epsilon^2)\)

    Nieznane funkcje xi(t) należy określić z wykorzystaniem równania ruchu i warunków początkowych. Wykorzystując regułę łańcuchową przechodzimy z równania różniczkowego o jednej zmiennej niezależnej t do równania różniczkowego o dwóch zmiennych niezależnych \(\tau\) i T

    \(\frac{d\tau}{dt} = 1\) \(\frac{dT}{dt} = \epsilon\)

    \[\frac{dx}{dt} = \frac{\partial x}{\partial \tau} \frac{\partial \tau}{\partial t} + \frac{\partial x}{\partial T} \frac{\partial T}{\partial t}\] \[\frac{dx}{dt} = \frac{\partial x}{\partial \tau} + \epsilon \frac{\partial x}{\partial T}\]

    \[\frac{d}{dt}\frac{dx}{dt} = \frac{d}{dt}\biggl(\frac{\partial x}{\partial \tau} + \epsilon \frac{\partial x}{\partial T}\biggr)\]

    \[\frac{d}{dt}\biggl(\frac{\partial x}{\partial \tau}\biggr)=\frac{\partial^2 x}{\partial \tau^2}\frac{d\tau}{dt} + \frac{\partial^2 x}{\partial \tau \partial T}\frac{dT}{dt} \] \[\frac{d}{dt}\biggl(\epsilon \frac{\partial x}{\partial T}\biggr)=\epsilon \biggl(\frac{\partial^2 x}{\partial T^2}\frac{dT}{dt} + \frac{\partial^2 x}{\partial T \partial \tau}\frac{d\tau}{dt}\biggr)\]

    \[\frac{d}{dt}\frac{dx}{dt} = \frac{\partial^2 x}{\partial \tau^2} + \epsilon\biggl(\frac{\partial^2 x}{\partial \tau \partial T} + \frac{\partial^2 x}{\partial T \partial \tau}\biggr) + \epsilon^2\biggl(\frac{\partial^2 x}{\partial T^2 }\biggr)\]

    Wykorzystując notację indeksu dolnego mamy: \[\dot x = \partial_{\tau}x + \epsilon \partial_{T}x\] \[\ddot x = \partial_{\tau \tau}x + 2\epsilon \partial_{\tau T}x + \epsilon^2 \partial_{TT}x\]

    Podstawiamy rozwinięcie w szereg potęgowy do równań o dwóch zmiennych niezależnych:

    \(\dot x = \partial_{\tau} (x_{0}(\tau,T) + \epsilon x_{1}(\tau,T)) + \epsilon \partial_{T}(x_{0}(\tau,T) + \epsilon x_{1}(\tau,T)) + O(\epsilon^2)\)

    \(\dot x = \partial_{\tau} x_{0} + \epsilon(\partial_{\tau}x_{1} + \partial_{T} x_{0}) + O(\epsilon^2)\)

    \[\ddot x = \partial_{\tau \tau}(x_{0} + \epsilon x_{1}) + 2\epsilon \partial_{\tau T}(x_{0} + \epsilon x_{1}) + \epsilon^2 \partial_{TT}(x_{0} + \epsilon x_{1})\]

    \[\ddot x = \partial_{\tau \tau}x_{0} + \epsilon(\partial_{\tau \tau} x_{1} + 2 \partial_{\tau T}x_{0}) + O(\epsilon^2)\]

  2. Przykład: przybliżone rozwiązanie dla tłumionego oscylatora liniowego

    \[\ddot x + 2\epsilon \dot x + x =0\]

    \[\partial_{\tau \tau}x_{0} + \epsilon(\partial_{\tau \tau} x_{1} + 2 \partial_{\tau T}x_{0}) +2\epsilon \partial_{\tau} x_{0} + 2\epsilon^2(\partial_{\tau}x_{1} + \partial_{T} x_{0}) + x_{0} + \epsilon x_{1}=0\]

    \[O(1): \partial_{\tau \tau}x_{0} + x_{0} = 0\] \[O(\epsilon): \partial_{\tau \tau} x_{1} + 2 \partial_{\tau T}x_{0} + 2\partial_{\tau} x_{0} +x_{1}=0\]

    \[\partial_{\tau \tau}x_{0}(\tau,T) + x_{0}(\tau,T) = 0\] \[\frac{\partial^2 x_{0}(\tau,T)}{\partial \tau^2} + x_{0}(\tau,T)=0\]

    Jeśli potraktujemy T jako stałą to wtedy:

    \[\ddot x_{0} + x_{0} = 0,\;\text{nietłumiony oscylator harmoniczny}\]

    Rozwiązaniem takiego układu jest:
    \(\dbinom{x_{0}}{v_{0}} = c_{1}\dbinom{\cos{t}}{ -\omega \sin{t}} + c_{2}\dbinom{\sin{t}}{\omega \cos{t}\)

    \(x_{0} = c_{1}\cos{t} + c_{2}\sin{t}\) \(x_{0} = c_{1}\cos{\tau} + c_{2}\sin{\tau}\)

    W nietłumionym oscylatorze współczynniki c1 i c2 są stałe. W dalszej analizie potraktujemy te współczynniki jak funkcje czasu wolnego: c1(T) i c2(T)

    Podstawiamy do równania:

    \[\partial_{\tau \tau} x_{1} + 2 \partial_{\tau T}(c_{1}\cos{\tau} + c_{2}\sin{\tau}) + 2\partial_{\tau} (c_{1}\cos{\tau} + c_{2}\sin{\tau}) +x_{1}=0\]

    \[\partial_{\tau \tau} x_{1} + x_{1} = 2c_{1}'\sin{\tau} - 2c_{2}'\cos{\tau} + 2c_{1}\sin{\tau} - 2c_{2}\cos{\tau}\] \[\partial_{\tau \tau} x_{1} + x_{1} = 2(c_{1}' + c_{1})\sin{\tau} - 2(c_{2}'+ c_{2})\cos{\tau}\]

    Naszym celem jest znalezienie przybliżonego rozwiązania x(τ,T). Kolejnym krokiem w metodzie podwójnego skali czasowej jest ustawienie współczynników czynnika rezonującego do zera:

    \(c_{1}' + c_{1} = 0,\;c_{2}' + c_{2} = 0\) \(\frac{dc_{1}}{dT} = - c_{1}\) \(\frac{dc_{2}}{dT} = - c_{2}\)

    \[-\ln{|c_{1}|}=T + C\] \[C = -\ln{|c_{01}|}\] \[-\ln{|c_{1}|}=T -\ln{|c_{01}|}\] \[-\ln{|c_{1}|}+\ln{|c_{01}|} =T\] \[\ln{\frac{|c_{01}|}{|c_{1}|} =T\] \[c_{1} = c_{01}e^{-T}\] \[c_{2} = c_{02}e^{-T}\]

    Warunki początkowe: \(x(0)=0, \dot x(0) = 1\)

    \(x(t,\epsilon) = x_{0}(\tau,T) + \epsilon x_{1}(\tau,T) + O(\epsilon^2)\) \(x(0) = 0 = x_{0} + \epsilon x_{1}\)

    \[\partial_{\tau \tau} x_{1} + x_{1} = 2(c_{1}' + c_{1})\sin{\tau} - 2(c_{2}'+ c_{2})\cos{\tau}\]

4.5 Bifurkacje

Bifurkacje są to przejścia pomiędzy stanami układów, które nie są sobie równoważne topologicznie tzn. pojawiają lub znikają punkty stałe, orbity oraz może ulec zmianie ich stabilność. Bifurkacje siodło-węzeł, transkrytyczne i widłowe są przykładami bifurkacji zerowej wartości własnej (ang. zero-eigen). Dodatkową cechą tych bifurkacji jest lokalność występowania. Są one związane z punktami stałymi. Bifurkacjami lokalnymi są także bifurkacje Hopfa, ale w przeciwieństwie do poprzednich bifurkacji, to w tej wartości własne są liczbami zespolonymi i w trakcie bifurkacji część rzeczywista się zeruje.

bif_saddlenode_eigen.png

Figure 123: Zmiana wartości własnych podczas bifuracji siodło-węzeł

bif_transcritical_eigen.png

Figure 124: Zmiana wartości własnych podczas bifuracji transkrytycznej

bif_pitchfork_eigen.png

Figure 125: Zmiana wartości własnych podczas bifuracji widłowej

bif_hopf_eigen.png

Figure 126: Zmiana wartości własnych podczas bifuracji Hopfa

4.5.1 Bfurkacja siodło-węzeł w układzie o zmiennych niesprzężonych

\[\dot x = \mu - x^2\] \[\dot y = - y\]

Punkt stały dla układu 1D \(\dot x = \mu - x^2\): \(x_{0} = \pm \sqrt{\mu},\;\mu>0\)

Punkt \(\mu=0\) jest punktem bifurkacji, gdyż zerowa wartość własna pojawia się w tym punkcie.

\[A = \begin{pmatrix} -2x & 0\\ 0 & 0\end{pmatrix}\] \[A(\sqrt{\mu},0) = \begin{pmatrix} -2\sqrt{\mu} & 0\\ 0 & 0\end{pmatrix},\;det(A)=0,\;\tau=-2\sqrt{\mu}\] \[A(\sqrt{\mu},0) = \begin{pmatrix} -2\sqrt{\mu} -\lambda & 0\\ 0 & 0-\lambda \end{pmatrix}=\lambda^2 + 2\sqrt{\mu}\lambda \] \[\Delta = 4\mu\] \[\lambda_{1,2} = \frac{-2\sqrt{\mu} \pm \sqrt{\Delta}}{2}\] \[\lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2-4det(A)}}{2}\] \[\lambda_{1} = 0,\;\lambda_{2} = -2\sqrt{\mu}\] Obliczanie wektorów własnych: \[\begin{pmatrix} -2\sqrt{\mu}-0 & 0\\ 0 & 0-0\end{pmatrix}\dbinom{v_{x}}{v_{y}}=0\] \[V=\dbinom{0}{1}\]

\[\begin{pmatrix} -2\sqrt{\mu}+2\sqrt{\mu} & 0\\ 0 & 0-2\sqrt{\mu}\end{pmatrix}\dbinom{w_{x}}{w_{y}}=0\] \[W=\dbinom{1}{0}\]

Ponieważ \(\lambda_{1}=0\) i V=(0,1) to układ wzdłuż osi Y się nie porusza co wynika z tego, że analizujemy układ 1D. Układ porusza się natomiast wzdłuż osi X z wartością własną \(\lambda_{2}\)

Dla drugiego punktu stałego:

\[A(-\sqrt{\mu},0),\;det(A)=0,\;\tau=2\sqrt{\mu}\] \[\lambda_{1} = 2\sqrt{\mu},\;\lambda_{2} = 0\]

Obliczanie wektorów własnych: \[\begin{pmatrix} -2\sqrt{\mu}-2\sqrt{\mu} & 0\\ 0 & 0-2\sqrt{\mu}\end{pmatrix}\dbinom{v_{x}}{v_{y}}=0\] \[V=\dbinom{1}{0}\]

\[\begin{pmatrix} -2\sqrt{\mu}-0 & 0\\ 0 & 0-0\end{pmatrix}\dbinom{w_{x}}{w_{y}}=0\] \[W=\dbinom{0}{1}\]

Punkt stały dla \(\dot y = -y\): \(y_{0} = 0\)

Punkt stały dla układu 2D: \(X_{*} = (0,0),\;\mu=0\)
\[A = \begin{pmatrix} -2x & 0\\ 0 & -1\end{pmatrix}\] \[A(0,0) = \begin{pmatrix} 0 & 0\\ 0 & -1\end{pmatrix},\;det(A) = 0,\;\tau=-1\] \[\lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2-4det(A)}}{2}\] \[\lambda_{1} = 0,\;\lambda_{2} = -1,\;\text{zero-eigen}\]

Obliczanie wektorów własnych: \[\begin{pmatrix} 0-0 & 0\\ 0 & -1-0\end{pmatrix}\dbinom{v_{x}}{v_{y}}=0\] \[V=\dbinom{1}{0}\]

\[\begin{pmatrix} 0+1 & 0\\ 0 & -1+1\end{pmatrix}\dbinom{w_{x}}{w_{y}}=0\] \[W=\dbinom{0}{1}\]

example81_bifurcation.gif

Figure 127: Bifurkacja siodło-węzeł

par mu=0.1

x'=mu - x^2
y'=-y

@ total=1000,dt=0.05,xp=x,yp=y,
@ xlo=-0.5, ylo=-0.5, xhi=0.5 yhi=0.5
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Dir.filed/flow —> Direct Field –> 16
# Nullcline —> New
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: 0.1 —> End: -0.1
# Kinescope —> Make AniGif

W układzie znajdują się punkty stałe, które są zależne od wartości parametru \(\mu\). Jeśli \(\mu=0\) to wtedy mamy jeden punkt stały \(x_{*}=(0,0)\). Gdy \(\mu>0\) to układ ma dwa punkty stałe; gdy \(\mu<0\) układ nie ma żadnego. Dla \(\mu < 0\) zachodzi zjawisko, które omawialiśmy podczas analizy jednowymiarowych układów - zjawisko ducha 43, 45, 46, które jest charakterystyczne dla układów z bifurkacją siodło-węzeł.

example81_ghost.gif

Figure 128: Duch

par mu=-0.01

x'=mu - x^2
y'=-y

@ total=49,dt=0.05,xp=t,yp=x,
@ xlo=0, ylo=-1.5, xhi=49 yhi=0.1
@ maxstor=2000000, bounds=100000,  nmesh=800
done

# GUI 
# Initialconds —> Range —> Range Over: mu —> Steps: 249 —> Start: -0.01 —> End: -0.0001
# Kinescope —> Make AniGif

Podczas omawiania oscylatora niejednorodnego powiedzieliśmy sobie o prawie skalowania pierwiastka kwadratowego, które dostrzegliśmy w równaniu na okres oscylacji: \(T= \frac{2\pi}{\sqrt{\omega^2-a^2}}\approxeq \sqrt{\frac{2}{\omega}}\frac{\pi}{\sqrt{\omega-a}}\)

Okres oscylacji jest skalowany przez czynnik \(\frac{1}{\sqrt{\omega-a}}\), który jest formalnym zapisem naszego prawa skalowania. Skalowanie zależy od odległości od punktu bifurkacji.

Prawo skalowania możemy wykorzystać teraz do oszacowania wzrostu czasu przejścia przez zwężenie szyjkowe. Rozważmy trzy przypadki: dla \(\mu_{1}=-0.01,\mu_{2}=-0.002,\mu_{3}=-0.0001\)

\[T = k(\mu - \mu_{c})^{-\frac{1}{2}}, \mu_{c}=x_{*}=0\] \[T_{1} = k\mu_{1}^{-\frac{1}{2}}=10k\] \[T_{2} = k\mu_{2}^{-\frac{1}{2}}=22.36k\] \[T_{3} = k\mu_{3}^{-\frac{1}{2}}=100k\]

\[\frac{T_{1}}{T_{2}}=\sqrt{\frac{\mu_{2}}{\mu_{1}}}=0.447\] \[\frac{T_{1}}{T_{3}}=\sqrt{\frac{\mu_{3}}{\mu_{1}}}=0.1\]

\[x(10,\mu_{1}) = -0.156\] \[x(10,\mu_{1})\cdot 0.447 = x(22.36,\mu_{2}) = -0.069\] \[x(10,\mu_{1}) \cdot 0.1 = x(100,\mu_{3}) = -0.0156\]

4.5.2 Bifurkacja siodło-węzeł w układzie o zmiennych sprzężonych

Bardziej złożony przykład pojawiania się ducha, ale mechanizm pozostaje taki sam: anihilacja punktów stałych w zależności od wartości parametrów:

Analiza izoklin zerowego wzrostu:
\[\dot x = -ax + y\] \[\dot y = \frac{x^2}{1+x^2} - by\] \[-ax+y=0,\;y=ax,\;\dot y = \frac{x^2}{1+x^2} - bax\] \[\frac{x^2}{1+x^2} - by=0,\;y=\frac{x^2}{b(1+x^2)},\;\dot x = -ax + \frac{x^2}{1+x^2} - by\]

example811_nullcline.png

Izoklina \(\dot x\) jest zaznaczona na czerwono, a izoklina \(\dot y\) na niebiesko. Układ zawsze będzie miał minimum jeden punkt stały \(X_{*}=(0,0)\) niezależnie od parametrów a i b, ale będą one miały wpływ na jego typ. Gdy \(a=\frac{1}{2b}\) to w punkcie \(X_{*}=(1,\frac{1}{2b})\) zachodzi bifurkacja siodło-węzeł - np. dla a = 1 i b = 0.5. Gdy dla takiego układu czerwoną linię obrócimy w lewo to punkt stały znika, a gdy obrócimy ją w prawo to pojawią się dwa punkty stałe. Łącznie będą wtedy w układzie trzy punkty stałe.

Dla \(a=\frac{1}{2b}\) izoklina \(\dot x\) jest styczna do izokliny \(\dot y\) w punkcie \(X_{*}=(1,\frac{1}{2b})\). Dla \(a < \frac{1}{2b}\) punkt znika, ale pojawia się duch tego punktu - układ zachowuje się tak, jakby ten punkt istniał, gdyż w okolicach tego punkt pochodne są prawie równe 0.

example811_ghost.png

Figure 130: Duch

par a=1
par b=0.5001

x'=-a*x + y
y'=x^2/(1+x^2) - b*y

init x=1.5
init y=1.5

@ total=1000,dt=0.05,xp=x,yp=y,
@ xlo=0, ylo=0, xhi=2 yhi=2
@ maxstor=2000000, bounds=1000000,  nmesh=800
done
#GUI
# Dir.filed/flow —> Flow –> 16
# Initialconds —> Go
# Nullcline —> New
# Makewindow —> Create
# Wybierz nowo powstałe okno
# Xi vs t —> Plot vs t: X —> ENTER
# Makewindow —> Create
# Wybierz nowo powstałe okno
# Xi vs t —> Plot vs t: Y —> ENTER

Sprawdzamy punkty stałe:
\[\frac{x^2}{1+x^2}-abx=0\] \[x^2 - abx(1+x^2)=0\] \[x(x - ab(1+x^2))=0\] \[x_{*} = 0\] \[x - ab(1+x^2)=0\] \[abx^2 - x + ab=0\] \[\Delta = 1 - 4a^2b^2\] \[\Delta = 0 \;-\; \text{jeden punkt stały},\;a=\pm \frac{1}{2b}\] \[\Delta > 0 \;-\; \text{dwa punkty stałe}\] \[\Delta < 0 \;-\; \text{brak punktów stałych}\]

Sprawdzamy jakiego typu są punkty stałe:

\[A = \begin{pmatrix} -a & 1\\ \frac{2x(1+x^2)-2x^3}{(1+x^2)^2} & -b\end{pmatrix}\] \[A = \begin{pmatrix} -a & 1\\ \frac{2x}{(1+x^2)^2} & -b\end{pmatrix}\] \[A(0,0) = \begin{pmatrix} -a & 1\\ 0 & -b\end{pmatrix},\;det(A)=ab,\;\tau=-a-b\] \[X_{*}=(0,0)- a,b > 0,\;\text{ściek}\] \[X_{*}=(0,0)- a,b < 0,\;\text{źródło}\] \[X_{*}=(0,0) - a \oplus b < 0,\;\text{siodło}\]

example811_null_all.png

Figure 131: Analiza izoklin zerowego rzędu

\[\Delta=0\] \[A(1,\frac{1}{2b}) = \begin{pmatrix} -a & 1\\ 0.5 & -b\end{pmatrix},\;det(A)=ab-0.5\] \[a=\frac{1}{2b},b>0\; A(1,\frac{1}{2b}),\;det(A)=0,\;\text{połowicznie stabilny}\] \[a=-\frac{1}{2b},b>0\; A(1,\frac{1}{2b}),\;det(A)=-1,\;\text{siodło}\] \[a=\frac{1}{2b},b<0\; A(1,\frac{1}{2b}),\;det(A)=0,\;\text{połowicznie stabilny}\] \[a=-\frac{1}{2b},b<0\; A(1,\frac{1}{2b}),\;det(A)=-1,\;\text{siodło}\]

Gdy \(\Delta=0\) w układzie mamy dwa punkty stałe:

  1. Jeśli a,b>0: \(X_{*}=(0,0)\) jest ściekiem a \(X_{**}=(1,\frac{1}{2b})\) jest połowicznie stabilny
  2. Jeśli a,b<0: \(X_{*}=(0,0)\) jest źródłem a \(X_{**}=(1,\frac{1}{2b})\) jest połowicznie stabilny
  3. Jeśli a>0 i b<0: \(X_{*}=(0,0)\) jest siodłem a \(X_{**}=(1,\frac{1}{2b})\) jest siodłem
  4. Jeśli a<0 i b>0: \(X_{*}=(0,0)\) jest siodłem a \(X_{**}=(1,\frac{1}{2b})\) jest siodłem

Dla a,b>0 \[\lambda_{1,2} = \frac{-a-b \pm \sqrt{(-a-b)^2}}{2}\] \[\lambda_{1} = -a-b,\;\lambda_{2} = 0,\;\text{zero-eigen}\]

Obliczanie wektorów własnych: \[\begin{pmatrix} -a+a+b & 1\\ 0.5 & -b+a+b\end{pmatrix}\dbinom{v_{x}}{v_{y}}=0\] \[bv_{x} + v_{y} = 0\] \[V=\dbinom{1}{-b}\]

\[\begin{pmatrix} -a-0 & 1\\ 0.5 & -b-0\end{pmatrix}\dbinom{w_{x}}{w_{y}}=0\] \[-aw_{x} + w_{y} = 0\] \[W=\dbinom{1}{a}\]

Zachowanie układu wzdłuż wektora W jest podobne to zachowania poprzedniego układu wzdłuż osi X. Wartość własna wzdłuż tych wektorów jest równa 0. To jest cechą charakterystyczną bifurkacji, podczas których są tworzone lub anihilowane punkty stałe. Do tych bifurkacji zalicza się siodło-węzeł, transkrytyczną i widłową.

\[\Delta>0\] \[1 - 4a^2b^2>0,\;a<\pm\frac{1}{2b},\;dla\; b=0.5,\;a=0.5\] \[\frac{1}{4}x^2 - x + \frac{1}{4}=0,\;x^2 - 4x + 1\] \[x_{1} = \frac{4-\sqrt{12}}{2}=2-\sqrt{3},\;x_{2}=2+\sqrt{3}\] \[-\frac{1}{2}x + y=0,\;y_{1}=1-\frac{\sqrt{3}}{2},\;y_{2}=1+\frac{\sqrt{3}}{2}\]

  1. a>0,b>0 (I ćwiartka układu współrzędnych kartezjańskich)

\[a<\frac{1}{2b},a=0.5,\;b=0.5\] \[A(2-\sqrt{3},1-\frac{\sqrt{3}}{2}) = \begin{pmatrix} -a & 1\\ 0.4665 & -b\end{pmatrix}=ab-0.4665\] \[a=0.5,b=0.5,b>0,a>0\; A=-0.2165,\;\text{siodło}\] \[A(2+\sqrt{3},1+\frac{\sqrt{3}}{2}) = \begin{pmatrix} -a & 1\\ 0.0335 & -b\end{pmatrix}=ab-0.0335\] \[a=0.5,b=0.5,b>0,a>0\; A=0.2165,\;\tau=-1,\;\text{ściek}\]

  1. a<0,b<0 (IV ćwiartka układu współrzędnych kartezjańskich)

\[-a<\frac{1}{2b},\;a>\frac{1}{2b},\;a=-0.5,\;b=-0.5\] \[x_{1} = \frac{4-\sqrt{12}}{2}=2-\sqrt{3},\;x_{2}=2+\sqrt{3}\] \[\frac{1}{2}x + y=0,\;y_{1}=-1+\frac{\sqrt{3}}{2},\;y_{2}=-1-\frac{\sqrt{3}}{2}\]

\[a=-0.5,b=-0.5,\; A_{1}=-0.2165,\;\text{siodło}\] \[a=-0.5,b=-0.5,\; A_{2}=0.2165,\;\tau=1,\;\text{źródło}\]

example811_null14.gif

Figure 132: Analiza izoklin zerowego rzędu dla I i IV ćwiartki układu współrzędnych kartezjańskich

par a=1
par b=0.5

x'=-a*x + y
y'=x^2/(1+x^2) - b*y

init x=0.3
init y=0.1

@ total=20,dt=0.05,xp=x,yp=y
@ xlo=0, ylo=-2, xhi=2 yhi=2
@ rangeover=a, rangestep=249, rangelow=1, rangehigh=-1
done
  1. a>0,b<0 (III ćwiartka układu współrzędnych kartezjańskich)

\[a<-\frac{1}{2b},a=0.5,\;b=-0.5\] \[-\frac{1}{4}x^2 - x - \frac{1}{4}=0,\;x^2 + 4x + 1\]

\[x_{1} = \frac{-4-\sqrt{12}}{2}=-2-\sqrt{3},\;x_{2}=-2+\sqrt{3}\] \[-\frac{1}{2}x + y=0,\;y_{1}=-1-\frac{\sqrt{3}}{2},\;y_{2}=-1+\frac{\sqrt{3}}{2}\]

\[a=0.5,b=-0.5,A_{1}=-0.2165,\;\text{siodło}\] \[a=0.5,b=-0.5,\; A_{2}=0.2165,\;\tau=0\;\text{środek?}\] \[dla\; a \in (-\frac{1}{2b},-b),\; \tau>0 - \text{stabilny wir}\] \[dla\; a \in (-b,0),\; \tau<0 - \text{niestabilny wir}\]

Gdy a i b są przeciwnego znaku i równe co do wartości bezwzględnej to pojawia się i znika oscylacja. To jest nowy typ bifurkacji, z którym nie mieliśmy do czynienia w układach jednowymiarowych.

  1. a<0,b>0 (II ćwiartka układu współrzędnych kartezjańskich)

\[-a<-\frac{1}{2b},\;a>-\frac{1}{2b},\;a=-0.5,\;b=0.5\] \[x_{1} = \frac{-4-\sqrt{12}}{2}=-2-\sqrt{3},\;x_{2}=-2+\sqrt{3}\] \[\frac{1}{2}x + y=0,\;y_{1}=1+\frac{\sqrt{3}}{2},\;y_{2}=1-\frac{\sqrt{3}}{2}\]

\[a=-0.5,b=0.5,A_{1}=-0.2165,\;\text{siodło}\] \[a=0.5,b=-0.5,\; A_{2}=0.2165,\;\tau=0\;\text{środek?}\] \[dla\; a \in (-\frac{1}{2b},-b),\; \tau>0 - \text{niestabilny wir}\] \[dla\; x \in (-b,0),\; \tau<0 - \text{stabilny wir}\]

example811_null23.gif

Figure 133: Analiza izoklin zerowego rzędu dla II i III ćwiartki układu współrzędnych kartezjańskich

par a=-1
par b=0.5

x'=-a*x + y
y'=x^2/(1+x^2) - b*y

init x=-0.3
init y=0.1

@ total=20,dt=0.05,xp=x,yp=y
@ xlo=0, ylo=-2, xhi=2 yhi=2
@ rangeover=a, rangestep=249, rangelow=-1, rangehigh=1
done

Jeśli minimum jeden parametr będzie równy zero to: \[x(x - ab(1+x^2))=0,\;x^2=0\]

W układzie będzie zawsze jeden punkt stały X* = (0,0) \[A(0,0) = \begin{pmatrix} -a & 1\\ 0 & -b\end{pmatrix},\;det(A)=ab,\;\tau=-a-b\] Na postawie analizy pola wektorowego możemy się przekonać, że ten punkt jest siodłem.

Podsumowanie: Gdy \(\Delta>0\) w układzie mamy trzy punkty stałe:

  1. Jeśli a,b>0:
    \(X_{0}=(0,0)\) jest ściekiem a z punktu połowicznie stabilnego \(X_{0}^*=(1,\frac{1}{2b})\) pojawiają się dwa punkty \(X_{0}^{**}=(2-\sqrt{3},1-\frac{\sqrt{3}}{2})\) - siodło i \(X_{0}^{***}=(2+\sqrt{3},1+\frac{\sqrt{3}}{2})\) - ściek. Rozmaitość niestabilna siodła z jednej strony łączy się z punktem X0 a z drugiej z punktem X0*.
  2. Jeśli a,b<0:
    \(X_{0}=(0,0)\) jest źródłem a z punktu połowicznie stabilnego \(X_{0}^*=(1,\frac{1}{2b})\) pojawiają się dwa punkty \(X_{0}^{**}=(2-\sqrt{3},-1+\frac{\sqrt{3}}{2})\) - siodło i \(X_{0}^{***}=(2+\sqrt{3},-1-\frac{\sqrt{3}}{2})\) - źródło. Rozłożenie rozmaitości siodła ulega zmianie (to już widać w punktach połowiczne stabilnych): teraz stabilna łączy się z punktem X0 a z drugiej strony z punktem X0*, gdyż X0* jest źródłem.
  3. Jeśli a>0 i b<0:
    \(X_{0}=(0,0)\) jest siodłem a z siodła \(X_{0}^*=(1,\frac{1}{2b})\) pojawiają się dwa punkty \(X_{0}^{**}=(-2-\sqrt{3},-1-\frac{\sqrt{3}}{2})\) - siodło i \(X_{0}^{***}=(-2+\sqrt{3},-1+\frac{\sqrt{3}}{2})\) - środek dla a=b, stabilny wir dla $a ∈ \((-\frac{1}{2b},-b)\) i niestabilny wir dla \(a \in (-b,0)\) .
  4. Jeśli a<0 i b>0: \(X_{0}=(0,0)\) jest siodłem a z siodła \(X_{0}^*=(1,\frac{1}{2b})\) pojawiają się dwa punkty \(X_{0}^{**}=(-2-\sqrt{3},1+\frac{\sqrt{3}}{2})\) - siodło i \(X_{0}^{***}=(-2+\sqrt{3},1-\frac{\sqrt{3}}{2})\) - środek dla a=b, stabilny wir dla $a ∈ \((-\frac{1}{2b},-b)\) i niestabilny wir dla \(a \in (-b,0)\) .

4.5.3 Bifurkacja transkrytyczna

\[\dot x = \mu x - x^2\] \[\dot y = - y\]

bif_transcritical_nullcline.png

Figure 134: Analiza bifurkacja na podstawie izloklin zerowego wzrostu

4.5.4 Bifurkacja widłowa nadkrytyczna

\[\dot x = \mu x - x^3\] \[\dot y = - y\]

bif_pitchfork_sup_nullcline.png

Figure 135: Analiza bifurkacja na podstawie izloklin zerowego wzrostu

4.5.5 Bifurkacja widłowa podkrytyczna

\[\dot x = \mu x + x^3\] \[\dot y = - y\]

bif_pitchfork_sub_nullcline.png

Figure 136: Analiza bifurkacja na podstawie izloklin zerowego wzrostu

4.5.6 Przykład: widłowa nadkrytyczna

\[\dot x = \mu x + y + \sin{x}\] \[\dot y = x - y\]

Pytanie: czy w punkcie \(X_{*}=(0,0)\) zachodzi bifurkacja widłowa nadktytyczna?
\[A = \begin{pmatrix} \mu + \cos{x} & 1\\ 1 & -1\end{pmatrix}\] \[A(0,0) = \begin{pmatrix} \mu + 1 & 1\\ 1 & -1\end{pmatrix},\;det(A)=-\mu-2,\;\tau=\mu\] \[-\mu-2 > 0,\; \mu < -2 \rightarrow \text{ściek}\] \[-\mu-2 < 0,\; \mu > -2 \rightarrow \text{siodło}\] \[-\mu-2 = 0,\; \mu = -2 \rightarrow \text{bifurkacja zero-eigen}\]

bif_pitchfork_sup.gif

Figure 137: Bifurkacja widłowa nadkrytyczna

par mu=-2.5

x'=mu*x + y + sin(x)
y'=x-y

@ total=10,dt=0.01,xp=x,yp=y, xlo=-3,ylo=-3, xhi=3 yhi=3
@ rangeover=mu, rangestep=249, rangelow=-2.5, rangehigh=-1.5
done
#GUI
# Dir.filed/flow —> Flow –> 16
# Nullcline —> New
# Initialconds —> Range —> Movie:Yes —> ENTER
# Kinescope —> Make AniGif

\[\lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2 - 4det(A)}}{2}\] \[\lambda_{1,2} = \frac{\mu \pm \sqrt{\mu^2 - 0}}{2}\] \[\lambda_{1} = \mu,\;\lambda_{2}=0\]

Obliczanie wektorów własnych: \[\begin{pmatrix} \mu+1-\mu & 1\\ 1 & -1 -\mu\end{pmatrix}\dbinom{v_{x}}{v_{y}}=0\] \[v_{x} + v_{y} = 0\] \[V=\dbinom{1}{-1}\]

Obliczanie wektorów własnych: \[\begin{pmatrix} \mu+1-0 & 1\\ 1 & -1 -0\end{pmatrix}\dbinom{w_{x}}{w_{y}}=0\] \[w_{x}(\mu+1) + w_{y} = 0\] \[V=\dbinom{1}{-\mu-1}\]

W punkcie bifurkacji kierunek ruchu jest wzdłuż wektora V, gdyż tylko dla tego wektora wartość własna jest różna od zera.

bif_pitchfork_sup_zeroeigen.png

Figure 138: Zachowanie układu w okolicach punktu stałego dla μ=-2

4.5.7 Bifurkacja Hopf nadkrytyczna

\[\dot r = \mu r - r^3\] \[\dot \theta = \omega + br\]

W układzie współrzędnych polarnych łatwiej jest dostrzec cykl graniczny. Podczas bifurkacji Hopfa nadkrytycznej zachodzi zmiana z wiru stabilnego na niestabilny i nad tym wirem niestabilnym pojawia się stabilny cykl graniczny. Może pojawić się także środek, ale wtedy mówimy, że bifurkacja jest zdegenerowana.

W tym przypadku widzimy, że mamy dwa punkty stałe: \(r_{0}=0,\; r_{0}^* = \sqrt{\mu}\)
Już na tym etapie możemy dostrzec orbitę zamkniętą dla \(r=\sqrt{\mu}\). Dla \(r > 0\) i \(r < \sqrt{\mu}\) pochodna jest dodatnia, a dla \(r > \sqrt{\mu}\) ujemna, co oznacza, że nasz cykl graniczny jest stabilny.

Aby zbadać, jak zmieniają się wartości własne, przejdźmy do układu współrzędnych kartezjańskich. \[\dot x = \dot r \cos{\phi} - \dot \theta r\sin{\phi}\] \[\dot y = \dot r \sin{\phi} + \dot \theta r\cos{\phi}\]

Podstawiając \(\dot r\) i \(\dot \theta\) do równań otrzymamy (końcowy rezultat): \[\dot x = x(\mu-x^2-y^2) - y(\omega +b\sqrt{x^2+y^2})\] \[\dot y = y(\mu-x^2-y^2) + x(\omega +b\sqrt{x^2+y^2})\]

\[A = \begin{pmatrix} \mu -3x^2-y^2+b\frac{x}{\sqrt{x^2+y^2}} & -2xy-\omega + b\frac{y}{\sqrt{x^2+y^2}}\\ -2xy+\omega - b\frac{x}{\sqrt{x^2+y^2}} & \mu-3y^2-x^2-b\frac{y}{\sqrt{x^2+y^2}}\end{pmatrix}\]

\[A(0,0) = \begin{pmatrix} \mu & -\omega\\ \omega & \mu\end{pmatrix},\;det(A)=\mu^2+\omega^2,\;\tau=2\mu\]

\[\lambda_{1} = \frac{2\mu + \sqrt{4\mu^2 - 4\mu^2 - 4\omega^2}}{2}\] \[\lambda_{2} = \frac{2\mu - \sqrt{4\mu^2 - 4\mu^2 - 4\omega^2}}{2}\]

\[\lambda_{1} = \mu + i\omega,\;\lambda_{2} = \mu - i\omega\]

Bifurkacja zachodzi, gdy część rzeczywista wartości własnych się zeruje, czyli dla \(\mu=0\). Jest to Bifurkacja Hopfa nadkrytyczna, gdyż cykl graniczny pojawia się nad wirem niestabilnym.

bif_hopf_sup.gif

Figure 139: Bifurkacja Hopfa nadkrytyczna

par mu=1
par omega=1
par b=0.1
x'=x*(mu-x^2-y^2) - y*(omega+b*sqrt(x^2+y^2))
y'=y*(mu-x^2-y^2) + x*(omega+b*sqrt(x^2+y^2))

init x=0.5, y=0

@ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2, yhi=2
@ rangeover=mu, rangestep=249, rangelow=-1, rangehigh=1
done
#GUI
# Dir.filed/flow —> Flow –> 16
# Nullcline —> New
# Initialconds —> Range —> Movie:Yes —> ENTER
# Kinescope —> Make AniGif

4.5.8 Bifurkacja Hopf podkrytyczna

\[\dot r = \mu r + r^3\] \[\dot \theta = \omega + br^2\]

Bifurkacja Hopfa nadkrytyczna zachodzi zmiana z wiru niestabilnego na stabilny i nad tym wirem stabilnym pojawia się niestabilny cykl graniczny.

W tym przypadku widzimy, że mamy dwa punkty stałe: \(r_{0}=0,\; r_{0}^* = \sqrt{-\mu}\)
Tak jak w poprzednim przykładzie, już na tym etapie możemy dostrzec orbitę zamkniętą dla \(r=\sqrt{-\mu},\mu<0\). Dla \(r > 0\) i \(r < \sqrt{-\mu}\) i \(\mu<0\) pochodna jest ujemna a dla \(r > \sqrt{\mu}\) jest dodatnia, co oznacza, że nasz cykl graniczny jest niestabilny.

Aby zbadać, jak zmieniają się wartości własne przejdźmy do układu współrzędnych kartezjańskich. \[\dot x = \dot r \cos{\phi} - \dot \theta r\sin{\phi}\] \[\dot y = \dot r \sin{\phi} + \dot \theta r\cos{\phi}\]

Podstawiając \(\dot r\) i \(\dot \theta\) do równań i zastosowaniu przekształceń otrzymamy: \[\dot x = x(\mu+x^2+y^2) - y(\omega +b(x^2+y^2))\] \[\dot y = y(\mu+x^2+y^2) + x(\omega +bx^2(x^2+y^2))\]

\[A = \begin{pmatrix} \mu + 3x^2+y^2-2byx & 2xy-\omega - bx^2y - 3by^2\\ 2xy+\omega + 5bx^4 + 3by^2x^2 & \mu + x^2 +3y^2+2bx^3y\end{pmatrix}\]

\[A(0,0) = \begin{pmatrix} \mu & -\omega\\ \omega & \mu\end{pmatrix}=\mu^2+\omega^2,\;\tau=2\mu\]

\[\lambda_{1} = \frac{2\mu + \sqrt{4\mu^2 - 4\mu^2 - 4\omega^2}}{2}\] \[\lambda_{2} = \frac{2\mu - \sqrt{4\mu^2 - 4\mu^2 - 4\omega^2}}{2}\]

\[\lambda_{1} = \mu + i\omega,\;\lambda_{2} = \mu - i\omega\]

Bifurkacja zachodzi gdy część rzeczywista wartości własnych się zeruje, czyli dla \(\mu=0\). Jest to Bifurkacja Hopfa podkrytyczna, gdyż cykl graniczny pojawia się nad wirem stabilnym.

bif_hopf_sub.gif

Figure 140: Bifurkacja Hopfa podkrytyczna

par mu=-0.001
par omega=1
par b=1

x' = x*(mu+x^2+y^2) - y*(omega +b*(x^2+y^2))
y' = y*(mu+x^2+y^2) + x*(omega +b*x^2*(x^2+y^2))

init x=0.05, y=0

@ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2, yhi=2
@ rangeover=mu, rangestep=249, rangelow=-0.01, rangehigh=0.01
done
#GUI
# Dir.filed/flow —> Flow –> 16
# Nullcline —> New
# Initialconds —> Range —> Movie:Yes —> ENTER
# Kinescope —> Make AniGif

4.5.9 Przykład: bifurkacja Hopfa nadkrytyczna

\[\dot x = \mu x-y+xy^2\] \[\dot y = x+ \mu y+y^3\]

Pytanie: czy w punkcie \(X_{*} = (0,0)\) zachodzi bifurkacja Hopfa nadkrytyczna?

\[A = \begin{pmatrix} \mu + y^2 & -1 + 2xy\\ 1 & \mu+3y^2\end{pmatrix}\] \[A(0,0) = \begin{pmatrix} \mu & -1\\ 1 & \mu\end{pmatrix},\;det(A)=\mu^2+1,\;\tau=2\mu\] \[\forall \mu, det(A(0,0)) > 0\] \[\mu>0,\; \text{niestabilny wir}\] \[\mu<0,\; \text{stabilny wir}\] \[\mu=0,\; \text{pukt bifurkacji Hopfa?}\]

\[\lambda_{1} = \frac{2\mu + \sqrt{4\mu^2 - 4\mu^2 - 4}}{2}\] \[\lambda_{2} = \frac{2\mu - \sqrt{4\mu^2 - 4\mu^2 - 4}}{2}\]

\[\lambda_{1} = \mu + i,\;\lambda_{2} = \mu - i\]

Na tym etapie analizy wiemy, że punkt μ=0 jest punktem bifurkacji podczas której zmienia się stabilność wira. Nie wiemy, czy nad stabilnym bądź niestabilnym wirem jest cykl graniczny.

W celu sprawdzenie, czy istnieje cykl graniczny przejdziemy do układu współrzędnych polarnych. \[x^2+y^2=r^2,\;\tan{\phi}=\frac{y}{x}\] \[x\dot x + y\dot y = r\dot r,\;\dot \theta=\frac{\dot y x - y \dot x}{x^2+y^2}\] Po podstawieniu i wykonaniu rachunków dostaniemy:
\[\dot r = \mu r + r^3 \sin^2{\phi},\;\dot \theta = r\]

Pytanie: czy na podstawie tych równań widać, że w układzie jest cykl graniczny?

Sprawdźmy, czy znajdziemy region pułapkowania:
\[\mu r_{max} + r_{max}^3 \sin^2{\phi} <0\] \[r_{max} < \sqrt{- \frac{\mu}{\sin^2{\phi}}}\]

Dla \(\mu >0\) to równanie nie ma rozwiązania. Ponieważ dla \(\mu>0\) mamy niestabilny wir to z poprzedniego faktu wynika, że nie jest możliwa bifurkacja Hopfa nadkrytyczna. Pochodna zawsze dodatnia \(\dot r = \mu r + r^3 \sin^2{\phi}>0\)

Nie jest to także bifurkacja Hopfa zdegenerowana, gdyż dla \(\mu = 0 \rightarrow r=0\), pochodna zawsze dodatnia \(\dot r = r^3 \sin^2{\phi}>0\)

Pozostaje nam do sprawdzenia, czy jest to Hopf podkrytyczny, czyli szukamy niestabilnego cyklu granicznego dla \(\mu<0\):
\[r_{max} < \sqrt{\frac{1}{\sin^2{\phi}}},\sin^2{\phi} \in (0,1)\] \[r_{min} > \sqrt{\frac{1}{\sin^2{\phi}}},\sin^2{\phi} \in (0,1)\] \[\phi \rightarrow 0,\; r_{min} \rightarrow \infty\]

Nasze równanie tego nie wyklucza, ale także nie dowodzi. Sprawdzamy numerycznie.

bif_hopf_sub.png

Figure 141: Bifurkacja Hopfa nadkrytyczna

5 XPPAut

5.1 Podstawowe elementy skryptów XPP

# Deklaracja stałych wspłóczynników
# współczynniki z możliwością modyfikacji w GUI
par a=1
par b=1
par c[1..2] # c1=0, c2=0

# współczynniki bez możliwości modyfikacji w GUI (ukryte) 
number z[1..2]=.123456 
number p01 = 1
number p10 = 2

# stała wielkość akutalizowana za każdym razem jak współczynniki zostaną zmienione 
# wartość tak zdefiniowanej wielkości można sprawdzić File--> Calculator--> 'nazwa_stałej' 
!pole=a*b*z1

# współczynniki prcesu Wienera (ruch Browna, proces stochastyczny)
# unormowany szum biały
wiener w[1..2]
# suma użytych w progamie stałe współczynników 
# jest pokazywana w konsoli po ruchommieniu programu:
# „Used X constatns“ 
#######################################################################

# zmienne: mogą być ciągłe lub markova

# zmienna markova; informacja w konsoli o liczbie tych zmiennych: nmark 
markov z 2
{0} {P01}
{P10} {0}

# zmienna ustalona; informacja w konsoli o liczbie tych zmiennych: nfix
# bez możliwości dostępu w GUI 
# tak zdefiniową zmienną możemy wykorzystac do zdefiniowania innych zmiennych
energy=0.5*y^2 -0.5*x^2 + 0.25*x^4

# zmienna pomocnicza; informacja w konsoli o liczbie tych zmiennych: naux
# dostęp w GUI
aux e1=energy
aux e2=0.5*y^2 -0.5*x^2 + 0.25*x^4

# definicja równań; informacja w konsoli o liczbie tych zmiennych: nvar
x'=-y
y'=x

# inicjacja zmiennych 
init x=0, y=0

# informacja w konsoli o liczbie równań: neq = nvar + naux + nmark
# W GUI informacja o równaniach jest dostepna pod zakładką „Eqns“
# informacja w konsoli o liczbie węzłów: node = nvar + naux + nfix

@ total=10,xp=x,yp=y, xlo=-2,ylo=-2, xhi=2
@ yhi=2, xmax=2,xmin=-2, ymax=2,ymin=-2
@ maxstor=2000000, bounds=100000,  nmesh=800 
done

5.2 Podstawowe opcje XPP modyfikowane za pomocą skryptów

# zakres numerycznego całkowania wraz z długością kroku tj. w tym
# przypadku liczymy total/dt = 200 punktów 
@ total=10, dt=0.05

# maksymalna liczba kroków przechowywana w pamięci oraz dopuszczalna wartość
# obliczeń. Gdy zostanie osiągnięta maksymalna wartość określona w ‚bounds‘ to
# dostaniemy komunikat: „Y out of bounds at t=a“; w konsoli informacje o
# ostatnich obliczonych wartościach zmiennych: np dla bounds=10 i init x=1, y=15
# variable  f(t-1)   f(t)
# X         1        1.67437 
# Y         15       14.661 
@ maxstor=2000000 {Intiger}, bounds=100000 {Intiger}
# maxstor=1 - maksymalna liczba kroków przechowywana w pamięci to 1: w zakładce ‚Data‘
# maxstor=20 - maksymalna liczba kroków przechowywana w pamięci to 20: w zakładce ‚Data‘  

# Metoda całkowania numerycznego
@ meth={discrete,euler,modeuler,rungekutta,adams,gear,volterra, backeul, qualrk,stiff,cvode,5dp,83dp,2rb,ymp} 

# całkowanie numeryczne od razu po uruchomieniu xpp (równoznaczne z
# Initialaconds—>Go)
@ runnow={0,1}

# dokładność graficzna izloklin zerowego rzędu. Izokliny są obliczane w taki
#  sposób, że przestrzeń fazowa jest dzielona na siatkę o dokładności określonej
#  przez ~nmesh~. W kążdym punkcie siatki jest oblicznane pole
#  wektorowe. Wyświetlane są miejsca w siatce, których wartosć jest równa zero.
@ nmesh=800 {Intiger}

# paramtery algorytmu (metoda Newtona) do znajdowania punktów zerowych funkcji
@ jac_eps=value
@ newt_tol=value
@ newt_iter=value

# Całkowanie dla określonego zakresu zmiannych lub parametrów
@ range=1 # do zastosowania w trybie wsadowym ~xppaut -silent~
@ rangeover=name, rangestep, rangelow, rangehigh

# oblicz średnią (1) lub wariancje (2) trajektori określonych w ‚rangeover‘
# równoważne z (1) Numerics—>Stochast—>Mean; (1) Numerics—>Stochast—>Variance
# do zastosowania w trybie wsadowym ~xppaut -silent~
@ stoch={1,2}  

#  Okres w toridalnej przestrzeni fazowej: domyślna wartość 2pi 
# równoznaczne Phasespace—>Choose—> Period: 2pi
@ tor_per=value

# wskazanie, która zmienna jest okresowa tells XPP that the variable ¡name¿ is to be considered
# modulo the period. You can repeat this for many variables.
# równoznaczne Phasespace—>Choose—>Period: 2pi: ENTER Fold with: name
@ fold=name

# Automatycznie odświeżanie wartości w tabeli po zmianie paramerów w AUTO
@ autoeval={0,1} 

# sets up a Poincare map for either sections of a variable or the extrema.
@ poimap={ section,maxmin} 
# sets the variable name whose section you are interested in finding.
@ poivar=name 
# is the value of the section; it is a floating point.
@ poipln=value 
# determines the direction of the section.
@ poisgn={ 1, -1, 0 } 
# means to stop the integration when the section is reached.
@ poistop=1 

# ustawiamy osie plota oraz ich zakres - otrzymamy wykres zależności pomiędzy x i y 
@ xp=x,yp=y, xlo=-2,ylo=-2, xhi=2, yhi=2

# Ilość wykresów na głównym ekranie oraz określenie nazw poszczególnych osi; maksymalnie 8 wyresów 
@ nplot=2 {Intiger}, xp2=t,yp2=v 

# tworzenie przycisku na górnym panelu: ‚fit‘ to nazwa a ‚wf‘ to działanie przycisku określone na
# podstawie skrótów: ‚wf‘ znaczy ‚Window‘ ‚Fit‘ 
@ but=fit:wf 

# zamiast do konosli zapisz komunikaty programu do pliku 
@ logfile=nazwa_pliku 

# wyłączenie/włączenie dzwięków do sygnalizacij wykonania poszczególnych zadań 
@ bell={0,1}

# wypisuj/nie wypisuj komunikaty programu do konsoli bądź pliku
@ quiet={0,1}

# wrażenie trójwymiarowości przycisków (wył/włącz)
@ grads={0,1}

5.3 Przykładowy plik konfiguracyjny .xpprc

@ maxstor=50000,bell=0,quiet=0
@ BIG=*adobe-helvetica-medium-r-normal*12*
@ SMALL=*adobe-helvetica-medium-r-normal*8*
@ grads=0,mwcolor=cecece,dwcolor=C6D6FF
@ backcolor=cecece,forecolor=000000
@ height=8, width=13

5.4 Podstawowe opcje programu AUTO, które można modyfikować za pomocą skryptów

#  This is the number of mesh intervals for discretization of periodic orbits. If
#  you are getting bad results or not converging, it helps to increase this. For
#  following period doubling bifurcations, this is automatically doubled so you
#  should reset it later.
@ ntst=15 {Intiger}

# Maksymalna liczba punktów obliczeniowych na gałęzi diagramu bifurkacji
@ nmax=200 {Intiger}

# Give complete info every Npr steps. Set this to a big number if you want to speed things along.
@ npr=50 {Intiger}

#  Inicjacja kroku zmiany parametru podczas obliczania bifurkcji. Znak przy
#  wartości określa kierunek zmian. Ds jest zmieniany przez program
#  automatycznie w granicach dsmin i dsmax
@ ds=0.02

# Minimalna i maksymalna długość kroku. 
@ dsmin=0.001, dsmax=0.5

# This is the left-hand/right-hand limit of the diagram for the principle parameter. The
#  calculation will stop if the parameter is less/greater than this.
@ parmin=0, parmax=2 

# The lower/upper bound for the L2 norm of the solution. If it is less/greater than this the
# calculation will stop.
@ normmin=0, normmax=1000 

# Number of collocation points used; 2<=ncol<=7
@ ncol=4

# Some sort of tolerances for AUTO. Keep 'em positive but make 'em smaller for
#  better accuracy. 1e-7 seems to be a popular choice!
@ epsu=0.0000001,epss=0.00001,epsl=0.0000001

# Zakres osi na plocie oraz nazwa osi Y;
@ autoxmin=value, autoxmax=value, autoymin=value, autoymax=value, autovar=variable

5.5 Komunikaty o błędach

  • Out of bounds: Some variable or auxiliary quantity went beyond the prescribed maximum. Fix this by increasing the Bounds in the nUmerics menu. But be careful, as some solutions to ODEs really do blow up.
  • Bad formula: The parser was unable to figure out what you meant. You have probably mistyped a name, forgotten a parenthesis, or used a name not known to XPPAUT.
  • Singular Jacobian In Newton's method, the matrix was not invertible: In BVPs, this sometimes means that you have not written the boundary conditions correctly. In integration, sometimes setting the time step to be smaller works.
  • Maximum iterates exceeded: Again, in Newton's method or in curve fitting, the specified tolerances could not be achieved in the given number of iterates. Either loosen tolerances or increase the maximal number of iterates.
  • Delay negative or too large: Increase the maximal allowable delay. If the delay is negative, you have an ill posed problem.
  • Could not converge to root: In solving a nonlinear equation, the Newton solver failed to converge. You should make a better guess.
  • Working too hard???: This error occurs when you have global flags. It is a difficult error to fix. It generally means that when condition A occurs it induces condition B which in turn induces condition A again, and so on. Clicking Escape many times will sometimes work, but occasionally you need to kill XPPAUT.
  • All curves used: When you freeze a plot, there are only 26 per graph window. Delete some.
  • Bad condition: Ignoring occurs in the computation of histograms. If the formula for the condition is bad, XPPAUT just ignores it.
  • Out of memory: Time to upgrade!
  • Out of film: In the Kinescope command, only so many screen shots can be saved.
  • No prior solution: You tried the Initialconds Last command without first computing a solution with new initial data.
  • No shooting data available: You have to find a fixed point to your system with a onedimensional invariant manifold before you can use this.
  • Incompatible parameters: You tried to load a saved set file that is not associated with the ODE file you are running.

Źródło: Ermentrout B., (2002), Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students, SIAM, ISBN: 978-0-89871-506-4

5.6 AUTO

Znaczenie skrótów

  • EP Endpoint of a branch
  • LP Limit point or turning point of a branch
  • TR Torus bifurcation from a periodic
  • PD Period doubling bifurcation
  • UZ User defined function
  • MX Failure to converge
  • BP Bifurcation or branch point

6 Metody numeryczne

6.1 Metoda Newtona: do znajdowania punktów zerowych funkcji

Metoda Newtona dla układu zdefiniowanego jak poniżej daje wynik x0=1.0003

x'=x^2-1
init x=.8
@ total=100, xlo=0,ylo=-2,xhi=5,yhi=2
done 

# GUI
# Initialconds—>Go
# Numerics—> sing pt ctrl—> Maximum interates: 2 ENTER Newton tolerance: 0.05 ENTER Jacobian essilon:0.001
# Sing Pts—>Go

Działanie algorytmu jest zobrazowane na poniższym rysunku:

set terminal pngcairo enhanced color dashed font "Alegreya, 14" rounded size 32.00cm, 19.20cm
f(x)= 4*x- 5 < 3.2 && 4*x-5 > 0 ? 4*x-5 : NaN 
g(x)= 5.0/2*x - 41.0/16 < 0.8 && 5.0/2*x - 41.0/16 > 0 ? 5.0/2*x - 41.0/16 : NaN 
h(x)= 41.0/20*x - 2.05 < 0.5 && 41.0/20*x - 2.05 > 0 ? 41.0/20*x - 2.05 : NaN 
set arrow 1 from 2,0 to 2,3 nohead dashtype '-' lc black lw 2
set arrow 2 from 5.0/4,0 to 5.0/4,0.6 nohead dashtype '-' lc black lw 2
set arrow 3 from 41.0/40,0 to 41.0/40,0.25 nohead dashtype '-' lc black lw 2
set arrow 4 from 1.0003,0 to 1.0003,0.25 nohead dashtype '-' lc black lw 2
set label 1 'x_{0}=2' at 1.95,-0.2 font "Alegreya, 10"
set label 2 'x_{1}=1.25' at 1.20,-0.2 font "Alegreya, 10"
set label 3 'x_{2}=1.025' at 1.020,-0.2 font "Alegreya, 10"
set label 4 'x_{3}=1.0003' at 0.9,0.25 font "Alegreya, 10"
plot [0.9:2.1] [-1:3.1] x**2 - 1 ls 1 notitle,\
   f(x) dashtype '-' lw 2 lc black notitle,\
   g(x) dashtype '-' lw 2 lc black notitle,\
   h(x) dashtype '-' lw 2 lc black notitle

newton_method.png

Figure 142: Metoda Newtona

Styczna do krzywej \(x^2 -1\) w punkcie xn:
\[f(x_{n+1})-f(x_{n}) = f'(x_n)(x_{n+1}-x_n)\] \[\text{Dla}\; x_{0} = 2:\] \[y-3 = 4(x_{1}-2),\;y=4x_{1}-5,\;4x_{1}-5=0,\;x_{1}=\frac{5}{4}=1.25\] \[y-\frac{9}{16} = \frac{5}{2}(x_{2}-\frac{5}{4}),\;y=\frac{5}{2}x_{2}-\frac{41}{16}=0,\;x_{2}=\frac{41}{40}=1.025\] \[y-0.05=\frac{41}{20}(x_{3}-\frac{41}{40}),\;y=\frac{41}{20}x_{2}-2.05=0,\;x_{3}=1.0003\]

NewtonIteration_Ani.gif

Figure 143: Metoda Newtona (źródło: Wikipedia)

7 Narzędzia

7.1 Konwersja gif to mp4:

ffmpeg -i anim.gif -movflags faststart -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" video.mp4

Źródło: https://unix.stackexchange.com/questions/40638/how-to-do-i-convert-an-animated-gif-to-an-mp4-or-mv4-on-the-command-line

7.2 Łączenie zdjęć

# Łączenie horyzontalne - od lewej do prawej
convert img1.png img2.png +append img_top.png
convert img3.png img4.png +append img_bottom.png
# łączenie wertykalne - od góry do dołu
convert img_top.png img_bottom.png -append img_all.png

Autor: Michał Pierzchalski

Created: 2021-05-25 wto 19:46

Validate