Runge-Kutta-Verfahren

In der numerischen Analyse werden die Runge-Kutta-Verfahren (engl: /ˈrʊŋəˈkʊtɑː/ (listen) RUUNG-ə-KUUT-tah) sind eine Familie von impliziten und expliziten iterativen Methoden, zu denen auch die Euler-Methode gehört, die bei der zeitlichen Diskretisierung für die Näherungslösungen von simultanen nichtlinearen Gleichungen verwendet werden. Diese Methoden wurden um 1900 von den deutschen Mathematikern Carl Runge und Wilhelm Kutta entwickelt. ⓘ

Vergleich der Runge-Kutta-Methoden für die Differentialgleichung $y'=\sin(t)^{2}\cdot y$ (rot ist die exakte Lösung) ⓘ

Das Runge-Kutta-Verfahren

Vom klassischen Runge-Kutta-Verfahren verwendete Steigungen ⓘ

Das bekannteste Mitglied der Runge-Kutta-Familie wird allgemein als "RK4", das "klassische Runge-Kutta-Verfahren" oder einfach als "Runge-Kutta-Verfahren" bezeichnet. ⓘ

Ein Anfangswertproblem sei wie folgt spezifiziert:

${\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}.$ ⓘ

Hier ist $y$ ist eine unbekannte Funktion (Skalar oder Vektor) der Zeit $t$, die wir approximieren möchten; man sagt uns, dass ${\frac {dy}{dt}}$die Geschwindigkeit, mit der $y$ sich ändert, eine Funktion ist von $t$ und von $y$ selbst ist. Zum Anfangszeitpunkt $t_{0}$ ist der entsprechende $y$ Wert gleich $y_{0}$. Die Funktion $f$ und die Anfangsbedingungen $t_{0}$, $y_{0}$ sind gegeben. ⓘ

Wählen Sie nun eine Schrittweite h > 0 und definieren Sie ⓘ

${\begin{aligned}y_{n+1}&=y_{n}+{\frac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right)h,\\t_{n+1}&=t_{n}+h\\\end{aligned}}$ ⓘ

für n = 0, 1, 2, 3, ..., mit

${\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}$

(Hinweis: Die obigen Gleichungen haben in verschiedenen Texten unterschiedliche, aber gleichwertige Definitionen). ⓘ

Hier ist $y_{n+1}$ ist die RK4-Approximation von $y(t_{n+1})$, und der nächste Wert ($y_{n+1}$) wird bestimmt durch den aktuellen Wert ($y_{n}$) plus dem gewichteten Durchschnitt von vier Inkrementen, wobei jedes Inkrement das Produkt aus der Größe des Intervalls h und einer geschätzten Steigung ist, die durch die Funktion f auf der rechten Seite der Differentialgleichung angegeben wird.

• $k_{1}$ ist die Steigung am Anfang des Intervalls unter Verwendung von $y$ (Eulersche Methode);
• $k_{2}$ ist die Steigung in der Mitte des Intervalls unter Verwendung von $y$ und $k_{1}$;
• $k_{3}$ ist wiederum die Steigung in der Mitte des Intervalls, aber jetzt unter Verwendung von $y$ und $k_{2}$;
• $k_{4}$ ist die Steigung am Ende des Intervalls, unter Verwendung von $y$ und $k_{3}$. ⓘ

Bei der Mittelung der vier Steigungen werden die Steigungen in der Mitte stärker gewichtet. Wenn $f$ unabhängig ist von $y$ist, so dass die Differentialgleichung einem einfachen Integral entspricht, dann ist RK4 die Simpsonsche Regel. ⓘ

Die RK4-Methode ist eine Methode vierter Ordnung, was bedeutet, dass der lokale Abschneidefehler in der Größenordnung von $O(h^{5})$liegt, während der gesamte akkumulierte Fehler in der Größenordnung von $O(h^{4})$. ⓘ

In vielen praktischen Anwendungen wird die Funktion $f$ unabhängig ist von $t$ (so genanntes autonomes System oder zeitinvariantes System, insbesondere in der Physik) und ihre Inkremente überhaupt nicht berechnet und nicht an die Funktion $f$übergeben, wobei nur die endgültige Formel für $t_{n+1}$ verwendet. ⓘ

Explizite Runge-Kutta-Verfahren

Die Familie der expliziten Runge-Kutta-Verfahren ist eine Verallgemeinerung des oben erwähnten RK4-Verfahrens. Sie ist gegeben durch

$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},$

wobei

${\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1})),\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2})),\\&\ \ \vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+h(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1})).\end{aligned}}$

(Hinweis: Die obigen Gleichungen können in einigen Texten unterschiedliche, aber gleichwertige Definitionen haben). ⓘ

Um eine bestimmte Methode zu spezifizieren, muss man die ganze Zahl s (die Anzahl der Stufen) und die Koeffizienten a_ij (für 1 ≤ j < i ≤ s), bi (für i = 1, 2, ..., s) und ci (für i = 2, 3, ..., s) angeben. Die Matrix [a_ij] wird als Runge-Kutta-Matrix bezeichnet, während bi und ci als die Gewichte und die Knoten bezeichnet werden. Diese Daten werden in der Regel in einer Gedächtnisstütze angeordnet, die als Butcher-Tableau (nach John C. Butcher) bezeichnet wird:

$0$
$c_{2}$	$a_{21}$
$c_{3}$	$a_{31}$	$a_{32}$
$\vdots$	$\vdots$		$\ddots$
$c_{s}$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
	$b_{1}$	$b_{2}$	$\cdots$	$b_{s-1}$	$b_{s}$

ⓘ

Eine Taylorreihenentwicklung zeigt, dass das Runge-Kutta-Verfahren nur dann konsistent ist, wenn

$\sum _{i=1}^{s}b_{i}=1.$

Es gibt auch begleitende Anforderungen, wenn man verlangt, dass die Methode eine bestimmte Ordnung p hat, was bedeutet, dass der lokale Abbruchfehler O(hp⁺¹) ist. Diese lassen sich aus der Definition des Trunkierungsfehlers selbst ableiten. Zum Beispiel hat eine zweistufige Methode die Ordnung 2, wenn b₁ + b₂ = 1, b₂c₂ = 1/2 und b₂a₂₁ = 1/2. Eine beliebte Bedingung für die Bestimmung der Koeffizienten ist

$\sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.$

Diese Bedingung allein ist jedoch weder ausreichend noch notwendig für die Konsistenz. ⓘ

Im Allgemeinen, wenn ein explizites $s$-Stufen-Runge-Kutta-Verfahren die Ordnung $p$hat, dann kann bewiesen werden, dass die Anzahl der Stufen folgende Bedingungen erfüllen muss $s\geq p$, und wenn $p\geq 5$, dann $s\geq p+1$. Es ist jedoch nicht bekannt, ob diese Schranken in allen Fällen scharf sind; zum Beispiel haben alle bekannten Methoden der Ordnung 8 mindestens 11 Stufen, obwohl es möglich ist, dass es Methoden mit weniger Stufen gibt. (Die obige Schranke legt nahe, dass es eine Methode mit 9 Stufen geben könnte; es könnte aber auch sein, dass die Schranke einfach nicht scharf ist.) In der Tat ist es ein offenes Problem was die genaue Mindestanzahl von Stufen $s$ für ein explizites Runge-Kutta-Verfahren mit der Ordnung $p$ in den Fällen, in denen noch keine Methoden gefunden wurden, die die obigen Schranken mit Gleichheit erfüllen. Einige Werte, die bekannt sind, sind:

${\begin{array}{c|cccccccc}p&1&2&3&4&5&6&7&8\\\hline \min s&1&2&3&4&6&7&9&11\end{array}}$

Die obigen beweisbaren Grenzen implizieren, dass wir keine Methoden der Ordnung $p=1,2,\ldots ,6$ zu finden, die weniger Stufen benötigen als die Methoden, die wir für diese Ordnungen bereits kennen. Es ist jedoch denkbar, dass wir eine Methode der Ordnung $p=7$ zu finden, die nur 8 Stufen hat, während die einzigen heute bekannten Verfahren mindestens 9 Stufen haben, wie in der Tabelle gezeigt. ⓘ

Beispiele

Die RK4-Methode fällt in diesen Rahmen. Ihr Tableau lautet ⓘ

0
1/2	1/2
1/2	0	1/2
1	0	0	1
	1/6	1/3	1/3	1/6

ⓘ

Eine leichte Abwandlung des Runge-Kutta-Verfahrens geht ebenfalls auf Kutta (1901) zurück und wird als 3/8-Regel bezeichnet. Der Hauptvorteil dieser Methode besteht darin, dass fast alle Fehlerkoeffizienten kleiner sind als bei der bekannten Methode, aber sie erfordert etwas mehr FLOPs (Gleitkommaoperationen) pro Zeitschritt. Das Butcher-Tableau lautet ⓘ

0
1/3	1/3
2/3	-1/3	1
1	1	-1	1
	1/8	3/8	3/8	1/8

ⓘ

Das einfachste Runge-Kutta-Verfahren ist jedoch das (Vorwärts-)Euler-Verfahren, das durch die Formel $y_{n+1}=y_{n}+hf(t_{n},y_{n})$. Dies ist das einzige konsistente explizite Runge-Kutta-Verfahren mit einer Stufe. Das entsprechende Tableau lautet ⓘ

0
	1

ⓘ

Verfahren zweiter Ordnung mit zwei Stufen

Ein Beispiel für ein zweistufiges Verfahren zweiter Ordnung ist die Mittelpunktsmethode:

$y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hf(t_{n},\ y_{n})\right).$

Das entsprechende Tableau lautet ⓘ

0
1/2	1/2
	0	1

ⓘ

Das Mittelpunktverfahren ist nicht das einzige Runge-Kutta-Verfahren zweiter Ordnung mit zwei Stufen; es gibt eine Familie solcher Verfahren, die durch α parametrisiert und durch die Formel ⓘ

$y_{n+1}=y_{n}+h{\bigl (}(1-{\tfrac {1}{2\alpha }})f(t_{n},y_{n})+{\tfrac {1}{2\alpha }}f(t_{n}+\alpha h,y_{n}+\alpha hf(t_{n},y_{n})){\bigr )}.$ ⓘ

Ihr Butcher-Tableau lautet ⓘ

0
$\alpha$	$\alpha$
	$(1-{\tfrac {1}{2\alpha }})$	${\tfrac {1}{2\alpha }}$

ⓘ

In dieser Familie, $\alpha ={\tfrac {1}{2}}$ die Mittelpunktsmethode, $\alpha =1$ die Heunsche Methode und $\alpha ={\tfrac {2}{3}}$ ist die Ralston-Methode. ⓘ

Verwenden Sie

Betrachten wir als Beispiel das zweistufige Runge-Kutta-Verfahren zweiter Ordnung mit α = 2/3, auch bekannt als Ralston-Verfahren. Es wird durch das Tableau ⓘ

	0 ⓘ
	2/3	2/3
		1/4	3/4

mit den entsprechenden Gleichungen ⓘ

${\begin{aligned}k_{1}&=f(t_{n},\ y_{n}),\\k_{2}&=f(t_{n}+{\tfrac {2}{3}}h,\ y_{n}+{\tfrac {2}{3}}hk_{1}),\\y_{n+1}&=y_{n}+h\left({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2}\right).\end{aligned}}$ ⓘ

Dieses Verfahren wird zur Lösung des Anfangswertproblems verwendet

${\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]$

mit der Schrittweite h = 0,025, so dass die Methode vier Schritte benötigt. ⓘ

Die Methode geht wie folgt vor:

$t_{0}=1\colon$		ⓘ
	$y_{0}=1$
$t_{1}=1.025\colon$
	$y_{0}=1$	$k_{1}=2.557407725$	$k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400$
	$y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}$
$t_{2}=1.05\colon$
	$y_{1}=1.066869388$	$k_{1}=2.813524695$	$k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})$
	$y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}$
$t_{3}=1.075\colon$
	$y_{2}=1.141332181$	$k_{1}=3.183536647$	$k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})$
	$y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}$
$t_{4}=1.1\colon$
	$y_{3}=1.227417567$	$k_{1}=3.796866512$	$k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})$
	$y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.$

Die numerischen Lösungen entsprechen den unterstrichenen Werten. ⓘ

Adaptive Runge-Kutta-Verfahren

Adaptive Methoden sind so konzipiert, dass sie eine Schätzung des lokalen Abbruchfehlers eines einzelnen Runge-Kutta-Schrittes liefern. Dies geschieht durch zwei Methoden, eine mit der Ordnung $p$ und eine mit der Ordnung $p-1$. Diese Methoden sind miteinander verwoben, d.h. sie haben gemeinsame Zwischenschritte. Dadurch ist der Rechenaufwand für die Schätzung des Fehlers im Vergleich zu einem Schritt mit der Methode höherer Ordnung gering oder vernachlässigbar. ⓘ

Während der Integration wird die Schrittweite so angepasst, dass der geschätzte Fehler unter einem vom Benutzer definierten Schwellenwert bleibt: Ist der Fehler zu groß, wird ein Schritt mit einer niedrigeren Schrittgröße wiederholt; ist der Fehler viel kleiner, wird die Schrittgröße erhöht, um Zeit zu sparen. Dies führt zu einer (fast) optimalen Schrittgröße, die Rechenzeit spart. Außerdem muss der Benutzer keine Zeit darauf verwenden, eine geeignete Schrittgröße zu finden. ⓘ

Der Schritt niedrigerer Ordnung ist gegeben durch

$y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i},$

wobei $k_{i}$ sind die gleichen wie bei der Methode höherer Ordnung. Dann ist der Fehler

$e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i},$

der ist $O(h^{p})$. Das Butcher-Tableau für diese Art von Verfahren wird erweitert, um die Werte von $b_{i}^{*}$:

	0 ⓘ
	$c_{2}$	$a_{21}$
	$c_{3}$	$a_{31}$	$a_{32}$
	$\vdots$	$\vdots$		$\ddots$
	$c_{s}$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_{1}$	$b_{2}$	$\cdots$	$b_{s-1}$	$b_{s}$
		$b_{1}^{*}$	$b_{2}^{*}$	$\cdots$	$b_{s-1}^{*}$	$b_{s}^{*}$

Das Runge-Kutta-Fehlberg-Verfahren hat zwei Verfahren der Ordnung 5 und 4. Ihr erweitertes Butcher-Tableau lautet:

	0 ⓘ
	1/4	1/4
	3/8	3/32	9/32
	12/13	1932/2197	−7200/2197	7296/2197
	1	439/216	−8	3680/513	-845/4104
	1/2	−8/27	2	−3544/2565	1859/4104	−11/40
		16/135	0	6656/12825	28561/56430	−9/50	2/55
		25/216	0	1408/2565	2197/4104	−1/5	0

Bei der einfachsten adaptiven Runge-Kutta-Methode wird jedoch die Heun-Methode der Ordnung 2 mit der Euler-Methode der Ordnung 1 kombiniert. Sein erweitertes Butcher-Tableau ist:

	0 ⓘ
	1	1
		1/2	1/2
		1	0

Weitere adaptive Runge-Kutta-Methoden sind die Bogacki-Shampine-Methode (Ordnung 3 und 2), die Cash-Karp-Methode und die Dormand-Prince-Methode (beide mit Ordnung 5 und 4). ⓘ

Nicht-konfluente Runge-Kutta-Verfahren

Ein Runge-Kutta-Verfahren wird als nicht-konfluent bezeichnet, wenn alle $c_{i},\,i=1,2,\ldots ,s$ verschieden sind. ⓘ

Runge-Kutta-Nyström-Verfahren

Runge-Kutta-Nyström-Methoden sind spezielle Runge-Kutta-Methoden, die für Differentialgleichungen zweiter Ordnung der folgenden Form optimiert sind:

${\frac {d^{2}y}{dt^{2}}}=f(y,{\dot {y}},t).$ ⓘ

Implizite Runge-Kutta-Verfahren

Alle bisher genannten Runge-Kutta-Verfahren sind explizite Verfahren. Explizite Runge-Kutta-Methoden sind im Allgemeinen für die Lösung von steifen Gleichungen ungeeignet, da ihr Bereich der absoluten Stabilität klein ist; insbesondere ist er begrenzt. Dieses Problem ist besonders wichtig bei der Lösung von partiellen Differentialgleichungen. ⓘ

Die Instabilität der expliziten Runge-Kutta-Methoden motiviert die Entwicklung impliziter Methoden. Ein implizites Runge-Kutta-Verfahren hat die Form ⓘ

$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},$

wobei

$k_{i}=f\left(t_{n}+c_{i}h,\ y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.$ ⓘ

Der Unterschied zu einem expliziten Verfahren ist, dass bei einem expliziten Verfahren die Summe über j nur bis i - 1 geht. Dies zeigt sich auch im Butcher-Tableau: Die Koeffizientenmatrix $a_{ij}$ einer expliziten Methode ist unterdreieckig. Bei einem impliziten Verfahren geht die Summe über j bis s, und die Koeffizientenmatrix ist nicht dreieckig, so dass sich ein Butcher-Tableau der Form ⓘ

${\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\&b_{1}^{*}&b_{2}^{*}&\dots &b_{s}^{*}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}$

Siehe oben unter Adaptive Runge-Kutta-Methoden für die Erklärung der $b^{*}$ Zeile. ⓘ

Dieser Unterschied hat zur Folge, dass bei jedem Schritt ein System algebraischer Gleichungen gelöst werden muss. Dies erhöht den Rechenaufwand erheblich. Wird ein Verfahren mit s Stufen zur Lösung einer Differentialgleichung mit m Komponenten verwendet, so hat das algebraische Gleichungssystem ms Komponenten. Dies kann mit impliziten linearen Mehrschrittverfahren (der anderen großen Familie von Verfahren für ODEs) verglichen werden: ein implizites lineares Mehrschrittverfahren mit s Schritten muss ein System algebraischer Gleichungen mit nur m Komponenten lösen, so dass die Größe des Systems nicht mit der Anzahl der Schritte zunimmt. ⓘ

Beispiele

Das einfachste Beispiel für ein implizites Runge-Kutta-Verfahren ist das Rückwärts-Euler-Verfahren:

$y_{n+1}=y_{n}+hf(t_{n}+h,\ y_{n+1}).\,$ ⓘ

Das Butcher-Tableau für dieses ist einfach:

${\begin{array}{c|c}1&1\\\hline &1\\\end{array}}$ ⓘ

Dieses Butcher-Tableau entspricht den Formeln ⓘ

$k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{and}}\quad y_{n+1}=y_{n}+hk_{1},$ ⓘ

die man umformulieren kann, um die oben genannte Formel für das Backward-Euler-Verfahren zu erhalten. ⓘ

Ein weiteres Beispiel für ein implizites Runge-Kutta-Verfahren ist die Trapezregel. Ihr Butcher-Tableau lautet:

${\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&1&0\\\end{array}}$ ⓘ

Die Trapezregel ist eine Kollokationsmethode (wie in diesem Artikel beschrieben). Alle Kollokationsverfahren sind implizite Runge-Kutta-Verfahren, aber nicht alle impliziten Runge-Kutta-Verfahren sind Kollokationsverfahren. ⓘ

Die Gauß-Legendre-Methoden bilden eine Familie von Kollokationsmethoden, die auf der Gauß-Quadratur basieren. Ein Gauß-Legendre-Verfahren mit s Stufen hat die Ordnung 2s (es können also Verfahren mit beliebig hoher Ordnung konstruiert werden). Die Methode mit zwei Stufen (und damit vierter Ordnung) hat ein Butcher-Tableau:

${\begin{array}{c|cc}{\frac {1}{2}}-{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {1}{6}}{\sqrt {3}}\\{\frac {1}{2}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&{\frac {1}{2}}+{\frac {1}{2}}{\sqrt {3}}&{\frac {1}{2}}-{\frac {1}{2}}{\sqrt {3}}\end{array}}$ ⓘ

Stabilität

Der Vorteil der impliziten Runge-Kutta-Methoden gegenüber den expliziten ist ihre größere Stabilität, insbesondere bei der Anwendung auf steife Gleichungen. Betrachten wir die lineare Testgleichung $y'=\lambda y$. Ein auf diese Gleichung angewandtes Runge-Kutta-Verfahren reduziert sich auf die Iteration $y_{n+1}=r(h\lambda )\,y_{n}$, wobei r gegeben ist durch ⓘ

$r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},$ ⓘ

wobei e für den Vektor der Einsen steht. Die Funktion r wird als Stabilitätsfunktion bezeichnet. Aus der Formel folgt, dass r der Quotient zweier Polynome vom Grad s ist, wenn das Verfahren s Stufen hat. Explizite Methoden haben eine streng untere Dreiecksmatrix A, was bedeutet, dass det(I - zA) = 1 ist und dass die Stabilitätsfunktion ein Polynom ist. ⓘ

Die numerische Lösung der linearen Testgleichung zerfällt zu Null, wenn | r(z) | < 1 mit z = hλ. Die Menge solcher z wird als der Bereich der absoluten Stabilität bezeichnet. Insbesondere wird die Methode als absolut stabil bezeichnet, wenn alle z mit Re(z) < 0 im Bereich der absoluten Stabilität liegen. Die Stabilitätsfunktion eines expliziten Runge-Kutta-Verfahrens ist ein Polynom, so dass explizite Runge-Kutta-Verfahren niemals A-stabil sein können. ⓘ

Wenn das Verfahren die Ordnung p hat, dann erfüllt die Stabilitätsfunktion $r(z)={\textrm {e}}^{z}+O(z^{p+1})$ als $z\to 0$. Daher ist es von Interesse, Quotienten von Polynomen bestimmten Grades zu untersuchen, die die Exponentialfunktion am besten approximieren. Diese sind als Padé-Approximanten bekannt. Eine Padé-Approximante mit einem Zähler vom Grad m und einem Nenner vom Grad n ist A-stabil, wenn und nur wenn m ≤ n ≤ m + 2. ⓘ

Das Gauß-Legendre-Verfahren mit s Stufen hat die Ordnung 2s, seine Stabilitätsfunktion ist also die Padé-Approximante mit m = n = s. Daraus folgt, dass das Verfahren A-stabil ist. Dies zeigt, dass das A-stabile Runge-Kutta-Verfahren eine beliebig hohe Ordnung haben kann. Im Gegensatz dazu kann die Ordnung von A-stabilen linearen Mehrschrittverfahren zwei nicht überschreiten. ⓘ

B-Stabilität

Das Konzept der A-Stabilität für die Lösung von Differentialgleichungen steht im Zusammenhang mit der linearen autonomen Gleichung $y'=\lambda y$. Dahlquist schlug die Untersuchung der Stabilität numerischer Verfahren bei der Anwendung auf nichtlineare Systeme vor, die eine Monotoniebedingung erfüllen. Die entsprechenden Konzepte wurden als G-Stabilität für Mehrschrittverfahren (und die damit verbundenen Einschrittverfahren) und als B-Stabilität (Butcher, 1975) für Runge-Kutta-Verfahren definiert. Ein Runge-Kutta-Verfahren, das auf ein nichtlineares System $y'=f(y)$angewendet, das $\langle f(y)-f(z),\ y-z\rangle <0$erfüllt, wird B-stabil genannt, wenn diese Bedingung $\|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|$ für zwei numerische Lösungen. ⓘ

Sei $B$, $M$ und $Q$ seien drei $s\times s$ Matrizen, definiert durch ⓘ

$B=\operatorname {diag} (b_{1},b_{2},\ldots ,b_{s}),\,M=BA+A^{T}B-bb^{T},\,Q=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-1}.$ ⓘ

Ein Runge-Kutta-Verfahren wird als algebraisch stabil bezeichnet, wenn die Matrizen $B$ und $M$ beide nicht-negativ definit sind. Eine hinreichende Bedingung für B-Stabilität ist: $B$ und $Q$ sind nicht-negativ definitiv. ⓘ

Herleitung des Runge-Kutta-Verfahrens vierter Ordnung

Im Allgemeinen kann ein Runge-Kutta-Verfahren der Ordnung $s$ kann wie folgt geschrieben werden:

$y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1}),$

wobei:

$k_{i}=y_{t}+h\cdot \sum _{j=1}^{s}\beta _{ij}f\left(k_{j},\ t_{n}+\alpha _{i}h\right)$

Inkremente sind, die durch die Auswertung der Ableitungen von $y_{t}$ in der $i$-ten Ordnung. ⓘ

Wir entwickeln die Ableitung für das Runge-Kutta-Verfahren vierter Ordnung unter Verwendung der allgemeinen Formel mit $s=4$ die, wie oben erläutert, am Anfangs-, Mittel- und Endpunkt eines beliebigen Intervalls $(t,\ t+h)$; daher wählen wir:

${\begin{aligned}&\alpha _{i}&&\beta _{ij}\\\alpha _{1}&=0&\beta _{21}&={\frac {1}{2}}\\\alpha _{2}&={\frac {1}{2}}&\beta _{32}&={\frac {1}{2}}\\\alpha _{3}&={\frac {1}{2}}&\beta _{43}&=1\\\alpha _{4}&=1&&\\\end{aligned}}$ ⓘ

und $\beta _{ij}=0$ sonst. Wir beginnen mit der Definition der folgenden Größen:

${\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},\ t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)\end{aligned}}$

wobei $y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}$ und $y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}$. Wenn wir definieren:

${\begin{aligned}k_{1}&=f(y_{t},\ t)\\k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hk_{3},\ t+h\right)\end{aligned}}$

und für die vorangegangenen Beziehungen können wir zeigen, dass die folgenden Gleichungen gelten, bis ${\mathcal {O}}(h^{2})$:

${\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},\ t\right),\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t},\ t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\right]\end{aligned}}$

wobei:

${\frac {d}{dt}}f(y_{t},\ t)={\frac {\partial }{\partial y}}f(y_{t},\ t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},\ t)=f_{y}(y_{t},\ t){\dot {y}}+f_{t}(y_{t},\ t):={\ddot {y}}_{t}$

ist die Gesamtableitung von $f$ nach der Zeit. ⓘ

Wenn wir nun die allgemeine Formel mit Hilfe der soeben abgeleiteten Formel ausdrücken, erhalten wir:

${\begin{aligned}y_{t+h}={}&y_{t}+h\left\lbrace a\cdot f(y_{t},\ t)+b\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right.+\\&{}+c\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]+\\&{}+d\cdot \left[f(y_{t},\ t)+h{\frac {d}{dt}}\left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f(y_{t},\ t)+\left.{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\={}&y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}$ ⓘ

und vergleichen dies mit der Taylorreihe von $y_{t+h}$ um $t$:

${\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{24}}y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},\ t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},\ t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},\ t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},\ t)\end{aligned}}$ ⓘ

erhalten wir ein System von Nebenbedingungen für die Koeffizienten:

${\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}$

das, wenn gelöst, folgendes ergibt $a={\frac {1}{6}},b={\frac {1}{3}},c={\frac {1}{3}},d={\frac {1}{6}}$ wie oben angegeben. ⓘ

Butcher-Tableau

Man kann die charakteristischen Koeffizienten $c_{j}$, $b_{j}$, $a_{jl}$ übersichtlich im Runge-Kutta-Tableau (auch Butcher-Schema, -Tableau oder engl. Butcher array genannt) anordnen. Hierbei ist die Matrix A bei einem expliziten Verfahren eine strikte untere Dreiecksmatrix (Nilpotente Dreiecksmatrix). ⓘ

${\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}={\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\\end{array}}\;,$  $c={\begin{bmatrix}c_{1}\\\vdots \\c_{j}\\\vdots \\c_{s}\end{bmatrix}}\;,\quad b={\begin{bmatrix}b_{1}\\\vdots \\b_{j}\\\vdots \\b_{s}\end{bmatrix}}\;,\quad A={\begin{bmatrix}a_{11}&\dots &a_{1l}&\dots &a_{1s}\\\vdots &\ddots &\vdots &\ddots &\vdots \\a_{j1}&\dots &a_{jl}&\dots &a_{js}\\\vdots &\ddots &\vdots &\ddots &\vdots \\a_{s1}&\dots &a_{sl}&\dots &a_{ss}\\\end{bmatrix}}\;.$ ⓘ

Konsistenzordnung und Konvergenzordnung

Eine wichtige Eigenschaft zum Vergleich von Verfahren ist die Konsistenzordnung, die auf dem Begriff des lokalen Diskretisierungsfehlers $\tau =y_{n+1}-y(t_{n+1})$ beruht. Dabei ist $y_{n+1}$ die numerische Lösung nach einem Schritt und $y(t_{n+1})$ die exakte Lösung. Ein Einschrittverfahren heißt konsistent von der Ordnung $p$ (hat Konsistenzordnung $p$), falls für den lokalen Diskretisierungsfehler gilt:

$\tau \leq {\mathcal {O}}(h^{p+1})$ (Zur Notation siehe Landau-Symbole). ⓘ

Die Konsistenzordnung kann durch Taylorentwicklung von $\tau$ oder der exakten und numerischen Lösung bestimmt werden. Allgemein gilt:

Konsistenzordnung $p$ und Stabilität $\Rightarrow$ Konvergenzordnung $p$ ⓘ

Bei Einschrittverfahren wie den Runge-Kutta-Verfahren gilt sogar, sofern $f$ und die Verfahrensvorschrift Lipschitz-stetig sind:

Konsistenzordnung $p$ $\Rightarrow$ Konvergenzordnung $p$ ⓘ

Aus der Konsistenzbedingung (z. B. soll das Verfahren Ordnung 4 haben) ergeben sich Konsistenzgleichungen (engl. conditions) für die Koeffizienten des Runge-Kutta-Verfahrens. Die Gleichungen und ihre Anzahl können mit Hilfe von Taylorentwicklung oder der Theorie der Butcher-Bäume ermittelt werden. Mit zunehmender Ordnung wächst die Zahl der zu lösenden nicht-linearen Konsistenzgleichungen schnell an. Das Aufstellen der Konsistenzgleichungen ist bereits nicht einfach, kann jedoch mit Hilfe der Butcher-Bäume von Computeralgebrasystemen erledigt werden. Das Lösen ist allerdings noch schwieriger und bedarf Erfahrung und Fingerspitzengefühl, um „gute“ Koeffizienten zu erhalten. ⓘ

Ein explizites $s$-stufiges Runge-Kutta-Verfahren hat höchstens Konvergenzordnung $s$, ein implizites dagegen bis zu $2s$. ⓘ

Um die Genauigkeit eines Ergebnisses zu verbessern, gibt es zwei Möglichkeiten:

1. Man kann die Schrittweite verkleinern, das heißt, man erhöht die Anzahl der Diskretisierungspunkte.
2. Man kann Verfahren höherer Konvergenzordnung wählen.

Welche Strategie die bessere ist, hängt von der konkreten Problemstellung ab, die Erhöhung der Konvergenzordnung ist allerdings nur bis zu einer bestimmten Grenze sinnvoll, da wegen der Butcher-Schranken die Stufenzahl $s$ schneller wächst als die Ordnung $p$. Für $s\geq 5$ existiert beispielsweise kein explizites $s$-stufiges RKV der Konvergenzordnung $q=s$. ⓘ

Zeitschrittweitensteuerung: Eingebettete Verfahren

Um die Effizienz der Verfahren zu erhöhen, ist es sinnvoll, den Zeitschritt einer Fehlertoleranz anzupassen. Runge-Kutta-Verfahren bieten hierzu über eingebettete Verfahren eine relativ einfache Möglichkeit. Diese bestehen aus einem zweiten Satz an Koeffizienten ${\hat {b}}$ für ein zweites Verfahren:

${\hat {y}}_{n+1}=y_{n}+h\sum _{j=1}^{s}{\hat {b}}_{j}k_{j},$

wobei die Koeffizienten so gewählt werden, dass sich ein schlechteres Verfahren, konkret ein Verfahren von niedrigerer Ordnung ergibt als das ursprüngliche. Dann ist die Differenz

$y_{n+1}-{\hat {y}}_{n+1}=h\sum _{j=1}^{s}(b_{j}-{\hat {b}}_{j})k_{j}$

eine Schätzung des lokalen Fehlers des ursprünglichen Verfahrens von der Ordnung wie das eingebettete Verfahren. Zur Berechnung sind keine neuen Funktionsauswertungen notwendig, sondern nur Linearkombinationen der bereits berechneten $k_{j}$. Die Bestimmung einer neuen Zeitschrittweite aus dem Fehlerschätzer kann über verschiedene Schrittweitensteuerungen erfolgen. ⓘ

Im expliziten Fall sind die bekanntesten eingebetteten Verfahren die Runge-Kutta-Fehlberg- sowie die Dormand-Prince-Formeln (DOPRI). ⓘ