Analiza układów dynamicznych

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

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.

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 x₀. Jeśli \(a > 0\) to układ jest niestabilny i punktem stałym jest \(x_{*}=0\) 2.

Rozwiązanie analityczne:

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

Zachowanie układu zależy od warunku początkowego x₀. 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\).

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 x₀ 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 x₀ = -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 x₀=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}\]

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.

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 x² - 1 jest stabilny lokalnie tzn. dla \(x \leq 1\)

\[\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})}}\]

• x₀ - wartość zmiennej niezależnej, dla której f(x₀) = 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))\]

Rozwiązanie analityczne:

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

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 N₀ od -2 do 4. Możemy się przekonać, że nasze obliczenia punktów stałych są poprawne.

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

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

Punkt N₁ = 0 jest punktem niestabilnym, gdyż wektor \(V_{x}=\frac{dx}{dt}\) jest skierowany od tego punktu; punkt N₂ = 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 N₂ = 2.

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:

using Plots
x = -2.2:0.1:2.2
p1 = plot(x,x.^2 .- 1)
p1 = hline!([0])
p2 = plot(x,x.*(1 .- x./2))
p2 = hline!([0])
p3 = plot(x,x.^2)
p3 = hline!([0])
p4 = plot(x,x.^3)
p4 = hline!([0])
plot(p1,p2,p3,p4,layout=(2,2),legend=false)
savefig("dyn_linear_analysis.png")

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ć:

using Plots
x = -3:0.1:3.5
s1 = x .* (1 .- 1/3 .* x.^(2))
s2 = x.^2 .* (x./(3*2)  .- 1/2)
p1 = plot(x,s1,yaxis=("\$V_1\$",(-2.5,2.5)))
p1 = hline!([0])
p2 = plot(x,s2,yaxis=("\$V_2\$",(-3,1)))
p2 = hline!([0])
plot(p1,p2,layout=(1,2),legend=false)
savefig("dyn_potential.png")

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")

for i in [1, 4, 0] 
    println(i) 
end

for i in [1, 4, 0] println(i) end

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] + k3 - (y[i-1,j] + 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")

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

using Plots
x = -5:0.1:5
p1 = plot(x,x.^2 .- 1)
p1 =   hline!([0])
p2 = plot(x,x.^2 .- 0)
p2 =   hline!([0])
plot(p1,p2,layout=(1,2),legend=false)
savefig("dyn_bifurcation.png")

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:

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\) (prototyp, postać normalna) ma zawsze bifurkacje siodło-węzeł. Na przykład, taką postać znajdziemy, szukając bifurkacji dla 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.

using Plots
x = -2:0.1:2
r = 1 #dla x=0 ta wartosc daje punkt stały
s1 = (r - 1) .- x.^2 ./ 2
s2 = (r - 1) .- x.^2 ./ 2 .+ 1/6 .* x.^3
s3 = r .- x .- exp.(-x)
p1 = plot(x,s1,label="T2")
p1 = plot!(x,s2,label="T3")
p1 = plot!(x,s3,label="f(x)")
plot(p1,layout=(1,1),legend=true)
savefig("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. Wizualizacja w Geogebrze: https://www.geogebra.org/calculator/duur4dgj

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,\text{x>0}\)
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 \in (-1,0)\); ujemna dla r < -1. Pierwsza pochodna badanej funkcji jest nieokreślona dla x=0.

using Plots
x = 0.1:0.01:3
x0=1.0
r = -1
s1 = -0.5*r.*x.^2 .+ 2*r.*x .- 1.5*r .+ x .- 1
s2 = r .* log.(x) .+ x .- 1
p1=plot(x,s1,label="T2")
p1=plot!(x,s2,label="f(x,r=-1)")

r = -2
s1 = -0.5*r.*x.^2 .+ 2*r.*x .- 1.5*r .+ x .- 1
s2 = r .* log.(x) .+ x .- 1
p2=plot(x,s1,label="T2")
p2=plot!(x,s2,label="f(x,r=-2)")

r = 1
s1 = -0.5*r.*x.^2 .+ 2*r.*x .- 1.5*r .+ x .- 1
s2 = r .* log.(x) .+ x .- 1
p3=plot(x,s1,label="T2")
p3=plot!(x,s2,label="f(x,r=1)")

plot(p1,p2,p3,layout=(1,3),legend=true)
savefig("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\). Wizualizacja w Geogebrze: https://www.geogebra.org/calculator/d98t4yeg

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

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.

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.

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.

using Plots
x = -2:0.01:2
x0=1.0
s1 = -1 .* x .- x.^3
s2 = 0 .* x .- x.^3
s3 = 1 .* x .- x.^3
s4 = 2 .* x .- x.^3
xh = sqrt(1/3) - sqrt(1/3)^3
p1=plot(x,s1,label="r=-1",xaxis=("x"),yaxis="dx/dt",title="r=-1")
p1 = hline!([0])
p1 = hline!([xh])
p2=plot(x,s2,label="r=0",xaxis=("x"),yaxis="dx/dt",title="r=0")
p2 = hline!([0])
p2 = hline!([xh])
p3=plot(x,s3,label="r=1",xaxis=("x"),yaxis="dx/dt",title="r=1")
p3 = hline!([0])
p3 = hline!([xh])
p4=plot(x,s4,label="r=2",xaxis=("x"),yaxis="dx/dt",title="r=2")
p4 = hline!([0])
p4 = hline!([xh])

plot(p1,p2,p3,p4,layout=(2,2),legend=false)
savefig("dyn_catastrophe.png")

Wizualizacja portretu fazowego: https://www.geogebra.org/calculator/ceetbwfj

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 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.

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?

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}\]

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.

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

init x=0.5

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

init x=-1.5

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.

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\)

\[\dot \theta = \omega\]

Jednorodność oscylatora wynika z faktu, że prędkość kątowa jest stała i w naszym przypadku jest równa częstotliwości kątowej (szybkości zmian kąta fazowego) \(\omega=2\pi f,\text{gdzie} f=\frac{1}{T}\) (jeden cykl odbywa się przez T). Dwa niezależne od 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ę dudnieniem (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 T₁ i T₂ 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 reprezentuje czas, kiedy jeden dzwon zdubluje w dzwonieniu drugi. Wtedy wielokrotność tej liczby reprezentuje czas, kiedy dzwony wybrzmią razem np.

• T₁ = 3 i T₂ = 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.
• T₁ = 3 i T₂ = 6 \(t = (\frac{1}{3} - \frac{1}{6})^{-1}=6\)

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_{2}\) powrócił do stanu początkowego, a przecięcia z osią Y oznaczają, że 4 razy oscylator \(\theta_{1}\) powrócił do stanu początkowego.

Jeśli okres dzwonienia będzie liczbą wymierną:

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

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

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?

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

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

using Plots, SymPy, LaTeXStrings;
@vars x real=true 
y1 = 2 - sin(x)
y2 = 2 - 2*sin(x)
y3 = 2 - 4*sin(x)
plot(y1,xlims=(-pi,pi),xaxis=(L"\theta"),yaxis=(L"\dot x"),label="a=1")
plot!(y2,label="a=2")
plot!(y3,label="a=4")
hline!([0])
savefig("ex432.png")

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. Gdy wartość parametru a zbliża się 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.

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

1. \(\omega = a,\, \max(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, których istnienie pokazano na poniższym rysunku dla \(\omega =1\) i a=2, to \(sin(\theta)=0.5\), \(\theta=\arcsin(0.5) = \frac{\pi}{6}\), ale również \(\theta = \pi-\frac{\pi}{6}\)

Stabilność: \(f'= -a \cdot cos(\theta)\), \(f'(\pi/6)=-2\cdot cos(\pi/6)<0\), \(f'(5\pi/6)=-2 \cdot cos(5\pi/6)>0\)

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}}\)

using Plots, SymPy;
@vars ω positive=true;
@vars a positive=true;
f(a) = 2*π/sqrt(2^2 - a^2)
plot(f,0,2,ylims=(0,20),xaxis="a",yaxis="T",legend=false)
vline!([2],line=:dot)
savefig("ex432_t.png")

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 (gdy \(a>\omega\)), to dwa punkty stałe łą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.

using Plots, SymPy
f(x) = x^2+0.01
plot(f,-1,1)
savefig("ex432_t2.png")

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\)

using Plots, SymPy, LaTeXStrings;
f(θ) = 1-sin(θ)
g(ϕ) = 1 - 1 + 1/2*ϕ^2
h(θ) = 1 - 1 + 1/2*(θ-π/2)^2
# h(θ) = 1 - 1 + 1/2*(θ^2 - θ*π + (π/2)^2)
label=map(string, styles)
plot(f,label=(L"1-\sin(\theta)"))
plot!(g,label=(L"1-1 + 0.5\theta^2"))
plot!(h,label=(L"1-1 + 0.5(\theta-0.5\pi)^2"))
savefig("ex432_taylor.png") 

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.

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. k₁ = -0.1 i k₂ = -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\)

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

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

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

\(\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})\)

f_i 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: f_i 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 X₀, 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\)

\(\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.

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 link1, link2) - 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 V₀. To jest analogiczne do rozwiązania jednowymiarowego układu dynamicznego. Pierwsza współrzędna v_x odpowiada osi X a v_y 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

\(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}\)

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.

\(\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).

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.

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. https://www.geogebra.org/calculator/zhmzutxw

\[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:

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\),

\[\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.

      

      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 x}{\dot y} = \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?

\[\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).
      

      

      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?

\(\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.
      

      

      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}\)?

\(\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\) i \(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:
      

      

      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

\[m\ddot x + kx =0\] \[\dot x = v\] \[\dot v = -\frac{k}{m}x = -\omega^2 x\] \[\omega=\sqrt{\frac{k}{m}},\text{częstość kątowa oscylatora}\],

\[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:
   

   

   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

\[m\ddot x + c\dot x + kx =0\] \[\dot x = v\] \[\dot v = -\frac{c}{m}v - \frac{k}{m} = -2\zeta\omega_{0} v -\omega_{0}^2 x\] \[\omega_{0}=\sqrt{\frac{k}{m}},\text{częstotść kątowa oscylatora nietłumionego}\], \[\zeta = \frac{c}{2\sqrt{mk}},\text{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?

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)\)
\(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:
   

   

   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

Pytanie: Jak zachowa się układ dla \(\zeta=1\)?

\(\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}{1-\sqrt{2}}\beta\)
\(3 = \alpha + \beta\), \(3 = \frac{\sqrt{2}+1}{1-\sqrt{2}}\beta + \beta\), \(3 = \frac{\sqrt{2}+1+1-\sqrt{2}}{1-\sqrt{2}}\beta\),
\(\beta = \frac{2}{1-\sqrt{2}}\), \(\alpha=\frac{2(\sqrt{2}+1)}{(1-\sqrt{2})^2}\)

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

\[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 Z_re i Z_im 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

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

\[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

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

\(\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

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

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 powtarzających się 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}\]

\[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

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

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

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ę.

Z analizy portretu fazowego 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

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. 

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?

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

\[\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)\]

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\)

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: \mathbb{R}^n \rightarrow \mathbb{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)\).

gdzie F(X) tworzy pole sił, to 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),\text{F to siła}\] \[\dot x = v\] \[\dot v = -\frac{1}{m}\nabla U(x)=-\frac{1}{m}\frac{dU}{dx}\]

Twierdzenie o nieliniowych środkach w układach zachowawczych:

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

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

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.

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.

\[\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 X₀ = (0,0) jest lokalnym minimum U i zgodnie z naszym twierdzeniem wszystkie trajektorie wystarczająco blisko tego punktu są zamknięte. Punkt X₀=(0,0) jest nieliniowym środkiem.

Ponieważ jest to nieliniowy środek w X₀=(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\]

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

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?

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 dodatni 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 pojedynczego punktu stałego. Jeśli \(I_{C}=-1\), to okrąża siodło; dla wszystkich pozostałych izolowanych punktów stałych \(I_{C}=1\)

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

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 krzywych są punktami stałymi.

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:

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

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

\[\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\)

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

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\]

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

\[\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?

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

Twierdzenie Liénarda:

Jeśli:

• funkcje f(x) i g(x) są C¹ 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

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 Rⁿ jest postaci układu równań różniczkowych:
\[\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.

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 R² 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\)

Funkcja Lapunowa jest C¹ 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; wszystkie trajektorie układu podążają w dół, w kierunku stabilnego punktu stałego.

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, gdyż mamy globalny i asymptotyczny punkt stabilny.

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

Sprawdzimy, czy istnieje funkcja Lapunowa, którą dobieramy sobie jako suma 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<0\]

W badanym układzie zatem nie ma orbity zamkniętej.

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 domkniętym i ograniczonym podzbiorem w przestrzeni fazowej
• pole wektorowe w R jest C¹
• 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\)

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

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, YC=0
# Initialconds —> New —> X=-2,V=3, XC=0, YC=0
# Initialconds —> New —> X=2,V=-3, XC=0, YC=0
# Parameters —> r ENTER 3 ENTER
# Initialconds —> New —> X=2,V=0, XC=0, YC=0
# Dir.filed/flow —> Scaled Dir.Fld –> 16

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

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

@ 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=x1, Y-axis=y1, Color=1
# Initialconds —> New —> R=2,THETA=0, XC=0, YC=0
# Parameters —> r ENTER 3 ENTER
# Initialconds —> New —> R=2,THETA=0, XC=0, YC=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^2\theta>0 \smallsetminus{\{\theta=\frac{\pi}{2}}\}\] \[cos(\frac{\pi}{2})=0, \; 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^2\theta>0 \smallsetminus{\{\theta=\frac{\pi}{2}}\}\] \[cos(\frac{\pi}{2})=0, \; r > \sqrt{\frac{1}{0}}\] \[r_{max} > \infty\]

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

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

\[\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 której krańcach wektory są skierowane do tej podprzestrzeni. W układzie współrzędnych polarnych określa się r_min i r_max, 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},\; \text{jeśli }\theta \in R, \mu \le 1 \text{ to} \; 1 - \mu cos\theta \ge 0, \; \max(cos\theta)=1\]

dla \(\mu < 1\) istnieje minimalny promień początku regionu pułapkowania np. 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},\; \text{jeśli }\theta \in R, \mu \le 1 \text{ to} \; 1 - \mu cos\theta \ge 0, \; \min(cos\theta)=-1\] dla \(\mu = 0.5,\; r_{max}>\frac{\sqrt{6}}{2}\)

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

\[\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\] Jeśli \(a + b^2 < 0\), to dostaniemy siodło, czyli dla takich warunków nie jest możliwa orbita zamknięta. Orbita zamknięta może wystąpić dla \(a + b^2 > 0\). Załóżmy, że a>0.

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

Z analizie wizualnej wykresu: https://www.geogebra.org/3d/xwjtzm3p 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.

\[0=\frac{b^2-a}{a+b^2} -(a + b^2)\] \[0=(b^2-a)-(a + b^2)(a+b^2)\] \[0=b^2-a-(a^2 + ab^2 + ab^2+b^4)\] \[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.

Orbita zamknięta może wystąpić, gdy \(\tau^2-4det(A)<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-4det(A) = -1.06555 <0\)

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

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

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

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óżny 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, dla \(4b^2> a\) \[\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}),m=-1\)

\[y = -(x - b) +\frac{b}{a}\]

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=-(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.

Przykładem jest oscylator van der Pola

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ą.

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

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

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

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

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 x₁ = 2 do x₂ = 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\]

Nasz układ w nowych zmiennych: \[\dot x = w-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\)

Faza szybka będzie zachodziła wzdłuż izokliny \(y=\frac{1}{\mu}F(x)\), więc interesuje nas \(\dot y\) dla tej fazy – aby obliczyć \(\dot y\), skorzystamy z reguły łańcuchowej: \[\frac{dy}{dt} = \frac{1}{\mu}\frac{dF(x)}{dx}\frac{dx}{dt}\] \[\frac{dy}{dt} = (x^2-1)\frac{dx}{dt}\]

Naszym celem jest uzyskać \(T\approx 2\int_{t1}^{t2}dt\), więc musimy znać 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.

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

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

\[\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 przestrzeń fazowa będzie podobna 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.

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

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

Punkt stały dla \(\dot x = \mu - x^2\): \(x_{0} = \pm \sqrt{\mu},\;\mu>0\)

Punkt stały dla \(\dot y = -y\): \(y_{0} = 0\)

Punkt stały i punkt bifurkacji 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}\]

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 44, 45, 46, które jest charakterystyczne dla układów z bifurkacją siodło-węzeł.

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

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\]

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.

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); ab > 0; \tau < 0 \;\text{ściek}; \tau > 0 \;\text{źródło}\] \[X_{*}=(0,0); ab < 0,\;\text{siodło}\]

\[\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, zero-eigen}\] \[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,zero-eigen}\] \[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}\]

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=-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}\]

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 X₀ a z drugiej z punktem X₀^*.
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 X₀ a z drugiej strony z punktem X₀^*, gdyż X₀^* 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)\) .

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

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

\[\dot x = \mu x + x^3\] \[\dot y = - y\]

\[\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}\]

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\] \[W=\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.

\[\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.

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

\[\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.

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

\[\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.