Runge-Kutta-Verfahren - Ingenieurmathematik (i939)

Man kann systematisch bessere Verfahren zum numerischen Lösen von Differentialgleichungen entwickeln. Diese Systematik stammt von Runge (1895) und Kutta (1901) und wurde von Butcher (1963) noch etwas vereinfacht. Geben hier nur einen groben Überblick und schauen uns an, wie die in diesem Kontext auftretenden Butcher-Tableaus zu verstehen sind.

6.6.1Ziel¶

Zunächst ist zu klären, wie man einem Verfahren ansieht, dass es “besser” ist als die Euler-Verfahren.

Grundsätzlich ist man für gegebene Schrittanzahl $n\in\bbN$ am Gesamtfehler

\max_{k=1,\ldots,n}|y(x_k)-y_k|

(1)

interessiert in Abhängigkeit von der Schrittweite $h$ . Für die beiden Euler-Verfahren gilt

\max_{k=1,\ldots,n}|y(x_k)-y_k|\leq c\,h

(2)

für kleine $h>0$ und eine Konstante $c>0$ . Besser ist

\max_{k=1,\ldots,n}|y(x_k)-y_k|\leq c\,h^p

(3)

mit $p>1$ , da zum Beispiel für $p=2$ das Halbieren der Schrittweite den Fehler auf ein Viertel reduziert. Man nennt den Exponent $p$ auch Konvergenzordnung des Verfahrens.

Je höher die Konvergenzordnung, desto größer kann $h$ gewählt werden und desto geringer also der Rechenaufwandt.

6.6.2Idee¶

Zur Erhöhung der Konvergenzordnung führt man Zwischenschritte ein, d.h. das Verfahren wertet $f$ an mehr Stellen aus und erhält so mehr Informationen über die Differentialgleichung. Dabei kann die Wahl der Stellen für die zusätzlichen Auswertungen auch von den Funktionswerten von $f$ abhängen.

Um geeignete Zwischenschritt zu finden, kann man aus der Forderung nach hoher Konvergenzordnung gewisse Gleichungen aufstellen. Das so entstehende nichtlineare Gleichungssystem muss dann symbolisch gelöst werden und liefert die unter behandelten Butcher-Tableaus. Man kann zeigen, dass beispielsweise für Konvergenzordnung $p=10$ mindestens 10 Zwischenschritte nötig sind, deren Berechnungsvorschriften aus der Lösung eines nichtlinearen Gleichungssystems mit 200 Gleichungen gewonnen werden.

6.6.3Heun-Verfahren¶

Als einfaches Verfahren mit einem Zwischenschritt betrachten wir kurz das Heun-Verfahren. Dieses führt zunächst einen Schritt des expliziten Euler-Verfahrens aus, um eine Näherung $\tilde{y}_{k+1}$ der Ableitung $y'(x_{k+1})$ zu erhalten. Anschließend führt es den eigentlichen Schritt aus, wobei der verwendete Anstieg gerade der Mittelwert aus $y_k$ und $\tilde{y}_{k+1}$ ist:

{\begin{align} \tilde{y}_{k+1}&:=y_k+h\,f(x_k,y_k),\\ y_{k+1}&:=y_k+h\,\frac{f(x_k,y_k)+f(x_{k+1},\tilde{y}_{k+1})}{2}. \end{align}}

(4)

Das Heun-Verfahren hat Konvergenzordnung $p=2$ und liefert damit bei gleicher Schrittweite genauere Ergebnisse.

6.6.4Allgemeine Formulierung¶

Die ein Runge-Kutta-Verfahren mit $s$ Zwischenschritten beschreibenden Zahlen werden üblicherweise als sogenanntes Butcher-Tableau geschrieben.

\begin{array}{c|ccccc} c_1&&&&&\\ c_2&a_{21}&&&&\\ c_3&a_{31}&a_{32}&&&\\ \vdots&\vdots&\vdots&\ddots&&\\ c_s&a_{s,1}&a_{s,2}&\cdots&a_{s,s-1}&\\\hline &b_1&b_2&\cdots&b_{s-1}&b_2 \end{array}

(5)

Die Iterationsvorschrift für das zugehörige Verfahren ist dann

{\begin{align} z_1&:=f\bigl(x_k+c_1\,h,y_k\bigr),\\ z_2&:=f\bigl(x_k+c_2\,h,y_k+h\,a_{21}\,z_1\bigr),\\ z_3&:=f\bigl(x_k+c_3\,h,y_k+h\,(a_{31}\,z_1+a_{32}\,z_2)\bigr),\\ &\;\vdots\\ z_s&:=f\bigl(x_k+c_s\,h,y_k+h\,(a_{s,1}\,z_1+\cdots+a_{s,s-1}\,z_{s-1})\bigr),\\ y_{k+1}&:=y_k+h\,\sum_{i=1}^s b_i\,z_i, \end{align}}

(6)

wobei die Werte $z_1,\ldots,z_s$ nach jedem Schritt verworfen werden können (und somit keinen Index $k$ tragen).

Ein Runge-Kutta-Verfahren mit $s$ Zwischenschritten wird auch als $s$ -stufiges Runge-Kutta-Verfahren bezeichnet.

6.6.5Implizite Runge-Kutta-Verfahren¶

Bei der hier formulierten Variante handelt es sich um explizite Runge-Kutta-Verfahren, da alle Schritt direkt berechnet werden können. Es gibt jedoch auch implizite Runge-Kutta-Verfahren, bei denen $z_1,\ldots,z_s$ als Lösung eines nichtlinearen Gleichungssystems entstehen. In diesem Fall hat das Butcher-Tableau keine Dreieckstruktur, sondern ist voll besetzt, also rechteckig.

Implizite Runge-Kutta-Verfahren haben bei gleichem $s$ üblicherweise eine höhere Konvergenzordnung als explizite Verfahren. Während bei expliziten Verfahren höchstens Kovergenzordnung $s$ möglich ist, kann man mit impliziten Verfahren die Konvergenzordnung $2\,s$ erreichen.

6.6.6Klassisches Runge-Kutta-Verfahren¶

Als klassisches Runge-Kutta-Verfahren oder als RK4-Verfahren wird das explizite Runge-Kutta-Verfahren mit dem Butcher-Tableau

\begin{array}{c|cccc} 0&&&&\\ \frac{1}{2}&\frac{1}{2}&&&\\ \frac{1}{2}&0&\frac{1}{2}&&\\ 1&0&0&1&\\\hline &\frac{1}{6}&\frac{1}{3}&\frac{1}{3}&\frac{1}{6} \end{array}

(10)

bezeichnet.

6.6.7Systeme¶

Analog zu den Euler-Verfahren kann das Runge-Kutta-Verfahren ohne Änderungen an den Formeln auf Systeme von Differentialgleichungen angewendet werden. Die Zwischenwerte $z_i$ sind dann Vektoren mit $n$ Komponenten, wobei $n$ die Anzahl der Gleichungen im System ist.

6.6.8Numerischer Vergleich¶

Die unterschiedliche Konvergenzordnung von Euler-, Heun- und klassischem Runge-Kutta-Verfahren kann man leicht numerisch nachvollziehen. Lösen dazu das Anfangswertproblem

y'=\sin x-\frac{y}{t},\qquad y(1)=1.

(11)

Wählen die Schrittweiten so, dass der Rechenaufwandt für die verschiedenen Verfahren in etwa gleich ist. Ist $h$ die Schrittweite für das Euler-Verfahren, so führen wir das Heun-Verfahren mit Schrittweite $2\,h$ und das klassische Runge-Kutta-Verfahren mit Schrittweite $4\,h$ aus.

import matplotlib.pyplot as plt
import numpy as np

f = lambda x, y: np.sin(x) - y / x
x0 = 1
y0 = 1

x_max = x0 + 8  # Intervall, auf dem die Lösung berechnet werden soll
h_min = 0.5  # kleinste Schrittweite

# Euler
h = h_min
x = [x0]
y = [y0]
while x[-1] < x_max:
    x_new = x[-1] + h
    y_new = y[-1] + h * f(x[-1], y[-1])
    x.append(x_new)
    y.append(y_new)
x_euler = x
y_euler = y

# Heun
h = 2 * h_min
x = [x0]
y = [y0]
while x[-1] < x_max:
    x_new = x[-1] + h
    y_tilde = y[-1] + h * f(x[-1], y[-1])
    y_new = y[-1] + 0.5 * h * (f(x[-1], y[-1]) + f(x_new, y_tilde))
    x.append(x_new)
    y.append(y_new)
x_heun = x
y_heun = y

# Runge-Kutta
h = 4 * h_min
x = [x0]
y = [y0]
while x[-1] < x_max:
    x_new = x[-1] + h
    z1 = f(x[-1], y[-1])
    z2 = f(x[-1] + 0.5 * h, y[-1] + 0.5 * h * z1)
    z3 = f(x[-1] + 0.5 * h, y[-1] + 0.5 * h * z2)
    z4 = f(x[-1] + h, y[-1] + h * z3)
    y_new = y[-1] + h * (1/6 * z1 + 1/3 * z2 + 1/3 * z3 + 1/6 * z4)
    x.append(x_new)
    y.append(y_new)
x_rk = x
y_rk = y

# Plot
fig, ax = plt.subplots(figsize=(10, 5))

sol_x = np.linspace(x0 - 0.5, x_max + 0.5, 100)
sol_y = (np.sin(sol_x) - sol_x * np.cos(sol_x) + 1 - np.sin(1) + np.cos(1)) / sol_x
ax.plot(sol_x, sol_y, '-r', label='exakte Lösung')

ax.plot(x_euler, y_euler, '-bo', label='Euler')
ax.plot(x_heun, y_heun, '-mo', label='Heun')
ax.plot(x_rk, y_rk, '-go', label='Runge-Kutta')

ax.plot(x0, y0, 'or', label='Anfangsbedingung')

ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid()
plt.show()