Einschrittverfahren

Einschrittverfahren sind in der numerischen Mathematik neben den Mehrschrittverfahren eine große Gruppe von Rechenverfahren zur Lösung von Anfangswertproblemen. Diese Aufgabenstellung, bei der eine gewöhnliche Differentialgleichung zusammen mit einer Startbedingung gegeben ist, spielt in allen Natur- und Ingenieurwissenschaften eine zentrale Rolle und gewinnt beispielsweise auch in den Wirtschafts- und Sozialwissenschaften immer mehr an Bedeutung. Anfangswertprobleme werden verwendet, um dynamische Vorgänge zu analysieren, zu simulieren oder vorherzusagen.

Die namensgebende Grundidee der Einschrittverfahren ist, dass sie ausgehend von dem gegebenen Anfangspunkt Schritt für Schritt entlang der gesuchten Lösung Näherungspunkte berechnen. Dabei verwenden sie jeweils nur die zuletzt bestimmte Näherung für den nächsten Schritt, im Gegensatz zu den Mehrschrittverfahren, die auch weiter zurückliegende Punkte in die Rechnung miteinbeziehen. Die Einschrittverfahren lassen sich grob in zwei Gruppen einteilen: in die expliziten Verfahren, die die neue Näherung direkt aus der alten berechnen, und in die impliziten Verfahren, bei denen dazu eine Gleichung gelöst werden muss. Letztere eignen sich auch für sogenannte steife Anfangswertprobleme.

Das einfachste und älteste Einschrittverfahren, das explizite Euler-Verfahren, wurde 1768 von Leonhard Euler veröffentlicht. Nachdem 1883 eine Gruppe von Mehrschrittverfahren vorgestellt worden war, entwickelten um 1900 Carl Runge, Karl Heun und Wilhelm Kutta deutliche Verbesserungen des eulerschen Verfahrens. Aus diesen ging die große Gruppe der Runge-Kutta-Verfahren hervor, die die wichtigste Klasse von Einschrittverfahren bildet. Weitere Entwicklungen des 20. Jahrhunderts sind beispielsweise die Idee der Extrapolation, vor allem aber Überlegungen zur Schrittweitensteuerung, also zur Wahl geeigneter Längen der einzelnen Schritte eines Verfahrens. Diese Konzepte bilden die Grundlage, um schwierige Anfangswertprobleme, wie sie in modernen Anwendungen auftreten, effizient und mit der benötigten Genauigkeit durch Computerprogramme lösen zu können.

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.^[1]

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.

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.^[2] 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}}$.^[3]

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.^[4] 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.

Mit den Überlegungen aus dem Einführungsabschnitt dieses Artikels kann der Begriff des Einschrittverfahrens wie folgt definiert werden: Gesucht sei die Lösung ${\displaystyle y}$ des Anfangswertproblems

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

Dabei werde vorausgesetzt, dass die Lösung

${\displaystyle y\colon I\to \mathbb {R} ^{d}}$

auf einem gegebenen Intervall ${\displaystyle I=[t_{0},T]}$ existiert und eindeutig bestimmt ist. Sind

${\displaystyle t_{0}<t_{1}<t_{2}<\ldots <t_{n}=T}$

Zwischenstellen im Intervall ${\displaystyle I}$ und ${\displaystyle h_{j}=t_{j+1}-t_{j}}$ die zugehörigen Schrittweiten, dann heißt das durch

${\displaystyle y_{j+1}=y_{j}+h_{j}\Phi (t_{j},y_{j},y_{j+1},h_{j})}$, ${\displaystyle \quad j=0,\dotsc ,n-1}$

gegebene Verfahren Einschrittverfahren mit Verfahrensfunktion ${\displaystyle \Phi }$. Wenn ${\displaystyle \Phi }$ nicht von ${\displaystyle y_{j+1}}$ abhängt, dann nennt man es explizites Einschrittverfahren. Anderenfalls muss in jedem Schritt ${\displaystyle j}$ eine Gleichung für ${\displaystyle y_{j+1}}$ gelöst werden, und das Verfahren wird implizit genannt.^[5]

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.^[6]

Die Konvergenzordnung ist die wichtigste Kenngröße für den Vergleich verschiedener Einschrittverfahren.^[7] 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.

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^[8]

${\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:^[9]

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.^[10]

Die Konvergenzordnung eines Verfahrens ist eine asymptotische Aussage, die das Verhalten der Näherungen beschreibt, wenn die Schrittweite gegen null konvergiert. Sie sagt jedoch nichts darüber aus, ob das Verfahren für eine gegebene feste Schrittweite auch tatsächlich eine brauchbare Näherung berechnet. Dass dies bei bestimmten Typen von Anfangswertproblemen tatsächlich ein großes Problem darstellen kann, beschrieben zuerst Charles Francis Curtiss und Joseph O. Hirschfelder 1952. Sie hatten beobachtet, dass bei manchen Differentialgleichungssystemen der chemischen Reaktionskinetik die Lösungen nicht mit expliziten numerischen Verfahren berechnet werden können, und nannten solche Anfangswertprobleme „steif“.^[11] Es existieren zahlreiche mathematische Kriterien, um für ein gegebenes Problem festzustellen, wie steif es ist. Anschaulich handelt es sich bei steifen Anfangswertproblemen meist um Differentialgleichungssysteme, bei denen einige Komponenten sehr schnell konstant werden, während andere Komponenten sich nur langsam ändern. Ein solches Verhalten tritt typischerweise bei der Modellierung chemischer Reaktionen auf. Die für die praktische Anwendung nützlichste Definition von Steifheit ist dabei jedoch: Ein Anfangswertproblem ist steif, wenn man bei seiner Lösung mit expliziten Einschrittverfahren die Schrittweite „zu klein“ wählen müsste, um eine brauchbare Lösung zu erhalten. Solche Probleme können also nur mit impliziten Verfahren gelöst werden.^[12]

Dieser Effekt lässt sich genauer darstellen, wenn man untersucht, wie die einzelnen Verfahren mit exponentiellem Zerfall zurechtkommen. Dazu betrachtet man nach dem schwedischen Mathematiker Germund Dahlquist die Testgleichung

${\displaystyle y'(t)=\lambda y(t),\quad y(0)=1}$

mit der für ${\displaystyle \lambda <0}$ exponentiell abfallenden Lösung ${\displaystyle y(t)=e^{\lambda t}}$. Die nebenstehende Grafik zeigt – exemplarisch für das explizite und das implizite Euler-Verfahren – das typische Verhalten dieser beiden Verfahrensgruppen bei diesem so einfach erscheinenden Anfangswertproblem: Verwendet man bei einem expliziten Verfahren eine zu große Schrittweite, dann ergeben sich stark oszillierende Werte, die sich im Laufe der Rechnung aufschaukeln und sich immer weiter von der exakten Lösung entfernen. Implizite Verfahren berechnen hingegen typischerweise die Lösung für beliebige Schrittweiten qualitativ richtig, nämlich als exponentiell fallende Folge von Näherungswerten.^[13]

Etwas allgemeiner betrachtet man die obige Testgleichung auch für komplexe Werte von ${\displaystyle \lambda }$. In diesem Fall sind die Lösungen Schwingungen, deren Amplitude genau dann beschränkt bleibt, wenn ${\displaystyle \operatorname {Re} (\lambda )\leq 0}$ gilt, also der Realteil von ${\displaystyle \lambda }$ kleiner oder gleich 0 ist. Damit lässt sich eine wünschenswerte Eigenschaft von Einschrittverfahren formulieren, die für steife Anfangswertprobleme eingesetzt werden sollen: die sogenannte A-Stabilität. Ein Verfahren heißt A-stabil, wenn es für beliebige Schrittweiten ${\displaystyle h>0}$ angewendet auf die Testgleichung für alle ${\displaystyle \lambda }$ mit ${\displaystyle \operatorname {Re} (\lambda )\leq 0}$ eine Folge von Näherungen ${\displaystyle y_{0},y_{1},y_{2},\dotsc }$ berechnet, die (wie die wahre Lösung) beschränkt bleibt. Das implizite Euler-Verfahren und das implizite Trapez-Verfahren sind die einfachsten Beispiele A-stabiler Einschrittverfahren. Andererseits lässt sich zeigen, dass ein explizites Verfahren niemals A-stabil sein kann.^[14][15]

Wie der französische Mathematiker Augustin-Louis Cauchy um 1820 bewies, besitzt das Euler-Verfahren die Konvergenzordnung 1.^[16] 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.^[17]

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.^[18] 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.^[19]

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.^[20] 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.^[21]

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.^[22] Einen frühen erfolgreichen Extrapolationsalgorithmus für Anfangswertprobleme veröffentlichten Roland Bulirsch und Josef Stoer im Jahr 1966.^[23]

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.^[24]

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.^[25]

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.^[26]

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.^[27] 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).^[28]

Für die in diesem Artikel überblicksartig dargestellten mathematischen Konzepte wurden zahlreiche Software-Implementierungen entwickelt, die dem Anwender die Möglichkeit geben, praktische Probleme auf einfache Weise numerisch zu lösen. Als konkretes Beispiel dazu soll nun eine Lösung der Lotka-Volterra-Gleichungen mit der verbreiteten numerischen Software Matlab berechnet werden. Die Lotka-Volterra-Gleichungen sind ein einfaches Modell aus der Biologie, das die Wechselwirkungen von Räuber- und Beutepopulationen beschreibt. Gegeben sei dazu das Differentialgleichungssystem

${\displaystyle {\begin{aligned}y_{1}'(t)&=ay_{1}(t)-by_{1}(t)y_{2}(t)\\y_{2}'(t)&=cy_{1}(t)y_{2}(t)-dy_{2}(t)\end{aligned}}}$

mit den Parametern ${\displaystyle a=1,b=2,c=1,d=1}$ und der Anfangsbedingung ${\displaystyle y_{1}(0)=3}$, ${\displaystyle y_{2}(0)=1}$. Hierbei entsprechen ${\displaystyle y_{1}}$ und ${\displaystyle y_{2}}$ der zeitlichen Entwicklung der Beute- bzw. der Räuberpopulation. Die Lösung soll auf dem Zeitintervall ${\displaystyle [0,20]}$ berechnet werden.

Für die Berechnung mithilfe von Matlab wird zunächst für die gegebenen Parameterwerte die Funktion ${\displaystyle f}$ auf der rechten Seite der Differentialgleichung ${\displaystyle y'=f(t,y)}$ definiert:

a = 1; b = 2; c = 1; d = 1;
f = @(t,y) [a*y(1) - b*y(1)*y(2); c*y(1)*y(2) - d*y(2)];

Außerdem werden das Zeitintervall und die Anfangswerte benötigt:

t_int = [0, 20];
y0 = [3; 1];

Anschließend kann die Lösung berechnet werden:

[t, y] = ode45(f, t_int, y0);

Die Matlab-Funktion ode45 implementiert ein Einschrittverfahren, das zwei eingebettete explizite Runge-Kutta-Verfahren mit den Konvergenzordnungen 4 und 5 zur Schrittweitensteuerung verwendet.^[29]

Die Lösung kann nun gezeichnet werden, ${\displaystyle y_{1}}$ als blaue Kurve und ${\displaystyle y_{2}}$ als rote; die berechneten Punkte werden durch kleine Kreise markiert:

figure(1)
plot(t, y(:,1), 'b-o', t, y(:,2), 'r-o')

Das Ergebnis ist unten im linken Bild dargestellt. Das rechte Bild zeigt die vom Verfahren verwendeten Schrittweiten und wurde erzeugt mit

figure(2)
plot(t(1:end-1), diff(t))

Lösungen der Lotka-Volterra-Gleichungen

Verwendete Schrittweiten

Dieses Beispiel kann ohne Änderungen auch mit der freien numerischen Software GNU Octave ausgeführt werden. Mit dem dort implementierten Verfahren ergibt sich allerdings eine etwas andere Schrittweitenfolge.

• John C. Butcher: Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, Chichester 2008, ISBN 978-0-470-72335-7. 
• Wolfgang Dahmen, Arnold Reusken: Numerik für Ingenieure und Naturwissenschaftler. 2. Auflage. Springer, Berlin/Heidelberg 2008, ISBN 978-3-540-76492-2, Kap. 11: Gewöhnliche Differentialgleichungen. 
• Peter Deuflhard, Folkmar Bornemann: Numerische Mathematik 2 – Gewöhnliche Differentialgleichungen. 3. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1. 
• David F. Griffiths, Desmond J. Higham: Numerical Methods for Ordinary Differential Equations – Initial Value Problems. Springer, London 2010, ISBN 978-0-85729-147-9. 
• Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, Kap. 7: Einschrittverfahren für Anfangswertprobleme. 
• Hans-Jürgen Reinhardt: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Walter de Gruyter, Berlin/Boston 2012, ISBN 978-3-11-028045-6. 
• Hans Rudolf Schwarz, Norbert Köckler: Numerische Mathematik. 8. Auflage. Vieweg+Teubner, Wiesbaden 2011, ISBN 978-3-8348-1551-4, Kap. 8: Anfangswertprobleme. 
• Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8. 

• Lars Grüne: Numerische Methoden für gewöhnliche Differentialgleichungen (Numerische Mathematik II). (PDF) 2008, abgerufen am 20. August 2018 (Vorlesungsskript der Universität Bayreuth). 
• Peter Spellucci: Numerik gewöhnlicher Differentialgleichungen. (PDF) 2007, abgerufen am 20. August 2018 (Vorlesungsskript der TU Darmstadt). 
• Hans U. Fuchs: Numerical Methods for Differential Equations. (PDF) 2007, abgerufen am 20. August 2018 (englisch, Vorlesungsskript der Zürcher Hochschule für Angewandte Wissenschaften). 
• Mathe Tutorial: Rechner für allgemeine Differentialgleichungen 1. Ordnung. In: Mathe Tutorial. Abgerufen am 20. August 2018. 

1. ↑ Thomas Sonar: 3000 Jahre Analysis. Springer, Berlin/Heidelberg 2011, ISBN 978-3-642-17203-8, S. 378–388 und 401–426. 
2. ↑ Jean-Luc Chabert u. a.: A History of Algorithms. Springer, Berlin/Heidelberg 1999, ISBN 978-3-540-63369-3, S. 374–378. 
3. ↑ Wolfgang Dahmen, Arnold Reusken: Numerik für Ingenieure und Naturwissenschaftler. 2. Auflage. Springer, Berlin/Heidelberg 2008, ISBN 978-3-540-76492-2, S. 386 f. 
4. ↑ Wolfgang Dahmen, Arnold Reusken: Numerik für Ingenieure und Naturwissenschaftler. 2. Auflage. Springer, Berlin/Heidelberg 2008, ISBN 978-3-540-76492-2, S. 386–392. 
5. ↑ Hans Rudolf Schwarz, Norbert Köckler: Numerische Mathematik. 8. Auflage. Vieweg+Teubner, Wiesbaden 2011, ISBN 978-3-8348-1551-4, S. 350 f. 
6. ↑ Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, S. 157. 
7. ↑ Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, S. 156. 
8. ↑ Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, S. 157. 
9. ↑ Hans-Jürgen Reinhardt: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Walter de Gruyter, Berlin/Boston 2012, ISBN 978-3-11-028045-6, S. 42 f. 
10. ↑ John C. Butcher: Numerical Methods for Ordinary Differential Equations. John Wiley & Sons, Chichester 2008, ISBN 978-0-470-72335-7, S. 95–100. 
11. ↑ J. C. Butcher: Numerical methods for ordinary differential equations in the 20th century. In: Journal of Computational and Applied Mathematics. Band 125, Nr. 1–2, 15. Dezember 2000, S. 21 f. (online). 
12. ↑ Peter Deuflhard, Folkmar Bornemann: Numerische Mathematik 2 – Gewöhnliche Differentialgleichungen. 3. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1, S. 228 f. 
13. ↑ Peter Deuflhard, Folkmar Bornemann: Numerische Mathematik 2 – Gewöhnliche Differentialgleichungen. 3. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1, S. 229–231. 
14. ↑ Wolfgang Dahmen, Arnold Reusken: Numerik für Ingenieure und Naturwissenschaftler. 2. Auflage. Springer, Berlin/Heidelberg 2008, ISBN 978-3-540-76492-2, S. 443 f. 
15. ↑ Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8, S. 258 f. 
16. ↑ Jean-Luc Chabert u. a.: A History of Algorithms. Springer, Berlin/Heidelberg 1999, ISBN 978-3-540-63369-3, S. 378 f. 
17. ↑ Jean-Luc Chabert u. a.: A History of Algorithms. Springer, Berlin/Heidelberg 1999, ISBN 978-3-540-63369-3, S. 381–388. 
18. ↑ Wolfgang Dahmen, Arnold Reusken: Numerik für Ingenieure und Naturwissenschaftler. 2. Auflage. Springer, Berlin/Heidelberg 2008, ISBN 978-3-540-76492-2, S. 406 f. 
19. ↑ J. C. Butcher: Numerical methods for ordinary differential equations in the 20th century. In: Journal of Computational and Applied Mathematics. Band 125, Nr. 1–2, 15. Dezember 2000, S. 4–6 (online). 
20. ↑ Peter Deuflhard, Folkmar Bornemann: Numerische Mathematik 2 – Gewöhnliche Differentialgleichungen. 3. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1, S. 160–162. 
21. ↑ Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8, S. 219–221. 
22. ↑ Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8, S. 79 ff. 
23. ↑ J. C. Butcher: Numerical methods for ordinary differential equations in the 20th century. In: Journal of Computational and Applied Mathematics. Band 125, Nr. 1–2, 15. Dezember 2000, S. 26 (online). 
24. ↑ Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, S. 171–173. 
25. ↑ Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8, S. 57–59. 
26. ↑ Peter Deuflhard, Folkmar Bornemann: Numerische Mathematik 2 – Gewöhnliche Differentialgleichungen. 3. Auflage. Walter de Gruyter, Berlin 2008, ISBN 978-3-11-020356-1, S. 199–204. 
27. ↑ Robert Plato: Numerische Mathematik kompakt. 4. Auflage. Vieweg+Teubner, Wiesbaden 2010, ISBN 978-3-8348-1018-2, Kap. 7: Einschrittverfahren für Anfangswertprobleme, S. 173–177. 
28. ↑ Karl Strehmel, Rüdiger Weiner, Helmut Podhaisky: Numerik gewöhnlicher Differentialgleichungen. 2. Auflage. Springer Spektrum, Wiesbaden 2012, ISBN 978-3-8348-1847-8, S. 64–70. 
29. ↑ ode45: Solve nonstiff differential equations — medium order method. MathWorks, abgerufen am 23. November 2017 (englisch).