Das mathematische Pendel

H.J. Lüdde

nach

N.J. Giordano, Computational Physics, Kapitel 3

Eines der einfachsten mechanischen Systeme ist das mathematische Pendel: an einem starren, masselosen 'Faden' der Länge l hängt ein Massepunkt der Masse m. Das System wird durch die einzige generalisierte Koordinate bestimmt, die Rückstellkraft ist nichtlinear. Diese Tatsache führt zu vielfältigen Bewegungsformen, insbesondere dann, wenn das Pendel durch eine äußere Kraft angeregt wird. Unter bestimmten Bedingungen wird die Pendelbewegung so kompliziert, dass sie nicht mehr voraussagbar erscheint. Trotzdem ist sie vollständig bestimmt: man nennt sie chaotisch.

In diesem Arbeitsblatt wird in Abschnitt 1) das freie mathematische Pendel mit seiner linearen Näherung (harmonischer Oszillator) verglichen. Abschnitt 2) befasst sich mit möglichen Bewegungsformen des getriebenen Pendels und Abschnitt 3) mit einer ersten qualitativen Analyse von Chaos.

1) Das freie mathematische Pendel

1.1) Das mathematische Pendel im Vergleich mit seiner linearen Näherung (harmonischer Oszillator)

> restart:

> with(plots):

Zunächst wollen wir uns einen Überblick darüber verschaffen, wie gut man die nichtlineare Rückstellkraft in einer Potenzreihe entwickeln kann:

> sp:=array(1..4):

> s:=taylor(sin(theta),theta=0,9);

> ii:=1:

> for i from 1 to 4 do

> ii:=ii+2:

> s:=taylor(sin(theta),theta=0,ii):

> sp[i]:=unapply(convert(s,polynom),theta):

> od:

> plot([seq(sp[i](theta),i=1..4),sin(theta)],theta=0..2*Pi,y=-1..1);

Für kleine Auslenkungen ( < 0.4) ist die lineare Näherung akzeptabel.

Wir definieren die DGL für das Pendel (ps1,ps2) bzw seine Linearisierung (ps1,px2):

> g:='g':

> l:='l':

> ps1:=diff(theta(t),t)=omega(t);

> ps2:=diff(omega(t),t)=-(g/l)*sin(theta(t));

> px2:=diff(omega(t),t)=-(g/l)*theta(t);

und lösen beide Differentialgleichungssysteme numerisch für g/l= =1:

> g:=9.81:

> l:=9.81:

> om0:=0:

> th0:=130:

> sols:=dsolve({ps1,ps2,theta(0)=(Pi/180)*th0,omega(0)=om0},{theta(t),omega(t)},type=numeric):

> solx:=dsolve({ps1,px2,theta(0)=(Pi/180)*th0,omega(0)=om0},{theta(t),omega(t)},type=numeric):

> fs:=odeplot(sols,[t,theta(t)],0..100,numpoints=500, color=green):

> fx:=odeplot(solx,[t,theta(t)],0..100,numpoints=500, color=blue):

> display(fs,fx,view=[0..20,-Pi..Pi],axes=boxed,labels=[t,theta], title=` Mathematisches Pendel (grün) vs harmonische Näherung (blau) `);

Man sieht einen deutlichen Unterschied zwischen den beiden Lösungen bei großen Auslenkungen. Das gleiche Bild erhält man im Phasenraum:

> fs:=odeplot(sols,[theta(t),omega(t)],0..100,numpoints=500, color=green):

> fx:=odeplot(solx,[theta(t),omega(t)],0..100,numpoints=500, color=blue):

> display(fs,fx,view=[-Pi..Pi,-Pi..Pi],axes=boxed,labels=[theta,omega], title=` Mathematisches Pendel (grün) vs harmonische Näherung (blau) `);

a) Interpretieren Sie das Ergebnis und variieren Sie die Anfangsbedingungen des Pendels.

1.2) Phasenraum Portrait des mathematischen Pendels

Überlegen wir einmal, welche Bewegungsformen wir bei großen Auslenkungen erwarten: da ist zunächst die gewohnte Oszillation um die Nullpunktlage für Anfangsauslenkungen < . Was passiert aber, wenn wir eine Anfangsauslenkung = wählen und dem Pendel gleichzeitig einen winzigen Stoß versetzen? Das Pendel wird um seinen Aufhängepunkt rotieren. Welche Bewegungsform realisiert wird hängt natürlich von der Anfangsbedingung ab, oder weil das System konservativ ist, von seiner Gesamtenergie. Das Phasenraum Portrait des Pendels zeigt die möglichen Trajektorien in Abhängigkeit von der Gesamtenergie.

> restart:

> with(plots):

Zunächst definiert man kinetische und potentialle Energie des Pendels:

> g:='g':

> m:='m':

> l:='l':

> T:=1/2*m*l^2*omega^2;

> U:=m*g*l*(1-cos(theta));

Für verschiedene Werte der Gesamtenergie E=T+U findet man in dem folgenden Contourplot die entsprechenden Trajektorien . Die Grenztrajektorie mit dem Wert E=2gl entspricht dem Fall, dass das Pendel bei startet und nach einem Umlauf (links bzw rechts) wieder bei stoppt. Für E>2gl rotiert das Pendel und für Energien E<2gl oszilliert es.

> g:=9.81:

> m:=1:

> l:=9.81: # Damit ist g/l=1

> contourplot(T+U,theta=-3/2*Pi..3/2*Pi,omega=-4..4,contours=[25,50,100,192.5,400,600,800],numpoints=1000,axes=boxed,filled=true, coloring=[red,black],title=`Phasenraum Portrait des mathematischen Pendels`);

b) Berechnen Sie für verschiedene Werte der Gesamtenergie .

2) Das getriebene mathematische Pendel - Chaotische Bewegungsform

Wir betrachten ein Pendel mit geschwindigkeitsabhängiger Reibung (Stokes) und äußerer periodischer Anregung (z.B. ein periodisches elektrisches Feld, das auf den geladenen Pendelkörper einwirkt). In Abhängigkeit von der Kopplungsstärke Fex=0, 0.5, 1.2 der äußeren Anregung beobachtet man verschiedene Bewegungsformen:

> restart:

> with(plots):

> g:='g':

> l:='l':

> k:='k':

> Omega:='Omega':

> Fex:='Fex':

> pf1:=diff(theta(t),t)=omega(t);

> pf2:=diff(omega(t),t)=-(g/l)*sin(theta(t))-k*omega(t)+Fex*sin(Omega*t);

> g:=9.81:

> l:=9.81:

> k:=0.5:

> Omega:=2/3:

> th0:=Pi/2:

> om0:=0:

> F:=array(1..3):

> F[1]:=0:

> F[2]:=0.5:

> F[3]:=1.2:

> f:=array(1..3):

> f1:=array(1..3):

> f2:=array(1..3):

Es folgt eine längere Rechnung. Bitte etwas Geduld!

> for i from 1 to 3 do

> Fex:=F[i]:

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

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

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

> sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure):

> f[i]:=odeplot(sol,[t,theta(t)],0..60,numpoints=500,color=COLOR(RGB, rf, gf, bf)):

> f1[i]:=odeplot(sol,[theta(t),omega(t)],0..60,numpoints=500,color=COLOR(RGB, rf, gf, bf)):

> f2[i]:=odeplot(sol,[theta(t),omega(t)],0..600,numpoints=2000,color=COLOR(RGB, rf, gf, bf)):

> od:

> display(f[1],f[2],f[3],view=[0..60,-4..4],axes=boxed,labels=[t,theta], title=` Getriebenes Pendel für Fex=0,0.5,1.2 `);

c) Diskutieren Sie die Ergebnisse für die drei Fälle. Wieso kommt es zu einer Oszillation mit konstanter Amplitude bei Fex=0.5?

> display(f1[1],f1[2],f1[3],view=[-4..4,-Pi..Pi],axes=boxed,labels=[theta,omega], title=` Getriebenes Pendel für Fex=0,0.5,1.2: Phasenraumdarstellung (t=0..60) `);

Man erkennt deutlich die unterschiedlichen Topologien der Lösungen: Fex=0 ergibt eine Trajektorie mit dem kritischen Punkt (0,0) als Fokus (positiver Attraktor), die Trajektorie zu Fex=0.5 zeigt asymptotisch ein Grenzzyklus Verhalten und die Trajektorie zu Fex=1.2 gibt uns zunächst ein Rätsel auf. Möglicherweise stabilisiert sich diese Trajektorie erst nach längerer Zeit?

> display(f2[1],f2[2],f2[3],view=[-5..70,-Pi..Pi],axes=boxed,labels=[theta,omega], title=` Getriebenes Pendel für Fex=0,0.5,1.2: Phasenraumdarstellung (t=0..600) `);

Offenbar nicht!! Diese Bewegungsform unterscheidet sich vollständig von den anderen. Sie erscheint vollkommen unvorhersehbar und trotzdem werden wir im nächsten Abschnitt erste Schritte zu ihrer Analyse kennen lernen.

d) Variieren Sie Anfangsbedingungen und Kopplungsstärke der äußeren Anregung.

3) Was ist Chaos?

3.1) Lyapunov Exponent

Wir werden in diesem Abschnitt die Bewegungsform für Fex=1.2 etwas genauer untersuchen. Insbesondere betrachten wir den Fall, was passiert, wenn wir die Anfangsbedingungen minimal verstimmen. Dazu berechnen wir nochmals die Lösung von Abschnitt 2) und parallel dazu ein Lösung deren Anfangsbedingung sich nur um 1ppm von der alten unterscheidet. Damit wir sehen, wie sich die Lösungen relativ zueinander verhalten, zeichnen wir die Differenz der Lösungen als Funktion von t.

> restart:

> with(plots):

> g:='g':

> l:='l':

> k:='k':

> Omega:='Omega':

> Fex:='Fex':

> pf1:=diff(theta(t),t)=omega(t);

> pf2:=diff(omega(t),t)=-(g/l)*sin(theta(t))-k*omega(t)+Fex*sin(Omega*t);

> g:=9.81:

> l:=9.81:

> k:=0.5:

> Omega:=2/3:

> th0:=Pi/2:

> om0:=0:

> F:=array(1..3):

> F[1]:=0:

> F[2]:=0.5:

> F[3]:=1.2:

> f3:=array(1..3):

Es folgt eine längere Rechnung. Bitte etwas Geduld!

> for i from 1 to 3 do

> Fex:=F[i]:

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

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

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

> sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure):

> th:=subs(sol,theta(t)):

> sol1:=dsolve({pf1,pf2,theta(0)=1.000001*th0,omega(0)=1.000001*om0},{theta(t),omega(t)},type=numeric,output=listprocedure):

> th1:=subs(sol1,theta(t)):

> f3[i]:=logplot('abs(th1(t)-th(t))',t=0..100,numpoints=500,color=COLOR(RGB, rf, gf, bf)):

> od:

> display(f3[1],f3[2],f3[3],view=[0..50,-10..-2],axes=boxed,labels=[t,dtheta], title=` Getriebenes Pendel für Fex=0,0.5,1.2: AB um 1ppm unterschieden `);

dtheta ist die Differenz der Lösungen, die sich in der Angangsbedingung minimal unterscheiden. Während für die beiden stabilen Fälle (Fex=0, 0.5) die Differenz im Mittel exponentiell gegen null strebt, wächst die Differenz für Fex=1.2 im Mittel exponentiell an (beachte die halblogarithmische Darstellung). Dabei handelt es sich offenbar um ein Kriterium für Chaos: Trajektorien zu verschiedenen Anfangsbedingungen entfernen sich bei chaotischen Bewegungsformen exponentiell voneinander. Die Steigung der mittleren Differenz dtheta(t) nennt man Ljapunovexponent. Er ist ein Maß für Chaos.

e) Variieren Sie die Anfangsbedingung in sol1.

f) Untersuchen Sie qualitativ den Ljapunovexponent für verschiedene Fex.

g) Machen Sie sich die komplizierte Bewegung des chaotischen Pendels klar, insbesondere mit Hilfe des Phasenraum Portraits.

3.2) Poincare Schnitte

Ein nützliches Hilfsmittel zur Analyse chaotischer Systeme sind Poincare Schnitte: man zeichnet nur solche Punkte des Phasenraumes, die in Phase mit der äußeren Anregung sind. Diese Bedingung wirkt in Analogie zu einem Stroboskop. Betrachten wir eine schnell rotierende Scheibe mit einer Aufschrift, so lässt sich der Text nicht lesen. Beleuchtet man die Scheibe mit einem Stroboskop in Phase (also immer dann wenn die Schrift aufrecht steht) so kann man den Text einwandfrei entziffern.

Im folgenden zeichnen wir alle Phasenraumpunkte der chaotischen Trajektorie, für die gilt .

Zunächst schreibt man die Pendelgleichung als autonomes System (siehe Vorlesung) und löst sie dann über einen langen Zeitraum (jj Teillösungen mit jeweils 50 Poincare Punkten). Damit das 'Punktesammeln' effektiver wird reduziert man die Pendelbewegung auf das Intervall , d.h. immer wenn setzt man (periodische Randbedingungen).

> restart:

> with(plots):

> g:='g':l:='l':k:='k':Fex:='Fex':Omega:='Omega':

> pf1:=diff(theta(t),t)=omega(t);

> pf2:=diff(omega(t),t)=-(g/l)*sin(theta(t))-k*omega(t)+Fex*sin(phi(t));

> pf3:=diff(phi(t),t)=Omega;

> g:=9.81:l:=9.81:k:=0.5:Fex:=1.2:Omega:=2/3: # Parameter für chaotische Lösung

> th0:=Pi/2:om0:=0:ph0:=0:jj:=200: # jj Lösungen a jeweils 50 Punkte (ACHTUNG: mit jj verlängert sich die Rechenzeit dramatisch)

> ic:=theta(0)=th0,omega(0)=om0,phi(0)=ph0: # Anfangsbedingungen

> for j from 1 to jj do # Schleife über kk Einzelkurven

> vars:=theta(t),omega(t),phi(t);

> sol:=dsolve({pf1,pf2,pf3,ic}, {vars},type=numeric,output=listprocedure);

> assign(sol);

> theta:=theta(t);

> omega:=omega(t);

> phi:=phi(t);

> npt:=50; # Zahl der Punkte in einer Kurve

> pts:=seq([theta(Pi*i/Omega)-trunc((theta(Pi*i/Omega)+sign(theta(Pi*i/Omega))*Pi)/(2*Pi))*2*Pi,omega(Pi*i/Omega)],i=1..npt); # (i) betrachte nur Punkte zu den 'Stroboskopzeiten' (ii) reduziere theta auf [-Pi..Pi]

> plt(j):=plot([pts],view=[-Pi..Pi,-2..1],color=blue):

> th0:=`th0`:om0:=`om0`:ph0:=`ph0`:

> th0:=theta(Pi*npt/Omega)-trunc((theta(Pi*npt/Omega)+sign(theta(Pi*npt/Omega))*Pi)/(2*Pi))*2*Pi:om0:=omega(Pi*npt/Omega):ph0:=phi(Pi*npt/Omega):

> unassign(theta,omega,phi);

> ic:=theta(0)=th0,omega(0)=om0,phi(0)=ph0; # neue Anfangsbedingungen

> od:

> display([seq(plt(j),j=1..jj)],symbol=POINT,style=point,view=[-Pi..Pi,-2..1],axes=boxed,labels=[theta,omega], title=` Poincare Schnitt für Fex=1.2 `);

Mit einer wachsenden Zahl jj von Teillösungen erhalten Sie ein deutlicheres Bild

Man sieht eine erstaunlich regelmäßige Struktur! Hätten man dieses Bild für einen der beiden stabilen Fälle Fex=0, 0.5 gezeichnet, so würde man genau EINEN Punkt sehen, da die Bewegungsformen in diesem Fall durch einen Attraktor (Fokus bzw Grenzzyklus) bestimmt sind. In Analogie dazu bezeichnet man die regelmäßige Struktur im Poincare Schnitt der chaotischen Bewegungsform als 'seltsamen Attraktor'.

h) Zeigen Sie, dass die Topologie des seltsamen Attraktors eine Eigenschaft des nichtlinearen Systems ist, also nicht von der Anfangsbedingung abhängt.

3.3) Periodenverdopplung - Wege zum Chaos

Bisher haben wir zwei typische Eigenschaften chaotischer Bewegungsformen diskutiert: (i) die extreme Abhängigkeit einer Lösung von der Anfangsbedingung und (ii) die charakteristische topologische Struktur ihrer Poincare Darstellung.

In diesem Kapitel wollen wir den Übergang von einer regelmäßigen in eine chaotische Bewegungsform untersuchen. Ist dieser Übergang singulär oder erkennt man eine kontinuierliche Veränderung der Bewegungsform in Abhängigkeit von einem Kontrollparameter?

Dazu betrachten wir noch einmal drei Lösungen des getriebenen Pendels für die Kopplungsstärken Fex=1.35, 1.44, 1.465:

> restart:

> with(plots):

> g:='g':

> l:='l':

> k:='k':

> Omega:='Omega':

> Fex:='Fex':

> pf1:=diff(theta(t),t)=omega(t):

> pf2:=diff(omega(t),t)=-(g/l)*sin(theta(t))-k*omega(t)+Fex*sin(Omega*t):

> g:=9.81:

> l:=9.81:

> k:=0.5:

> Omega:=2/3:

> th0:=0.1:

> om0:=0:

> F:=array(1..3):

> F[1]:=1.35:

> F[2]:=1.44:

> F[3]:=1.465:

Es folgt eine längere Rechnung. Bitte etwas Geduld!

> for j from 1 to 3 do

> Fex:=F[j]:

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

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

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

> sol:=dsolve({pf1,pf2,theta(0)=th0,omega(0)=om0},{theta(t),omega(t)},type=numeric,output=listprocedure):

> assign(sol);

> theta:=theta(t);

> npt:=200; # Zahl der Punkte in einer Kurve

> pts:=seq([i/2,theta(i/2)-trunc((theta(i/2)+sign(theta(i/2))*Pi)/(2*Pi))*2*Pi-(j-1)*2*Pi],i=1..npt); # reduziere theta auf [-Pi..Pi]

> f(j):=plot([pts],color=COLOR(RGB, rf, gf, bf)):

> unassign(theta,omega);

> od:

> display([seq(f(j),j=1..3)],view=[0..100,-5*Pi..Pi],axes=boxed,labels=[t,theta], title=` Getriebenes Pendel für Fex=1.35, 1.44, 1.465 `);

>

Nach einem Einschwingvorgang (t=20) oszilliert das Pendel für Fex=1.35 mit der Periode der äußeren Anregung. Für Fex=1.44 hat sich jedoch die Periode verdoppelt und für Fex=1.465 vervierfacht (dies lässt sich nur für eine große Punktzahl npt verifizieren (Achtung: lange Rechenzeiten)). Trägt man die (auf das Intervall reduzierten Maxima in als Funktion von Fex auf, so erhält man ein Bifurkationsdiagramm (Bifurkation: Gabelung), das sich in immer kürzeren Abständen Fex immer weiter verzweigt. D.h. immer kleinere Änderungen im Kontrollparameter Fex führen zu weiteren Periodenverdopplungen (Verdopplung der Anzahl der verschiedenen Maxima) bis für einen Grenzwert F die Zahl der verschiedenen Maxima unendlich ist, und somit die Pendelbewegung nicht mehr periodisch sondern chaotisch wird. Im Phasenraum sieht man dies am Übergang von einer geschlossenen in eine offene Trajektorie.