Einschrittverfahren

Zweck und Grundprinzip

Einschrittverfahren sind numerische Methoden zur Lösung von Anfangswertproblemen gewöhnlicher Differentialgleichungen. Solche Probleme beschreiben dynamische Vorgänge und bestehen aus einer Differentialgleichung sowie einem bekannten Startwert. Gesucht ist eine Funktion y mit

y'(t)=f(t,y(t)),\quad y(t_0)=y_0.

Die Lösung kann skalare Werte oder Vektoren y(t)=(y_1(t),\dotsc,y_d(t)) haben; im zweiten Fall liegt ein d-dimensionales Differentialgleichungssystem vor. Nicht für jedes Problem lässt sich die Lösung mit Methoden der Analysis ausdrücklich berechnen. Einschrittverfahren liefern dann Näherungen an einzelnen Stellen t_0<t_1<\dotsc<t_n.

Sie berechnen nacheinander y_0,y_1,\dotsc als Näherungen für y(t_0),y(t_1),\dotsc. Entscheidend ist: Für den nächsten Wert y_{j+1} verwenden sie nur die aktuelle Näherung y_j. Mehrschrittverfahren beziehen dagegen zusätzlich ältere Näherungen ein. Allgemein hat ein Einschrittverfahren die Form

y_{j+1}=y_j+h_j\Phi(t_j,y_j,y_{j+1},h_j),\quad h_j=t_{j+1}-t_j.

\Phi heißt Verfahrensfunktion. Ist sie unabhängig von y_{j+1}, ist das Verfahren explizit und der neue Wert wird direkt berechnet. Hängt sie von y_{j+1} ab, ist das Verfahren implizit; dann muss in jedem Schritt eine Gleichung gelöst werden.

Euler-Verfahren und einfache Verbesserungen

Das explizite Euler-Verfahren nähert die Lösung in jedem Schritt durch ihre Tangente am aktuellen Näherungspunkt an. Seine Vorschrift lautet

y_{j+1}=y_j+h_jf(t_j,y_j).

Es besitzt Konvergenzordnung 1. Bei exponentiellem Wachstum y'(t)=\lambda y(t) mit y(0)=y_0 lautet die exakte Lösung y(t)=y_0e^{\lambda t}.

Das implizite Euler-Verfahren verwendet dagegen die Steigung am Endpunkt:

y_{j+1}=y_j+h_jf(t_{j+1},y_{j+1}).

Das implizite Trapez-Verfahren mittelt die Steigungen an beiden Endpunkten und hat Konvergenzordnung 2:

y_{j+1}=y_j+\frac{h}{2}\Big(f(t_j,y_j)+f(t_{j+1},y_{j+1})\Big).

Das explizite Heun-Verfahren ersetzt darin den unbekannten Endwert durch einen Euler-Schritt und hat ebenfalls Ordnung 2:

y_{j+1}=y_j+\frac{h}{2}\Big(f(t_j,y_j)+f(t_{j+1},y_j+hf(t_j,y_j))\Big).

Auch das verbesserte Euler-Verfahren ist explizit und von Ordnung 2. Es verwendet eine angenäherte Steigung in der Schrittmitte:

y_{j+1}=y_j+hf\Big(t_j+\frac{h}{2},y_j+\frac{h}{2}f(t_j,y_j)\Big).

Fehler, Konsistenz und Konvergenz

Der globale Fehler an einer Stelle t_j ist |y_j-y(t_j)|; eine Vektornorm misst dabei den Abstand von Näherung und exakter Lösung. Mit h als größter verwendeter Schrittweite besitzt ein Verfahren Konvergenzordnung p\geq1, wenn für hinreichend kleines h gilt:

\max_{j=0,\dotsc,n}|y_j-y(t_j)|\leq Ch^p,

wobei C>0 nicht von h abhängt. Eine höhere Ordnung bedeutet im Allgemeinen einen kleineren Gesamtfehler bei gleicher Schrittweite. Bei p=1 halbiert sich der Fehler bei Halbierung von h ungefähr; bei p=4 verringert er sich ungefähr um den Faktor (1/2)^4=1/16.

Der lokale Abschneidefehler betrachtet nur einen einzelnen Schritt, der an der exakten Lösung startet. Für ein explizites Verfahren mit konstanter Schrittweite h ist er

\eta(t,h)=y(t)+h\Phi(t,y(t),h)-y(t+h).

Ein Verfahren hat Konsistenzordnung p, wenn |\eta(t,h)|\leq Ch^{p+1} gilt. Der globale Fehler enthält zusätzlich die Fehler aus früheren Schritten. Ist \Phi Lipschitz-stetig, folgt aus Konsistenzordnung p auch Konvergenzordnung p. Vereinfacht gilt in der allgemeinen Theorie: Ein Verfahren ist genau dann konvergent, wenn es konsistent und stabil ist.

Steifheit und A-Stabilität

Die Konvergenzordnung beschreibt nur den Grenzfall h gegen null. Bei einer festen Schrittweite kann ein Verfahren dennoch ungeeignet sein, besonders bei steifen Anfangswertproblemen. Diese haben meist Komponenten, die sehr schnell konstant werden, während andere sich langsam ändern. Praktisch heißt steif: Ein explizites Einschrittverfahren müsste eine unbrauchbar kleine Schrittweite verwenden, um eine brauchbare Lösung zu liefern; dafür werden implizite Verfahren benötigt.

Zur Untersuchung dient die Testgleichung

y'(t)=\lambda y(t),\quad y(0)=1.

Für \lambda<0 ist die exakte Lösung y(t)=e^{\lambda t} exponentiell fallend. Bei einer zu großen Schrittweite können explizite Verfahren stark oszillierende und anwachsende Näherungswerte erzeugen. Implizite Verfahren liefern typischerweise auch bei beliebigen Schrittweiten qualitativ richtig fallende Werte.

Ein Verfahren heißt A-stabil, wenn es für jedes h>0 bei der Testgleichung für alle komplexen \lambda mit \operatorname{Re}(\lambda)\leq0 eine beschränkte Folge y_0,y_1,y_2,\dotsc erzeugt. Das implizite Euler- und das implizite Trapez-Verfahren sind A-stabil. Kein explizites Verfahren kann A-stabil sein.

Runge-Kutta, Extrapolation und variable Schrittweiten

Runge-Kutta-Verfahren bilden die wichtigste Klasse der Einschrittverfahren. Sie berechnen pro Schritt mehrere Hilfssteigungen k_1,\dotsc,k_s und verwenden anschließend ein gewichtetes Mittel als Verfahrenssteigung. Bei expliziten Runge-Kutta-Verfahren werden die Hilfssteigungen nacheinander direkt berechnet; bei impliziten entstehen sie aus einem Gleichungssystem.

Das klassische explizite Runge-Kutta-Verfahren der Ordnung 4 verwendet

k_1=f(t_j,y_j), k_2=f(t_j+\frac{h}{2},y_j+\frac{h}{2}k_1), k_3=f(t_j+\frac{h}{2},y_j+\frac{h}{2}k_2), k_4=f(t_j+h,y_j+hk_3)

und die Steigung \frac16k_1+\frac13k_2+\frac13k_3+\frac16k_4. Für explizite Verfahren der Ordnung 5 sind mindestens sechs Stufen nötig, für Ordnung 8 mindestens 11. Bei impliziten Runge-Kutta-Verfahren gibt es für jede Stufenzahl s ein Verfahren maximaler Ordnung p=2s.

Extrapolation kombiniert Ergebnisse bei mehreren Schrittweiten. Für ein Verfahren der Ordnung p werden zwei Näherungen y_{h_1} und y_{h_2} mit h_2=\frac12h_1 genutzt. Der extrapolierte Wert lautet

y_{h_2}+\frac{y_{h_2}-y_{h_1}}{2^p-1}.

Er ist im Allgemeinen genauer; die so erhaltene Ordnung ist mindestens p+1.

Bei der Schrittweitensteuerung sollen Fehlertoleranzen eingehalten und zugleich möglichst wenige Schritte verwendet werden. Kleine Schritte sind dort nötig, wo sich die Lösung stark ändert, große bei nahezu konstantem Verlauf. Bei der Schrittweitenhalbierung wird ein Schritt mit zwei Halb-Schritten verglichen; ist die Fehlerschätzung zu groß, wird der Schritt verworfen und mit kleinerem h wiederholt. Moderne Programme verwenden meist eingebettete Verfahren: Zwei Verfahren verschiedener Ordnung teilen möglichst viele Berechnungen. Eingebettete Runge-Kutta-Verfahren nutzen dieselben Hilfssteigungen, mitteln sie aber unterschiedlich.

Beispiel mit numerischer Software

Als Beispiel werden die Lotka-Volterra-Gleichungen für Beute- und Räuberpopulationen verwendet:

y_1'(t)=ay_1(t)-by_1(t)y_2(t), y_2'(t)=cy_1(t)y_2(t)-dy_2(t).

Für a=1, b=2, c=1, d=1 sowie y_1(0)=3 und y_2(0)=1 wird die Lösung auf dem Intervall [0,20] berechnet. Dabei beschreibt y_1 die Beute- und y_2 die Räuberpopulation.

In Matlab berechnet ode45 die Lösung. Die Funktion verwendet zwei eingebettete explizite Runge-Kutta-Verfahren mit den Konvergenzordnungen 4 und 5 zur Schrittweitensteuerung. Das Beispiel kann auch mit GNU Octave ausgeführt werden; dort ergibt sich jedoch eine etwas andere Folge verwendeter Schrittweiten.