Gewöhnliche Differentialgleichungen

Die Entwicklung der Differential- und Integralrechnung durch den englischen Physiker und Mathematiker Isaac Newton und unabhängig davon durch den deutschen Universalgelehrten Gottfried Wilhelm Leibniz im letzten Drittel des 17. Jahrhunderts war ein wesentlicher Impuls für die Mathematisierung der Wissenschaft in der frühen Neuzeit. Diese Methoden bildeten den Startpunkt des mathematischen Teilgebiets der Analysis und sind in allen Natur- und Ingenieurwissenschaften von zentraler Bedeutung. Während Leibniz von dem geometrischen Problem, Tangenten an gegebene Kurven zu bestimmen, zur Differentialrechnung geführt wurde, ging Newton von der Fragestellung aus, wie sich Änderungen einer physikalischen Größe zu einem bestimmten Zeitpunkt bestimmen lassen.

Zum Beispiel ergibt sich bei der Bewegung eines Körpers dessen Durchschnittsgeschwindigkeit einfach als die zurückgelegte Strecke dividiert durch die dafür benötigte Zeit. Um jedoch die Momentangeschwindigkeit ${\displaystyle v(t)}$ des Körpers zu einem bestimmten Zeitpunkt ${\displaystyle t}$ mathematisch zu formulieren, ist ein Grenzübergang notwendig: Man betrachtet kurze Zeitspannen der Länge ${\displaystyle \Delta t}$, die dabei zurückgelegten Wegstrecken ${\displaystyle \Delta x}$ und die zugehörigen Durchschnittsgeschwindigkeiten ${\displaystyle {\tfrac {\Delta x}{\Delta t}}}$. Wenn man nun die Zeitspanne ${\displaystyle \Delta t}$ gegen null konvergieren lässt und wenn sich dabei die Durchschnittsgeschwindigkeiten ebenfalls einem festen Wert annähern, dann wird dieser Wert die (Momentan-)Geschwindigkeit ${\displaystyle v(t)}$ zu dem gegebenen Zeitpunkt ${\displaystyle t}$ genannt. Bezeichnet ${\displaystyle x(t)}$ die Position des Körpers zur Zeit ${\displaystyle t}$, dann schreibt man ${\displaystyle v(t)=x'(t)}$ und nennt ${\displaystyle v}$ die Ableitung von ${\displaystyle x}$.

Der entscheidende Schritt in die Richtung der Differentialgleichungsmodelle ist nun die umgekehrte Fragestellung: Im Beispiel des bewegten Körpers sei also zu jedem Zeitpunkt ${\displaystyle t}$ die Geschwindigkeit ${\displaystyle v(t)}$ bekannt und daraus soll seine Position ${\displaystyle x(t)}$ bestimmt werden. Es ist anschaulich klar, dass zusätzlich die Anfangsposition des Körpers zu einem Zeitpunkt ${\displaystyle t_{0}}$ bekannt sein muss, um dieses Problem eindeutig lösen zu können. Es ist also eine Funktion ${\displaystyle x(t)}$ mit ${\displaystyle x'(t)=v(t)}$ gesucht, die die Anfangsbedingung ${\displaystyle x(t_{0})=x_{0}}$ mit gegebenen Werten ${\displaystyle t_{0}}$ und ${\displaystyle x_{0}}$ erfüllt.

Im Beispiel der Bestimmung der Position ${\displaystyle x}$ eines Körpers aus seiner Geschwindigkeit ist die Ableitung der gesuchten Funktion explizit gegeben. Meist liegt jedoch der wichtige allgemeine Fall gewöhnlicher Differentialgleichungen für eine gesuchte Größe ${\displaystyle y}$ vor: Aufgrund der Naturgesetze oder der Modellannahmen ist ein Funktionszusammenhang bekannt, der angibt, wie die Ableitung ${\displaystyle y'(t)}$ der zu bestimmenden Funktion aus ${\displaystyle t}$ und aus dem (unbekannten) Wert ${\displaystyle y(t)}$ berechnet werden kann. Zusätzlich muss wieder eine Anfangsbedingung gegeben sein, die beispielsweise aus einer Messung der gesuchten Größe zu einem fest gewählten Zeitpunkt erhalten werden kann. Zusammengefasst liegt also der folgende allgemeine Aufgabentyp vor: Man finde die Funktion ${\displaystyle y}$, die die Gleichungen

${\displaystyle y'(t)=f(t,y(t)),\quad y(t_{0})=y_{0}}$

erfüllt, wobei ${\displaystyle f}$ eine gegebene Funktion ist.

Ein einfaches Beispiel ist eine Größe ${\displaystyle y}$, die exponentiell wächst. Das bedeutet, dass die momentane Änderung, also die Ableitung ${\displaystyle y'(t)}$ proportional zu ${\displaystyle y(t)}$ selbst ist. Es gilt also ${\displaystyle y'(t)=\lambda y(t)}$ mit einer Wachstumsrate ${\displaystyle \lambda }$ und beispielsweise einer Anfangsbedingung ${\displaystyle y(0)=y_{0}}$. Die gesuchte Lösung ${\displaystyle y}$ lässt sich in diesem Fall bereits mit elementarer Differentialrechnung finden und mithilfe der Exponentialfunktion angeben: Es gilt ${\displaystyle y(t)=y_{0}e^{\lambda t}}$.

Die gesuchte Funktion ${\displaystyle y}$ in einer Differentialgleichung kann vektorwertig sein, das heißt, für jedes ${\displaystyle t}$ kann ${\displaystyle y(t)=(y_{1}(t),\dotsc ,y_{d}(t))}$ ein Vektor mit ${\displaystyle d}$ Komponenten sein. Man spricht dann auch von einem ${\displaystyle d}$-dimensionalen Differentialgleichungssystem. Im Anschauungsfall eines bewegten Körpers ist dann ${\displaystyle y(t)}$ seine Position im ${\displaystyle d}$-dimensionalen euklidischen Raum und ${\displaystyle y'(t)}$ seine Geschwindigkeit zur Zeit ${\displaystyle t}$. Durch die Differentialgleichung ist also in jedem Zeit- und Raumpunkt die Geschwindigkeit der gesuchten Bahnkurve mit Richtung und Betrag vorgegeben. Daraus soll die Bahn selbst berechnet werden.

Grundidee der Einschrittverfahren

Bei der oben als Beispiel betrachteten einfachen Differentialgleichung des exponentiellen Wachstums ließ sich die Lösungsfunktion direkt angeben. Das ist bei komplizierteren Problemen im Allgemeinen nicht mehr möglich. Man kann dann zwar unter gewissen Zusatzvoraussetzungen an die Funktion ${\displaystyle f}$ zeigen, dass eine eindeutig bestimmte Lösung des Anfangswertproblems existiert; diese kann aber dann nicht mehr durch Lösungsverfahren der Analysis (wie beispielsweise Trennung der Variablen, einen Exponentialansatz oder Variation der Konstanten) explizit berechnet werden. In diesem Fall können numerische Verfahren verwendet werden, um Näherungen für die gesuchte Lösung zu bestimmen.

Die Verfahren zur numerischen Lösung von Anfangswertproblemen gewöhnlicher Differentialgleichungen lassen sich grob in zwei große Gruppen einteilen: die Einschritt- und die Mehrschrittverfahren. Beiden Gruppen ist gemeinsam, dass sie schrittweise Näherungen ${\displaystyle y_{0},y_{1},y_{2},\dotsc }$ für die gesuchten Funktionswerte ${\displaystyle y(t_{0}),y(t_{1}),y(t_{2}),\dotsc }$ an Stellen ${\displaystyle t_{0}<t_{1}<t_{2}<\ldots }$ berechnen. Die definierende Eigenschaft der Einschrittverfahren ist dabei, dass zur Bestimmung der folgenden Näherung ${\displaystyle y_{j+1}}$ nur die „aktuelle“ Näherung ${\displaystyle y_{j}}$ verwendet wird. Bei Mehrschrittverfahren gehen im Gegensatz dazu zusätzlich bereits zuvor berechnete Näherungen mit ein; ein Dreischrittverfahren würde also beispielsweise außer ${\displaystyle y_{j}}$ auch noch ${\displaystyle y_{j-1}}$ und ${\displaystyle y_{j-2}}$ zur Bestimmung der neuen Näherung ${\displaystyle y_{j+1}}$ verwenden.

Das einfachste und grundlegendste Einschrittverfahren ist das explizite Euler-Verfahren, das der Schweizer Mathematiker und Physiker Leonhard Euler 1768 in seinem Lehrbuch Institutiones Calculi Integralis vorstellte. Die Idee dieser Methode ist es, die gesuchte Lösung durch eine stückweise lineare Funktion anzunähern, bei der in jedem Schritt von der Stelle ${\displaystyle t_{j}}$ zur Stelle ${\displaystyle t_{j+1}}$ die Steigung des Geradenstücks durch ${\displaystyle f(t_{j},y_{j})}$ gegeben ist. Genauer betrachtet: Durch die Problemstellung ist bereits ein Wert der gesuchten Funktion gegeben, nämlich ${\displaystyle y(t_{0})=y_{0}}$. Aber auch die Ableitung an dieser Stelle ist bekannt, denn es gilt ${\displaystyle y'(t_{0})=f(t_{0},y_{0})}$. Damit kann die Tangente an den Graphen der Lösungsfunktion bestimmt und als Näherung verwendet werden. An der Stelle ${\displaystyle t_{1}>t_{0}}$ ergibt sich mit der Schrittweite ${\displaystyle h_{0}:=t_{1}-t_{0}}$

${\displaystyle y(t_{1})\approx y_{0}+h_{0}f(t_{0},y_{0})=:y_{1}}$.

Dieses Vorgehen kann nun in den folgenden Schritten fortgesetzt werden. Insgesamt ergibt sich damit für das explizite Euler-Verfahren die Rechenvorschrift

${\displaystyle y_{j+1}=y_{j}+h_{j}f(t_{j},y_{j}),\quad j=0,1,2,\dotsc }$

mit den Schrittweiten ${\displaystyle h_{j}=t_{j+1}-t_{j}}$.

Das explizite Euler-Verfahren ist der Ausgangspunkt für zahlreiche Verallgemeinerungen, bei denen die Steigung ${\displaystyle f(t_{j},y_{j})}$, durch Steigungen ersetzt wird, die das Verhalten der Lösung zwischen den Stellen ${\displaystyle t_{j}}$ und ${\displaystyle t_{j+1}}$ genauer annähern. Eine zusätzliche Idee für Einschrittverfahren bringt das implizite Eulerverfahren, das ${\displaystyle f(t_{j+1},y_{j+1})}$ als Steigung verwendet. Diese Wahl erscheint auf den ersten Blick wenig geeignet, da ja ${\displaystyle y_{j+1}}$ unbekannt ist. Als Verfahrensschritt erhält man aber nun die Gleichung

${\displaystyle y_{j+1}=y_{j}+h_{j}f(t_{j+1},y_{j+1})}$

aus der ${\displaystyle y_{j+1}}$ (gegebenenfalls durch ein numerisches Verfahren) berechnet werden kann. Wählt man beispielsweise als Steigung das arithmetische Mittel aus den Steigungen des expliziten und des impliziten Euler-Verfahrens, so erhält man das implizite Trapez-Verfahren. Aus diesem lässt sich wiederum ein explizites Verfahren gewinnen, wenn man zum Beispiel das unbekannte ${\displaystyle y_{j+1}}$ auf der rechten Seite der Gleichung durch die Näherung des expliziten Euler-Verfahrens nähert, das sogenannte Heun-Verfahren. All diesen Verfahren und allen weiteren Verallgemeinerungen ist die Grundidee der Einschrittverfahren gemeinsam: der Schritt

${\displaystyle y_{j+1}=y_{j}+h_{j}\Phi }$

mit einer Steigung ${\displaystyle \Phi }$, die von ${\displaystyle t_{j}}$, ${\displaystyle y_{j}}$ und ${\displaystyle h_{j}}$ sowie (bei impliziten Verfahren) von ${\displaystyle y_{j+1}}$ abhängen kann.

Konvergenzordnung

Für ein praxistaugliches Einschrittverfahren sollen die berechneten ${\displaystyle y_{j}}$ gute Näherungen für die Werte ${\displaystyle y(t_{j})}$ der exakten Lösung ${\displaystyle y}$ an der Stelle ${\displaystyle t_{j}}$ sein. Da im Allgemeinen die Größen ${\displaystyle d}$-dimensionale Vektoren sind, misst man die Güte dieser Näherung mit einer Vektornorm als ${\displaystyle \|y_{j}-y(t_{j})\|}$, dem Fehler an der Stelle ${\displaystyle t_{j}}$. Es ist wünschenswert, dass diese Fehler für alle ${\displaystyle j}$ schnell gegen null konvergieren, falls man die Schrittweiten gegen null konvergieren lässt. Um auch den Fall nicht konstanter Schrittweiten zu erfassen, definiert man dazu genauer ${\displaystyle h}$ als das Maximum der verwendeten Schrittweiten und betrachtet das Verhalten des maximalen Fehlers an allen Stellen ${\displaystyle t_{j}}$ im Vergleich zu Potenzen von ${\displaystyle h}$. Man sagt, das Einschrittverfahren zur Lösung des gegebenen Anfangswertproblems habe die Konvergenzordnung ${\displaystyle p\geq 1}$, wenn die Abschätzung

${\displaystyle \max _{j=0,\dotsc ,n}\|y_{j}-y(t_{j})\|\leq Ch^{p}}$

für alle hinreichend kleinen ${\displaystyle h}$ mit einer von ${\displaystyle h}$ unabhängigen Konstante ${\displaystyle C>0}$ gilt.

Die Konvergenzordnung ist die wichtigste Kenngröße für den Vergleich verschiedener Einschrittverfahren. Ein Verfahren mit höherer Konvergenzordnung ${\displaystyle p}$ liefert im Allgemeinen bei vorgegebener Schrittweite einen kleineren Gesamtfehler bzw. umgekehrt sind weniger Schritte nötig, um eine vorgegebene Genauigkeit zu erreichen. Bei einem Verfahren mit ${\displaystyle p=1}$ ist zu erwarten, dass sich bei einer Halbierung der Schrittweite auch der Fehler nur ungefähr halbiert. Bei einem Verfahren der Konvergenzordnung ${\displaystyle p=4}$ kann man hingegen davon ausgehen, dass sich dabei der Fehler ungefähr um den Faktor ${\displaystyle {\bigl (}{\tfrac {1}{2}}{\bigr )}^{4}={\tfrac {1}{16}}}$ verringert.

Globaler und lokaler Fehler

Die in der Definition der Konvergenzordnung betrachteten Fehler ${\displaystyle \|y_{j}-y(t_{j})\|}$ setzen sich auf eine zunächst kompliziert erscheinende Weise aus zwei Einzelkomponenten zusammen: Natürlich hängen sie zum einen von dem Fehler ab, den das Verfahren in einem einzelnen Schritt macht, indem es die unbekannte Steigung der gesuchten Funktion durch die Verfahrensfunktion annähert. Zum anderen ist aber zusätzlich zu berücksichtigen, dass bereits der Startpunkt ${\displaystyle (t_{j},y_{j})}$ eines Schrittes im Allgemeinen nicht mit dem exakten Startpunkt ${\displaystyle (t_{j},y(t_{j}))}$ übereinstimmt; der Fehler nach diesem Schritt hängt also auch von allen Fehlern ab, die bereits in den vorangegangenen Schritten gemacht wurden. Aufgrund der einheitlichen Definition der Einschrittverfahren, die sich nur in der Wahl der Verfahrensfunktion ${\displaystyle \Phi }$ unterscheiden, lässt sich aber beweisen, dass man (unter gewissen technischen Voraussetzungen an ${\displaystyle \Phi }$) direkt von der Fehlerordnung in einem einzelnen Schritt, der sogenannten Konsistenzordnung, auf die Konvergenzordnung schließen kann.

Der Begriff der Konsistenz ist ein allgemeines und zentrales Konzept der modernen numerischen Mathematik. Während man bei der Konvergenz eines Verfahrens untersucht, wie gut die numerischen Näherungen zur exakten Lösung passen, stellt man sich bei der Konsistenz vereinfacht gesprochen die „umgekehrte“ Frage: Wie gut erfüllt die exakte Lösung die Verfahrensvorschrift? In dieser allgemeinen Theorie gilt, dass ein Verfahren genau dann konvergent ist, wenn es konsistent und stabil ist. Um die Notation zu vereinfachen, soll in der folgenden Überlegung angenommen werden, dass ein explizites Einschrittverfahren

${\displaystyle y_{j+1}=y_{j}+h\Phi (t_{j},y_{j},h)}$

mit konstanter Schrittweite ${\displaystyle h}$ vorliegt. Mit der wahren Lösung ${\displaystyle t\mapsto y(t)}$ definiert man den lokalen Abschneidefehler (auch lokaler Verfahrensfehler genannt) ${\displaystyle \eta }$ als

${\displaystyle \eta (t,h)=y(t)+h\Phi (t,y(t),h)-y(t+h)}$.

Man nimmt also an, dass die exakte Lösung bekannt ist, startet einen Verfahrensschritt an dem Punkt ${\displaystyle (t,y(t))}$ und bildet die Differenz zur exakten Lösung an der Stelle ${\displaystyle t+h}$. Damit definiert man: Ein Einschrittverfahren hat die Konsistenzordnung ${\displaystyle p\geq 1}$, wenn die Abschätzung

${\displaystyle \|\eta (t,h)\|\leq Ch^{p+1}}$

für alle hinreichend kleinen ${\displaystyle h}$ mit einer von ${\displaystyle h}$ unabhängigen Konstante ${\displaystyle C>0}$ gilt.

Der auffällige Unterschied zwischen den Definitionen der Konsistenzordnung und der Konvergenzordnung ist die Potenz ${\displaystyle h^{p+1}}$ anstelle von ${\displaystyle h^{p}}$. Das lässt sich anschaulich so deuten, dass beim Übergang vom lokalen zum globalen Fehler eine Potenz der Schrittweite „verloren geht“. Es gilt nämlich der folgende, für die Theorie der Einschrittverfahren zentrale Satz:

Ist die Verfahrensfunktion ${\displaystyle \Phi }$ Lipschitz-stetig und hat das zugehörige Einschrittverfahren die Konsistenzordnung ${\displaystyle p}$, dann hat es auch die Konvergenzordnung ${\displaystyle p}$.

Die Lipschitz-Stetigkeit der Verfahrensfunktion als Zusatzvoraussetzung für die Stabilität ist im Allgemeinen immer dann erfüllt, wenn die Funktion ${\displaystyle f}$ aus der Differentialgleichung selbst Lipschitz-stetig ist. Diese Forderung muss für die meisten Anwendungen sowieso vorausgesetzt werden, um die eindeutige Lösbarkeit des Anfangswertproblems zu garantieren. Nach dem Satz genügt es also, die Konsistenzordnung eines Einschrittverfahrens zu bestimmen. Das lässt sich prinzipiell durch Taylor-Entwicklung von ${\displaystyle \eta (t,h)}$ nach Potenzen von ${\displaystyle h}$ erreichen. In der Praxis werden die entstehenden Formeln für höhere Ordnungen sehr kompliziert und unübersichtlich, sodass zusätzliche Konzepte und Notationen benötigt werden.

Einfache Verfahren der Ordnung 1 und 2

Wie der französische Mathematiker Augustin-Louis Cauchy um 1820 bewies, besitzt das Euler-Verfahren die Konvergenzordnung 1. Wenn man die Steigungen ${\displaystyle f(t_{j},y_{j})}$ des expliziten Euler-Verfahrens und ${\displaystyle f(t_{j+1},y_{j+1})}$ des impliziten Euler-Verfahrens, wie sie an den beiden Endpunkten eines Schritts vorliegen, mittelt, kann man hoffen, eine bessere Näherung über das ganze Intervall zu erhalten. Tatsächlich lässt sich beweisen, dass das so erhaltene implizite Trapez-Verfahren

${\displaystyle y_{j+1}=y_{j}+{\frac {h}{2}}{\Big (}f(t_{j},y_{j})+f(t_{j+1},y_{j+1}){\Big )}}$

die Konvergenzordnung 2 hat. Dieses Verfahren weist sehr gute Stabilitätseigenschaften auf, ist allerdings implizit, sodass in jedem Schritt eine Gleichung für ${\displaystyle y_{j+1}}$ gelöst werden muss. Nähert man diese Größe auf der rechten Seite der Gleichung durch das explizite Euler-Verfahren an, so entsteht das explizite Verfahren von Heun

${\displaystyle y_{j+1}=y_{j}+{\frac {h}{2}}{\Big (}f(t_{j},y_{j})+f{\big (}t_{j+1},y_{j}+hf(t_{j},y_{j}){\big )}{\Big )}}$,

das ebenfalls die Konvergenzordnung 2 besitzt. Ein weiteres einfaches explizites Verfahren der Ordnung 2, das verbesserte Euler-Verfahren, erhält man durch folgende Überlegung: Eine „mittlere“ Steigung im Verfahrensschritt wäre die Steigung der Lösung ${\displaystyle y}$ in der Mitte des Schritts, also an der Stelle ${\displaystyle t_{j}+{\tfrac {h}{2}}}$. Da die Lösung aber unbekannt ist, nähert man sie durch einen expliziten Euler-Schritt mit halber Schrittweite an. Es ergibt sich die Verfahrensvorschrift

${\displaystyle y_{j+1}=y_{j}+hf{\big (}t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}f(t_{j},y_{j}){\big )}}$.

Diese Einschrittverfahren der Ordnung 2 wurden als Verbesserungen des Euler-Verfahrens alle 1895 von dem deutschen Mathematiker Carl Runge veröffentlicht.

Runge-Kutta-Verfahren

Die erwähnten Ideen für einfache Einschrittverfahren führen bei weiterer Verallgemeinerung zur wichtigen Klasse der Runge-Kutta-Verfahren. Zum Beispiel lässt sich das Verfahren von Heun übersichtlicher so präsentieren: Zuerst wird eine Hilfssteigung ${\displaystyle k_{1}=f(t_{j},y_{j})}$ berechnet, nämlich die Steigung des expliziten Euler-Verfahrens. Damit wird eine weitere Hilfssteigung bestimmt, hier ${\displaystyle k_{2}=f(t_{j}+h,y_{j}+hk_{1})}$. Die tatsächlich verwendete Verfahrenssteigung ${\displaystyle \Phi }$ ergibt sich anschließend als ein gewichtetes Mittel der Hilfssteigungen, im Verfahren von Heun also ${\displaystyle {\tfrac {1}{2}}k_{1}+{\tfrac {1}{2}}k_{2}}$. Dieses Vorgehen lässt sich auf mehr als zwei Hilfssteigungen verallgemeinern. Ein ${\displaystyle s}$-stufiges Runge-Kutta-Verfahren berechnet zunächst Hilfssteigungen ${\displaystyle k_{1},\dotsc ,k_{s}}$ durch Auswertung von ${\displaystyle f}$ an geeigneten Stellen und anschließend ${\displaystyle \Phi }$ als gewichtetes Mittel. Bei einem expliziten Runge-Kutta-Verfahren werden die Hilfssteigungen ${\displaystyle k_{1},k_{2},k_{3},\dotsc }$ der Reihe nach direkt berechnet, bei einem impliziten ergeben sie sich als Lösungen eines Gleichungssystems. Als typisches Beispiel sei das explizite klassische Runge-Kutta-Verfahren der Ordnung 4 angeführt, das mitunter einfach als das Runge-Kutta-Verfahren bezeichnet wird: Dabei werden zunächst die vier Hilfssteigungen

${\displaystyle {\begin{aligned}k_{1}&=f(t_{j},y_{j})\\k_{2}&=f(t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}k_{1})\\k_{3}&=f(t_{j}+{\tfrac {h}{2}},y_{j}+{\tfrac {h}{2}}k_{2})\\k_{4}&=f(t_{j}+h,y_{j}+hk_{3})\end{aligned}}}$

berechnet und dann als Verfahrenssteigung das gewichtete Mittel

${\displaystyle {\tfrac {1}{6}}k_{1}+{\tfrac {1}{3}}k_{2}+{\tfrac {1}{3}}k_{3}+{\tfrac {1}{6}}k_{4}}$

verwendet. Dieses bekannte Verfahren veröffentlichte der deutsche Mathematiker Wilhelm Kutta im Jahr 1901, nachdem ein Jahr zuvor Karl Heun ein dreistufiges Einschrittverfahren der Ordnung 3 gefunden hatte.

Die Konstruktion von expliziten Verfahren noch höherer Ordnung mit möglichst kleiner Stufenzahl ist ein mathematisch recht anspruchsvolles Problem. Wie John C. Butcher 1965 zeigen konnte, gibt es zum Beispiel für die Ordnung 5 nur minimal sechsstufige Verfahren; ein explizites Runge-Kutta-Verfahren der Ordnung 8 benötigt mindestens 11 Stufen. 1978 fand der österreichische Mathematiker Ernst Hairer ein Verfahren der Ordnung 10 mit 17 Stufen. Die Koeffizienten für ein solches Verfahren müssen 1205 Bestimmungsgleichungen erfüllen. Bei impliziten Runge-Kutta-Verfahren ist die Situation einfacher und übersichtlicher: Für jede Stufenzahl ${\displaystyle s}$ existiert ein Verfahren der Ordnung ${\displaystyle p=2s}$; das ist zugleich die maximal erreichbare Ordnung.

Extrapolationsverfahren

Die Idee der Extrapolation ist nicht auf die Lösung von Anfangswertproblemen mit Einschrittverfahren beschränkt, sondern lässt sich sinngemäß auf alle numerischen Verfahren anwenden, die das zu lösende Problem mit einer Schrittweite ${\displaystyle h}$ diskretisieren. Ein bekanntes Beispiel eines Extrapolationsverfahrens ist etwa die Romberg-Integration zur numerischen Berechnung von Integralen. Sei daher ganz allgemein ${\displaystyle v}$ ein Wert, der numerisch bestimmt werden soll, im Fall dieses Artikels etwa der Wert der Lösungsfunktion eines Anfangswertproblems an einer gegebenen Stelle. Ein numerisches Verfahren, beispielsweise ein Einschrittverfahren, berechnet dafür einen Näherungswert ${\displaystyle {\tilde {v}}(h)}$, der von der Wahl der Schrittweite ${\displaystyle h>0}$ abhängt. Dabei sei angenommen, dass das Verfahren konvergent ist, also, dass ${\displaystyle {\tilde {v}}(h)}$ gegen ${\displaystyle v}$ konvergiert, wenn ${\displaystyle h}$ gegen null konvergiert. Diese Konvergenz ist jedoch nur eine rein theoretische Aussage, da bei der realen Anwendung des Verfahrens zwar Näherungswerte ${\displaystyle {\tilde {v}}(h_{1}),{\tilde {v}}(h_{2}),\dotsc ,{\tilde {v}}(h_{m})}$ für endlich viele verschiedene Schrittweiten ${\displaystyle h_{1}>h_{2}>\ldots >h_{m}}$ berechnet werden können, man aber selbstverständlich nicht die Schrittweite „gegen null konvergieren“ lassen kann. Die berechneten Näherungen für verschiedene Schrittweiten lassen sich jedoch als Information über die (unbekannte) Funktion ${\displaystyle {\tilde {v}}}$ auffassen: Bei den Extrapolationsverfahren wird dabei ${\displaystyle {\tilde {v}}}$ durch ein Interpolationspolynom angenähert, also durch ein Polynom ${\displaystyle P}$ mit

${\displaystyle P(h_{k})={\tilde {v}}(h_{k})}$

für ${\displaystyle k=1,2,\dotsc ,m}$. Der Wert ${\displaystyle P(0)}$ des Polynoms an der Stelle ${\displaystyle h=0}$ wird dann als berechenbare Näherung für den nicht berechenbaren Grenzwert von ${\displaystyle {\tilde {v}}(h)}$ für ${\displaystyle h}$ gegen null verwendet. Einen frühen erfolgreichen Extrapolationsalgorithmus für Anfangswertprobleme veröffentlichten Roland Bulirsch und Josef Stoer im Jahr 1966.

Ein konkretes Beispiel im Falle eines Einschrittverfahrens der Ordnung ${\displaystyle p}$ kann das allgemeine Vorgehen der Extrapolation verständlich machen. Bei einem solchen Verfahren lässt sich die berechnete Näherung für kleine Schrittweiten ${\displaystyle h}$ gut durch ein Polynom der Form

${\displaystyle P(h)=a+bh^{p}}$

mit zunächst unbekannten Parametern ${\displaystyle a}$ und ${\displaystyle b}$ annähern. Berechnet man nun mit dem Verfahren für eine Schrittweite ${\displaystyle h_{1}}$ und für die halbe Schrittweite ${\displaystyle h_{2}={\tfrac {1}{2}}h_{1}}$ zwei Näherungen ${\displaystyle y_{h_{1}}}$ und ${\displaystyle y_{h_{2}}}$, erhält man aus den Interpolationsbedingungen ${\displaystyle P(h_{1})=y_{h_{1}}}$ und ${\displaystyle P(h_{2})=y_{h_{2}}}$ zwei lineare Gleichungen für die Unbekannten ${\displaystyle a}$ und ${\displaystyle b}$. Der auf ${\displaystyle h=0}$ extrapolierte Wert

${\displaystyle P(0)=a=y_{h_{2}}+{\frac {y_{h_{2}}-y_{h_{1}}}{2^{p}-1}}}$

stellt dann im Allgemeinen eine deutlich bessere Näherung dar als die beiden zunächst berechneten Werte. Es lässt sich zeigen, dass die Ordnung des so erhaltenen Einschrittverfahrens mindestens ${\displaystyle p+1}$ ist, also um mindestens 1 größer als beim ursprünglichen Verfahren ist.

Verfahren mit Schrittweitensteuerung

Ein Vorteil der Einschrittverfahren ist, dass in jedem Schritt ${\displaystyle j}$ unabhängig von den anderen Schritten eine beliebige Schrittweite ${\displaystyle h_{j}}$ verwendet werden kann. In der Praxis stellt sich dabei offensichtlich die Frage, wie ${\displaystyle h_{j}}$ gewählt werden soll. In realen Anwendungen wird stets eine Fehlertoleranz gegeben sein, mit der die Lösung eines Anfangswertproblems berechnet werden soll; zum Beispiel wäre es sinnlos, eine numerische Näherung zu bestimmen, die wesentlich „genauer“ ist als die mit Messabweichungen behafteten Daten für Startwerte und Parameter des gegebenen Problems. Das Ziel wird also sein, die Schrittweiten so zu wählen, dass einerseits die vorgegebenen Fehlertoleranzen eingehalten werden, andererseits aber auch möglichst wenige Schritte zu verwenden, um den Rechenaufwand klein zu halten. Dies wird sich im Allgemeinen nur dann erreichen lassen, wenn die Schrittweiten an den Verlauf der Lösung angepasst werden: kleine Schritte, wo sich die Lösung stark ändert, große Schritte, wo sie nahezu konstant verläuft.

Bei gut konditionierten Anfangswertproblemen lässt sich zeigen, dass der globale Verfahrensfehler ungefähr gleich der Summe der lokalen Abschneidefehler ${\displaystyle \eta _{j}:=\|\eta (t_{j},h_{j})\|}$ in den einzelnen Schritten ist. Daher sollte als Schrittweite ein möglichst großes ${\displaystyle h_{j}}$ gewählt werden, für das ${\displaystyle \eta _{j}}$ unter einer gewählten Toleranzschwelle liegt. Das Problem dabei ist, dass ${\displaystyle \eta _{j}}$ nicht direkt berechnet werden kann, da es ja von der unbekannten exakten Lösung ${\displaystyle y(t_{j})}$ des Anfangswertproblems an der Stelle ${\displaystyle t_{j}}$ abhängt. Die Grundidee der Schrittweitensteuerung ist es daher, ${\displaystyle y(t_{j})}$ mit einem Verfahren anzunähern, das genauer ist als das zugrundeliegende Basisverfahren.

Zwei grundlegende Ideen zur Schrittweitensteuerung sind die Schrittweitenhalbierung und die eingebetteten Verfahren. Bei der Schrittweitenhalbierung wird zusätzlich zum eigentlichen Verfahrensschritt als Vergleichswert das Ergebnis für zwei Schritte mit der halben Schrittweite berechnet. Aus beiden Werten wird dann durch Extrapolation eine genauere Näherung für ${\displaystyle y(t_{j})}$ bestimmt und damit der lokale Fehler ${\displaystyle \eta _{j}}$ geschätzt. Ist dieser zu groß, wird dieser Schritt verworfen und mit einer kleineren Schrittweite wiederholt. Ist er deutlich kleiner als die vorgegebene Toleranz, kann im nächsten Schritt die Schrittweite vergrößert werden. Der zusätzliche Rechenaufwand für dieses Schrittweitenhalbierungsverfahren ist relativ groß; deshalb verwenden moderne Implementierungen meist sogenannte eingebettete Verfahren zur Schrittweitensteuerung. Die Grundidee ist dabei, in jedem Schritt zwei Näherungen für ${\displaystyle y(t_{j})}$ mit zwei Einschrittverfahren zu berechnen, die verschiedene Konvergenzordnungen haben, und damit den lokalen Fehler zu schätzen. Um den Rechenaufwand zu optimieren, sollten die beiden Verfahren möglichst viele Rechenschritte gemeinsam haben: Sie sollten „ineinander eingebettet“ sein. Eingebettete Runge-Kutta-Verfahren verwenden beispielsweise die gleichen Hilfssteigungen und unterscheiden sich nur darin, wie sie diese mitteln. Bekannte eingebettete Verfahren sind unter anderem die Runge-Kutta-Fehlberg-Verfahren (Erwin Fehlberg, 1969) und die Dormand-Prince-Verfahren (J. R. Dormand und P. J. Prince, 1980).