Doppelpendel - Matlab (mit Kontrollen)
|
Achtung, hier wird geschummelt! | |||||
|
Ein Doppelpendel wird definiert durch die beiden Pendelmassen m1 und m2, die auf die jeweiligen Schwerpunkte bezogenen Massenträgheitsmomente JS1 und JS2, die Schwerpunktabstände von den Drehpunkten s1 und s2 und den Abstand l1 der beiden Drehpunkte voneinander.
Die Bewegung soll durch die Funktionen φ1(t) und φ2(t) beschrieben werden, die für das Zeitintervall t = 0 ... 10 s zu berechnen sind.
Die Aufgabe wird in den Kapiteln "Prinzipien der Mechanik" und "Verifizieren von Computerrechnungen" behandelt.
Die gegebenen Werte gelten für zwei schlanke Stäbe gleicher Masse und gleicher Länge. Sie sollen aus der nebenstehend skizzierten Anfangslage ohne Anfangsgeschwindigkeiten freigelassen werden, so dass folgende Anfangsbedingungen gelten:
Unter Verwendung der in der Aufgabenstellung skizzierten Koordinaten gelten folgende Bewegungs-Differenzialgleichungen (die ausführlich kommentierte Herleitung findet man im Kapitel "Prinzipien der Mechanik"):
Das Anfangswertproblem soll mit zwei verschiedenen Verfahren gelöst werden: Matlab-Standard-Solver ode45 (Dormand-Prince-Verfahren mit variabler Schrittweite) und klassisches Runge-Kutta-Verfahren 4. Ordnung mit konstanter Schrittweite (Function rk4, die mit dem im Kapitel "Kinetik des Massenpunktes" gelisteten Formelsatz arbeitet).
Beide Verfahren erwarten ein Differenzialgleichungssystem 1. Ordnung. Deshalb werden die beiden neuen Variablen ω1 und ω2 eingeführt:
Damit kann das Anfangswertproblem so formuliert werden, wie man es dem Matlab-ode45-Solver anbieten kann:
Die Function rk4 erwartet das nach den ersten Ableitungen aufgelöste System, das mit den oben bereits angegebenen Abkürzungen a11, a12, a22, b1 und b2 zum Beispiel so formuliert werden kann (auch die Anfangsbedingungen gelten wie oben formuliert):
Die Berechnung wird mit dem Standard-Solver von Matlab (ode45) ausgeführt. Dieser steuert die Schrittweite automatisch. Anschließend wird die Runge-Kutta-Rechnung mit der gleichen Anzahl von Schritten ausgeführt, die ode45 berechnet hat. Die Rechnungen werden mit folgenden Strategien kontrolliert:
Weil das Bewegungs-Differenzialgleichungssystem in den Beschleunigungsgliedern gekoppelt ist, entsteht nach dem Umschreiben auf ein System 1. Ordnung auf der linken Seite des Differenzialgleichungssystems eine "Massenmatrix" M. Die Matlab-Solver können damit umgehen, man muss nur über die "options", die dem Solver (hier: ode45) übergeben werden, ankündigen, dass eine gesonderte Function diese Matrix beschreibt (hier wird für diese Function der Name Massenmatrix gewählt).
Die beiden Functions, mit denen das Differenzialgleichungssystem für den Solver ode45 definiert wird (hier willkürlich Massenmatrix und RechteSeite genannt) empfangen die gleichen Werte: Zeit t und Vektor der Funktionswerte mit φ1(t), φ2(t), ω1(t), ω2(t), so dass gemeinsam mit den gegebenen Parametern der Aufgabenstellung alle Informationen verfügbar sind, um die aktuelle Massenmatrix und die aktuelle rechte Seite des Differenzialgleichungssystems zu berechnen und abzuliefern.
Für Bewegungs-Differenzialgleichungen sind die Standard-Einstellungen in Matlab in der Regel nicht ausreichend. Deshalb sollte immer die "MaxStep"-Option (Voreinstellung: 1) verwendet werden, um die Schrittweite zu verkleinern. Hier wird für eine erste Rechnung der Wert 0.1 eingestellt.
Der in der Datei rk4.m realisierte Runge-Kutta-Algorithmus erwartet nur eine Function, die ihm das Differenzialgleichungssystem beschreibt. Hier wurde dafür die Funtion DopendelDgln eingefügt.
Das Matlab-Script liefert die in der nachstehenden Graphik zu sehenden Bewegungsgesetze φ1(t) und φ2(t). Nach relativ kurzer Zeit gehen die Kurven, die mit den beiden Verfahren berechnet wurden, stark auseinander. Das ist natürlich ein sicheres Indiz dafür, dass mindestens eine der beiden Kurven falsch ist.
Die nachfolgend beschriebenen genaueren Rechnungen zeigen, dass beide falsch sind.
Rechnung mit MaxStep-Parameter 0.1:
Die Ergebnisse können nur verbessert werden, wenn man über den MaxStep-Parameter feinere Schrittweiten erzwingt. Es wurden drei Rechnungen mit den MaxStep-Werten 0.1, 0.01 und 0.001 durchgeführt. Im Command Window oben rechts sieht man die Protokolle dieser Rechnungen. Es wurden Ergebnisse für 512, 4016 bzw. 40012 Zeitschritte berechnet. Oben links und nachfolgend sieht man die graphischen Darstellungen der Bewegungsgesetze und der Kontrollfunktion.
Die zur Kontrolle ausgegebenen Werte für φ1 und φ2 am Ende des Integrationsintervalls sind ein sehr schönes Indiz dafür, ob mit ausreichend feiner Schrittweite gerechnet wurde. Bei der Rechnung mit der feinsten Schrittweite stimmen alle Werte für beide Verfahren im Rahmen der ausgegebenen Stellen überein, auch die Maxima und Minima der Kontrollfunktion entsprechen in diesem Rahmen genau den Erwartungen.
In der graphischen Darstellung der Kontrollfunktion Tges(t) sieht man allerdings in jedem Fall noch Zacken, was natürlich auf die automatische Maßstabseinstellung zurückzuführen ist. Wenn man zum Beispiel auf der vertikalen Achse einen Bereich von -1 bis +1 einstellt, dann zeigt sich die Kontrollfunktion als horizontale Gerade mit dem konstanten Funktionswert -0.5.
Rechnung mit MaxStep-Parameter 0.01: |
Rechnung mit MaxStep-Parameter 0.001: |