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 t0 gegeben, und man möchte u(t) für t > t0 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
du

dt
=f(t,u)
mit einer zugehörigen Anfangsbedingung
u(t0)=u(0)
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
d2u

dt2
=f
t,u, du

dt

gegeben, so führt man die zwei neuen abhängigen Variablen
u0
=
u
u1
=
du

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

dt
=
f(t,u0,u1)
du0

dt
=
u1
ergibt. Dieser Vorgangsweise entspricht z.B. bei der Diskussion der Newtonschen Bewegungsgleichungen der Mechanik die Einführung der neuen Variablen "Geschwindigkeit", wodurch aus m d2x/dt2=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
t0, t1, t2, t3,
zu bestimmen. Häufig sind diese Stützstellen äquidistant, also Vielfache tn=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
+ Dt2

2
  d2u(t)

dt2
+ Dt3

6
  d3u(t)

dt3
+ O(Dt4)
u(t-Dt)
=
u(t) - Dt  du(t)

dt
+ Dt2

2
  d2u(t)

dt2
- Dt3

6
  d3u(t)

dt3
+ O(Dt4)
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)

Dt
= du(t)

dt
+ 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)

Dt
= du(t)

dt
+ 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)

dt
+ O(Dt2)
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 Dt2d3u/dt3).
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)

Dt2
= d2u(t)

dt2
+ O(Dt2)
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(tn), den Wert der unbekannten Funktion zum nächsten Zeitpunkt berechnet. Es genügt dann ja, diese Vorgangsweise sukzessive auf u(t0), u(t1), u(t2), ... anzuwenden, um die Lösung der Differentialgleichung beliebig weit in die Zukunft fortzusetzen. Der Startwert u(t0) dieses Verfahrens ist durch die Anfangsbedingung vorgeben.
Integriert man die Differentialgleichung
du

dt
=f(t,u)
formal zwischen den Grenzen tn und tn+1, so erhält man

tn+1

tn 
dt  du

dt
= u(tn+1) - u(tn) =
tn+1

tn 
dt f(t,u)
oder, mit der Bezeichnung u(n)=u(tn),
u(n+1) = u(n) +
tn+1

tn 
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 tn und tn+1 kennen-insbesondere natürlich zum Zeitpunkt tn+1 selbst! Die verschiedenen numerischen Integrationsverfahren unterscheiden sich daher hauptsächlich dadurch, wie dieses Integral in der Praxis approximiert wird.

 onestep.png 
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(n+1)
=
u(n) +
tn+1

tn 
dt 

f(t,u)|tn+ (t-tn df(t,u)

dt



tn 
+ (t-tn)2

2
  d2f(t,u)

dt2



tn 
+

=
u(n) + Dt f(t,u)|tn+ Dt2

2
  df(t,u)

dt



tn 
+ Dt3

6
  d2f(t,u)

dt2



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

Explizites Eulerverfahren



 eulerexpl.png 
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(n) unmittelbar berechnen kann) als Höhe der Rechtecksfläche zu approximieren. Das ergibt
u(n+1) = u(n) + Dt f(t,u)|tn
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(n) berechnen.

Implizites Eulerverfahren



 eulerimpl.png 
Wird als Höhe der Rechtecksfläche der Funktionswert am rechten Rand des Intervalls verwendet, so ergibt sich
u(n+1) = u(n) + Dt f(t,u)|tn+1
Wie man durch Einsetzen der Taylorentwicklung von f(t,u)|tn+1 um tn feststellt, ist auch dieses Verfahren von erster Ordnung. Es ist aber implizit, d.h. u(n+1) 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(n+1) 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



 implicit2o.png 
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(n+1) = u(n) + Dt

2
 [f(t,u)|tn + f(t,u)|tn+1]
(Beweis wieder durch Einsetzen der Taylorentwicklung von f(t,u)|tn+1). Da u(n+1) 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 [tn,tn+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.

 auxstep2o.png 
Um z.B. das Integral mit der "Mittelpunktsregel"
u(n+1) = u(n) + Dt f(t,u)|tn+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(tn+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(n) + Dt

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

2
,u(n)+ Dt

2
f(tn,u(n))

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
f1
=
f(tn,u(n))
f2
=
f
tn+ Dt

2
,u(n)+ Dt

2
 f1
f3
=
f
tn+ Dt

2
,u(n)+ Dt

2
 f2
f4
=
f(tn+Dt,u(n)+Dt f3)
Der eigentliche Integrationsschritt lautet dann
u(n+1) = u(n) + Dt

6
 [f1+2f2+2f3+f4]
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 tn bekannte) Funktionswerte zu verwenden.

 leapfrog.png 
Beim Leapfrog Verfahren integriert man f(t,u) nicht zwischen tn und tn+1, sondern zwischen tn-1 und tn+1
u(n+1) = u(n-1) +
tn+1

tn-1 
dt f(t,u)
und approximiert das Integral durch die Mittelpunktsregel
u(n+1) = u(n-1) + 2Dt f(t,u)|tn
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(0) gegeben ist, muß man sich die Lösung zum ersten [oder "minus ersten"] Zeitpunkt, u(1) [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(n), 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) uP für u(n+1) ermittelt, diesen beim impliziten Verfahren an Stelle von u(n+1) auf der rechten Seite der Gleichung einsetzt und daraus einen verbesserten Wert (Korrektor) uC für u(n+1) berechnet.
Am Beispiel des Eulerverfahrens könnte das z.B. so aussehen, daß man für den Prädiktor das explizite Verfahren verwendet
uP = u(n)+Dt f(tn,u(n))
und diesen in das implizite Verfahren einsetzt
u(n+1) = uC = u(n)+Dt f(tn+1,uP)
wobei nach dem Korrektorschritt gleich u(n+1)=uC 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(n) und u(n-1) dazu benützt, die Funktion f(t,u) für t > tn linear zu extrapolieren und von tn bis tn+1 zu integrieren
u(n+1)
=
u(n) +
tn+1

tn 
dt
f(tn,u(n))+(t-tn f(tn,u(n))-f(tn-1,u(n-1))

Dt

=
u(n) + Dt
3

2
 f(tn,u(n))- 1

2
 f(tn-1,u(n-1))
Der auf diese Weise gewonnene, vorläufige Wert von u(n+1) 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(n) und u(n+1) zwischen tn und tn+1 linear interpoliert
u(n+1)
=
u(n) +
tn+1

tn 
dt
f(tn,u(n))+(t-tn f(tn+1,u(n+1))-f(tn,u(n))

Dt

=
u(n) + Dt
1

2
 f(tn+1,u(n+1))+ 1

2
 f(tn,u(n))
(Das ist nichts anderes als das schon besprochene implizite Verfahren zweiter Ordnung.)

 abam2o.png 
Beim kombinierten Adams-Bashforth/Adams-Moulton Verfahren zweiter Ordnung wird also zuerst der Prädiktor uP aus
uP = u(n) + Dt

2
[3f(tn,u(n))-f(tn-1,u(n-1))]
berechnet und dann im Korrektorschritt an Stelle von u(n+1) auf der rechten Seite der impliziten Gleichung eingesetzt
u(n+1) = u(n) + Dt

2
[f(tn+1,uP)+f(tn,u(n))]
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(n) bzw. (beim Korrektor) u(n-2), u(n-1), u(n) und u(n+1), betrachten und bei der Integration zwischen tn und tn+1 durch Polynome dritten Grades extra- bzw. interpolieren
uP
=
u(n) + Dt

24
[55f(tn,u(n))-59f(tn-1,u(n-1))+37f(tn-2,u(n-2))-9f(tn-3,u(n-3))]
u(n+1)
=
u(n) + Dt

24
[9f(tn+1,uP)+19f(tn,u(n))-5f(tn-1,u(n-1))+f(tn-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 Dt0 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
m d2x

dt2
= - m w2x
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)

w
 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 = m

2
 v2 + mw2

2
 x2
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
dv

dt
=
- w2x
dx

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

dt
= d

dt



v
x



=


- w2x
v



angewendet.
Für das explizite Eulerverfahren ergibt sich
v(n+1)
=
v(n) - Dt w2 x(n)
x(n+1)
=
x(n) + Dt v(n)
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(n+1)
=
v(n) - Dt w2 x(n+1)
x(n+1)
=
x(n) + Dt v(n+1)
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(n+1)
=
v(n) - Dt w2 x(n)
x(n+1)
=
x(n) + Dt v(n+1)
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 t1/2, t3/2, t5/2, ... für die Geschwindigkeit und t0, t1, t2, ... für den Ort aus:

 lfgrd.png 
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 w2 x(n)
x(n+1)
=
x(n) + 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(n) durch
v(n) 1

2
 (v(n-1/2)+v(n+1/2))
approximieren. Den fehlenden Startwert v(-1/2) kann man aus den Anfangsbedingungen x(0) und v(0) z.B. durch Taylorentwicklung
v
t0- Dt

2

= v(t0)- Dt

2
dv

dt



t0 
+ 1

2

Dt

2

2

 
d2v

dt2



t0 
+
und Einsetzen von dv/dt=-w2x und d2v/dt2=-w2dx/dt=-w2v (aus der Differentialgleichung) gewinnen
v-1/2 v(0) + Dt

2
 w2x(0)- 1

2

Dt

2

2

 
 w2v(0)



File translated from TEX by TTH, version 3.80.
On 14 Apr 2008, 20:24.