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 höchste vorkommende Ableitung nennt man den Grad der Differentialgleichung. Meist ist man an der Lösung eines Anfangswertproblems interessiert, d.h. der Wert der unbekannten Funktion u ist zu einem 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

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. in der Mechanik die Einführung der neuen Variablen ``Geschwindigkeit'' in die Newtonschen Bewegungsgleichungen, 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 darauf beschränken, 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 t_n=nDt eines konstanten Zeitschritts Dt.

Sind oder wären die Funktionswerte an den Stützstellen bekannt, ergibt sich umgekehrt das Problem, die Ableitungen von u an diesen 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 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) 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

Das Anfangswertproblem 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, denn man muß dieses Verfahren nur sukzessive auf u(t₀), u(t₁), u(t₂), ... anwenden, um die Lösung der Differentialgleichung beliebig weit in die Zukunft zu propagieren. Der Startwert u(t₀) ist durch die Anfangsbedingung gegeben.

Integriert man die Differentialgleichung

du dt =f(t,u)

formal zwischen den Grenzen t_n und t_n+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⁽ⁿ⁾=u(t_n) usw.,

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 Werte von t zwischen t_n und t_n+1 kennen! Die verschiedenen numerischen Integrationsverfahren unterscheiden sich dadurch, wie dieses Integral in der Praxis approximiert wird.

  

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

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  + ¼

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(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⁽ⁿ⁾ berechnen.

Implizites Eulerverfahren

  

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 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 im allgemeinen aber durch größere Stabilität auszeichnen.

Implizites Verfahren zweiter Ordnung

  

Approximiert man das Integral statt durch eine Rechtecks- genauer durch eine 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⁽ⁿ⁺¹⁾ zur Berechnung der rechten Seite schon bekannte 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(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(t_n+Dt/2), der dann in f(t,u) eingesetzt wird. Das Verfahren 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))

Es 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 explizit 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 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(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 aufeinander folgenden Zeitpunkten. Da bei einer Differentialgleichung erster Ordnung durch die Anfangsbedingung nur u⁽⁰⁾ gegeben ist, muß man sich daher die Lösung zum ersten Zeitpunkt, u⁽¹⁾, mit Hilfe eines anderen Verfahrens oder durch Taylorentwicklung aus der Differentialgleichung selbst bestimmen. Man bezeichnet daher das Leapfrog-Verfahren als ``nicht selbst-startend''.

Kriterien

Es gibt leider keinen Universal-Algorithmus 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

m d2x dt2 = - m w2x

ist nicht nur eines der wenigen exakt behandelbaren physikalischen Modellsysteme, 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 dient hier neben dem direkten Vergleich der erzeugten Trajektorien vor allem die Erhaltung der Gesamtenergie

E = m 2  v2 + mw2 2  x2

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

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 kommt aber langfristig zum Stillstand.

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 t_1/2, t_3/2, t_5/2, ... für die Geschwindigkeit bzw. 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 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⁽ⁿ⁾ 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⁽⁰⁾ und v⁽⁰⁾ 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=-w²x und d²v/dt²=-w²dx/dt=-w²v (aus der Differentialgleichung) gewinnen

v-1/2 » v(0) + Dt 2  w2x(0)- 1 2 æ ç è Dt 2 ö ÷ ø 2    w2v(0)