13. Kapitel

Explizites, implizites oder symplektisches Verfahren?

Im folgenden sollen anhand der EULER-CAUCHYschen Lösungsmethode die wesentlichen Unterschiede und Eigenschaften expliziter, impliziter bzw. symplektischer Lösungsansätze erläutert werden.

Explizite Lösungsverfahren

Das im Beispiel behandelte Lösungsverfahren () ist ein explizites Verfahren. Explizit deshalb, weil die Berechnung aller Funktionswerte von jeweils nur einem bekannten Datensatz beginnend mit s₀ und v₀ zur Zeit t₀ über s_i und v_i zur Zeit t_i erfolgt.
Im streng mathematischen Sinn bedeutet das für die Lösung von Differentialgleichungen 2. Ordnung, dass die Berechnung sowohl des Weges s als auch der Geschwindigkeit v mit einem Datensatz, der von ein und demselben Zeitpunkt t stammt, erfolgt. An dem in gezeigten Programmabschnitt kann dies nachvollzogen werden. Da Weg s und Geschwindigkeit v vom selben Zeitpunkt stammen, spielt die Integrationsreihenfolge keine Rolle. Um sicher zu stellen, dass beide Integrationen mit einem vom selben Zeitpunkt stammenden Wertepaar arbeiten, wird der Wert für die Geschwindigkeit v vor der Aktualisierung der Geschwindigkeit zwischengespeichert (v_) und für die folgende Wegberechnung verwendet. Im Programmcode wird anstelle des Ausdrucks Δt der Ausdruck dt verwendet.

Abb. Implementierung des expliziten EC-Lösungsverfahrens

Wie das Beispielprogramm () demonstriert, hat das explizite Lösungsverfahren die Neigung Energie in das System zu pumpen. Es akkumuliert die Prinzip bedingten Fehler und führt so zu instabilen Lösungen.

Explizite Lösungen neigen zur Instabilität und verstoßen gegen die Erhaltungssätze von Energie und Impuls.

Implizite Lösungsverfahren

Eine verbesserte Lösung des Integrals in sollte erreicht werden, wenn weitere Stützstellen der gesuchten Funktion in die Integration einbezogen werden. Nicht nur wie bisher am linken Rand (y_i) sondern zusätzlich auch die Stützstelle y_i+1 am rechten Rand des Integrationsintervalls werden herangezogen ().

Nun taucht natürlich sofort die Frage auf, woher nehmen, wenn nicht stehlen? v_i+1 ist ja noch gar nicht bekannt, das soll ja erst berechnet werden! Und genau das ist das charakteristische an den impliziten Lösungsverfahren!
Für eine Differentialgleichung 2. Ordnung heißt das, dass sowohl der bekannte Ort s_i als auch die bekannte Geschwindigkeit v_i als Ausgangswerte für die Berechnung der nächsten Werte v_i+1 bzw. s_i+1 herangezogen werden (a und b). Aber dann! Die Zuwächse an Geschwindigkeit f(s_i+1,v_i+1)·Δt werden bereits mit den Ergebnissen dieses Schrittes berechnet:

(a)

Abb. implizites EC-Verfahren

Analog wird die Berechnung der neuen Ortsvariablen ausgeführt:

(b)

Statt die Steigung im Punkt v_i zur Grundlage der Berechnung des nächsten Schrittes zu verwenden, wird nun die Steigung im Punkt v_i+1; s_i+1 verwendet. Damit bekommt der Algorithmus eine "Ahnung", wie sich die gesuchte Funktion weiterhin verhalten wird. Damit liegt diese Schätzung (schwarzer Punkt) näher am wirklichen Wert (roter Punkt) als bei der expliziten Lösungsmethode. Daher verhält sich die Abweichung so, dass dem System keine zusätzliche Energie zugeführt wird. Die Lösung ist immer stabil!
Es wird also eine Geradengleichung zur Gewinnung des Schätzwertes

{\tilde{v}}_{i + 1}

angewendet, deren Ausgangspunkt der i-te Wert ist, aber die Steigung aus dem Wert von i+1 bezieht.

Wollen wir uns das mal genauer am Beispiel des aufgehängten Feder-Masse-Systems ansehen. Zunächst wird die 1. Integration ausgeführt, um so auf der Basis der noch nicht bekannten Größen v_i+1 und y_i+1 den neuen Geschwindigkeitswert v_i+1 zu berechnen:

(a)

Der neue Ortswert y_i+1 wird wie beim expliziten EC-Verfahren auf der Grundlage des alten Geschwindigkeitswertes v_i ermittelt:

(b)

Betrachten wir zunächst die Abschätzung der Geschwindigkeit. Wir sehen, dass v_i+1 auf beiden Seiten der a) vorkommt. Glücklicher Weise sogar in einer Form, die die Separierung von v_i+1 leicht macht. Beachten wir noch, dass

y_{i + 1} = y_{i} + v_{i+1} \cdot Δ t

ist (v_i+1 deshalb, weil wir ja des Ergebnis des i+1-ten Schrittes für das implizite Verfahren verwenden wollen). Wir lösen a) nach v_i+1 auf:

()

und bekommen so eine Lösung für den neuen Geschwindigkeitswert v_i+1 in expliziter Form. Allerdings muss darauf hingewiesen werden, dass das gewählte Beispiel nicht besonders gut zur Demonstration der Eigenschaften des impliziten EC-Algorithmus geeignet ist! Denn sehen wir uns a) an, stellen wir fest, dass der Einfluss des noch zu berechnenden Geschwindigkeitswertes v_i+1 stark von der Dämpfung δ abhängig ist. In einem ungedämpften System gibt es diesen Einfluss gar nicht! Nur über den Federort gibt es eine unbedingt implizite Abhängigkeit.

Leider sind nicht alle impliziten Gleichungen in der Praxis so einfach in eine explizite Form zu überführen. In den meisten Fällen sind die impliziten Zusammenhänge nichtlinear (z.B. tritt beim einfachen Pendel ein sinus-Abhängigkeit auf) und die lassen sich nicht so leicht umstellen.

zeigt den wesentlichen Auszug der Implementierung des impliziten EC-Lösungsverfahrens und deren ausführbares Beispiel.

Abb. Implementierung des impliziten EC-Lösungsverfahrens

Das Beispielprogramm zeigt das Verhalten eines schwingenden Systems, wenn das implizite EC-Lösungsverfahren angewendet wird. Mit den Schiebereglern können Eigenfrequenz ω₀ und Dämpfung δ des Systems eingestellt werden. Mit dem Button push kann die Schrittweite der Integration gewählt werden.

Bitte einen Augenblick Geduld
während das Programm geladen wird!

Abb. Implizites Lösungsverfahren

Im Experiment erkennen wir den Unterschied zur expliziten Lösung. Die Schwingungsamplitude bleibt konstant oder verringert sich. Das bedeutet, das System nimmt mit der Zeit keine Energie auf!
Dennoch gibt es ein ABER! Der graue Graf in zeigt eine Vergleichsschwingung mit genau der Resonanzfrequenz ω₀, die auch der impliziten Lösung zugrunde liegt. Es sollte also erwartet werden, dass beide Schwingungen identische Periodendauern aufweisen. Das ist aber nicht der Fall! Erst bei sehr geringen Schrittweiten verlaufen beide Schwingungen deckungsgleich. Sonst ist ein Voraneilen der impliziten Lösung zu beobachten.
So, wie unser Beispiel der hängenden Feder-Masse-Anordnung nahe legt, scheint mit dem impliziten Verfahren, außer einer stabilen Lösung, nicht viel gewonnen zu sein. Die erreichte Genauigkeit lässt auch zu wünschen übrig. Eine Verbesserung der Genauigkeit wird mit dem HEUN-Verfahren erreicht. Das implizite EC-Verfahren ist der Ursprung vieler weiterer Lösungsmethoden mit verbesserter Genauigkeit.

Implizite Lösungen sind stabil und verstoßen trotzdem gegen die Erhaltungssätze von Energie und Impuls.

Symplektische Lösungsverfahren

Im Gegensatz zu den expliziten bzw. impliziten Lösungsmethoden, die auf Differentialgleichungen beliebiger Ordnung anwendbar sind, werden die symplektischen Verfahren nur zur Lösung von Differentialgleichung 2. Ordnung herangezogen. Der Begriff symplektisch stammt aus der Geometrie und bedeutet so viel wie "schief symmetrisch". Das Prinzip beruht darauf, dass die Ortsberechnung zu einem anderen Zeitpunkt als die Geschwindigkeitsberechnung erfolgt (). Die Berechnung des neuen Ortes wird mit dem alten Ort und der neuen Geschwindigkeit berechnet:

()

Dies ist die intuitive Lösungsmethode, die aus der analytischen Lösung von Differentialgleichungen 2. Ordnung folgt. Die Lösung erfolgt ja in zwei Schritten:

Im ersten Schritt wird aus dem Kräftegleichgewicht (der Beschleunigung) und dem Anfangswert für die Geschwindigkeit eine Funktion der Geschwindigkeit und
im zweiten Schritt aus dieser Geschwindigkeitsfunktion mit dem Anfangswert für den Ort eine Funktion des Ortes berechnet.

Und genau so wird bei der symplektischen numerischen Lösung vorgegangen (das ist auch die Lösungsmethode, die wir bisher auch meistens angewendet haben!).

Nehmen wir wieder unser Beispiel des hängenden Feder-Masse-Systems. Die zu lösende Aufgabe lautet jetzt:

()

Kommt Dir der Programmauszug () für die symplektische Lösung bekannt vor? Na, klar! Das ist das Muster, nach dem 99% der aufgetretenen Differentialgleichungen gelöst wurden.

Abb. Implementierung des symplektischen EC-Lösungsverfahrens

Auch in diesem Beispiel wird die Parametertreue überprüft. Was verstehe ich unter Parametertreue? Nun ein schwingungsfähiges System hat eine Eigenschwingung, die durch die Parameter des Systems vorgegeben ist. So z.B. wird die Eigenschwingung ω₀ eines Feder-Masse-Systems durch die Federsteifigkeit n und seine Masse m bestimmt. Das Quadrat der Eigenkreisschwingung ergibt sich dann aus ω₀² = n/m.
Nun ist offensichtlich, dass die erzeugte Schwingung eines zu untersuchenden Lösungsverfahrens auch diesem Maßstab entsprechen muss. Also nicht allein der Energieerhalt, sondern auch die Parametertreue zeichnet die Qualität eines Lösungsalgorithmus aus.
Und wie verhält sich nun die Lösung des symplektischen EC-Verfahrens? Urteile selbst:

Bitte einen Augenblick Geduld
während das Programm geladen wird!

Abb. Symplektisches Lösungsverfahren

Die erhaltene Zeitfunktion deckt sich mit unseren Erwartungen! Trotz verfahrensbedingter Ungenauigkeiten bleibt die Schwingungsamplitude konstant. Aber auch in diesem Fall stimmen die Frequenzen der erwarteten und die der errechnete Schwingung nicht überein! Unseren Anwendungen im Spielebereich tut dies keinen Abbruch! In den meisten Anwendungsfällen ist das symplektische EC-Verfahren das Verfahren der Wahl.

Das symplektische Verfahren ist in jedem Fall stabil trotz einfacher Implementation. Seiner Vorteile wegen ist dieses Verfahren Ursprung vieler weiterer Lösungsmethoden mit verbesserter Genauigkeit.

Symplektische Lösungen sind trotz Genauigkeitsdefiziten stabil und verstoßen nicht gegen die Erhaltungssätze von Energie und Impuls. Die Parametertreue ist hinreichend.

In diesem Programm werden alle drei EULER-CAUCHY-Verfahren gemeinsam dargestellt und mit einander verglichen. Mit Hilfe einer Energie- bzw. Frequenzmessung werden die Qualitätsmerkmale vergleichbar gemacht.

Das energetische Verhalten eine durch numerische Lösung einer Differentialgleichung gewonnene Bewegungsgleichung sagt sehr viel über die Stabilität der Lösung aus. Insbesondere wird dies deutlich, wenn ein Verlust freies System beobachtet wird. Dann nämlich muss die Summe aller Energien konstant bleiben! Also nicht steiegn, aber auch nicht fallen:

()

Zur Berechnung eines möglichen Energiedefizits wird die Differenz zwischen der Energie der Lage W_niveau und der in der Feder gespeicherten Energie W_spring plus der Bewegungsenergie W_kin gebildet.
In unsrem Beispiel bedeutet das:

()

Daher

()

Beziehen wir noch auf die Masse m, erhalten wir eine gut zu implementierende Form:

()

Nun ist offensichtlich, dass die erzeugte Schwingung eines zu bewertenden Lösungsverfahrens auch in Bezug auf die Frequenz der Schwingung untersucht werden muss. Also nicht allein der Energieerhalt, sondern auch die Parametertreue zeichnet die Qualität eines Lösungsalgorithmus aus.
Hinsichtlich Stabilität unterscheiden sich die implizite und die symplektische Löschung (bei einer Dämpfung δ = 0) nicht. Anders die Parametertreue. Da ist die implizite der symplektischen Lösung überlegen, wenn die Dämpfung δ > 0 ist. Diese Aussage trifft auch für kleinere Schrittweiten zu. Der im beigefügten Programmbeispiel verwendete Frequenzmesser hat eine Beobachtungsdauer von jeweils 4 s und gibt dabei die Frequenz der zu messenden Schwingung mit einer Genauigkeit von 10^-6 aus.

Symbol	Bedeutung
s	Strecke, Weg [m]
v	Geschwindigkeit [m/s]
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]
δ	Dämpfung [1/s] ()
ω₀	Kreisfrequenz der Eigenschwingung [rad/s] ()
W_pot	potentielle Energie [Nm] oder [Ws]
W_kin	kinetische Energie [Nm] oder [Ws]
t	Zeit [s]
Δt	Zeitquant, Bildwechseldauer [s]