Anfangswertproblem, Lösung mit MATLAB

Die skizzierte Masse m ist durch eine Feder gefesselt und kann sich auf der horizontalen Führung reibungsfrei bewegen. Sie ist durch eine (lineare) Feder gefesselt, die im entspannten Zustand die Länge b hat. Die Masse wird um x_anf ausgelenkt und zum Zeitpunkt t = 0 ohne Anfangsgeschwindigkeit freigelassen. Für das Zeitintervall 0 ≤ t ≤ 30 s sollen das Weg-Zeit-Gesetz und das Geschwindigkeits-Zeit-Gesetz berechnet und graphisch dargestellt werden. Ort und Geschwindigkeit der Masse am Ende des Zeitintervalls sind anzugeben.

Geg.: m = 10 kg ;    c = 5 N/m ;    a = 0,4 m ;    b = 0,8 m .

Anfangsbedingungen:      x_anf,1 = - 2000 mm   ,   x_anf,2 = - 1000 mm   ,   x_anf,3 = - 1131,5 mm.

• Man berechne zunächst die Lösungen unter Verwendung eines Runge-Kutta-Verfahrens 4. Ordnung mit konstanter Schrittweite und verändere die Schrittweite, bis man sicher ist, dass die Ergebnisse korrekt sind:

xanf,1 = - 2 m Zeitschritte n 100 200 500 1000   xend [m] 1,5226 1,5233 1,5233 1,5233   vend [m/s] 0,68439 0,68413 0,68412 0,68412

xanf,1 = - 1 m Zeitschritte n 100 200 500 1000   xend [m] - 0,65790 - 0,65793 - 0,65793     vend [m/s] - 0,19473 - 0,19474 - 0,19474

xanf,1 = - 1,1315 m Zeitschritte n 100 200 500 1000 2000 xend [m] - 0,0014782 - 0,0019625 - 0,0019809 - 0,0019811 - 0,0019811 vend [m/s] - 0,0068817 - 0,0071102 - 0,0071186 - 0,0071187 - 0,0071187

• Danach ermittle man die Lösungen unter Verwendung einer “MATLAB-ode”-Funktion (Empfehlung: ode45), die mit automatischer Schrittweitensteuerung arbeitet, vergleiche die Ergebnisse mit denen der Runge-Kutta-Rechnung und entscheide, ob die “MaxStep-Option” gesetzt werden muss. Die Anzahl der mit der MATLAB-Funktion berechneten Zeitpunkte ist als zusätzliches Ergebnis auszugeben.

xanf,1 = - 2 m MaxStep Ohne Option 0,5 0,2 0,1   Zeitschritte n 156 276 636 1232   xend [m] 1,5182 1,5236 1,5233 1,5233   vend [m/s] 0,68673 0,68405 0,68412 0,68412

xanf,1 = - 1 m MaxStep Ohne Option 0,5 0,2 0,1   Zeitschritte n 120 272 632 1232   xend [m] - 0,66005 - 0,65794 - 0,65793 - 0,65793   vend [m/s] -0,19468 - 0,19474 - 0,19474 - 0,19474

xanf,1 = - 1,1315 m MaxStep Ohne Option 0,5 0,2 0,1 0,05 Zeitschritte n 132 272 632 1228 2428 xend [m] 0,0067547 - 0,0016937 - 0,0019789 - 0,0019811 - 0,0019811 vend [m/s] 0,0078380 - 0,0069873 - 0,0071176 - 0,0071187 - 0,0071187

• Weil keine Energieverluste (durch Reibung oder andere Bewegungswiderstände) während der Bewegung auftreten, muss die Gesamtenergie

(Summe aus der kinetischen Energie und der jeweils in der Feder gespeicherten Energie) konstant sein. Man gebe zusätzlich diese Funktion als Kontrolle in ein Graphikfenster aus.

Die folgenden Rechnungen zeigen in den beiden oberen Fenstern die Bewegungsgesetze nach beiden Verfahren, im unteren Fenster den Verlauf von T_ges (wenn nur eine Kurve zu sehen ist, dann sind die Ergebnisse nach beiden Verfahren annähernd identisch). Dabei wurden stets etwa gleiche Anzahl von Zeitschritten für beide Verfahren gewählt (für die ode45 -Rechnung ergibt sich die Schrittanzahl durch die automatische Schrittweitensteuerung, die mit der MaxStep-Option beeinflusst werden kann). Allerdings ist die Anzahl der mit ode45 ausgeführten Zeitschritte nur ein Viertel der abgelieferten Werte, die zusätzlichen Werte werden durch Interpolation ermittelt.

Rechnung ohne MaxStep-Option mit   x_anf = - 2 m:

Rechnung mit MaxStep-Option 0,2 mit  x_anf = - 1 m:

Rechnung mit MaxStep-Option 0,1 mit  x_anf = - 1,1315 mm:

Fazit der Rechnugen: Die Verläufe dürfen bei geeigneter Wahl der Schrittweite bzw. der MaxStep-Option im betrachteten Intervall als ausreichend genau angesehen werden. Dass der ode45-Algorithmus stets etwas ungenauere Ergebnisse liefert als der Runge-Kutta-Algorithmus 4. Ordnung, liegt an oben bereits beschriebenen Tatsache, dass die Anzahl der tatsächlich berechneten Integrationsschritte nur ein Viertel der ausgewiesenen Ergebnisse beträgt.

Diese Aussage wird auch durch die T_ges-Kurven- verdeutlicht, auch wenn man beachten muss, dass z. B. bei der letzten Rechnung die mit 5 Stellen angezeigten T_ges-Werte auf der vertikalen Skala sich nicht ändern (tatsächlich ergibt sich eine Änderung erst in der 9. Stelle).

Die oben beschriebenen Berechnungen wurden mit dem Script FedMassHoriz.m ausgeführt.

Ein Script, das die Bewegung des Feder-Masse-Systems verdeutlicht (nebenstehendes Bild), steht als FedMassHorizAni.m zur Verfügung. Dieses benötigt zusätzlich die beiden Function Files SpringCoords.m und MassCoords.m .

Lösung der Zusatzaufgabe: Man kann mit einer Energiebetrachtung die Anfangsbedingung exakt ermitteln, bei der theoretisch der Stillstand bei  x = 0 eintreten müsste. In beiden Zuständen ist nur Energie in der Feder gespeichert. Gleichsetzen der Energien entsprechend

liefert neben der trivialen Lösung x₀ = 0 zwei weitere Lösungen:

Mit den gegebenen Werten für a und b erhält man also z. B.:  x₀ = - 1,13137084989848 m .

Auch wenn man die Rechnung mit dem Anfangswert xanf = -2*b*sqrt(1-a/b) startet, bleibt die Masse selbst bei feinster Einteilung des Integrationsintervalls bzw. sehr kleinem MaxStep -Wert nicht länger als etwa bis t = 50 s in Ruhe, was nicht den Integrationsverfahren anzulasten ist, sondern der internen Darstellung der Zahlen mit begrenzter Länge.

Die nebenstehende Graphik zeigt die Lösung mit 24000 Runge-Kutta-Schritten und MaxStep = 0.01 für die ode45-Lösung..