Das Lorenz Modell

H.J. Lüdde

Das Lorenz Modell basiert auf den Navier-Stokes Gleichungen und stellt ein stark vereinfachtes Modell der Dynamik in einer Flüssigkeitsschicht dar (Rayleigh-Benard Konvektion). Ein einfaches Beispiel ist das langsame Erwärmen einer ca 2-3cm hohen Wasserschicht in einem Gefäß mit großer Grundfläche. Zu Beginn ist die Temperaturdifferenz zwischen unterer und oberer Wasserschicht Null. Mit der Zeit erhöht sich diese Temperaturdifferenz und in der ersten Phase findet ein Temperaturausgleich durch Wärmeleitung innerhalb der Wasserschichten statt: die Flüssigkeit ist in Ruhe. Mit Überschreiten einer ersten kritischen Temperaturdifferenz ändert sich das Verhalten der Flüssigkeit spontan: der Wärmetransport erfolgt in dieser Phase durch Konvektion (warme Wasserschichten steigen auf, kalte Wasserschichten sinken ab). D.h. die Temperatur- und Geschwindigkeitsverteilungen werden zeitabhängig und es stellt sich mit der Zeit ein dynamisches Gleichgewicht ein. In dieser Phase bilden sich sichtbare Strukturen (kreisförmige Rollen, Hexagonalstrukturen) aus (Selbstorganisation). Bei weiterer Erwärmung erreicht die Temperaturdifferenz zwischen den Kontaktschichten einen weiteren kritischen Wert: die Flüssigkeit beginnt zu sieden und kurz darauf zu kochen. Der Temperaturausgleich innerhalb der Flüssigkeit ist jetzt charakterisiert durch Turbulenzen, also völlig unregelmäßige Strukturen: mathematisch markiert der zweite kritische Wert den Beginn von Chaos.

Stark vereinfacht lässt sich diese Beobachtung mit Hilfe der Lorenzgleichung studieren:

Dabei bedeuten: x die Geschwindigkeitsverteilung innerhalb der Flüssigkeit senkrecht zu den Grenzflächen, y die Temperaturdifferenz zwischen auf- und absteigender Flüssigkeit und z der mittlere konvektive Wärmestrom. s,b,r sind Kontrollparameter, welche die Flüssigkeit charakterisieren (s,b) bzw die Randbedingungen, gegeben durch die Temperaturdifferenz r zwischen den Kontaktschichten (die Temperatur der oberen Schicht ist gegeben durch die Umgebungstemperatur (Dampfphase oberhalb der Flüssigkeit), die Temperatur der unteren Kontaktschicht durch die Temperatur des Topfbodens).

Wir wollen die Lorenzgleichungen diskutieren für s=10, b=8/3 und in Abhängigkeit von r [1].

Referenz

[1] E.N. Lorenz: Deterministic Nonperiodic Flow, J. Atmos. Sc. 20 , 130 (1963)

1) Topologische Analyse

In diesem Abschnitt untersuchen wir die Fixpunktlösungen der Lorenzgleichungen, um einen Eindruck von den Stabilitätseigenschaften der Lösungen für verschiedene Kontrollparameter zu erhalten. In der klassischen Arbeit von Lorenz werden s und b konstant gewählt und r variiert [2].

Referenz

[2] E-Skript zu Differentialgleichungen: http://www.th.physik.uni-frankfurt.de/~luedde/E-Skript/

> restart:

> with(plots):

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

Definition der Gleichungen für die kritischen Punkt Lösungen

> eqncp:=s*(y-x)=0,r*x-y-x*z=0,x*y-b*z=0;

Lösung des Gleichungssystems

> solcp:=solve({eqncp},{x,y,z});

> solcp2:=allvalues(solcp[2]); # Auswertung der Wurzeln

Charakteristische Gleichung des linearisierten Lorenzsystems

> A:=matrix(3,3,[-s-p,s,0,r-z,-1-p,-x,y,x,-b-p]);

> DA:=unapply(det(A),x,y,z);

Setze ersten kritischen Punkt ein

> assign(solcp[1]):

> x,y,z;

> DA1:=DA(x,y,z);

...und bestimme die drei Eigenwerte p

> solDA1:=[solve(DA1,p)];

Alle drei Eigenwerte sind reell. Sind sie darüberhinaus alle negativ, so handelt es sich bei dem kritischen Punkt solcp[1] um einen stabilen Knoten. Offenbar hängt diese Stabilitätseigenschaft von r ab:

> r:=1:

> solDA1; # kritische Lösung für r=1

> for i to nops(solDA1) do

> simplify(solDA1[i],assume=positive):

> od;

Somit haben wir folgende Situation: für r<1 ist solcp[1] ein attraktiver Knoten, bei r>1 ein Sattel. D.h. bei r=1 ändern sich die Stabilitätseigenschaften der Lorenzgleichung.

Im nächsten Schritt untersuchen wir die Eigenschaften des zweiten kritischen Punktes solcp[2]:

> x:='x':y:='y':z:='z':

> s:='s':b:='b':r:='r':

> assign(solcp2[1]):

> x,y,z;

> DA2:=DA(x,y,z);

> simplify(%); # fasse Terme zusammen

> DA2:=collect(%,p); # sortiere nach Potenzen von p (Eigenwert)

Für r=1 erhält man

> r:=1:

> DA2;

> solDA2:=[solve(DA2,p)];

Für r=1 fallen die drei Punktlösungen solcp2[1,2] und solcp[1] zusammen und bilden einen Sattel. Für r>1 bleibt solcp[1] ein Sattel, während solcp[1,2] ihre Stabilitätseigenschaften ändern.

Wir versuchen die Eigenwerte p für r>1 zu berechnen:

> r:='r':

> solDA2:=[solve(DA2,p)]:

> for i to nops(solDA2) do

> simplify(solDA2[i]):

> od;

Ooobs...... Jetzt haben wir zwar eine Lösung für die Eigenwerte für r>1, die Ausdrücke sind aber unlesbar!!

Um eine Aussage über die Stabilität der Punktlösung zu machen, benötigen wir zum Glück nicht die explizite Form der Eigenwerte: es ist ausreichend zu wissen, dass alle Realteile negativ sind, damit der kritische Punkt ein stabiler Attraktor ist. Umgekehrt reicht es aus zu zeigen, dass nur einer der Eigenwerte einen positiven Realteil besitzt, mit der Folge, dass der kritische Punkt instabil ist. Wenn wir wie bisher s und b konstant halten, stellt sich die Frage, ob ein kritischer Wert für r existiert, an dem der Realteil mindestens eines Eigenwertes sein Vorzeichen ändert. Man nennt diesen Fall eine Hopfbifurkation.

Das Hurwitz Kriterium:

Alle Wurzeln der Gleichung

...

haben dann und nur dann negative Realteile, wenn alle Determinanten

, , ...

mit für positiv sind.

In unserem Fall gilt , d.h. D1>0. Wir suchen eine Bedingung, für die D2=0:

> D2:=det(matrix(2,2,[b*(s+r),2*b*s*(r-1),1,s+1+b]))=0;

> rc:=unapply(solve(%,r),s,b);

Man kann überprüfen, dass diese Bedingung auch D3 verschwinden lässt. D.h. - für ist rc im erlaubten Parameterbereich - für ist der kritische Punkt ein stabiler Fokus, für wird er instabil.

Für die Standardwerte s=10 und b=8/3 erhält man somit

> rc(10,8/3);

> evalf(%);

>

2) Numerische Lösungen

Mit Kenntnis der Stabilität der Punktlösungen für verschiedene Werte von r werden die Lorenzgleichungen für die klassische Parameterkombination s=10, b=8/3 gelöst.

> restart;

> with(plots): with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> lorenzeqs:=diff(x(t),t)=s*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t);

> vars:={x(t),y(t),z(t)}:

> s:=10:b:=8/3:

> rho[1]:=0.8:rho[2]:=1.2:rho[3]:=3:rho[4]:=10:rho[5]:=22:rho[6]:=28:

> ics:=x(0)=1,y(0)=1,z(0)=1:

> for i from 1 to 6 do

> rf:=rand()/10^12:gf:=rand()/10^12:bf:=rand()/10^12:

> r:=rho[i]:

> lorenzsol:= dsolve({lorenzeqs, ics}, vars, type=numeric, method=classical[rk4],output=listprocedure):

> lorenzx[i]:=odeplot(lorenzsol,[t,x(t)],0..100,numpoints=500, color=COLOR(RGB,rf,gf,bf)):

> if i=6 then

> lorenzxyz:=odeplot(lorenzsol,[x(t),y(t),z(t)],0..100, numpoints=5000):

> fi:

> od:

> display(seq(lorenzx[i],i=1..3),view=[0..6,0..3],axes=boxed,labels=[t,x], title=` Verschiedene Lösungen der Lorenzgleichungen für r=0.8,1.2,3 `);

Für r<1 befindet sich die Flüssigkeit in ruhendem Zustand: alle Zustandsvariablen konvergieren gegen die Punktlösung (0,0,0) (attraktiver Knoten). Für r>1 ist (0,0,0) keine stabile Punktlösung mehr (Sattel), die Flüssigkeit tritt in die konvektive Phase. Lösungen konvergieren gegen den zweiten Attraktor ( ) einen stabilen Fokus.

> display(seq(lorenzx[i],i=4..6),view=[0..100,-20..20],axes=boxed,labels=[t,x], title=` Verschiedene Lösungen der Lorenzgleichungen für r=10,20,28 `);

Mit wachsenden Werten für r wird die Einschwingphase (transiente Phase) immer ausgeprägter, bis der attraktive Fokus für r>rc über eine Hopfbifurkation in eine instabile Punktlösung übergeht: der Bereich chaotischer Lösungen ist erreicht (anschaulich vergleichbar mit den Turbulenzen, die in einer kochenden Flüssigkeiten auftreten).

> display(lorenzxyz,axes=boxed,labels=[x,y,z],orientation=[100,82], title=` Chaotische Lösungen der Lorenzgleichungen für r=28 `,shading=xyz,thickness=0,style=line);

Die Phasenraumdarstellung einer chaotischen Lösung der Lorenzgleichungen: die Zustandsvariablen springen willkürlich zwischen zwei Attraktorbecken.

3) Poincare Abbildung (Lorenz Plot)

Analog zur Stroboskop Abbildung des getriebenen Pendels, versucht man für ein autonomes System eine 'Wiederkehrabbildung' zu definieren, indem man aufeinanderfolgende Extrema einer Lösungskomponente (hier z(n+1) vs z(n)) gegeneinander aufträgt. Diese spezielle Form der Poincare Abbildung nennt man Lorenz Plot (siehe Kap 2.1 in Abschnitt II).

> restart;

> with(plots): with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> lorenzeqs:=diff(x(t),t)=s*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t):

> vars:={x(t),y(t),z(t)}:

> ics:=x(0)=1,y(0)=1,z(0)=1:

> s:=10:b:=8/3:tmax:=200:dt:=0.02:zahl:=0:r:=28:

> for k from 1 to tmax/4 do

> lorenzsol:=dsolve({lorenzeqs,ics},vars,numeric,output=listprocedure);

> assign(lorenzsol);

> X:=x(t):

> Y:=y(t):

> Z:=z(t):

> for i from 0 to 4 by dt do

> dz1:=(Z(i)-Z(i+dt));dz2:=(Z(i+dt)-Z(i+2*dt));zi:=Z(i):xi:=X(i):yi:=Y(i):

> if dz1*dz2< 0 then zahl:=zahl+1:pts[zahl]:=Z(i+dt): fi; # Bedingung für Maxima bzw Minima

> od:

> unassign(x,y,z);

> ics:=x(0)=xi,y(0)=yi,z(0)=zi;

> od:

> plot([seq([pts[j],pts[j+1]],j=0..(zahl-1))],style=point,symbol=POINT,color=blue,view=[0..50,0..50],font=[TIMES,ROMAN,12],title=` Lorenzplot für r=28`,axes=boxed,labels=[`Z(n)`,`Z(n+1)`]);

Lorenz Plot einer chaotischen Lösung des Lorenz Modells. Offenbar existiert eine eindimensionale diskrete Abbildung (Differenzengleichung) z(n+1)=f(z(n)), welche die Extremwerte der Koordinate z(t) miteinander verknüpft. Eine entsprechende Abbildung erhält man für die anderen Koordinaten: wie sehen sie aus?

4) Intermittenz

Die folgenden 3 Abbildungen zeigen einen möglichen Weg von einer stabil oszillierenden zu einer chaotischen Lösung: die Intermittenz. Eine zunächst stabile Lösung wird durch immer häufiger auftretende 'Störsignale' unterbrochen, bis der Abstand zwischen zwei unregelmäßigen Ausschlägen so klein ist, dass die Lösung irregulär wird: der Intermittenzpfad zum Chaos.

> restart;

> with(plots): with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> lorenzeqs:=diff(x(t),t)=s*(y(t)-x(t)),diff(y(t),t)=r*x(t)-y(t)-x(t)*z(t),diff(z(t),t)=-b*z(t)+x(t)*y(t):

> s:=10:b:=8/3:

> vars:={x(t),y(t),z(t)}:

> rho[1]:=166:rho[2]:=166.2:rho[3]:=166.4:

> ics:=x(0)=1,y(0)=0,z(0)=0:

> for i from 1 to 3 do

> rf:=rand()/10^12:gf:=rand()/10^12:bf:=rand()/10^12:

> r:=rho[i]:

> lorenzsol:= dsolve({lorenzeqs, ics}, vars, type=numeric, method=classical[rk4],output=listprocedure):

> lorenzf[i]:=odeplot(lorenzsol,[t,z(t)],0..60, numpoints=20000, color=COLOR(RGB,rf,gf,bf)):

> od:

> for i from 1 to 3 do

> display(lorenzf[i],view=[30..60,60..250], axes=boxed,labels=[t,z], title=` Lösungen der Lorenzgleichungen im Intermittenzbereich `);

> od;

>

5) Übung

Diskutieren Sie analog zu den Lorenzgleichungen das Rössler-Modell:

für a=b=0.2 in Abhängigkeit von c. Das System besitzt einen stabilen Grenzzyklus und durchläuft mit wachsenden Werten für c eine Folge von Periodenverdopplungen. Für c>4.2 ist die Dynamik chaotisch. Beginnen Sie mit der Anfangsbedingung x(0)=y(0)=z(0)=1.