Integration gewöhnlicher Differentialgleichungen

Anfangswertprobleme

Die Integration von Differentialgleichungen bildet eine der Grundaufgaben der Physik und verwandter naturwissenschaftlicher Disziplinen. Man denke nur an die Integration von Bewegungsgleichungen oder an die Tatsache, daß Feld- oder Kontinuumstheorien meist mit Hilfe von partiellen Differentialgleichungen formuliert werden. Da es nur in Ausnahmefällen gelingt, die Lösungen von Differentialgleichungen geschlossen anzugeben, ist man vor allem bei praktischen Problemen in der Regel auf numerische Verfahren angewiesen.

Im einfacheren Fall von gewöhnlichen (im Gegensatz zu partiellen) Differentialgleichungen hängt die gesuchte Funktion (bzw. hängen im Fall eines Systems von gewöhnlichen Differentialgleichungen die Funktionen) u(t) nur von einer einzigen unabhängigen Variablen-meist "Zeit" t genannt-ab, und die Differentialgleichung ist eine Beziehung zwischen u, den Ableitungen von u nach t und der unabhängigen Variablen t. Die Ordnung der höchsten vorkommenden Ableitung nennt man den Grad der Differentialgleichung. Meist ist man an der Lösung eines Anfangswertproblems interessiert, d.h. die Werte der unbekannten Funktion u sowie (im Fall einer Differentialgleichung höherer Ordnung) einer Reihe von Ableitungen von u sind zum Zeitpunkt t₀ gegeben, und man möchte u(t) für t > t₀ bestimmen. Im Gegensatz dazu sind bei einem Randwertproblem die Werte von u oder der Ableitungen von u an den Enden eines Intervalls gegeben und die Lösung im Inneren gesucht.

Wir schreiben das Anfangswertproblem abstrakt meist in der Form einer Differentialgleichung oder eines Systems von Differentialgleichungen erster Ordnung

=f(t,u)

mit einer zugehörigen Anfangsbedingung

u(t₀)=u⁽⁰⁾

Wenn es sich um ein System von Differentialgleichungen handelt, steht u als Abkürzung für einen Vektor von Unbekannten; ebenso steht auf der rechten Seite f für einen ganzen Vektor von Funktionen. Hat die unabhängige Variable nicht die Bedeutung einer Zeit-, sondern einer Ortskoordinate, so kann man statt t natürlich auch x schreiben.

Gewöhnliche Differentialgleichungen höherer Ordnung kann man formal immer auf Systeme erster Ordnung zurückführen. Ist z.B. die Differentialgleichung zweiter Ordnung

d²u

dt²

æ
è

t,u,

ö
ø

gegeben, so führt man die zwei neuen abhängigen Variablen

u₀

u₁

ein, für die sich dann das System von Differentialgleichungen

du₁

f(t,u₀,u₁)

du₀

u₁

ergibt. Dieser Vorgangsweise entspricht z.B. bei der Diskussion der Newtonschen Bewegungsgleichungen der Mechanik die Einführung der neuen Variablen "Geschwindigkeit", wodurch aus m d²x/dt²=F die Differentialgleichung erster Ordnung m dv/dt=F mit der Definition von v=dx/dt als zweiter Gleichung wird.

Approximation von Ableitungen durch Differenzenquotienten

Da Computer nur endlich viele Zahlen speichern können und man in endlicher Zeit auch nur endlich viele Rechenschritte durchführen kann, wird bei der numerischen Behandlung von Differentialgleichungen die unabhängige Variable t diskretisiert, d.h. man muß sich damit zufriedengeben, die unbekannte Funktion nicht für alle Werte von t, sondern nur auf einem diskreten Gitter von Stützstellen

t₀, t₁, t₂, t₃, ¼

zu bestimmen. Häufig sind diese Stützstellen äquidistant, also Vielfache t_n=nDt eines konstanten Zeitschritts Dt.

Schränkt man sich auf die Kenntnis der Funktion u an diskreten Stützstellen ein, ergibt sich umgekehrt das Problem, aus diesem Satz von Funktionswerten die Ableitungen von u an den Stützstellen zu berechnen. Einen einfachen Satz von Näherungsformeln erhält man, indem man die Taylorentwicklung von u um einen beliebigen Zeitpunkt t betrachtet

u(t+Dt)

u(t) + Dt

du(t)

Dt²

d²u(t)

dt²

Dt³

d³u(t)

dt³

+ O(Dt⁴)

u(t-Dt)

u(t) - Dt

du(t)

Dt²

d²u(t)

dt²

Dt³

d³u(t)

dt³

+ O(Dt⁴)

Bringt man in der ersten Gleichung u(t) auf die linke Seite und dividiert die ganze Gleichung durch Dt, so erhält man die sogenannte Vorwärtsdifferenz

u(t+Dt)-u(t)

du(t)

+ O(Dt)

Die erste Ableitung von u läßt sich also aus den Funktionswerten am Stützpunkt t und an seinem rechten Nachbarpunkt approximieren. Der Fehler dieser Näherung geht für kleine Zeitschritte linear mit Dt gegen Null, und man sagt daher, daß die Formel "von erster Ordnung" ist. Analog ergibt sich aus der zweiten Gleichung die Rückwärtsdifferenz

u(t)-u(t-Dt)

du(t)

+ O(Dt)

die ebenfalls von erster Ordnung ist. Subtrahiert man die zweite von der ersten Gleichung, so erhält man nach Division durch 2Dt die Zentraldifferenz

u(t+Dt)-u(t-Dt)

2Dt

du(t)

+ O(Dt²)

Da die geraden Ableitungen von u in den beiden Taylorentwicklungen einander kompensieren, ist diese Näherung von zweiter Ordnung (der führende Korrekturterm ist proportional zu Dt²d³u/dt³).

Addiert man die Taylorentwicklungen für u(t+Dt) und u(t-Dt), so fallen die ungeraden Ableitungen von u weg, und man erhält eine symmetrische Formel für die zweite Ableitung

u(t+Dt)-2u(t)+u(t-Dt)

Dt²

d²u(t)

dt²

+ O(Dt²)

die ebenfalls von zweiter Ordnung ist.

Durch Hinzunahme weiterer Stützpunkte kann man auf ähnliche Weise auch Approximationen höherer Ordnung gewinnen; die letzten beiden Formeln sind aber trotz ihrer Einfachheit sehr genau und werden daher gerne verwendet.

Einschrittverfahren

Grundproblem

Eine Differentialgleichung vom Typ des Anfangswertproblems kann als gelöst betrachtet werden, wenn es gelingt, eine Methode anzugeben, wie man, ausgehend von u(t_n), den Wert der unbekannten Funktion zum nächsten Zeitpunkt berechnet. Es genügt dann ja, diese Vorgangsweise sukzessive auf u(t₀), u(t₁), u(t₂), ... anzuwenden, um die Lösung der Differentialgleichung beliebig weit in die Zukunft fortzusetzen. Der Startwert u(t₀) dieses Verfahrens ist durch die Anfangsbedingung vorgeben.

Integriert man die Differentialgleichung

=f(t,u)

formal zwischen den Grenzen t_n und t_n+1, so erhält man

ó
õ

t_n+1

t_n

= u(t_n+1) - u(t_n) =

ó
õ

t_n+1

t_n

dt f(t,u)

oder, mit der Bezeichnung u⁽ⁿ⁾=u(t_n),

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ +

ó
õ

t_n+1

t_n

dt f[t,u(t)]

Diese Gleichung und die folgende Abbildung zeigen auch schon das Grundproblem aller Integrationsverfahren: Um den Funktionswert zum neuen Zeitpunkt bestimmen zu können, müßte man zur Berechnung des Integrals auf der rechten Seite die Lösung u(t) schon für alle Zeiten zwischen t_n und t_n+1 kennen-insbesondere natürlich zum Zeitpunkt t_n+1 selbst! Die verschiedenen numerischen Integrationsverfahren unterscheiden sich daher hauptsächlich dadurch, wie dieses Integral in der Praxis approximiert wird.

Die Güte von Integrationsverfahren beurteilt man u.a. durch den Vergleich der numerischen Approximation mit der formalen Taylorentwicklung der exakten Lösung, die durch

u⁽ⁿ⁺¹⁾

u⁽ⁿ⁾ +

ó
õ

t_n+1

t_n

ì
í
î

f(t,u)|_{t_n}+ (t-t_n)

df(t,u)

ê
ê

t_n

(t-t_n)²

d²f(t,u)

dt²

ê
ê

t_n

+ ¼

ü
ý
þ

u⁽ⁿ⁾ + Dt f(t,u)|_{t_n}+

Dt²

df(t,u)

ê
ê

t_n

Dt³

d²f(t,u)

dt²

ê
ê

t_n

+ ¼

gegeben ist. Ein numerisches Verfahren, das bis zum Term proportional Dt^p mit dieser Entwicklung übereinstimmt, nennt man ein Verfahren p-ter Ordnung.

Explizites Eulerverfahren

Die einfachste Näherung besteht darin, das Integral über f(t,u) durch eine "Riemannsumme" mit dem Funktionswert am linken Stützpunkt (dem einzigen Wert der Kurve f[t,u(t)], den man mit Hilfe von u⁽ⁿ⁾ unmittelbar berechnen kann) als Höhe der Rechtecksfläche zu approximieren. Das ergibt

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ + Dt f(t,u)|_{t_n}

Wie ein Vergleich mit der Taylorentwicklung der exakten Lösung zeigt, ist das ein Verfahren erster Ordnung. Es ist außerdem explizit, d.h. der Wert von u zum neuen Zeitpunkt kommt nur auf der linken Seite der Gleichung vor, und die rechte Seite läßt sich unmittelbar durch Einsetzen des alten Wertes u⁽ⁿ⁾ berechnen.

Implizites Eulerverfahren

Wird als Höhe der Rechtecksfläche der Funktionswert am rechten Rand des Intervalls verwendet, so ergibt sich

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ + Dt f(t,u)|_{t_n+1}

Wie man durch Einsetzen der Taylorentwicklung von f(t,u)|_{t_n+1} um t_n feststellt, ist auch dieses Verfahren von erster Ordnung. Es ist aber implizit, d.h. u⁽ⁿ⁺¹⁾ kommt sowohl auf der linken als auch der rechten Seite vor, und es hängt von der speziellen Gestalt der Gleichung ab, ob und wie leicht man sie nach u⁽ⁿ⁺¹⁾ auflösen kann. Dies ist der prinzipielle Nachteil von impliziten Verfahren, die sich aber gegenüber den expliziten im allgemeinen durch größere Stabilität auszeichnen.

Implizites Verfahren zweiter Ordnung

Approximiert man das Integral nicht durch eine Rechtecks-, sondern durch die in der Regel genauere Trapezfläche, so erhält man das Verfahren zweiter Ordnung

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ +

[f(t,u)|_{t_n} + f(t,u)|_{t_n+1}]

(Beweis wieder durch Einsetzen der Taylorentwicklung von f(t,u)|_{t_n+1}). Da u⁽ⁿ⁺¹⁾ zur Berechnung der rechten Seite schon bekannt sein müßte, handelt es sich auch hier wieder um ein implizites Verfahren.

Mehrschrittverfahren

Verfahren zweiter Ordnung mit Hilfsschritt

Eine Möglichkeit, den Verlauf des Integranden f[t,u(t)] im Intervall [t_n,t_n+1] besser zu beschreiben, besteht darin, zusätzliche Schätzwerte für diese Funktion im Inneren oder am Rand des Intervalls zu ermitteln und in eine Integrationsformel höherer Ordnung einzusetzen. Allerdings ist das mit einem größeren Rechenaufwand verbunden, da die vorkommenden Funktionen pro Zeitschritt öfter ausgewertet werden müssen.

Um z.B. das Integral mit der "Mittelpunktsregel"

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ + Dt f(t,u)|_{t_n+1/2}

(Stützpunkt in der Mitte des Intervalls) zu berechnen, verschafft man sich zunächst mit Hilfe des expliziten Eulerverfahrens in einem Hilfsschritt der Schrittweite Dt/2 einen Schätzwert für u^(n+1/2)=u(t_n+Dt/2), der dann in f(t,u) eingesetzt wird. Der Algorithmus besteht also aus einem Hilfs- und dem eigentlichen Integrationsschritt:

u^(n+1/2)

u⁽ⁿ⁾ +

f(t_n,u⁽ⁿ⁾)

u⁽ⁿ⁺¹⁾

u⁽ⁿ⁾ + Dt f(t_n+1/2,u^(n+1/2))

Dieses Verfahren ist explizit und, wie man durch Taylorentwicklung aller Größen um t_n feststellt, von zweiter Ordnung. Eine andere Schreibweise, bei der man die Hilfsgröße gleich in f(t,u) einsetzt, ist

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ + Dt f

é
ë

t_n+

,u⁽ⁿ⁾+

f(t_n,u⁽ⁿ⁾)

ù
û

Runge-Kutta Verfahren vierter Ordnung

Das obige Verfahren ist der einfachste Vertreter einer Klasse von expliziten Integrationsalgorithmen, bei denen durch geschickte Wahl der Hilfspunkte eine möglichst hohe Ordnung erreicht wird. Beim Runge-Kutta Standardverfahren vierter Ordnung müssen pro Zeitschritt vier Schätzwerte für f in der Mitte und am rechten Rand des Intervalls berechnet werden

f₁

f(t_n,u⁽ⁿ⁾)

f₂

æ
è

t_n+

,u⁽ⁿ⁾+

f₁

ö
ø

f₃

æ
è

t_n+

,u⁽ⁿ⁾+

f₂

ö
ø

f₄

f(t_n+Dt,u⁽ⁿ⁾+Dt f₃)

Der eigentliche Integrationsschritt lautet dann

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ +

[f₁+2f₂+2f₃+f₄]

Ist ein Verfahren von hoher Ordnung, so besagt dies im Prinzip nur, daß es für kleine Zeitschritte sehr genau wird. Bei größeren Zeitschritten kann aber durchaus ein Verfahren niedrigerer Ordnung das "bessere" sein.

Leapfrog Verfahren

Eine andere Möglichkeit, die Ordnung eines Verfahrens zu erhöhen, besteht darin, weiter in der Vergangenheit zurückliegende (und daher zum Zeitpunkt t_n bekannte) Funktionswerte zu verwenden.

Beim Leapfrog Verfahren integriert man f(t,u) nicht zwischen t_n und t_n+1, sondern zwischen t_n-1 und t_n+1

u⁽ⁿ⁺¹⁾ = u^(n-1) +

ó
õ

t_n+1

t_n-1

dt f(t,u)

und approximiert das Integral durch die Mittelpunktsregel

u⁽ⁿ⁺¹⁾ = u^(n-1) + 2Dt f(t,u)|_{t_n}

Dieses Verfahren ist explizit, wegen der Verwendung der Mittelpunktsregel von zweiter Ordnung und für gewisse Arten von Differentialgleichungen sehr gut geeignet. Zum Unterschied von Einschrittverfahren verknüpft es Funktionswerte zu drei aufeinanderfolgenden Zeitpunkten. Da bei einer Differentialgleichung erster Ordnung durch die Anfangsbedingung nur u⁽⁰⁾ gegeben ist, muß man sich die Lösung zum ersten [oder "minus ersten"] Zeitpunkt, u⁽¹⁾ [oder u^(-1)], mit Hilfe eines anderen Verfahrens oder durch Taylorentwicklung aus der ursprünglichen Differentialgleichung bestimmen. Man bezeichnet daher das Leapfrog Verfahren als "nicht selbst-startend".

Prädiktor-Korrektor Verfahren

Adams-Bashforth/Adams-Moulton Verfahren

Implizite Verfahren sind im allgemeinen stabiler als explizite, haben aber den Nachteil, daß im Algorithmus u⁽ⁿ⁾, der Wert der Unbekannten zum neuen Zeitschritt, auf beiden Seiten einer Gleichung vorkommt und man die Lösung eigentlich schon kennen müßte, um sie berechnen zu können. Das kann man zu umgehen versuchen, indem man zunächst mit einem expliziten Verfahren einen Schätzwert (Prädiktor) u^P für u⁽ⁿ⁺¹⁾ ermittelt, diesen beim impliziten Verfahren an Stelle von u⁽ⁿ⁺¹⁾ auf der rechten Seite der Gleichung einsetzt und daraus einen verbesserten Wert (Korrektor) u^C für u⁽ⁿ⁺¹⁾ berechnet.

Am Beispiel des Eulerverfahrens könnte das z.B. so aussehen, daß man für den Prädiktor das explizite Verfahren verwendet

u^P = u⁽ⁿ⁾+Dt f(t_n,u⁽ⁿ⁾)

und diesen in das implizite Verfahren einsetzt

u⁽ⁿ⁺¹⁾ = u^C = u⁽ⁿ⁾+Dt f(t_n+1,u^P)

wobei nach dem Korrektorschritt gleich u⁽ⁿ⁺¹⁾=u^C gesetzt wurde. Im Prinzip könnte der Korrektorschritt auch mehrere Male wiederholt werden. Dies würde auf eine iterative Lösung der impliziten Gleichung hinauslaufen, liefert in der Regel aber keine wesentliche Verbesserung, da sich dadurch die Ordung des Verfahrens (hier erste Ordnung) nicht ändert.

Eine Erhöhung der Ordnung (und damit der Genauigkeit) läßt sich zwar durch Runge-Kutta-artige Verfahren mit Hilfsschritten erreichen, doch wird das um den Preis eines vermehrten Rechenaufwandes erkauft, da bei jedem Zeitschritt die Funktion f(t,u) mehrere Male ausgewertet werden muß. Wie beim Leapfrog-Schema kann man einen Algorithmus höherer Ordnung jedoch auch dadurch gewinnen, daß man statt der Hilfsschritte auf schon berechnete, frühere Werte der Lösung, u^(n-1), u^(n-2), ..., zurückgreift. Das kostet zwar-ebenso wie Hilfsschritte-zusätzlichen Speicherplatz, ist aber kaum mit rechnerischem Mehraufwand verbunden. Dieses Prinzip soll an Hand des Adams-Bashforth/Adams-Moulton Verfahrens zweiter Ordnung erläutert werden.

Beim Prädiktorschritt, dem (expliziten) Adams-Bashforth Verfahren zweiter Ordnung, werden u⁽ⁿ⁾ und u^(n-1) dazu benützt, die Funktion f(t,u) für t > t_n linear zu extrapolieren und von t_n bis t_n+1 zu integrieren

u⁽ⁿ⁺¹⁾

u⁽ⁿ⁾ +

ó
õ

t_n+1

t_n

é
ë

f(t_n,u⁽ⁿ⁾)+(t-t_n)

f(t_n,u⁽ⁿ⁾)-f(t_n-1,u^(n-1))

ù
û

u⁽ⁿ⁾ + Dt

é
ë

f(t_n,u⁽ⁿ⁾)-

f(t_n-1,u^(n-1))

ù
û

Der auf diese Weise gewonnene, vorläufige Wert von u⁽ⁿ⁺¹⁾ wird nun auf der rechten Seite des (impliziten) Adams-Moulton Verfahrens zweiter Ordnung eingesetzt, das sich ergibt, wenn man f(t,u) mit Hilfe von u⁽ⁿ⁾ und u⁽ⁿ⁺¹⁾ zwischen t_n und t_n+1 linear interpoliert

u⁽ⁿ⁺¹⁾

u⁽ⁿ⁾ +

ó
õ

t_n+1

t_n

é
ë

f(t_n,u⁽ⁿ⁾)+(t-t_n)

f(t_n+1,u⁽ⁿ⁺¹⁾)-f(t_n,u⁽ⁿ⁾)

ù
û

u⁽ⁿ⁾ + Dt

é
ë

f(t_n+1,u⁽ⁿ⁺¹⁾)+

f(t_n,u⁽ⁿ⁾)

ù
û

(Das ist nichts anderes als das schon besprochene implizite Verfahren zweiter Ordnung.)

Beim kombinierten Adams-Bashforth/Adams-Moulton Verfahren zweiter Ordnung wird also zuerst der Prädiktor u^P aus

u^P = u⁽ⁿ⁾ +

[3f(t_n,u⁽ⁿ⁾)-f(t_n-1,u^(n-1))]

berechnet und dann im Korrektorschritt an Stelle von u⁽ⁿ⁺¹⁾ auf der rechten Seite der impliziten Gleichung eingesetzt

u⁽ⁿ⁺¹⁾ = u⁽ⁿ⁾ +

[f(t_n+1,u^P)+f(t_n,u⁽ⁿ⁾)]

In analoger Weise läßt sich z.B. auch ein Prädiktor-Korrektor Verfahren vierter Ordnung herleiten. Man muß nur, statt zwei Zeitpunkte zu benützen und linear zu extra- bzw. interpolieren, jeweils vier aufeinanderfolgende Zeitpunkte, nämlich (beim Prädiktor) u^(n-3), u^(n-2), u^(n-1) und u⁽ⁿ⁾ bzw. (beim Korrektor) u^(n-2), u^(n-1), u⁽ⁿ⁾ und u⁽ⁿ⁺¹⁾, betrachten und bei der Integration zwischen t_n und t_n+1 durch Polynome dritten Grades extra- bzw. interpolieren

u^P

u⁽ⁿ⁾ +

[55f(t_n,u⁽ⁿ⁾)-59f(t_n-1,u^(n-1))+37f(t_n-2,u^(n-2))-9f(t_n-3,u^(n-3))]

u⁽ⁿ⁺¹⁾

u⁽ⁿ⁾ +

[9f(t_n+1,u^P)+19f(t_n,u⁽ⁿ⁾)-5f(t_n-1,u^(n-1))+f(t_n-2,u^(n-2))]

Ein offensichtlicher Nachteil der Adams-Bashforth/Adams-Moulton Verfahren besteht darin, daß sie-wie Leapfrog-nicht selbst-startend sind und man die Lösung für die ersten paar Zeitschritte mit einem anderen Verfahren oder durch Taylorentwicklung aus der Differentialgleichung und den gebenen Anfangsbedingungen berechnen muß.

Kriterien

Es gibt leider keinen Universalalgorithmus zur numerischen Lösung von Differentialgleichungen, sondern eine Vielzahl von speziellen Verfahren, die jeweils nur für bestimmte Arten von Gleichungen geeignet sind. Welches Lösungsverfahren man im konkreten Fall verwendet, wird von verschiedenen Faktoren abhängen. Auf jeden Fall muß das gewählte Differenzenverfahren konsistent mit der Differentialgleichung sein, d.h. es muß im Limes Dt®0 gegen die vorgegebene Differentialgleichung streben. Weiters muß das Verfahren stabil sein. Das heißt, daß sich von einem Zeitschritt zum nächsten Rundungsfehler nicht verstärken dürfen, sondern beschränkt bleiben müssen. Auch Stabilität garantiert aber noch nicht Genauigkeit bzw. Konvergenz, d.h. daß die numerische Lösung bei Verkleinerung der Schrittweite tatsächlich gegen die exakte Lösung der Differentialgleichung strebt.

Kann man die Eigenschaften eines Integrationsverfahrens nicht theoretisch analysieren, so empfiehlt es sich, den Algorithmus wenigstens an einem exakt lösbaren Problem ähnlicher Art zu überprüfen.

Beispiele

Die Bewegungsgleichung des Harmonischen Oszillators in einer Dimension

d²x

dt²

= - m w²x

ist nicht nur eines der wenigen exakt behandelbaren physikalischen Modelle, sondern auch der Prototyp für eine Differentialgleichung zweiter Ordnung mit oszillierender, beschränkter (nicht zunehmender) Lösung. Die analytische Lösung lautet

x(t) = x(0) coswt +

v(0)

sinwt

wobei x(0) und v(0) durch die Anfangsbedingungen gegeben sind. Als Kriterium für die Brauchbarkeit eines numerischen Lösungsverfahrens kann hier neben dem direkten Vergleich der erzeugten Trajektorien vor allem die Erhaltung der Gesamtenergie

E =

v² +

mw²

x²

dienen.

Die Differentialgleichung zweiter Ordnung wird zunächst auf eine System von Gleichungen erster Ordnung umgeschrieben, indem man für die erste Ableitung dx/dt den neuen "Namen" v einführt

- w²x

Die Unbekannten v und x werden zu einem Vektor u=(v,x) zusammengefaßt und die Formeln der vorigen Abschnitte komponentenweise auf

æ
ç
ç
è

ö
÷
÷
ø

æ
ç
ç
è

- w²x

ö
÷
÷
ø

angewendet.

Für das explizite Eulerverfahren ergibt sich

v⁽ⁿ⁺¹⁾

v⁽ⁿ⁾ - Dt w² x⁽ⁿ⁾

x⁽ⁿ⁺¹⁾

x⁽ⁿ⁾ + Dt v⁽ⁿ⁾

Das Verfahren ist jedoch für den Harmonischen Oszillator instabil: die Gesamtenergie nimmt im Lauf der Zeit systematisch zu, und die Lösung explodiert.

Da die Differentialgleichung linear ist, führt das implizite Eulerverfahren auf ein System von linearen Gleichungen

v⁽ⁿ⁺¹⁾

v⁽ⁿ⁾ - Dt w² x⁽ⁿ⁺¹⁾

x⁽ⁿ⁺¹⁾

x⁽ⁿ⁾ + Dt v⁽ⁿ⁺¹⁾

das leicht nach den Werten von v und x zum neuen Zeitschritt auflösbar ist. Das Verfahren ist stabil, die Gesamtenergie ist jedoch nicht erhalten, sondern nimmt systematisch ab, d.h. Rundungsfehler werden zwar nicht verstärkt, der Oszillator wird aber langfristig sogar zum Stillstand kommen.

Daß kleine Änderungen im Integrationsschema entscheidend für die Brauchbarkeit eines Verfahrens sein können, zeigt das Cromer-Verfahren: Es unterscheidet sich von den beiden Eulerverfahren nur dadurch, daß die Gleichung für die Geschwindigkeit explizit, jene für den Ort aber implizit behandelt wird

v⁽ⁿ⁺¹⁾

v⁽ⁿ⁾ - Dt w² x⁽ⁿ⁾

x⁽ⁿ⁺¹⁾

x⁽ⁿ⁾ + Dt v⁽ⁿ⁺¹⁾

d.h. man setzt den in der ersten Gleichung berechneten neuen (aktuellen) Wert von v sofort in die zweite Gleichung ein. Dieses Verfahren ist beim Harmonischen Oszillator für hinreichend kleine Zeitschritte stabil und erhält langfristig die Gesamtenergie.

Für Systeme mit geschwindigkeitsunabhängigen Kräften geht man bei der Implementation des Leapfrog Verfahrens meist von versetzten Gittern von der Form t_1/2, t_3/2, t_5/2, ... für die Geschwindigkeit und t₀, t₁, t₂, ... für den Ort aus:

Formuliert man das Leapfrog Verfahren nämlich für eine Diskretisierung mit der Schrittweite Dt/2, so erhält man die Gleichungen

v^(n+1/2)

v^(n-1/2) - Dt w² x⁽ⁿ⁾

x⁽ⁿ⁺¹⁾

x⁽ⁿ⁾ + Dt v^(n+1/2)

d.h. daß man zur Berechnung der Geschwindigkeiten zu den halbzahligen Zeitpunkten nur die Orte zu den ganzzahligen Zeitpunkten braucht und umgekehrt. (Geschwindigkeiten zu ganzzahligen und Orte zu halbzahligen Zeitpunkten werden einfach nicht berechnet.) Das Leapfrog Verfahren ist beim Harmonischen Oszillator ebenfalls für hinreichend kleine Zeitschritte stabil, und die Gesamtenergie ist erhalten. Ein Nachteil besteht darin, daß Orte und Geschwindigkeiten nicht zu denselben Zeitpunkten bekannt sind, doch kann man, konsistent mit der Ordnung des Verfahrens, v⁽ⁿ⁾ durch

v⁽ⁿ⁾ »

(v^(n-1/2)+v^(n+1/2))

approximieren. Den fehlenden Startwert v^(-1/2) kann man aus den Anfangsbedingungen x⁽⁰⁾ und v⁽⁰⁾ z.B. durch Taylorentwicklung

æ
è

t₀-

ö
ø

= v(t₀)-

ê
ê

t₀

æ
è

ö
ø

d²v

dt²

ê
ê

t₀

+ ¼

und Einsetzen von dv/dt=-w²x und d²v/dt²=-w²dx/dt=-w²v (aus der Differentialgleichung) gewinnen

v^-1/2 » v⁽⁰⁾ +

w²x⁽⁰⁾-

æ
è

ö
ø

w²v⁽⁰⁾

File translated from T_EX by T_TH, version 3.80.
On 14 Apr 2008, 20:24.