Badanie chaosu w układzie Lorenza z polem magnetycznym

Jerzy Pryga
Astrofizyka i kosmologia
11.05.2022

Plan prezentacji

- Czym jest chaos?
- Metody rozpoznawania zachowań chaotycznych w układach dynamicznych.
- Klasyczny układ Lorenza.
- Bardziej ogólny układ Lorenza-Macka.
- Punkty krytyczne układu i ich stabilność.
- Część konserwatywna i dyssypatywna układu.
- Zastosowania astrofizyczne.
- Podsumowanie.

Chaos  - definicja

Gdy mamy ograniczony (ang. bounded) układ dynamiczny:

((d(Overscript[x, ->]))/dt)=f((Overscript[x, ->]),t)

Mówimy, że jest on chaotyczny (deterministycznie), jeśli jest aperiodyczny i wrażliwy na warunki początkowe.

Definicja Devaney’a:
Niech X będzie przestrzenią metryczną. Ciągła mapa f: X -> X jest chaotyczna jeśli:
1. f jest topologicznie tranzytywna.
2. Punkty okresowe są gęste w X.
3. f jest wrażliwe na warunki początkowe (co tak na prawdę wynika z powyższych).

Jak rozpoznać chaos?

Istnieją różne narzędzia, które potrafią określić czy układ jest chaotyczny czy też nie.

Najpopularniejszym jest chyba tak zwany wykładnik Lyapunova:

(Subscript[λ, L,i])=(Underscript[\[Limit], t->∞])(1/t)ln((((Subscript[δx, i])(t))/((Subscript[(Overscript[δx,   ]), i])(0))))

Gdzie (Subscript[δx, i])(t) jest różnicą pomiędzy dwoma orbitami. Z wykładnikiem Lyapunova łączy się tak zwany czas Lyapunova, który określa najdłuższy czas w którym zachowanie systemu jest przewidywalne. Czas Lyapunova jest odwrotnością największego wykładnika.

Innym sposobem na określenie czy i kiedy układ jest chaotyczny jest analiza stabilności jego punktów krytycznych.

Twierdzenie Sharkovskiego:
Niech ciągłe f będzie miało n-okresową orbitę. Jeśli n ◁ m, wtedy f posiada też m-okresową orbitę, gdzie ◁ to uporządkowanie Sharkovskiego.
1. Jeśli n ≠ (2^k), to f ma nieskończenie wiele orbit.
3. Jeśli f posiada 3-okresową orbitę, to posiada wszystkie orbity.

Układ Lorenza - zarys wyprowadzenia1

Graphics[{{Thick,Line[{{0,0},{1,0}}],Line[{{1,(1/2)},{0,(1/2)}}]},{Arrowheads[{-0,.05}],Arrow[{{(1/2),0},{(1/2),(3/4)}}],Arrow[{{(1/2),0},{(4/3),0}}]},{Blue,Arrowheads[{-.03,.03}],Arrow[{{(3/4),0},{(3/4),(1/2)}}],Text[Style["(Subscript[δT, 0])",Large],{(3/4)+0.05,(1/4)}],Text[Style["(Subscript[T, 0])",Large],{(3/4),(1/2)+0.05}],Text[Style["(Subscript[T, 0])+(Subscript[δT, 0])",Large],{(3/4),0-0.05}]},{Orange,Arrowheads[{-.03,.03}],Arrow[{{(1/4),0},{(1/4),(1/2)}}],Text[Style["h",Large],{(1/4)+0.05,(1/4)}]},Text[Style["z",Large],{(1/2)+0.05,(3/4)-0.05}],Text[Style["x",Large],{(4/3)-0.05,0-0.05}]}]

Równanie dyfuzji termicznej opisujące temperaturę T przyjmuje postać:

((∂T)/(∂t))(Overscript[+v·, ->])∇T=κ(∇^2)T

Równania Navier’a-Stokes’a, opisujące zasadę zachowania pędu w poruszającym się płynie, w rozważanej 2-wymiarowej geometrii przyjmują następującą postać postać (ρ to gęstość płynu, p to ciśnienie, (Overscript[v, ->]) to prędkości, a (Overscript[f, ->]) to przyspieszenia związane z zewnętrznymi siłami, w tym przypadku grawitacyjne)2:

((d(Overscript[v, ->]))/dt)=-(1/ρ)∇p+(μ/ρ)(∇^2)(Overscript[v, ->])+(Overscript[f, ->])

Postać skalarna:

ρ((∂(Subscript[v, z]))/(∂t)) +ρ (Overscript[v·, ->]) ∇(Subscript[v, z])=-ρ g -((∂p)/(∂z))+μ(∇^2)(Subscript[v, z])
ρ((∂(Subscript[v, x]))/(∂t)) +ρ (Overscript[v·, ->]) ∇(Subscript[v, x])= -((∂p)/(∂x))+μ(∇^2)(Subscript[v, x])

Funkcja dewiacji temperatury θ(x,z,t) opisująca jak zmiana temperatury odchyla się od zależności liniowej:

θ(x,z,t) =T(x,z,t)-(Subscript[T, 0])-(Subscript[δT, 0])(1-(z/h))

Funkcja strumienia (Overscript[Ψ, ->]):

(Overscript[v, ->])=∇ (Overscript[× Ψ,        ->])
∇·(Overscript[v, ->])=0

Zakładamy że szukana funkcja będzie postaci iloczynu funkcji zależnych wyłącznie od jednego z trzech parametrów x, z, t. Ponieważ mamy do czynienia z geometrią prostokątną oczekujemy, że powinno istnieć ogólne rozwiązanie w postaci szeregu Fouriera ((Subscript[λ, m,n]) to długości fal modów przestrzennych, a (Subscript[ω, m,n]) to odpowiadające im częstotliwości):

(Underscript[Σ, m,n]) (e^((Subscript[ω, m,n])t))((Subscript[A, m]) cos( (Subscript[λ, m])z) + (Subscript[B, m])sin ((Subscript[λ, m])z))((Subscript[C, n]) cos( (Subscript[λ, n])x) + (Subscript[D, n])sin ((Subscript[λ, n])x))

Dla modelu Lorenza, należy przyjąć pewne warunki brzegowe, które musi spełniać funkcja strumienia i funkcja dewiacji temperatury. Dla funkcji θ, warunki brzegowe są proste, ponieważ temperatury na górnym i dolnym brzegu komórki płynu są ustalone, to:

θ=0 dla z=0,h

Zakładamy, ze dla górnej i dolnej granicy wertykalna składowa prędkości musi być równa 0. Warunek ten przyjmuje zatem postać ∂(Subscript[v, x])/∂z =0 dla z = 0 i z = 1. W modelu Lorenza te warunki są spełnione przez następujące funkcje, będące rozwiązaniem szczególnym i zarazem najprostszym rozwiązania ogólnego (parametr a jest do wyznaczenia):

Ψ(x,z,t)=ψ(t)sin(π z)sin(a x)
θ(x,z,t)=(Subscript[T, 1])(t)sin(π z)cos(a x)-(Subscript[T, 2])(t)sin(2π z)

Przestrzenna część funkcji strumienia odpowiada za konwektywne walce obserwowane kiedy w płynie zaczyna się konwekcja. Funkcja dewiacji temperatury ma dwie składowe, pierwsza, (Subscript[T, 1]), daje różnicę temperatury pomiędzy górną a dolną częścią konwektywnej komórki. Druga, (Subscript[T, 2]), daje odchylenie od liniowej wariacji w centrum konwektywnej komórki jako funkcję z.

Aby otrzymać standardową formę równań Lorenza, dokonujemy pewnych bezpośrednich zamian zmiennych. Dla czasu dokonujemy zamiany t ->(1/( ((π^2) + (a^2))))t, a następnie bierzemy następujące podstawienia:

X(t)=((a π)/(((π^2)+(a^2))(Sqrt[2])))ψ(t)
Y(t)=((r π)/(Sqrt[2]))(Subscript[T, 1])(t)
Z(t)=π r (Subscript[T, 2])(t)

Gdzie r to tak zwana zredukowana liczba Rayleigh’a r = ((a^2)/((((a^2)+(π^2)))^3))R. Wprowadzamy również nowy parametr b = ((4 (π^2))/((a^2)+(π^2))). Z tymi wszystkimi podstawieniami otrzymaliśmy standardową postać równań układu Lorenz’a:

(OverDot[X])=σ(Y-X)
(OverDot[Y])=r X-X Z-Y
(OverDot[Z])=X Z-b Z

Znaczenie fizyczne parametrów

(OverDot[X])=-σ X+σ Y
(OverDot[Y])=r X-X Z-Y
(OverDot[Z])=X Z-b Z

σ = ν/κ - stosunek lepkości kinematycznej do przewodności cieplnej.
r ~ R ~ (Subscript[δT, 0]) - określa gradient temperatury.
b - geometria układu (relacja szerokości komórki konwektywnej i wysokości warstwy cieczy).

X ~ Ψ - określa przepływ cieczy.
Y ~ (Subscript[T, 1]) - proporcjonalne do różnicy temperatur pomiędzy prądem wznoszącym cieczy i opadającym.
Z ~ (Subscript[T, 2]) - proporcjonalne do odchylenia profilu temperatury od liniowości.

Przybliżenia przyjęte w czasie wyprowadzenia

Dla wyrazów bez sił zewnętrznych przyjęto, że płyn jest nieściśliwy:

ρ(T) = (Subscript[ρ, 0])

W pozostałych przypadkach zmiany gęstości w zależności od temperatury ograniczono do wyrazu linowego:

ρ(T) = (Subscript[ρ, 0])[1-β(T-(Subscript[T, 0]))]

Opuszczono wyrazy z (sin^3)((π/h)z).

Układ Lorenza - właściwości

Kładziemy wartości σ = 10 i b = (8/3), czyli takie jak w oryginalnej pracy Lorenza, i patrzymy jak zmieniają się przykładowe wartości X, Y i Z w zależności od wartości r i czasu t:

Lorenz1={x'[t]==σ(y[t]-x[t]),
              y'[t]==-x[t]z[t]+R1 x[t]-y[t],
              z'[t]==x[t]y[t]-b z[t]};
σ=10;
b=(8/3);
Manipulate[
warunkipoczatkowe1={x[0]==x1s,y[0]==x2s,z[0]==x3s};
ZagadnieniePoczatkowe1=Thread[{Lorenz1/.R1-> r1,warunkipoczatkowe1}];
solution1=NDSolve[Flatten[ZagadnieniePoczatkowe1],{x,y,z},{t,0,100},Method->Automatic];

ParametricPlot[{{Evaluate[{t,x[t]}/.solution1],Evaluate[{t,y[t]}/.solution1],Evaluate[{t,z[t]}/.solution1]}},{t,0,tmax},PlotPoints->10000,PlotRange->All,PlotLegends->{"X","Y","Z"},AspectRatio->2/3, AxesLabel->{"t","X(t)|Y(t)|Z(t)"}],{{x1s,3},-10,10},{{x2s,-1},-10,10},{{x3s,5},-10,10},{{tmax,1},1,100},{{r1,0.1},0,30}]

Możemy odtworzyć coś na kształt mapy logistycznej, żeby zobaczyć mechanizm składania funkcji (r = 28).

Weźmy przykładowe współrzędne początkowe ((Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0])) = (3, 5, -1). I sprawdźmy jaka będzie trajektoria układu w przestrzeni fazowej w zależności od wartości r.

LorenzMacek0={x'[t]==σ(y[t]-x[t]),
              y'[t]==-x[t]z[t]+R0 x[t]-y[t],
              z'[t]==x[t]y[t]-b z[t]};
σ=10;
b=(8/3);

Manipulate[
warunkipoczatkowe0={x[0]==3,y[0]==5,z[0]==-1};
ZagadnieniePoczatkowe0=Thread[{LorenzMacek0/.R0-> r,warunkipoczatkowe0}];
solution0=NDSolve[Flatten[ZagadnieniePoczatkowe0],{x,y,z,w},{t,0,100},Method->Automatic];

ParametricPlot3D[{{Evaluate[{x[t],y[t],z[t]}/.solution0]}},{t,0,tmax},PlotPoints->10000,PlotRange->{{-40,40},{-40,40},{0,70}},AspectRatio->2/3, AxesLabel->{"X","Y","Z"}],{{tmax,1},1,100},{{r,0.1},0,30}]

Dla pewnych konfiguracji występują również zamknięte orbity okresowe.

Widać zatem, że w zależności od wartości parametru r można wyróżnić trzy podstawowe stany układu:
- małe r (od 0 do 1) -> brak ruchów cieczy (punkt krytyczny (Subscript[C, 0])),
- większe r  (od 1 do ≈ 24.74) -> stabilna konwekcja (punkty krytyczne (C^±)),

Graphics[{{Thick,Line[{{0,0},{1,0}}],Line[{{1,(1/2)},{0,(1/2)}}]},Table[{Red,Arrow[{{(i/7),(1/15)},{(i/7),(7/15)}}]},{i,1,6}],{Arrowheads[{-0,.05}],Arrow[{{(1/2),0},{(1/2),(3/4)}}],Arrow[{{(1/2),0},{(4/3),0}}]},{Blue,Arrowheads[{-.03,.03}],Arrow[{{(3/4),0},{(3/4),(1/2)}}],Text[Style["(Subscript[δT, 0])",Large],{(3/4)+0.05,(1/4)}],Text[Style["(Subscript[T, 0])",Large],{(3/4),(1/2)+0.05}],Text[Style["(Subscript[T, 0])+(Subscript[δT, 0])",Large],{(3/4),0-0.05}]},{Orange,Arrowheads[{-.03,.03}],Arrow[{{(1/4),0},{(1/4),(1/2)}}],Text[Style["h",Large],{(1/4)+0.05,(1/4)}]},Text[Style["z",Large],{(1/2)+0.05,(3/4)-0.05}],Text[Style["x",Large],{(4/3)-0.05,0-0.05}]}]

Graphics[{{Thick,Line[{{0,0},{1,0}}],Line[{{1,(1/2)},{0,(1/2)}}]},Table[{Red,Dashed,Circle[{(i/7),(1/4)},{(1/15),(1/5)}]},{i,1,6}],{Arrowheads[{-0,.05}],Arrow[{{(1/2),0},{(1/2),(3/4)}}],Arrow[{{(1/2),0},{(4/3),0}}]},{Blue,Arrowheads[{-.03,.03}],Arrow[{{(3/4),0},{(3/4),(1/2)}}],Text[Style["(Subscript[δT, 0])",Large],{(3/4)+0.05,(1/4)}],Text[Style["(Subscript[T, 0])",Large],{(3/4),(1/2)+0.05}],Text[Style["(Subscript[T, 0])+(Subscript[δT, 0])",Large],{(3/4),0-0.05}]},{Orange,Arrowheads[{-.03,.03}],Arrow[{{(1/4),0},{(1/4),(1/2)}}],Text[Style["h",Large],{(1/4)+0.05,(1/4)}]},Text[Style["z",Large],{(1/2)+0.05,(3/4)-0.05}],Text[Style["x",Large],{(4/3)-0.05,0-0.05}]}]

- jeszcze większe r -> chaos.

Graphics[{{Thick,Line[{{0,0},{1,0}}],Line[{{1,(1/2)},{0,(1/2)}}]},Table[{Red,Dashed,Arrow[{{((RandomInteger[{1,15}])/15),((RandomInteger[{1,15}])/32)},{((RandomInteger[{1,15}])/15),((RandomInteger[{1,15}])/32)}}]},{i,1,30}],{Arrowheads[{-0,.05}],Arrow[{{(1/2),0},{(1/2),(3/4)}}],Arrow[{{(1/2),0},{(4/3),0}}]},{Blue,Arrowheads[{-.03,.03}],Arrow[{{(3/4),0},{(3/4),(1/2)}}],Text[Style["(Subscript[δT, 0])",Large],{(3/4)+0.05,(1/4)}],Text[Style["(Subscript[T, 0])",Large],{(3/4),(1/2)+0.05}],Text[Style["(Subscript[T, 0])+(Subscript[δT, 0])",Large],{(3/4),0-0.05}]},{Orange,Arrowheads[{-.03,.03}],Arrow[{{(1/4),0},{(1/4),(1/2)}}],Text[Style["h",Large],{(1/4)+0.05,(1/4)}]},Text[Style["z",Large],{(1/2)+0.05,(3/4)-0.05}],Text[Style["x",Large],{(4/3)-0.05,0-0.05}]}]

Układ Lorenza z polem magnetycznym - zarys wyprowadzenia3

Graphics[{{Thick,Line[{{0,0},{1,0}}],Line[{{1,(1/2)},{0,(1/2)}}]},{Arrowheads[{-0,.05}],Arrow[{{(1/2),0},{(1/2),(3/4)}}],Arrow[{{(1/2),0},{(4/3),0}}]},{Blue,Arrowheads[{-.03,.03}],Arrow[{{(3/4),0},{(3/4),(1/2)}}],Text[Style["(Subscript[δT, 0])",Large],{(3/4)+0.05,(1/4)}],Text[Style["(Subscript[T, 0])",Large],{(3/4),(1/2)+0.05}],Text[Style["(Subscript[T, 0])+(Subscript[δT, 0])",Large],{(3/4),0-0.05}]},{Orange,Arrowheads[{-.03,.03}],Arrow[{{(1/4),0},{(1/4),(1/2)}}],Text[Style["h",Large],{(1/4)+0.05,(1/4)}]},Text[Style["z",Large],{(1/2)+0.05,(3/4)-0.05}],Text[Style["x",Large],{(4/3)-0.05,0-0.05}],{Purple,Dashed,Arrow[{{(1/8),(1/8)},{(7/8),(1/8)}}],Arrow[{{(1/8),(3/8)},{(7/8),(3/8)}}],Arrow[{{(1/8),(5/8)},{(7/8),(5/8)}}],Text[Style["(Overscript[B, ->])",Large],{(5/8),(1/8)+0.05}]}}]

Równania Navier’a-Stokes’a przyjmuje teraz postać:

((d(Overscript[v, ->]))/dt) =-(1/ρ)∇(p+((B^2)/(2(Subscript[μ, 0]))))+((((Overscript[B, ->])·∇)(Overscript[B, ->]))/((Subscript[μ, 0])ρ))+ν(∇^2)(Overscript[v, ->])+(Overscript[f, ->])

Oprócz powyższego, potrzebne jest również równanie opisujące dyfuzję magnetyczną (ang. magnetic advection-diffusion), gdzie η jest rezystancją magnetyczną:

((d(Overscript[B, ->]))/dt) =((Overscript[B, ->])·∇)(Overscript[v, ->])+η(∇^2)(Overscript[B, ->])

Oprócz funkcji strumienia (Overscript[Ψ, ->]), wprowadzamy jeszcze potencjał wektorowy dla pola magnetycznego(Overscript[ A , ->])= (0, α(x,z,t), 0):

(Overscript[v, ->])=∇ (Overscript[× Ψ,          ->])
∇·(Overscript[v, ->])=0

(Overscript[B, ->])=∇ (Overscript[× A,          ->])
∇·(Overscript[B, ->])=0

Następnie wyprowadzenie przebiega analogicznie jak w przypadku klasycznego układu, tylko tutaj wprowadzamy jeszcze jedną zmodyfikowaną zmienną W (gdzie α’(t) jest nową bezwymiarową zmienną):

W(t)=((1/(κ(Sqrt[2]))) (a/(1+(a^2))))α'(t)

W ten sposób otrzymujemy układ 4 równań:

(OverDot[X])=-σ X+σ Y-(Subscript[ω, 0])W(t)
(OverDot[Y])=-X Z+r X-Y
(OverDot[Z])=XY-b Z
(OverDot[W])= (Subscript[ω, 0])X-(Subscript[σ, m])W

Znaczenie fizyczne parametrów

(OverDot[X])=-σ X+σ Y-(Subscript[ω, 0])W(t)
(OverDot[Y])=r X-X Z-Y
(OverDot[Z])=XY-b Z
(OverDot[W])= (Subscript[ω, 0])X-(Subscript[σ, m])W

σ = ν/κ - stosunek lepkości kinematycznej do przewodności cieplnej.
r ~ R ~ (Subscript[δT, 0]) - określa gradient temperatury.
b - geometria układu (relacja szerokości komórki konwektywnej i wysokości warstwy cieczy).
(Subscript[ω, 0]) ~ B - określa natężenie przyłożonego pola magnetycznego.
(Subscript[σ, m])= η/κ - stosunek oporności i przewodności.

X ~ Ψ - określa intensywność przepływu cieczy.
Y ~ (Subscript[T, 1]) - proporcjonalne do różnicy temperatur pomiędzy prądem wznoszącym cieczy i opadającym.
Z ~ (Subscript[T, 2]) - proporcjonalne do odchylenia profilu temperatury od liniowości.
W ~ A - określa natężenie indukowanego pola magnetycznego.

Układ Lorenza Macka - punkty krytyczne

Analizę układu zaczynamy od szukania punktów krytycznych, czyli takich gdzie (d/dt)(Overscript[X, ->]) = 0. Jeden jest trywialny, w punkcie (Subscript[C, 0]) = (0,0,0,0), a pozostałe dwa, nazwijmy je (C^±), zależą od parametrów układu i są symetryczne względem transformacji odbicia:

t->t,X->-X,Y->-Y,Z->Z,W->-W

Chcemy również aby punkty (C^±) były rzeczywiste, co nakłada nam ograniczenie na wartości parametrów:

r≥((σ (Subscript[σ, m])+(Subscript[ω, 0]^2))/(σ (Subscript[σ, m])))

Istnieje wiele sposobów badania stabilności punktów krytycznych np: szukanie funkcji Lapunova, twierdzenie o układzie pierwszego przybliżenia, kryterium Hurtwitza.

Stabilność punktów krytycznych

Zlinearyzowany układ ((Subscript[C, x]),(Subscript[C, y]),(Subscript[C, z]) - poszczególne współrzędne punktów krytycznych):

(OverDot[X])=-σ X+σ Y-(Subscript[ω, 0])W(t)
(OverDot[Y])=(r-(Subscript[C, x])-(Subscript[C, Z])) X-(Subscript[C, x]) Z-Y
(OverDot[Z])=2(Subscript[C, Y])X-b Z
(OverDot[W])= (Subscript[ω, 0])X-(Subscript[σ, m])W

Twierdzenie o układzie pierwszego przybliżenia - jeśli wszystkie współczynniki zlinearyzowanego układu równań są rzeczywiste i wszystkie rozwiązania równania charakterystycznego (Subscript[λ, i]) mają części rzeczywiste ujemne, to punkt krytyczny (Subscript[C, i]) jest stabilny i jest też stabilnym punktem układu ogólnego.

Zobaczmy to dla puntu (Subscript[C, 0]).Kładziemy σ = 1, b = 8/3 i sprawdzamy zależność od pozostałych parametrów.

λ00=(1/6) (-2 (1+c1+c5)-(2 (2^(1/3)) (1+(c1^2)-3 (c4^2)+c1 (-1+3 c2-c5)-c5+(c5^2)))/(((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3))-(2^(2/3)) (((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3)));


λ11=(1/12) (-4 (1+c1+c5)+(2 (2^(1/3)) (1+I (Sqrt[3])) (1+(c1^2)-3 (c4^2)+c1 (-1+3 c2-c5)-c5+(c5^2)))/(((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3))+(2^(2/3)) (1-I (Sqrt[3])) (((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3)));


λ22=(1/12) (-4 (1+c1+c5)+(2 (2^(1/3)) (1-I (Sqrt[3])) (1+(c1^2)-3 (c4^2)+c1 (-1+3 c2-c5)-c5+(c5^2)))/(((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3))+(2^(2/3)) (1+I (Sqrt[3])) (((2-3 c1-3 (c1^2)+2 (c1^3)+9 c1 c2+9 (c1^2) c2+18 (c4^2)-9 c1 (c4^2)-3 c5+12 c1 c5-3 (c1^2) c5-18 c1 c2 c5-9 (c4^2) c5-3 (c5^2)-3 c1 (c5^2)+2 (c5^3)+√(4 (((-1-(c1^2)+3 (c4^2)+c5-(c5^2)+c1 (1-3 c2+c5)))^3)+(((2 (c1^3)+(-2+c5) (-1-9 (c4^2)+c5+2 (c5^2))+(c1^2) (9 c2-3 (1+c5))-3 c1 (1+3 (c4^2)-4 c5+(c5^2)+c2 (-3+6 c5))))^2))))^(1/3)));
Manipulate[Plot[{-c3/.c3->(8/3),Re[λ00]/.c1-> 1/.c3-> (8/3)/.c5->σM/.c2->r,Re[λ11]/.c1-> 10/.c3-> (8/3)/.c5->σM/.c2->r,Re[λ22]/.c1-> 10/.c3-> (8/3)/.c5->σM/.c2->r},{c4,0,30},PlotRange->All,AxesLabel-> {"(Subscript[ω, 0])","(Subscript[λ, i])"},PlotLegends->{"(Subscript[λ, 0])","(Subscript[λ, 1])","(Subscript[λ, 2])","(Subscript[λ, 3])"},PlotStyle->{Blue,Red,Green,Purple}],{{r,1},0,10},{{σM,1},0,10}]

Zobaczmy to teraz dla puntu (SuperPlus[C]).

DetermM=(c4^2) (c3+(x10^2)+λ+c3 λ+(λ^2))+(-c5-λ) (-x10 (c1 (x10-2 x20)+x10 λ)-(c3+λ) (λ+(λ^2)+c1 (1-c2+x10+x30+λ)));

Sol=Solve[DetermM==0,λ];

pX02={x10->(((Sqrt[c3]) (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(Sqrt[(c4^2)+c1 c5])),x20-> (((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]) (c4+((c1 c5)/c4)))/(c1 c5 (Sqrt[(c4^2)+c1 c5]))),x30-> ((c4 (-(c4^2)-c1 c5+c1 c2 c5) (c4+((c1 c5)/c4)))/(c1 c5 ((c4^2)+c1 c5))),x40-> (((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(c5 (Sqrt[(c4^2)+c1 c5])))};
λ0=Sol[[1,1,2]]/.pX02/.c1-> 1/.c3-> (8/3);
λ1=Sol[[2,1,2]]/.pX02/.c1-> 1/.c3-> (8/3);
λ2=Sol[[3,1,2]]/.pX02/.c1-> 1/.c3-> (8/3);
λ3=Sol[[4,1,2]]/.pX02/.c1-> 1/.c3-> (8/3);

Manipulate[Plot[{Re[λ0]/.c2->r/.c5-> σM,Re[λ1]/.c2->r/.c5-> σM,Re[λ2]/.c2->r/.c5-> σM,Re[λ3]/.c2->r/.c5-> σM},{c4,0,5},PlotRange->All,AxesLabel-> {"(Subscript[ω, 0])","(Subscript[λ, i])"},PlotLegends->{"(Subscript[λ, 0])","(Subscript[λ, 1])","(Subscript[λ, 2])","(Subscript[λ, 3])"},PlotStyle->{Blue,Red,Green,Purple}],{{r,1.5},0,30},{{σM,1},0,10}]

Kryterium Hurwitz’a - podaje warunki na ujemne części rzeczywiste pierwiastków równania charakterystycznego. Najpierw tworzy się macierz Hurwitz’a, a następnie liczy się jej wyznaczniki główne. Jeśli wszystkie wyznaczniki są dodatnie, to wszystkie części rzeczywiste (Subscript[λ, i]) są ujemne.

Zobaczmy to dla punktu (Subscript[C, 0]):

D1=(1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30));
D2 = (1+c3) (c4^4)+(((1+c3))^2) (c4^2) c5+c5 (c3+(x10^2)) (c3+c5+c3 c5+(x10^2))+(c1^2) (1-c2+c3+x10+x30) ((c5^2)+x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))+c1 ((c5^2)-c2 (c5^2)+(c5^2) x10+2 c5 (x10^2)-2 c2 c5 (x10^2)+2 c5 (x10^3)+(x10^4)+2 (c5^2) x10 x20-2 (x10^3) x20+(c5^2) x30+2 c5 (x10^2) x30+(c3^2) (1-c2+(c4^2)+2 c5+(c5^2)+x10+x30)+(c4^2) (1-c2 (1+c5)+x10-2 x10 x20+x30+c5 (2+x10+x30))+c3 (2 (c5^2)+2 (c4^2) (1-c2+c5+x10+x30)+2 c5 (1-c2+x10+(x10^2)+x30)+x10 ((x10^2)-2 x20+x10 (2-c2+x30))));
D3=-((((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))))^2)+(1+c1+c3+c5) (-((1+c1+c3+c5) (x10 ((c4^2) x10+c1 c5 (x10-2 x20))+c3 ((c4^2)+c1 c5 (1-c2+x10+x30))))+(c3+(c4^2)+c5+c3 c5+(x10^2)+c1 (1-c2+c3+c5+x10+x30)) ((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))));
D4=-((((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))))^2)+(1+c1+c3+c5) (-((1+c1+c3+c5) (x10 ((c4^2) x10+c1 c5 (x10-2 x20))+c3 ((c4^2)+c1 c5 (1-c2+x10+x30))))+(c3+(c4^2)+c5+c3 c5+(x10^2)+c1 (1-c2+c3+c5+x10+x30)) ((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))));
Manipulate[ListPlot[{{Re[D1],Im[D1]},{Re[D2],Im[D2]},{Re[D4],Im[D4]}}/.x10-> 0/.x20-> 0/.x30-> 0/.x40-> 0/.c1-> σ/.c2-> r/.c3-> b/.c4-> ω0/.c5-> σM,PlotRange->{{-Rmax,Rmax},{-Rmax,Rmax}},AxesLabel->{"Re((Subscript[Det, i]))","Im((Subscript[Det, i]))"},PlotMarkers->{Automatic,PointSize[0.05]}],{{Rmax,40},0,100},{{σ,0.926},0,2},{{r,0.563},0,2},{{b,0.97},0,2},{{ω0,0.46},0,2},{{σM,0.95},0,2}]

Zobaczmy to dla punktu (C^±):

D1=(1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30));
D2 = (1+c3) (c4^4)+(((1+c3))^2) (c4^2) c5+c5 (c3+(x10^2)) (c3+c5+c3 c5+(x10^2))+(c1^2) (1-c2+c3+x10+x30) ((c5^2)+x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))+c1 ((c5^2)-c2 (c5^2)+(c5^2) x10+2 c5 (x10^2)-2 c2 c5 (x10^2)+2 c5 (x10^3)+(x10^4)+2 (c5^2) x10 x20-2 (x10^3) x20+(c5^2) x30+2 c5 (x10^2) x30+(c3^2) (1-c2+(c4^2)+2 c5+(c5^2)+x10+x30)+(c4^2) (1-c2 (1+c5)+x10-2 x10 x20+x30+c5 (2+x10+x30))+c3 (2 (c5^2)+2 (c4^2) (1-c2+c5+x10+x30)+2 c5 (1-c2+x10+(x10^2)+x30)+x10 ((x10^2)-2 x20+x10 (2-c2+x30))));
D3=-((((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))))^2)+(1+c1+c3+c5) (-((1+c1+c3+c5) (x10 ((c4^2) x10+c1 c5 (x10-2 x20))+c3 ((c4^2)+c1 c5 (1-c2+x10+x30))))+(c3+(c4^2)+c5+c3 c5+(x10^2)+c1 (1-c2+c3+c5+x10+x30)) ((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))));
D4=-((((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))))^2)+(1+c1+c3+c5) (-((1+c1+c3+c5) (x10 ((c4^2) x10+c1 c5 (x10-2 x20))+c3 ((c4^2)+c1 c5 (1-c2+x10+x30))))+(c3+(c4^2)+c5+c3 c5+(x10^2)+c1 (1-c2+c3+c5+x10+x30)) ((1+c3) (c4^2)+c5 (c3+(x10^2))+c1 (x10 (x10-2 x20)+c5 (1-c2+x10+x30)+c3 (1-c2+c5+x10+x30))));
Manipulate[ListPlot[{{{Re[D1],Im[D1]},{Re[D2],Im[D2]},{Re[D4],Im[D4]}}/.x10-> -(((Sqrt[c3]) (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(Sqrt[(c4^2)+c1 c5]))/.x20-> -(((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]) (c4+((c1 c5)/c4)))/(c1 c5 (Sqrt[(c4^2)+c1 c5])))/.x30-> ((c4 (-(c4^2)-c1 c5+c1 c2 c5) (c4+((c1 c5)/c4)))/(c1 c5 ((c4^2)+c1 c5)))/.x40-> -(((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(c5 (Sqrt[(c4^2)+c1 c5])))/.c1->σ/.c2-> r/.c3-> b/.c4-> ω0/.c5-> σM,{{Re[D1],Im[D1]},{Re[D2],Im[D2]},{Re[D4],Im[D4]}}/.x10-> (((Sqrt[c3]) (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(Sqrt[(c4^2)+c1 c5]))/.x20-> (((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]) (c4+((c1 c5)/c4)))/(c1 c5 (Sqrt[(c4^2)+c1 c5])))/.x30-> ((c4 (-(c4^2)-c1 c5+c1 c2 c5) (c4+((c1 c5)/c4)))/(c1 c5 ((c4^2)+c1 c5)))/.x40-> (((Sqrt[c3]) c4 (Sqrt[-(c4^2)-c1 c5+c1 c2 c5]))/(c5 (Sqrt[(c4^2)+c1 c5])))/.c1->σ/.c2-> r/.c3-> b/.c4-> ω0/.c5-> σM},PlotRange->{{-Rmax,Rmax},{-Rmax,Rmax}},AxesLabel->{"Re((Subscript[Det, i]))","Im((Subscript[Det, i]))"},PlotMarkers->{Automatic,PointSize[0.05]}],{{Rmax,40},0,100},{{σ,0.926},0,2},{{r,0.563},0,2},{{b,0.97},0,2},{{ω0,0.46},0,2},{{σM,0.95},0,2}]

Ewolucja układu w czasie

Kładziemy wartości σ = 1, b = (8/3) i (Subscript[σ, M]) = 1, czyli takie jak w pracy Macka, i patrzymy jak zmieniają się przykładowe wartości X, Y i Z w zależności od wartości r i czasu t:

LorenzMacek1={x'[t]==σ(y[t]-x[t])-ω0 w[t],
              y'[t]==-x[t]z[t]+RM1 x[t]-y[t],
              z'[t]==x[t]y[t]-b z[t],
w'[t]==ω0 x[t]-σm w[t]};
σ=10;
b=(8/3);
σm = 1;
Manipulate[
warunkipoczatkoweM1={x[0]==x1s,y[0]==x2s,z[0]==x3s,w[0]==x4s };
ZagadnieniePoczatkoweM1=Thread[{LorenzMacek1/.RM1-> r1/.ω0-> ω,warunkipoczatkoweM1}];
solutionM1=NDSolve[Flatten[ZagadnieniePoczatkoweM1],{x,y,z,w},{t,0,100},Method->Automatic];

ParametricPlot[{{Evaluate[{t,x[t]}/.solutionM1],Evaluate[{t,y[t]}/.solutionM1],Evaluate[{t,z[t]}/.solutionM1],Evaluate[{t,w[t]}/.solutionM1]}},{t,0,tmax},PlotPoints->10000,PlotRange->All,PlotLegends->{"X","Y","Z","W"},AspectRatio->2/3, AxesLabel->{"t","X(t)|Y(t)|Z(t)|W(t)"}],{{x1s,3},-10,10},{{x2s,-1},-10,10},{{x3s,5},-10,10},{{x4s,5},-10,10},{{tmax,1},1,100},{{r1,0.1},0,30},{{ω,0.1},0,30}]

Ponieważ układ jest 4-wymiarowy (więc wykres przestrzeni fazowej też musiałby taki być) zobaczmy jego rzut na 3-wymiarową przestrzeń. Weźmy przykładowe współrzędne początkowe ((Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0])) = (3, 5, -1,-2) i sprawdźmy jak to wygląda w zależności od wartości r i (Subscript[ω, 0]).

LorenzMacek2={x'[t]==σ(y[t]-x[t])-ω0 w[t],
              y'[t]==-x[t]z[t]+RM2 x[t]-y[t],
              z'[t]==x[t]y[t]-b z[t],
w'[t]==ω0 x[t]-σm w[t]};
σ=10;
b=(8/3);
σm = 1;
Manipulate[
warunkipoczatkoweM2={x[0]==3,y[0]==5,z[0]==-1,w[0]==-2 };
ZagadnieniePoczatkoweM2=Thread[{LorenzMacek2/.RM2-> r/.ω0-> ω,warunkipoczatkoweM2}];
solutionM2=NDSolve[Flatten[ZagadnieniePoczatkoweM2],{x,y,z,w},{t,0,100},Method->Automatic];

ParametricPlot3D[{{Evaluate[{x[t],y[t],z[t]}/.solutionM2]}},{t,0,tmax},PlotPoints->10000,PlotRange->{{-30,30},{-30,30},{0,60}},AspectRatio->2/3, AxesLabel->{"X","Y","Z"}],{{tmax,1},1,100},{{r,0.1},0,30},{{ω,0.1},0,30}]

Dziwny atraktor

Można pokazać, że objętość przestrzeni fazowej dąży do 0, co jest cechą charakterystyczną dziwnego atraktora.

1. Zakładamy, że (Subscript[f, t])((Subscript[x, 0])) jest rozwiązaniem układu dynamicznego z warunkiem początkowym (Overscript[X, ->])(0)=(Subscript[(Overscript[X, ->]), 0]).
2. Dodajemy stronami układ dynamiczny i porządkujemy wyrazy.

(1/2)(d/dt)((X^2)+(Y^2)+(((Z-σ-r))^2)+(W^2))=-(σ (X^2)+(Y^2)+((b(Z-((σ+r)/2)))^2)(Subscript[σ, M]) (W^2))+((b(((σ+r)/2)))^2)

3. W przestrzeni (R^4) rozważamy teraz kulę
	(Subscript[K, 0])={((x(t))^2)+((y(t))^2)+(((z(t)-σ-r))^2)+((w(t))^2)≤(c^2)}.
4. Jeżeli kula (Subscript[K, 0]) zawiera w swoim wnętrzu elipsoidę:
	σ ((x(t))^2)+((y(t))^2)+((b(z(t)-((σ+r)/2)))^2)+σm ((w(t))^2)=b (((((σ+r)/2)))^2)
to na brzegu kuli (Subscript[K, 0]) (Taka sytuacja jest możliwa, gdyż możemy zawsze wybrać dostatecznie duże c.) lewa strona równania z fioletowym tłem jest ujemna! Oznacza to, że odwzorowanie (Subscript[f, 1]) przekształca kulę (Subscript[K, 0]) w kulę (Subscript[K, 1]). Stąd mamy zstępujący ciąg zbiorów:
	(Subscript[K, 1])=(Subscript[f, 1])((Subscript[K, 0])), (Subscript[K, 2])=(Subscript[f, 1])((Subscript[K, 1]))=(Subscript[f, 2])((Subscript[K, 0])) ...
5. Ostatecznie otrzymamy, że objętość (Subscript[V, k]) zbioru (Subscript[R, k]) jest równa
	(Subscript[V, k])=(Subscript[V, 0])(e^(-k(1+σ+b+(Subscript[σ, m])))).

Zobaczmy również, jak objawia się wrażliwość na warunki początkowe (dla wartości parametrów, przy których występuje zachowanie chaotyczne). Weźmy nasz układ (σ = 1, b = (8/3), r = 28, (Subscript[ω, 0])=3 i (Subscript[σ, M]) = 1) z warunkami początkowymi (Subscript[X, 0]) = 19,  (Subscript[Y, 0]) = 20,  (Subscript[Z, 0]) = 50 i  (Subscript[W, 0]) = 10.

parametry:={σ->10,r->28,b->8/3,ω0->3,σm->1}
roz=NDSolve[{x'[t]==σ (y[t]-x[t])-ω0 w[t],y'[t]==x[t] (r-z[t])-y[t],z'[t]== x[t] y[t] -b z[t],w'[t]==ω0 x[t]-σm w[t],x[0]== 19,y[0]== 20,z[0]== 50,w[0]==10}/.parametry,{x,y,z,w},{t,0,100},MaxSteps-> Infinity,Method->{"ExplicitRungeKutta","DifferenceOrder"->4}];
rozz[t_]:=z[t]/.roz
plotz=Plot[rozz[t],{t,0,30},PlotStyle->{Black},AxesLabel->{"t","Z"}]

Zaburzmy teraz (Subscript[Z, 0]) o ΔZ = (10^(-6)).

δ:=(10^(-3))
rozδ3=NDSolve[{x'[t]==σ (y[t]-x[t])-ω0 w[t],y'[t]==x[t] (r-z[t])-y[t],z'[t]== x[t] y[t] -b z[t],w'[t]==ω0 x[t]-σm w[t],x[0]== 19,y[0]== 20,z[0]== 50+δ,w[0]==10}/.parametry,{x,y,z,w},{t,0,100},MaxSteps-> Infinity,Method->{"ExplicitRungeKutta","DifferenceOrder"->4}];
rozzδ3[t_]:=z[t]/.rozδ3
plotzδ3=Plot[rozz[t],{t,0,30},PlotStyle->{Red},AxesLabel->{"t","Z"}]

Plot[rozz[t]-rozzδ3[t],{t,0,30},PlotStyle->{Blue},AxesLabel->{"t","ΔZ"}]

A teraz zaburzmy (Subscript[Z, 0]) o ΔZ = (10^(-3)).

δ:=(10^(-6))
rozδ6=NDSolve[{x'[t]==σ (y[t]-x[t])-ω0 w[t],y'[t]==x[t] (r-z[t])-y[t],z'[t]== x[t] y[t] -b z[t],w'[t]==ω0 x[t]-σm w[t],x[0]== 19,y[0]== 20,z[0]== 50+δ,w[0]==10}/.parametry,{x,y,z,w},{t,0,100},MaxSteps-> Infinity,Method->{"ExplicitRungeKutta","DifferenceOrder"->4}];
rozzδ6[t_]:=z[t]/.rozδ6
plotzδ6=Plot[rozz[t],{t,0,30},PlotStyle->{Orange},AxesLabel->{"t","Z"}]

Plot[rozz[t]-rozzδ6[t],{t,0,30},PlotStyle->{Magenta},AxesLabel->{"t","ΔZ"}]

Część konserwatywna i dyssypatywna

Zapiszmy sobie nasz układ równań w innej formie:
- (Overscript[X, ->])(t) - wektor (X, Y, Z, W).
- M - macierz współczynników.
- F - prawa strona układu, pole wektorowe.

Wtedy nasz układ przyjmuje postać następującą:

(Overscript[X, ·])=F=M×(Overscript[X, ->])

Można policzyć dywergencję pola F:

∇·F= -1-σ-b-(Subscript[σ, M])

Rozbijamy M na część symetryczną R i antysymetryczną J:

R=(1/2)(M+(M^T))
J=(1/2)(M-(M^T))

Zapisując to jawnie:

R=(({{-σ, ((σ+r)/2), 0, 0}, {((σ+r)/2), -1, 0, 0}, {0, 0, -b, 0}, {0, 0, 0, -(Subscript[σ, M])}}))
J=(({{0, ((σ-r)/2), 0, -(Subscript[ω, 0])}, {(1/2)(r-σ), 0, -X, 0}, {0, X, 0, 0}, {(Subscript[ω, 0]), 0, 0, 0}}))

Następnie przechodzimy do rozkładu pola wektorowego F na część dyssypatywną (Subscript[F, d]) = R ×X i część konserwatywną (Subscript[F, c])= J × X.

(Overscript[X, ·])=(Subscript[F, c])+(Subscript[F, d])

Dywergencje poszczególnych części:

∇·(Subscript[F, c])=0
∇·(Subscript[F, d])= -1-σ-b-(Subscript[σ, M])

Można zauważyć, że układ dynamiczny dyssypatywny (Subscript[F, d]) jest liniowy, zatem ścisłe rozwiązanie jest znane.

Wprowadzamy potencjał pola wektorowego ϕ taki że:

(Subscript[F, d])=-∇ϕ
∇×(Subscript[F, d])=0

ϕ=(1/2)σ (X^2)-(1/2)(σ+r) X Y -(1/2)(Y^2)-(1/2)b (Z^2)+(1/2)(Subscript[σ, M])W

Warto zauważyć, że:
1. Znika rotacja pola wektorowego (Subscript[F, d]).
2. Pole (Subscript[F, d]) niesie pełną dywergencję pola F.
3. Kształt potencjału w zależności od parametrów jest zbadany i są to hiperkwadryki.

Zajmijmy się teraz konserwatywną częścią pola F, czyli (Subscript[F, c]), dla której dywergencja równa jest zero. Oznacza to że można zdefiniować potencjał wektorowy A, taki że (Subscript[F, c]) = ∇×A.

Da się znaleźć wszystkie niezależne całki pierwsze układu dynamicznego (Subscript[F, c]):

(Subscript[I, 1])=(Y^2)+σ Z-r Z+(Z^2)
(Subscript[I, 2])=(X^2)-σ Z+r Z+(W^2)

Warunek na otrzymanie trzeciej całki jest następujący:

(Subscript[I, 1])+(1/4)((((Subscript[c, 2])-(Subscript[c, 1])))^2)>0

(Subscript[I, 3])=-(Subscript[ω, 0])arctan(((σ-r+2 Z)/(2 Y)))+W

Wnioski z tego są następujące:
1. Układ dynamiczny (Subscript[F, c]) jest kompletnie całkowalny.
2. Można jeszcze wyznaczyć czwartą całkę pierwszą w przestrzeni {X,Y,Z,W,t}.

Warto zauważyć, że przez analogię do układów trzeciego rzędu można wprowadzić energię Helmholtza w następujący sposób:

H=(1/2)((Subscript[I, 1])+(Subscript[I, 2]))

Jej zmianę w czasie definiujemy następująco:

(dH/dt)=∇H·(({{((Subscript[dx, 1])/dt)}, {((Subscript[dx, 2])/dt)}, {((Subscript[dx, 3])/dt)}, {((Subscript[dx, 4])/dt)}}))=∇H·F=∇H·((Subscript[F, c])+(Subscript[F, d]))=∇H·(Subscript[F, c])+∇H·(Subscript[F, d])=∇H·(Subscript[F, d])

Zatem dla naszego układu można policzyć, że:

(Overscript[H, ·])=-(Subscript[c, 1])((Subscript[x, 1])^2)+((Subscript[c, 1])+(Subscript[c, 2]))(Subscript[x, 1])(Subscript[x, 2])-(Subscript[x, 2])-(Subscript[c, 3])(Subscript[x, 3])-(Subscript[c, 5])((Subscript[x, 4])^2)

Z powyższych rachunków wynika, że energia Helmholtza ma ciekawe właściwości:
1. Jest oczywiście całką układu dynamicznego.
2. Gradient H jest ortogonalny do pola (Subscript[F, c]).
3. (OverDot[H])=(∇H)·(Subscript[F, d])=-(Subscript[c, 1]) ((Subscript[x, 1])^2)+((Subscript[c, 1])+(Subscript[c, 2])) (Subscript[x, 1]) (Subscript[x, 2])-(Subscript[x, 2])-(Subscript[c, 3]) (Subscript[x, 3])-(Subscript[c, 5]) ((Subscript[x, 4])^2).
4. Na podstawie punktu 3 widać, że energia Helmholtza sprzęga pole wektorowe (Subscript[F, c]) z polem wektorowym (Subscript[F, d]).

Co ciekawego można z tym zrobić?

Można sprawdzić jaki wpływ ma część dyssypatywna na zachowanie układu. Mnożymy (Subscript[F, d]) przez parametr ϵ z przedziału od 0 do 1.

(Overscript[X, ·])=(Subscript[F, c])+(Subscript[ϵF, d])

Będziemy sterować wartością tego parametru “zaburzając” część konserwatywną i sprawdzać jak to wpłynie na właściwości układu, na przykład badając jak będzie zmieniać się (OverDot[H]).

Przykłady możliwych zastosowań astrofizycznych

- Plazma wewnątrz Słońca i w koronie słonecznej4.

- Płynne jądro Ziemi i innych planet5.

- Ośrodek międzygwiazdowy i międzygalaktyczny.
- Dyski akrecyjne.
I wiele wiele innych nieastrofizycznych przykładów...

Wnioski

- Tego typu układ ma 3 podstawowe sposoby zachowania: stacjonarny płyn, konwekcja, chaos.
- Pole magnetyczne może stabilizować układ lub wręcz przeciwnie.
- Badanie układu LM może pomóc zrozumieć i modelować niektóre układy astrofizyczne.

Bibliografia

1	Macek, Wiesław M., and Marek Strumik. "Model for hydromagnetic convection in a magnetized fluid." Physical Review E 82.2 (2010): 027301.

2	Saltzman, Barry. "Finite amplitude free convection as an initial value problem—I." Journal of atmospheric sciences 19.4 (1962): 329-341.

3	Macek, Wiesław M. "Nonlinear dynamics and complexity in the generalized Lorenz system." Nonlinear Dynamics 94.4 (2018): 2957-2968.

4	Boffetta, Guido, et al. "Power laws in solar flares: self-organized criticality or turbulence?." Physical review letters 83.22 (1999): 4662.

5	Finlay, Christopher C., and Andrew Jackson. "Equatorially dominated magnetic field change at the surface of Earth's core." Science 300.5628 (2003): 2084-2086.