Advanced Games Physics
Inhaltsverzeichnis zum
Inhalt

Lösung von nicht elementar lösbaren Aufgaben mittels numerischer Methoden

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 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 ΔWrel korreliert mit dem Wert (Δt) p. Mit der Ordnung steigt aber auch der Rechenaufwand!


Vergleiche mit Hilfe des Beispielprogramms die relativen Energiedifferenzen ΔWrel in Abhängigkeit von der gewählten Schrittweite dt. Der Zusammenhang ist linear wegen p = 1 ()!
run program
Fehlerentwicklung beim EC-Verfahren (Ordnung 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)
Formel
(b)
Formel

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
()
Formel
in der zweiten Zeile der Mittelwert aus altem f(si,vi) und neuen Funktionswert f(si+1,vi+1) zur Berechnung der neuen Geschwindigkeit vi+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:
()
Formel


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 vi+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 ():
()
Formel
wobei der Ausdruck vi* 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 vi+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 vi+1 durch die Zwischengeschwindigkeit vi*. 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 ():
()
Formel
... 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 vi+1. Daher auch die gleiche Kur zur Lösung des Problems:
()
Formel
... 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.
download processing
download p5.js
run program
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.
Fehlerentwicklung über der Schrittweite

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 hat zwei Paare von Parametern: eine Radial­komponente (r Radius, v Radial­geschwindig­keit) und eine Tangen­tialkompo­nente (φ Winkel, ω Winkelgeschwindigkeit). Für beide Variablen­paare gibt es je eine Differential­gleichung a) und b):
(a)
Formel
(b)
Formel
Kometenbahn

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 Differential­gleichungen. Das macht die numerische Lösung der Differential­gleichungen etwas schwieriger, als die Lösung einfacher DGln. Versuchen wir es also zunächst mit dem symplektischen EC-Algorithmus.
(a)
Formel
(b)
Formel
(c)
Formel
(d)
Formel

Beim EC-Algorithmus sieht die Lösung noch sehr überschaubar aus. Jeder neue Geschwindig­keits­wert vi+1 oder ωi+1 wird aus einem Satz alter Werte berechnet. Die Ortswerte ri+1 bzw. φi+1 werden symplektisch aus den neuen Geschwindig­keits­werten berechnet.

Eine bessere, d.h. genauere, Lösungsvariante ergibt sich durch die Anwendung des VERLET-Verfahrens. Mit den Gradientenfunktionen f1 und f2 kann das VERLET-Verfahren auf beide Differentialgleichungen separat angewendet werden:
(a)
Formel
(b)
Formel
Geistesblitz
on/off

Zunächst auf die Radial-Komponente...
(a)
Formel
... und dann auf die Tangential-Komponente:
(b)
Formel
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 vi+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ösungs­genauigkeit bedeutet dies natürlich auch eine deutliche Erhöhung der Schrittweite Δt. Im Beispiel­programm kann die Schrittweite wieder per Button gewählt werden. Beide Lösungen nach (symplektischer EC) und (velocity VERLET) wurden implementiert.
download processing
download p5.js
run program
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 Wkin und der aktuellen potentiellen Energie Wpot, 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 Spiele­entwicklung 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.!).