Was du hier lernst
- Differentialgleichungen
- Lagrange-Formalismus
- Runge-Kutta-Verfahren
- Energieerhaltung als Kontrolle
- Deterministisches Chaos
- Logarithmische Skala
Wenn ein Pendel nicht reicht
In der Danksagung im Vorwort bedanke ich mich bei Tino Wagner für eine sehr schöne Implementierung eines Doppelpendels. Sie hat mich damals dazu verleitet, im Buch etwas zur Chaostheorie zu bringen und selbst ein einfaches Pendel zu programmieren – daraus wurde das Kapitel über die Pendelkette. Das Doppelpendel selbst stand nie im Buch. Auf dieser Website gibt es jetzt eine eigene, neue Umsetzung, die im Browser läuft.
Ein Doppelpendel ist schnell gebaut: An die Kugel eines Pendels hängt man ein zweites Pendel. Die Bewegung, die dabei herauskommt, ist dagegen alles andere als einfach. Das obere Pendel schwingt, reißt das untere mit, das untere überschlägt sich, bremst das obere – und kein Schwung gleicht dem vorigen. Das Doppelpendel ist das wohl bekannteste Beispiel für deterministisches Chaos: Die Bewegung folgt strengen Gesetzen, und trotzdem lässt sie sich nicht vorhersagen.
Keine Formel mehr
Bei der Pendelkette gab es eine Formel, die den Winkel für jeden Zeitpunkt direkt liefert – eine Näherung, aber eine gute. Beim Doppelpendel gibt es so eine Formel nicht, auch keine ungenaue. Schon im Pendelkapitel des Buchs steht der Ausweg: ein numerisches Verfahren wie das Runge-Kutta-Verfahren. Statt die Bewegung auf einmal auszurechnen, gehen wir in winzigen Zeitschritten vor. Aus dem jetzigen Zustand berechnen wir, wie er sich im nächsten Augenblick ändert, dann im übernächsten und so weiter.
Der Zustand eines Doppelpendels besteht aus vier Zahlen: den Winkeln θ₁ und θ₂ der beiden Stäbe, gemessen von der Senkrechten, und ihren Winkelgeschwindigkeiten ω₁ und ω₂. Was wir brauchen, sind die Winkelbeschleunigungen – also die Antwort auf die Frage, wie sich die Geschwindigkeiten gerade ändern.
Die Bewegungsgleichungen
Die Kräfte am Doppelpendel direkt aufzuschreiben, ist mühsam: Die Stäbe ziehen aneinander, und die Richtungen ändern sich ständig. Eleganter geht es über die Energie, mit dem Lagrange-Formalismus. Man schreibt nur zwei Größen auf, die Bewegungsenergie T und die Lageenergie V:
T = ½·(m₁ + m₂)·l₁²·ω₁² + ½·m₂·l₂²·ω₂² + m₂·l₁·l₂·ω₁·ω₂·cos(θ₁ − θ₂)
V = −(m₁ + m₂)·g·l₁·cos θ₁ − m₂·g·l₂·cos θ₂
Daraus bildet man L = T − V. Die Euler-Lagrange-Gleichungen, d/dt (∂L/∂ωᵢ) − ∂L/∂θᵢ = 0, liefern für jeden Winkel eine Gleichung. Die beiden Gleichungen hängen voneinander ab, lassen sich aber nach den beiden Winkelbeschleunigungen auflösen. Das Ergebnis ist länglich, aber für den Computer kein Problem:
def derivatives(state):
# die Bewegungsgleichungen aus dem Lagrange-Formalismus: aus Winkeln und
# Winkelgeschwindigkeiten werden die Winkelbeschleunigungen
theta1, omega1, theta2, omega2 = state
delta = theta1 - theta2
den = 2 * mass1 + mass2 - mass2 * math.cos(2 * delta)
alpha1 = (-gravity * (2 * mass1 + mass2) * math.sin(theta1)
- mass2 * gravity * math.sin(theta1 - 2 * theta2)
- 2 * math.sin(delta) * mass2
* (omega2 ** 2 * length2 + omega1 ** 2 * length1 * math.cos(delta))) / (length1 * den)
alpha2 = (2 * math.sin(delta)
* (omega1 ** 2 * length1 * (mass1 + mass2)
+ gravity * (mass1 + mass2) * math.cos(theta1)
+ omega2 ** 2 * length2 * mass2 * math.cos(delta))) / (length2 * den)
return (omega1, alpha1, omega2, alpha2)
Die Funktion bekommt den Zustand und liefert zurück, wie schnell sich jede seiner vier Zahlen ändert: Die Winkel ändern sich mit ω₁ und ω₂, die Winkelgeschwindigkeiten mit den gerade berechneten Beschleunigungen α₁ und α₂. Beide Stäbe sind hier einen Meter lang, beide Kugeln wiegen ein Kilogramm, und die Stäbe selbst gelten als masselos.
Schritt für Schritt: Runge-Kutta
Das einfachste Verfahren, aus solchen Änderungsraten eine Bewegung zu machen, ist das Euler-Verfahren: Man nimmt die Steigung am Anfang eines Zeitschritts und geht damit geradeaus weiter. Genau so rechnet die Mondlandung. Das Problem: Eine Bahn ist fast nie gerade, und mit jedem Schritt entsteht ein kleiner Fehler. Beim Doppelpendel summieren sich diese Fehler so schnell, dass das Pendel nach wenigen Sekunden Energie aus dem Nichts gewinnt.
Das klassische Runge-Kutta-Verfahren vierter Ordnung ist schlauer. Es probiert die Steigung an vier Stellen aus – am Anfang des Schritts, zweimal in der Mitte und am Ende – und bildet daraus einen gewichteten Mittelwert:
def rungeKutta(state, h):
# das klassische Runge-Kutta-Verfahren vierter Ordnung:
# vier Steigungen ausprobieren und gewichtet mitteln
k1 = derivatives(state)
k2 = derivatives([s + h / 2 * k for s, k in zip(state, k1)])
k3 = derivatives([s + h / 2 * k for s, k in zip(state, k2)])
k4 = derivatives([s + h * k for s, k in zip(state, k3)])
return [s + h / 6 * (a + 2 * b + 2 * c + d) for s, a, b, c, d in zip(state, k1, k2, k3, k4)]
k2 ist die Steigung in der Mitte des Schritts, wenn man mit k1 dorthin gegangen wäre; k3 dieselbe Mitte, aber mit der besseren Steigung k2 erreicht; k4 die Steigung am Ende. Die Mitten zählen doppelt. Der Fehler pro Schritt schrumpft dadurch mit der fünften Potenz der Schrittweite: Halbierst du den Schritt, wird er 32-mal kleiner. Das Programm macht pro Bild vier solcher Schritte, bei 60 Bildern pro Sekunde also 240 Schritte in der Sekunde.
Energie als Kontrolle
Woher weiß man, dass eine Simulation richtig rechnet? Beim reibungsfreien Pendel gibt es eine eingebaute Probe: Die Gesamtenergie aus Bewegungs- und Lageenergie muss gleich bleiben. Die Funktion energy berechnet genau die T + V von oben, und oben links auf der Leinwand steht der Wert des roten Pendels neben dem Startwert. Mit Runge-Kutta ändert er sich in 40 Sekunden um weniger als zwei Zehntausendstel Joule – bei gut 22 Joule Gesamtenergie. In der Anzeige mit drei Nachkommastellen bleibt er stehen. Mit dem Euler-Verfahren wären es nach zehn Sekunden über sechs Joule mehr. Probier es unter »Probier mal« aus.
Genau hingeschaut
Dass die Energie stimmt, heißt nicht, dass die Bahn stimmt. Energie ist nur eine Zahl, die Bahn sind vier. Beim Doppelpendel kann selbst ein perfektes Verfahren die Bahn nicht lange genau berechnen – warum, zeigt der nächste Abschnitt.
Zwei Pendel, ein Tausendstel Grad
Jetzt kommt das eigentliche Experiment. Das Programm startet zwei Doppelpendel, ein rotes und ein blaues, beide aus der Klasse DoublePendulum. Sie sind in allem gleich, nur ist beim blauen der untere Stab um ein Tausendstel Grad weiter ausgelenkt:
return (DoublePendulum(angle1, angle2, red),
DoublePendulum(angle1, angle2 + difference, blue))
Die untere Kugel des blauen Pendels liegt damit anfangs gerade einmal 17 Mikrometer neben der roten – weniger als die Dicke eines Haares. Auf dem Bildschirm liegen beide Pendel genau übereinander, du siehst nur das blaue. Doch dann wächst der Abstand, erst unmerklich, dann immer schneller. Nach etwa drei Sekunden sind es Millimeter, nach sechs Sekunden zehn Zentimeter, und kurz danach schwingen die Pendel, als hätten sie nie etwas miteinander zu tun gehabt.
Das Diagramm unten zeigt den Abstand der unteren Kugeln auf einer logarithmischen Skala: Jede Linie steht für das Tausendfache der darunter. Auf dieser Skala wird exponentielles Wachstum zu einer ansteigenden Geraden. Im Mittel verdoppelt sich der Abstand etwa alle halbe Sekunde, bis er so groß ist wie das Pendel selbst. Genau das ist Chaos: Kleine Unterschiede werden nicht einfach größer, sie werden in gleichen Zeitabständen um denselben Faktor größer.
Es ist dasselbe Prinzip wie bei der Zahlenfolge im Chaos-Abschnitt des Buchs, die du auf der Seite Julia-Menge und Chaos findest: Dort verdoppelt sich ein winziger Unterschied bei jedem Schritt, bis er alles bestimmt. Und es hat eine bittere Folge für jede Vorhersage. Machst du den Unterschied am Start eine Million Mal kleiner, gewinnst du nicht eine Million Mal mehr Zeit, sondern nur ein paar Sekunden – die Variante »Eine Million Mal genauer« zeigt es. Deshalb reicht der Wetterbericht nur einige Tage in die Zukunft, egal wie gut die Messgeräte werden.
Spuren zeichnen
Die geschwungenen Linien sind die Spuren der unteren Kugeln. Jedes Pendel merkt sich seine letzten 240 Positionen in einer deque mit fester Länge: Kommt ein neuer Punkt hinzu, fällt der älteste von selbst heraus.
self.trail.append(toScreen(positions(self.state)[1]))
Gezeichnet wird die Spur in Stücken von 20 Punkten, und jedes ältere Stück bekommt eine geringere Deckkraft. So verblasst die Spur nach hinten, und du siehst, woher die Kugel gerade kommt. Ein Klick ins Bild startet das Experiment mit neuen, zufälligen Winkeln zwischen 90 und 180 Grad – also immer mit so viel Schwung, dass das Chaos eine Chance hat.
Probier mal
Jede Variante ändert den Code oben. Ein Klick auf »Ausprobieren« übernimmt die Änderung in den Editor und startet das Programm. »Zurücksetzen« holt das Original zurück.
-
01Euler statt Runge-Kutta
Das einfachste Verfahren für Differentialgleichungen ist das Euler-Verfahren: Steigung am Anfang des Schritts nehmen und geradeaus weitergehen. Genau so rechnet übrigens die Mondlandung. Beim Doppelpendel geht das schief: Schau auf die Energie oben links. Sie wächst mit jeder Sekunde, nach zehn Sekunden hat das Pendel über ein Viertel mehr Energie als am Start – es beschleunigt sich aus dem Nichts.
- class DoublePendulum(object): + def euler(state, h): + # das Euler-Verfahren: nur die Steigung am Anfang des Schritts + k = derivatives(state) + return [s + h * d for s, d in zip(state, k)] + + + class DoublePendulum(object): - self.state = rungeKutta(self.state, h) + self.state = euler(self.state, h) -
02Eine Million Mal genauer
Vielleicht war ein Tausendstel Grad einfach zu viel? Hier unterscheiden sich die Pendel nur noch um ein Milliardstel Grad – eine Million Mal weniger. Trotzdem trennen sie sich, nur etwa neun Sekunden später: nach gut 15 statt nach 6 Sekunden. Im Diagramm siehst du, warum. Der Abstand wächst jede Sekunde um denselben Faktor, eine Million mehr Genauigkeit kauft also nur ein paar Sekunden Vorhersage. Genau deshalb hält der Wetterbericht nicht ewig.
- difference = 0.001 # so viel Grad ist das blaue Pendel unten weiter ausgelenkt + difference = 0.000000001 # ein Milliardstel Grad -
03Kleine Auslenkung – kein Chaos
Chaos braucht Energie. Starten beide Stäbe mit nur 10 Grad Ausschlag, schwingt das Doppelpendel brav und regelmäßig, und der Abstand der Kugeln bleibt winzig – die Linie im Diagramm wächst nicht. Für kleine Winkel gilt wieder die Kleinwinkelnäherung aus der Pendelkette. Ein Klick startet übrigens wieder mit großen, zufälligen Winkeln.
- startAngle1 = 130.0 + startAngle1 = 10.0 - startAngle2 = 170.0 + startAngle2 = 10.0 -
04Mit Luftreibung
Nach jedem Rechenschritt verlieren beide Kugeln ein wenig Winkelgeschwindigkeit, wie durch Luftwiderstand. Anfangs ändert das wenig: Die Pendel trennen sich wie gewohnt. Dann sinkt die Energie oben links, die Ausschläge werden kleiner, und im Diagramm kommen sich die Kugeln langsam wieder näher. Nach ein paar Minuten hängen beide still. Mit Reibung ist das Ziel vorhersagbar, nur der Weg dorthin nicht.
- self.state = rungeKutta(self.state, h) + self.state = rungeKutta(self.state, h) + self.state[1] *= 1 - 0.1 * h # Luftreibung bremst beide Kugeln + self.state[3] *= 1 - 0.1 * h
Hinter den Kulissen
Das Doppelpendel, für das ich mich im Vorwort bei Tino Wagner bedanke, lag den Quelltexten zum Buch als Beispiel bei. Es zeichnete mit pyglet und OpenGL in ein eigenes Fenster und brachte ein eigenes Modul zum Lösen von Differentialgleichungen mit. Das Programm auf dieser Seite ist davon unabhängig neu geschrieben: in reinem Python, ohne Zusatzbibliotheken, mit der Zeichen-API c4f für die Leinwand im Browser.
Aufbau. Die Physik steckt in vier Funktionen: derivatives (die Bewegungsgleichungen), rungeKutta (ein Zeitschritt), positions (Winkel in Kugelpositionen) und energy (die Probe). Die Klasse DoublePendulum hält Zustand, Farbe und Spur eines Pendels zusammen – so wie die Klasse Pendulum in der Pendelkette. Alle Werte sind in SI-Einheiten: Meter, Kilogramm, Sekunden. Erst toScreen rechnet mit 115 Pixeln pro Meter in Leinwandkoordinaten um.
Zeit. Das Programm rechnet mit fester Schrittweite: 60 Bilder pro Sekunde, vier Runge-Kutta-Schritte pro Bild, also h = 1/240 Sekunde. await screen.frame(fps) hält die Bildrate so gleichmäßig wie möglich. Kommt der Browser einmal nicht hinterher, läuft die Simulation etwas langsamer, aber nie ungenauer – die Schrittweite bleibt gleich.
Warum kein NumPy? Für vier Zahlen pro Pendel lohnt sich keine Bibliothek. Die Listen-Ausdrücke in rungeKutta sind kurz genug, und zwei Pendel mit je 240 Schritten pro Sekunde rechnet Python auch im Browser mühelos.
Grenzen des Modells. Die Stäbe sind masselos und starr, die Kugeln punktförmig, und ohne die Variante »Mit Luftreibung« gibt es keine Reibung. Ein echtes Doppelpendel aus Holz oder Metall verhält sich deshalb im Detail anders – chaotisch ist es aber genauso. Und auch der Computer hat eine Grenze: Gleitkommazahlen haben etwa 16 gültige Stellen. Selbst ohne absichtlichen Unterschied würden zwei Rechnungen mit minimal anderer Reihenfolge der Operationen irgendwann auseinanderlaufen.
