Nichtlineare diskrete Abbildungen

H.J. Lüdde

nach

R. Mahnke Kapitel 8

Diskrete Abbildungen nennt man Iterationen der Form

Motivation:

(i) Wir haben bei Poincare Schnitten gesehen, dass die 'Stroboskopaufnahmen' einer chaotischen Phasenraumtrajektorie regelmäßige Strukturen aufweist (seltsame Attraktoren): Ein Poincare Schnitt diskretisiert die Trajektorie. (ii) Geht man davon aus, dass eine Zeitreihe x(t) deterministischer Natur ist, so sollte bezogen auf ein festes Intervall dt ein diskreter Zusammenhang bestehen zwischen x(n+1)=x(t+(n+1)dt) und x(n). (iii) Für komplexe Systeme sind Zeitreihen oft nur diskret verfügbar: z.B. die Populationen zu aufeinanderfolgenden Generationen (logistische Abbildung).

Es gibt verschiedene in der Literatur diskutierte diskrete Abbildungen: (i) Logistische Abbildung, (ii) Standardabbildung, (iii) Henon Abbildung.

Die logistische Abbildung

wurde erstmals von Verhulst zur Beschreibung einer einfachen Populationskinetik mit nichtlinearer Dämpfung vorgeschlagen.

Für s=0 entnimmt man der Iteration folgende Eigenschaften:

r<1: die Population stirbt,

r=1: die Population bleibt stabil,

r>1: die Population wächst exponentiell.

Die logistische Abbildung kann mit Hilfe der Skalentransformation x(n)=s/r y(n) auf Normalform gebracht werden

Wir wollen diese Abbildung im Einzelnen diskutieren:

> restart;

> with(plots):

1) Die logistische Abbildung für verschiedene Kontrollparameter r

In diesem Abschnitt betrachten wir numerische Lösungen der logistischen Abbildung für verschiedene Werte des Kontrollparameters r.

> r:=array(1..8):x0:=array(1..8):x:=array(1..401):

> r[1]:=0.5:r[2]:=0.95:r[3]:=1:r[4]:=1.2:r[5]:=2:r[6]:=3.2:r[7]:=3.5:r[8]:=4:

> x0[1]:=0.5:x0[2]:=0.5:x0[3]:=0.5:x0[4]:=0.001:x0[5]:=0.001:x0[6]:=0.001:x0[7]:=0.001:x0[8]:=0.001:

> nmax:=400:imax:=8:

> for i from 1 to imax do

> ri:=r[i]:x[1]:=x0[i]:

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

> for n from 1 to nmax do

> x[n+1]:=ri*x[n]*(1-x[n]):

> od:

> ff.i:=plot(ri*y*(1-y),y=0..1,color=COLOR(RGB, rf, gf, bf)):

> xf.i:=plot([seq([n-1,x[n]+i-1],n=1..nmax)],color=COLOR(RGB, rf, gf, bf)):

> od:

> display(seq(ff.i,i=1..imax),view=[0..1,0..1],axes=boxed,title=` Logistische Abbildung für verschiedene Kontrollparameter r `);

> display(seq(xf.i,i=1..imax),view=[0..100,0..8.5],axes=boxed,title=` Lösungen der logistische Abbildung für verschiedene Kontrollparameter r `);

(Beachte, dass die Nulllinie um jeweils 1 nach größeren Werten verschoben wird). Man sieht die unterschiedlichen Lösungstopologien: Für r<1 konvergiert die Lösung gegen 0 (r=1 unendlich langsame Konvergenz), für 1<r<3.2 konvergiert die Lösung ohne Oszillationen gegen einen r-abhängigen Wert, bei r=3.2 beobachtet man eine Periode 1 Oszillation, bei r=3.5 eine Periode 2 Oszillation und bei r=4 schließlich eine chaotische Lösung.

>

2) Feigenbaum Diagramme der logistischen Abbildung

Ein Feigenbaum Diagramm stellt die asymptotische Lösung einer diskreten Abbildung als Funktion des Kontrollparameters dar. Man kann anhand dieser Abbildung die verschiedenen Topologien möglicher Lösungen untersuchen.

Iteration der logistischen Gleichung

> LogisticMap := proc(r)

> local n,x;

> x:=array(1..501):x[1]:=0.001:

> # Transienter Bereich

> for n from 1 to 299 do

> x[n+1]:=r*x[n]*(1-x[n]):

> od:

> # Beiträge zum Feigenbaum Diagramm

> for n from 300 to 500 do

> x[n+1]:=r*x[n]*(1-x[n]):

> od:

> seq([r,x[n]],n=300..500);

> end:

Prozedur für Feigenbaum Diagramm

> Feigenbaum := proc(rmin,rmax,imax)

> local dr,ri,xf,i,j;

> dr:=(rmax-rmin)/imax:ri:=rmin:

> for i from 1 to imax while ri<=4 do

> xf.i:=plot([LogisticMap(ri)],color=blue):

> ri:=rmin+i*dr:

> od:

> display(seq(xf.j,j=1..i-1),symbol=POINT,style=point, view=[rmin..rmax,0..1],axes=boxed,title=` Feigenbaum Diagramm der logistischen Abbildung `):

> end:

> Feigenbaum(1,4,200); # ACHTUNG: lange Rechenzeit

Man sieht deutlich das Einsetzen von Periodenverdopplungen bei r = 3, 3.45, 3.54. Um detailliertere Informationen zu erhalten, muss man das Bild vergrößern

> Feigenbaum(3.52,3.62,200); # ACHTUNG: lange Rechenzeit

und erhält die neuen Bifurkationswerte r = 3.56, 3.569....

Untersuchen Sie das Feigenbaum Diagramm und machen Sie sich die Strukturen klar.

>

3) Analytische Untersuchung der Fixpunkte

Ein Fixpunkt einer diskreten Abbildung ist definiert durch x=f(x), d.h. die Iteration bildet einen Fixpunkt in sich selbst ab.

> restart:

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> f:=unapply(r*x*(1-x),x,r); # rhs der logistischen Abbildung

> Fixpunkte:=solve(x=f(x,r),x); # Bedingung für Fixpunkte

Zur Bestimmung der Stabilität eines Fixpunktes, untersucht man eine Lösung aus der Umgebung dieser Punktlösung. Sei also x ein Fixpunkt und y(n)=x+dx(n) eine Lösung der Iteration mit einem Anfangswert aus der Umgebung von x, so gilt:

y(n+1) = f(y(n)) = f(x+dx(n)) = f(x) + f '(x) dx(n) + .....

y(n+1) = x + f '(x) dx(n) + ...

dx(n+1) = f '(x) dx(n) + ...

und bis zur ersten Ordnung:

dx(n) = dx(0) = exp( n) dx(0)

= ln(f '(x)).

Ist <0, konvergiert jede Lösung aus der Umgebung gegen den Fixpunkt: Man nennt den Fixpunkt stabil. Andernfalls ist der Fixpunkt instabil. heißt Lyapunovexponent.

> fs:=unapply(diff(f(x,r),x),x,r);

Lyapunovexponent:

> lambda:=unapply(ln(abs(fs(x,r))),x,r);

Für den ersten Fixpunkt erhält man (beachte r>0):

> rc1:=solve(lambda(Fixpunkte[1],r)=0,r);

Stabilität für r<1 und Instabilität für r>1.

Für den zweiten Fixpunkt erhält man nur Lösungen für r>1 und

> rc2:=solve(lambda(Fixpunkte[2],r)=0,r);

also im Bereich 1<r<3.

Für größere Werte von r beobachtet man ein neues Phänomen: Oszillationen. Dabei wiederholt sich für Periode k der Wert der Iteration in jedem -ten Iterationsschritt. Somit erhält man als Bedingung einer Periode k=1 Oszillation

> eqs2:=x=f(f(x,r),r);

> Periode1:=solve(eqs2,x);

Achtung: der dritte Ausdruck liefert nur dann eine reelle Lösung der logistischen Gleichung, wenn der Radikand positiv, d.h.

> rc3:=solve(r^2-2*r-3,r);

r>3 ist (negative Werte für r sind nicht erlaubt). Somit haben wir für r=3 die erste Bifurkation identifiziert, die Aufgabelung von einer Periode 0 Lösung (keine Oszillation, ein asymptotisch stabiler Wert der Iteration) in eine Periode 1 Lösung (Oszillation mit zwei Werten).

Den nächsten kritischen Wert des Kontrollparameters r erhalten wir beim Übergang von Periode 1 nach Periode 2 Lösungen:

> eqs4:=0=x-f(f(f(f(x,r),r),r),r);

> sol:=solve(eqs4,x);

Die entscheidende Lösung erhält man aus einem Polynom 12. Ordnung. Sie lässt sich nicht mehr einfach analytisch bestimmen. Betrachtet man das Feigenbaumdiagramm, so kann man für den nächsten Bifurkationspunkt ablesen

> rc4:=1+sqrt(6);

a) Vegleichen Sie die obigen Aussagen mit den numerischen Ergebnissen.

b) Verifizieren Sie die folgende Aussage: Trägt man die Kontrollparameter der Bifurkationspunkte gegen die Ordnung auf, so erhält man eine geometrische Reihe der Form

r(k) = -

Dabei ist derjenige Wert für r, bei dem Chaos einsetzt, und q sind unbestimmte Konstanten.

Im nächsten Schritt soll bestimmt werden.

> rck:=unapply(rho-alpha*q^k,alpha,k,q,rho);

> Fk:=(rck(alpha,k+1,Q,rho)-rck(alpha,k,Q,rho))/(rck(alpha,k+2,Q,rho)-rck(alpha,k+1,Q,rho));

> simplify(Fk,{Q=q}); # Achtung Trick: der Umweg über Q=q ist notwendig, damit MAPLE sich 'bequemt', den Bruch zu kürzen

Im Grenzfall k -> entspricht Fk der Feigenbaumkonstante . D.h.

> q:=1/delta;

bestimmt sich aus dem ersten Bifurkationswert:

> solalpha:=solve(rck(alpha,1,q,rho)=rc3[2],alpha);

und kann man schließlich aus dem zweiten Bifurkationswert ermitteln:

> solrho:=solve(rck(solalpha,2,q,rho)=rc4,rho);

Als Zahlenwert erhält man mit

> delta:=4.6692016091;

> rho:=evalf(solrho);

Nun können Sie alle Bifurkationspunkte berechnen. Versuchen Sie es und überprüfen Sie Ihre Ergebnisse mit Hilfe der numerischen Lösungen des 2. Abschnittes.