13. Kapitel

Lösungsverfahren höherer Genauigkeit

Wir wollen hier nicht die ganze Theorie der numerischen Integrationsmethoden diskutieren, dafür gibt es hinreichend und gute Literatur, wie z.B.

J. Dankert EinfNumIntAWP.pdf
J. Werner ode.pdf

Viel mehr wollen wir uns um die Anwendungsseite kümmern. Also Fragen wie

Stabilität der Lösung,
erzielbare Genauigkeit,
Rechenaufwand und
Komplexität und Implementierbarkeit der Lösung

sollen hier besprochen werden.

Definitionen

Um ein Verfahren nach den genannten Gesichtspunkten beurteilen zu können, bedarf es einiger Begriffe und das Wissen um ihre Bedeutung:

Schritt: ein Schritt (eigentlich Integrationsschritt) bedeutet die Berechnung eines Wertes i+1 auf dem Wert i basierend. Wird dabei nur auf den Wert i zurück gegriffen, handelt es sich um ein Einschrittverfahren.

Alle EULER-CAUCHY-Verfahren sind Einschrittverfahren. Die Anzahl der Schritte nimmt Einfluss auf die Genauigkeit und auf den Rechenaufwand.

Stufe: ein Verfahren wird s-stufig bezeichnet, wenn innerhalb eines Schrittes der Wert der Gradientenfunktion f(x,y) s-mal berechnet werden muss. EULER-CAUCHY-Verfahren sind einstufig.

Das einfache EULER-CAUCHY-Verfahren ist ein 1-stufiges Verfahren. Anzahl der Stufen nimmt starken Einfluss auf die Genauigkeit und auf den Rechenaufwand. Der Implementierungsaufwand steigt nur mäßig.

Ordnung: ein Verfahren ist von der Ordnung (eigentlich Konsistenz-Ordnung) p, wenn der Restfehler kleiner als (Δt)^p ist.

Das einfache EULER-CAUCHY-Verfahren ist von der Ordnung p = 1. Verfahren niedriger Ordnung sind nicht geeignet für Lösungen, die über längere Zeiten hohe Genauigkeit verlangen. Mit einer verkleinerten Schrittweite allein ist das nicht zu schaffen!

Die Ordnung ist ein direkter Ausdruck der Genauigkeit. Unsere Genauigkeitskriterium ΔW_rel korreliert mit dem Wert (Δt)^p. Mit der Ordnung steigt aber auch der Rechenaufwand!

Vergleiche mit Hilfe des Beispielprogramms die relativen Energiedifferenzen ΔW_rel in Abhängigkeit von der gewählten Schrittweite dt. Der Zusammenhang ist linear wegen p = 1 ()!

Abb. Fehlerentwicklung beim EC-Verfahren (Ordnung p = 1)

Leapfrog- und VERLET-Verfahren

Die einfachsten Lösungsverfahren höherer Genauigkeit sind das Leapfrog- und das VERLET-Verfahren. Beide wurden ursprünglich für konservative Systeme, d.h. also reibungsfreie Anordnungen, entwickelt. Sie gehen auf die selbe Grundidee zurück: durch lineare Interpolation von Zwischenwerten von Ort s sowie Geschwindigkeit v und Anwendung des Energie erhaltenden symplektischen EULER-CAUCHY-Verfahrens werden Ort und Geschwindigkeit stets zu unterschiedlichen Zeitpunkten berechnet. Der Zwischenschritt wird erzeugt, indem das Zeitintervall Δt halbiert wird. Zu je einem halbierten Intervall i -> i+1/2 bzw. i+1/2 -> i+1 der Dauer Δt/2 erfolgt je eine symplektische Berechnung von Ort und Geschwindigkeit ( a und b):

(a)

(b)

Zunächst sieht es so aus, als würde nur die Schrittweite verringert und dadurch ein Genauigkeitsvorteil gewonnen. Das ist aber nicht die ganze Wahrheit. Werden nämlich a) und b) in einander überführt, zeigt sich, dass

()

in der zweiten Zeile der Mittelwert aus altem f(s_i,v_i) und neuen Funktionswert f(s_i+1,v_i+1) zur Berechnung der neuen Geschwindigkeit v_i+1 herangezogen werden und dass die Ortsberechnung in der ersten Zeile mit der Ordnung p = 2 erfolgt, wodurch die Genauigkeit um eine 10er-Potenz erhöht wird.
Wenn wir nun noch etwas aufräumen, erhalten wir mit den VERLET-Algorithmus, um genau zu sein, den velocity VERLET-Algorithmus:

()

Ich erwähnte schon, dass beide Verfahren für Systeme ohne Reibungseinfluss wie z.B. die Planetenbewegung oder die Molekularbewegung entwickelt worden sind. Dort sind die Differentialgleichungen nicht von der Geschwindigkeit v abhängig und f(s,v) wird zu f(s).

Warum erwähne ich das? verlangt nämlich bei der Geschwindigkeitsberechnung den Wert v_i+1 auf beiden Seiten der Gleichung. Damit hätten wir hier eine implizite Gleichung vorliegen, die oftmals nicht elementar lösbar ist. Tun wir also zunächst so, als würden wir hier ohne Reibungseinfluss arbeiten. Wie wir die Reibung doch noch unter bekommen, besprechen wir anhand der Implementierung beider Verfahren.

Das Leapfrog-Verfahren nach a) und b) kann, weil die 2. und die 3. Formel in einen Schritt zusammengefasst werden können, vereinfacht auch so geschrieben werden ():

()

wobei der Ausdruck v_i* nur ein Zwischenergebnis darstellt, das in der weiteren Berechnung mehrfach benötigt wird.

Wenden wir jetzt auf das schon bekannte Problem des hängend angeordneten Feder-Masse-Systems an. Die ersten beiden Zeilen bereiten keine Probleme, aber in der dritten Zeile tritt nun die implizite Abhängigkeit der Geschwindigkeit v_i+1 auf. Da uns dieser Wert aktuell nicht zur Verfügung steht, dieser aber "nur" auf die Reibungswirkung Einfluss nimmt, erstzen wir hier die Geschwindigkeit v_i+1 durch die Zwischengeschwindigkeit v_i*. Das ist nicht korrekt, aber zu verschmerzen. Wissen wir doch aus den vorangegangenen Experimenten (explizites EC-Verfahren), dass ein fehlerhaftes Zeitverhalten durch mehr oder weniger Reibungseinfluss kompensiert werden kann. Das werden wir nun auch tun ():

()

... und hier das zugehörige Listing:

Abb. Listing des Leapfrog-Algorithmus

Die Implementierung des VERLET-Algorithmus nach stößt auf das selbe Problem der implizit auftretenden Geschwindigkeit v_i+1. Daher auch die gleiche Kur zur Lösung des Problems:

()

... auch hier das zugehörige Listing:

Abb. Listing des VERLET-Algorithmus

... und das Beispielprogramm, das beide Verfahren zeigt:

Das Beispielprogramm zeigt das Verhalten eines schwingenden Systems, wobei zwei Lösungsmethoden - VERLET und leapfrog - gegenüber gestellt werden. Mit den Schiebereglern können wieder Eigenfrequenz und Dämpfung des Systems eingestellt werden. Mit dem Button Width kann die Schrittweite der Integration gewählt werden. Zum Vergleich ist auch das implizite EULER-CAUCHY in die Betrachtung einbezogen.

Zur Kontrolle der Lösungsgenauigkeit werden sowohl die Abweichungen der Gesamtenergien wie auch der Resonanzfrequenzen im Vergleich zur theoretisch vorhergesagten angezeigt.

Achtung! die Mittlungsdauer für die Fehlerermittlung ist relativ groß, so dass Du etwa 10 s warten musst, bevor gültige Ergebnisse angezeigt werden!

Fazit: Sowohl das leapfrog- als auch das VERLET-Verfahren sind einstufige Einschritt-Verfahren der 2. Konsistenzordnung. Wie zeigt, ist der energetische Fehler, den beide Verfahren aufweisen, um eine 10er-Potenz kleiner als bei den einfachen EC-Verfahren. Zudem weist die quadratische Abhängigkeit des Energiefehlers auf die Ordnung p = 2 hin.

Abb. Fehlerentwicklung über der Schrittweite

Warum höhere Genauigkeit?

Warum eigentlich benötigen wir Lösungsverfahren mit einer höheren Genauigkeit? Unsere bisherigen Erfahrungen, ob bei translatorischen Aufgaben oder bei der Rotation, mit dem einfachen Symplektische Lösungsverfahren waren doch alle hinreichend gut. Hinreichend natürlich im Sinne der Spielephilosophie! So lange es sich um einfache Differentialgleichungen 2. Ordnung handelt, müsste hinzugefügt werden. Die eigentliche Herausforderung bilden nämlich die verkoppelten Differentialgleichungen, die teilweise sogar nichtlinear sind. Verkoppelte Differentialgleichungen beschreiben zusammegesetzte Systeme, bei denen eine wechselseitige Abhängigkeit der Variablen der beteiligten Einzelsysteme bestehen. Dabei gibt es Fälle von extremen Abhängigkeiten, den sog. steifen Differentialgleichungen. Der Begriff "steif" ist nirgendwo exakt definiert und rührt offenbar von Feder-Masse-Systemen mit einer steifen Feder her. "Steife" Feder bedeutet, dass die Federkonstante n sehr groß ist, was zu einer hohen Eigenfrequenz des Systems führt. Dies wiederum hat schnelle Wechsel der Geschwindigkeits- bzw. Ortsfunktion zur Folge. Ist dann die Schrittweite des Lösungsverfahrens nicht klein genug, kann es zu großen Abweichungen von der realen Schwingung kommen und das System wird schnell instabil.

Um die Auswirkung mangelnder Genauigkeit auf verkoppelte Systeme zu untersuchen, betrachten wir im folgenden einen Fall ohne Reibungseinfluss. Ein solcher Fall ist die Bewegung eines Planeten oder Kometen um ein Zentralgestirn - in der Tat gänzlich reibungsfrei.

zeigt die Bahn eines Kometen um die Sonne. Wie bereits im Kapitel Bewegungsgleichungen für Planeten und Kometen besprochen, hat die Bahn zwei Paare von Parametern: eine Radialkomponente (r Radius, v Radialgeschwindigkeit) und eine Tangentialkomponente (φ Winkel, ω Winkelgeschwindigkeit). Für beide Variablenpaare gibt es je eine Differentialgleichung a) und b):

(a)

(b)

Abb. Kometenbahn

Wie an a) und b) zu sehen ist, kommen beide Variablen bzw. ihre Ableitungen in beiden Gleichungen vor, es handelt sich also um ein System verkoppelter Differentialgleichungen. Das macht die numerische Lösung der Differentialgleichungen etwas schwieriger, als die Lösung einfacher DGln. Versuchen wir es also zunächst mit dem symplektischen EC-Algorithmus.

(a)

(b)

(c)

(d)

Beim EC-Algorithmus sieht die Lösung noch sehr überschaubar aus. Jeder neue Geschwindigkeitswert v_i+1 oder ω_i+1 wird aus einem Satz alter Werte berechnet. Die Ortswerte r_i+1 bzw. φ_i+1 werden symplektisch aus den neuen Geschwindigkeitswerten berechnet.

Unter Gradientenfunktion verstehe ich die rechte Seite der Normalform einer Differentialgleichung 1. Ordnung, bei der auf der linken Seite die erste Ableitung der gesuchten Variablen nach der Zeit steht, also z.B.:

hat die Gradientenfunktion

Eine bessere, d.h. genauere, Lösungsvariante ergibt sich durch die Anwendung des VERLET-Verfahrens. Mit den Gradientenfunktionen f₁ und f₂ kann das VERLET-Verfahren auf beide Differentialgleichungen separat angewendet werden:

(a)

(b)

on/off

Zunächst auf die Radial-Komponente...

(a)

... und dann auf die Tangential-Komponente:

(b)

Allerdings, wir sehen es sofort, hat diese Implementierung einen Haken! Nämlich, die letzten Zeilen in den Gleichungen a) und b) sind implizite Funktionen der neu zu berechnenden Werte. Jetzt können wir diesen Fakt aber nicht wie in der obigen Erklärung dargelegt vernachlässigen! Jetzt muss eine andere Lösung her.
Diese Lösung besteht nun darin, dass wir mit Hilfe des einfachen EC-Algorithmus aus a) und c) eine Schätzung der Werte v_i+1 und ω_i+1 gewinnen und diese in den rechten Seiten der untersten Gleichungen von a) und b) einsetzen. Das ist zwar nicht das Ideal, aber besser als die Verwendung der alten Funktionswerte (dieses Verfahren wird übrigens bei allen impliziten Lösungsmethoden angewendet, sofern die Gradientengleichungen eine einfache Auflösung nicht erlauben).

Im nebenstehenden Programmbeispiel lassen wir einen Kometen oder Planeten um die Sonne kreisen. Selbstverständlich ist das ein Vorgang der Jahre dauern kann, es sei denn, wir betrachten den Umlauf im Zeitraffer. Im Beispiel wurde die Zeit um einen Faktor timeScale = 500 gestaucht. Für die Lösungsgenauigkeit bedeutet dies natürlich auch eine deutliche Erhöhung der Schrittweite Δt. Im Beispielprogramm kann die Schrittweite wieder per Button gewählt werden. Beide Lösungen nach (symplektischer EC) und (velocity VERLET) wurden implementiert.

Als Qualitätsmaß der Lösungsmethoden findet wieder das energetische Defizit Verwendung. Diesmal in prozentualen Angaben, bezogen auf die Gesamtenergie zu Beginn des Experiments. Die aktuelle Energie errechnet sich aus der Summe der aktuellen kinetischen Energie W_kin und der aktuellen potentiellen Energie W_pot, die gemäß berechnet wird.

Ergebnisdiskussion: Obwohl das eigentlich für energetischen Erhalt stehende symplektische EULER-CAUCHY-Verfahren im ersten Fall angewendet wird, reicht bei einer Schrittweite von Δt ≈ 80 s die Genauigkeit nicht aus, um eine wie auch immer geartete Umlaufbahn zu erzeugen. Selbst bei noch kleineren Schrittweiten sind die erzeugten Umläufe so falsch, dass eine geschlossene Umlaufbahn nicht zu Stande kommt. Das ist auch für die Spieleentwicklung nicht zun tolerieren! Hingegen erreicht der VERLET-Algorithmus dieses Ziel. Fehler die dennoch auftreten, treten in der Nähe des Zentralgestirns besonders deutlich in Erscheinung. Sie sind in jedem Fall deutlich geringer als beim EC-Verfahren.
Die Tatsache des vergrößerten Energiefehlers in der Nähe des Zentralgestirns deutet darauf hin, dass die Grandientenfunktionen dort besonders schnellen Änderungen unterliegen (siehe steife DGl.!).

Symbol	Bedeutung
p	Konsistenz-Ordnung
s	Anzahl der Stufen (nicht Strecke!)
s_i	Strecke, Weg [m], Rechenergebnis im i-ten Zeitquant
v_i	Geschwindigkeit [m/s], Rechenergebnis im i-ten Zeitquant
g	Erdbeschleunigungskonstante [m/s²]
s₀	Startort [m]
v₀	Anfangsgeschwindigkeit [m/s]
m	Masse [kg]
r	Reibzahl bei STOKESscher Reibung [kg/s]
n	Federkonstante [N/m]
l₀	Ruhefederlänge [m]
δ	Dämpfung [1/s] ()
ω₀	Kreisfrequenz der Eigenschwingung [rad/s] ()
φ	Winkel [rad]
ω	Drehgeschwindigkeit [rad/s] ()
r	Abstand, Radius einer Planetenbahn [m]
γ	Gravitationskonstante [m³/kgs²]
t	Zeit [s]
Δt	Zeitquant, Bildwechseldauer [s]