13. Kapitel
Numerische Lösung von Differentialgleichungen 2. Ordnung
Einführung in die numerischen Lösungsmethoden für Differentialgleichungen 2. Ordnung
Das Verfahren von EULER-CAUCHY
Die Spiele-Physik beruht in den meisten Fällen auf der Lösung von Differentialgleichungen. Glücklicher Weise sind dies Differentialgleichungen maximal 2. Ordnung. D.h. es kommen höchstens 2te Ableitungen vor. Dies hat seinen Grund im 2. NEWTONschen Axioms - Kraft ist gleich Masse mal Beschleunigung! Darum haben wir es prinzipiell mit diesem Gleichungstyp zu tun:
()
Differentialgleichungen mit konstanten Koeffizienten werden gewöhnliche>
Differentialgleichungen - im Englischen ordinary differential equation (ODE)
- genannt. Numerische Lösungsmethoden für diesen Typ Differentialgleichung werden
deshalb auch als ODE-Solver bezeichnet.
Dem Kräftegleichgewicht entsprechend sagt aus, dass die Beschleunigung (zweite Ableitung der Ortskoordinate nach der Zeit) abhängig ist vom Ort (z.B. Dehnung von Federn), von der Geschwindigkeit (erste Ableitung der Ortskoordinate nach der Zeit - z.B. Strömungswiderstand) und von der Zeit (z.B. bei periodischer Anregung eines Feder-Masse-Systems). Da aber solche Zeitabhängigkeiten meist nicht analytisch vorliegen (z.B. Interaktion des Nutzers) und zumeist sehr langsam verlaufen, legen wir hier fest, dass für unsere Belange alle Parameter der zu lösenden Differentialgleichung als zeitinvariant angesehen werden können. Damit vereinfacht sich die Differentialgleichung zu :
()
Wollen wir nun diese Differentialgleichung lösen, so ist es zweckmäßig, dies in zwei aufeinander folgenden Schritten zu tun (so haben wir das auch in den ersten Kapiteln "Beschleunigte translatorische Bewegung" und "Reibung" getan). Im ersten Schritt erfolgt die Reduktion der Ordnung der DGl. von 2. Ordnung auf 1. Ordnung, indem die erste Ableitung des Ortes ds/dt durch die Geschwindigkeit v substituiert wird. Dadurch erhalten wir nunmehr zwei DGln. erster Ordnung:
()
Mit den Anfangsbedingungen ()
()
()
()
()
Und wie verhält sich diese Lösungsmethode in Bezug auf die Interaktionsfähigkeit? und machen den Unterschied zu einer analytischen Lösung einer Differentialgleichung deutlich. Während die analytische Lösung von einmal gesetzten Anfangsbedingungen ausgeht und den Bewegungsablauf für alle Zeiten im Voraus beschreibt, geht die numerische Lösung schrittweise vor: auf der Grundlage eines jeden vorangegangenen Rechenschrittes wird ein neues Ergebnis berechnet. Dadurch können während des Übergangs von einem zum nächsten Schritt Parameter geändert werden, ohne dass die in der Vergangenheit berechneten Werte ihre Gültigkeit verlieren. So wie das in experimentell erprobt werden kann.
zeigt in einem Programmauszug, welche Unterschiede zwischen der analytischen und der numerischen Lösungsmethoden im Code bestehen.
Abb. Programmauszug: analytische versus numerische Lösungsmethoden
Genauigkeit numerischer Lösungen
zeigt am Beispiel der Geschwindigkeit, wie die Berechnung
eines neuen Funktionswertes erfolgt. Hier wird aber auch der Mangel des Verfahrens deutlich.
Der Schätzwert
weicht vom wirklichen Wert v, den wir allerdings nicht kennen, ab. Die
Berechnung des nächsten geschätzten Funktionswertes
,
geht zwangsläufig von dem fehlerbehafteten Wert des ersten Schrittes
aus und setzt zu dem bereits vorhanden Fehler noch einen weiteren Fehler hinzu, der
aus dem folgenden 2. Schritt her rührt. Auf diese Weise erfolgt eine Fehler-Fortpflanzung über alle
vorangegangenen Fehler, so dass letztendlich eine ziemliche Abweichung von der
Realität aufgebaut wird ().
Wenn wir über die Berechnungsgenauigkeit von Bewegungskennwerten sprechen, müssen wir zwischen Rechengenauigkeit und Verfahrensgenauigkeit unterscheiden. Die Rechengenauigkeit wird durch die Hardware des Rechners bestimmt. Wieviele Bytes werden für die Zahlendarstellung bereit gestellt. Bei in processing geschriebenen Programmen ist nach Möglichkeit die Genauigkeitsstufe double zu wählen, um eine ausreichende Rechengenauigkeit zu erreichen. In p5.js ist die Genauigkeit von Variablen nicht wählbar, aber mit einer Auflösung von 253 ausreichend.
Wenn wir über die Berechnungsgenauigkeit von Bewegungskennwerten sprechen, müssen wir zwischen Rechengenauigkeit und Verfahrensgenauigkeit unterscheiden. Die Rechengenauigkeit wird durch die Hardware des Rechners bestimmt. Wieviele Bytes werden für die Zahlendarstellung bereit gestellt. Bei in processing geschriebenen Programmen ist nach Möglichkeit die Genauigkeitsstufe double zu wählen, um eine ausreichende Rechengenauigkeit zu erreichen. In p5.js ist die Genauigkeit von Variablen nicht wählbar, aber mit einer Auflösung von 253 ausreichend.
Abb. Fehlerfortpflanzung
Modellversuch
Welche Auswirkungen die Fehlerfortpflanzung hat, kann an verfolgt werden. Ein hängendes Feder-Masse-System, wie schon im Kapitel "Feder-Masse unter Gravitationseinfluss" behandelt, wird durch die Differentialgleichung () beschrieben:
()
()
Nach Betätigung des Startbuttons fällt die Masse zunächst aus ihrer Ruhelage, die durch die Ruhefederlänge l0 bestimmt ist. Damit beginnt das Feder-Masse-System zu schwingen.
Bitte einen Augenblick Geduld
während das Programm geladen wird!
während das Programm geladen wird!
Abb. Fehlerfortpflanzung am Beispiel einer Feder-Masse-Anordnung
Der für die Dämpfung zuständige Regler ist zunächst auf den Wert δ = 0 eingestellt. Das bedeutet, dass die Anordnung verlustfrei - also konservativ - arbeitet, die Schwingung also ewig und unveränderlich fortwährt. Und was beobachten wir? Die Amplitude der Schwingung wird immer größer, bis schließlich ihr Wert über alle Grenzen steigt. Das Verhalten ist instabil und widerspricht allen Erwartungen. Was passiert? Die fehlerhafte Integration akkumuliert offenbar die Einzelfehler, die bei jedem Schritt auftreten. Im Effekt scheint dem System Energie wie von Geisterhand zugeführt worden zu sein. Und noch etwas Unangenehmes ist zu beobachten: je höher die Resonanzfrequenz ist, um so stärker tritt die Instabilität in Erscheinung! Erst wenn der Dämpfungsregler auf eine höhere Dämpfung eingestellt wird, stabilisiert sich die Schwingung, weil so die fehlerhaft "zugeführte" Energie in Reibungswärme umgesetzt wird.
Um eine stabile Schwingung zu erreichen, muss also die Dämpfung erhöht werden. Das ist aber keine saubere Lösung, weil die Abhängigkeit der Stabilität z.B. von der Resonanzfrequenz ω0 stets andere Einstellungen für die Dämpfung δ erfordert.
Noch eine weiterer Effekt kann beobachtet werden: Die Schwingungsperiode der nach dem EC-Verfahren berechneten Zeitfunktion stimmt nicht mit der tatsächlichen (grauer Graf) überein. Dieser Fehler ist aber im Vergleich zum instabilen Verhalten der Funktion fast nicht zu beobachten.
Einfluss der Schrittweite auf die Genauigkeit
Die Verfahrensgenauigkeit hingegen wird einzig (noch wissen wir es nicht besser) durch die Wahl der Schrittweite bestimmt. Nun könntest Du denken, je kleiner die Schrittweite ist, desto genauer die Rechnung. Aber das stimmt nicht. Abgesehen davon, dass dabei die Rechenleistung enorm in die Höhe getrieben wird (allein eine Halbierung der Schrittweite für zu einer Verdopplung des Rechenaufwandes!), gibt es für die Schrittweite ein Optimum. Jede Multiplikation führt zunächst zu einer Verdopplung der Wortlänge. Damit das Produkt wieder in die normale Wordgröße kommt, muss gerundet werden. Wird nun die Schrittweite sehr klein, wird auch das Produkt sehr klein. Und nun wird gerundet. Dabei fallen die für die Genauigkeit wichtigen Stellen unter den Tisch. Damit ist das ganze Gegenteil vom beabsichtigten eingetreten.Gewisse Funktionen (z.B. eine Kometenbahn in der Nähe der Sonne) zeigen Phasen großer Änderungsgeschwindigkeit ihrer Kenngrößen (Ort und Geschwindigkeit), in denen sinnvoller Weise mit kleinerer Schrittweite gearbeitet werden sollte, um die Feinheiten der Bahn besser abbilden zu können. Andrerseits gibt es aber auch Phasen, in denen relativ wenig passiert, also könnte die Schrittweite ohne Schaden größer gewählt werden, um Rechenleistung zu sparen. In diese Richtung stoßen wir allerdings an eine Grenze, die maximale Schrittweite sollte durch die Bildwechselfrequenz frameRate vorgegeben werden.
Und wie kann ich trotzdem die Schrittweite in diesem vorgegebenen Rahmen verändern? Nun ganz einfach: das durch die Bildwechselfrequenz vorgegebne Zeitraster wird durch einen Integerwert geteilt und so verkleinert. Damit aber der Zeitmaßstab unverändert bleibt, muss nun die Berechnung genau so oft wie das der Integerwert vorgibt, wiederholt werden (siehe : push-Button "Schrittweite"). Dies ist notwendig, auch wenn die so gewonnenen Zwischenwerte gar nicht dargestellt werden können, liegen sie doch innerhalb der Zeitspanne des Bildwechsels. Diese Zwischenwerte werden dennoch nicht umsonst berechnet! Sie liefern ja stets die Ausgangsbasis für die nächste Iteration.
Das alles legt die Idee einer adaptiven Schrittweitensteuerung nahe. Wir werden im Kapitel Numerische Methoden/Dormand-Prince-Verfahren genauer auf dieses Thema eingehen.
Nun ist es allerdings so, dass in der Spielephysik die Genauigkeit gar nicht an erster Stelle steht. Viel wichtiger ist, dass die Bewegungsabläufe stabil und erwartungsgemäß sind. Darum müssen wir die Frage beantworten, gibt es andere Wege oder Verfahren die Stabilität zu erhöhen. Damit sollen sich die folgenden Abschnitte befassen.
