Doppelpendel, wie lange stimmt die Rechnung?

Die Aufgabe

Das Beispiel stammt aus dem Buch von Cleve Moler: "Numerical Computing with MATLAB" (das Buch steht komplett zum Download zur Verfügung unter http://www.mathworks.com/moler/chapters.html).

Ein einfaches Doppelpendel besteht aus zwei Punktmassen, die mit masselosen starren Stangen verbunden sind, keine Reibung, kein Luftwiderstand). Die Differenzialgleichungen werden am Ende des Kapitels 7 hergeleitet:

Die Masse m₂ wird an einen durch die Koordinaten im x-y-System zu beschreibenden Punkt geführt, und das System wird danach ohne Anfangsgeschwindigkeiten freigelassen.

Das Function-Script swinger.m

Die in dem Buch von Cleve Moler beschriebenen Matlab-Scripts stehen unter http://www.mathworks.com/moler/ncmfilelist.html) zum Download zur Verfügung. Für das hier betrachtete Doppelpendel heißt die Datei swinger.m, ein ansehenswertes Script, das u. a. zeigt, wie man parallel zur Abarbeitung der ode-Functions die Zwischenergebnisse verarbeiten kann, in dem man z. B. eine Animation des Bewegungsablaufs erzeugt.

Im Script swinger.m erhalten alle Problemparameter (m₁, m₂, l₁, l₂, g) den Wert 1, so das sich für t keine der vertrauten Dimensionen für die Zeit ergibt, aber das ist für die Betrachtungen hier unerheblich. Das t-Intervall, für das die Lösung erzeugt wird, ist praktisch unendlich groß (obere Grenze: 10⁶), so dass die Bewegung solange läuft, bis der Benutzer den Stop-Button anklickt.

Die Anfangslage wird durch die Koordinaten der Masse m₂ festgelegt, hier soll eine Bewegung betrachtet werden, bei der das Doppelpendel eine horizontale Anfangslage einnimmt (nebenstehende Graphik links), die Koordinaten der Masse m₂ sind dafür x₂ = 2 ; y₂ = 0.

Das Script kann mit Angabe dieser beiden Koordinaten aus dem Command Window mit

swinger(2,0)

gestartet werden, und man sieht im Graphik-Window die Animation der Bewegung (rechts ein Schnappschuss bei t = 38). Die Animation läuft bis zum Abbruch durch den Benutzer. Die Frage ist, wie lange die dargestellte Bewegung tatsächlich das in der Aufgabenstellung formulierte Problem abbildet.

Wie lange stimmt die Rechnung?

Zur Untersuchung dieser Frage wird (unter Beibehaltung aller wesentlichen Passagen des Scripts swinger.m) ein vereinfachtes Script swinger1.m erzeugt (nebenstehend sind die wesentlichen Passagen dieses Scripts zu sehen). Dieses verzichtet auf die Animation zugunsten der graphischen Darstellung der Bewegungsgesetze φ₁(t) und  φ₂(t). Außerdem wird zweimal gerechnet, um den Einfluss des MaxStep-Parameters auf die Ergebnisse zu zeigen:

• In den Zeilen 12 und 13 werden die Einstellungen aus dem Script swinger.m übernommen (dort wird mit dem MaxStep-Wert 0,05 der Tatsache Rechnung getragen, dass für konplizierte Bewegungen dieser Art die Default-Einstellungen nicht ausreichend sind).
   
• In den Zeilen 16 und 17 wird die Rechnung mit dem deutlich kleineren Wert 0,005 für MaxStep wiederholt, der mit wesentlich mehr Aufwand (Integrationsschritten) ein deutlich genaueres Ergebnis erwarten lässt.

In den Zeilen 20 und 21 werden die Funktionen φ₁(t) und  φ₂(t) für beide Rechnungen in jeweils einem Graphik-Fenster dargestellt (die "Original-swinger-Rechnung" in rot, die "genauere Rechnung" in blau).

Die nebenstehende Graphik zeigt das Ergebnis. Am Anfang überdecken sich die beiden Kurven, aber ab etwa t = 45  ergeben sich deutliche Unterschiede, die dann sehr schnell dazu führen, dass die ungenauere Rechnung mit der korrekten Lösung nichts mehr zu tun (weil auch weitere Testrechnungen mit anderen ode-Functions und anderen Programmen durchgeführt wurden, darf man der Aussage vertrauen, dass die blauen Kurven die korrekte Lösung in diesem Integrationsbereich darstellen).

Die Abweichungen sind von den gewählten Anfangsbedingungen abhängig, bei einigen Anfangsbedingungen liefen die beiden Kurven schon wesentlich eher auseinander, bei anderen (kleinen) Anfangsausschlägen überdeckten sie sich im gesamten hier betrachteten Bereich.

Fazit: Auch die Vorsichtsmaßnahme mit dem MaxStep-Wert 0,05 in swinger.m genügt nicht, um über einen längeren Zeitbereich brauchbare Lösungen zu erzeugen. Man beachte vor allen Dingen, dass die Rechnung weder abstürzt noch irgendwelche Anzeichen dafür liefert, dass die Ergebnisse nichts mehr mit der Aufgabenstellung zu tun haben. Die Animation in swinger.m würde sogar bis  t = 10⁶ hübsch anzusehen sein (wie übrigens auch eine korrekte Lösung, wenn man sie denn erzeugen könnte).

Vorstellbare Parameter, Zeit in Sekunden

Weil das swinger-Beispiel mit Einheitswerten arbeitet, soll hier noch eine Vorstellung davon gegeben werden, wie groß der Zeitbereich tatsächlich ist, in dem brauchbare Lösungen geliefert werden. Dazu muss nur ein Parameter geändert werden, denn man kann sich für die Einheitswerte der Massen und der Längen vorstellen, dass sie in kg bzw. m zu interpretieren sind. Wenn man dann für g den Wert der Erdbeschleunigung  g = 9,81 m/s² annimmt, ergibt sich die Zeit t in s.

Man muss also nur die Zeile 32 in der Function swinger1 in

g = 9.81 ;

ändern, um ein Doppelpendel mit folgenden Parametern zu berechnen:

m₁ = m₂ = 1 kg ; l₁ = l₂ = 1 m ; g = 9,81 m/s² .

Rechts sieht man das Ergebnis dieser Rechnung: Schon nach etwa 10 Sekunden laufen die beiden Kurven auseinander, was nicht einmal heißt, dass die genauere Rechung richtig ist. Tatsächlich sind ab etwa t = 12 s beide Lösungen falsch.

Physikalisches Doppelpendel

Die ausführliche Analyse eines physikalischen Doppelpendels (zwei starre Körper, keine Punktmassen, Differenzialgleichungen gleichen Typs, nur etwas komplizierter) findet man in den Kapiteln "Prinzipien der Mechanik" und "Verifizieren von Computerrechnungen". Die Berechnung mit unterschiedlichen Verfahren und unterschiedlichen Programmen findet man über die Seite "Doppelpendel mit großen Ausschlägen".

Auf der Seite "Doppelpendel - Grenzen der Simulation" werden für dieses Problem die zeitlichen Grenzen ausgelotet, für die eine Simulation sinnvoll ist. Dort wird auch gezeigt, warum ein "Dauerläufer", wie er nebenstehend zu sehen ist, ganz gewiss nicht auf einer seriösen Berechnung beruhen kann.