Runge-Kutta-Verfahren

Differentialgleichungen – Kapitel 12

M. Ladecký, L. Pastewka

IMTEK, Universität Freiburg

Lernziele

Die Studierenden sollen…

• … die Grundidee der Runge-Kutta-Verfahren beschreiben können.
• … den Unterschied zwischen verschiedenen Ordnungen von Runge-Kutta-Verfahren erklären können.
• … das Konzept der Konvergenzordnung erklären können.
• … Runge-Kutta-Verfahren zur numerischen Lösung von Differentialgleichungen anwenden können.
• … die Funktion scipy.integrate.solve_ivp verstehen und anwenden können.

Motivation

• Euler-Verfahren: Annahme konstanter Steigung → Verfahren erster Ordnung → sehr kleine Zeitschritte für akzeptable Genauigkeit nötig
• Runge-Kutta-Verfahren: Familie von Verfahren höherer Ordnung
• Idee: Steigung nicht nur am Intervallanfang, sondern auch an Zwischenpunkten auswerten
• Geeignete Gewichtung der Steigungen reduziert den Fehler systematisch — ohne höhere Ableitungen
• Grundlage moderner ODE-Löser wie scipy.integrate.solve_ivp

Ausgangspunkt: das Euler-Verfahren

\[ \mathbf{y}(t + \Delta t) = \mathbf{y}(t) + \Delta t \, \mathbf{g}(\mathbf{y}(t), t) + \mathcal{O}(\Delta t^2) \]

• \(\mathcal{O}(\Delta t^2)\): lokaler Abbruchfehler (Fehler pro Schritt)
• Integration über Zeit \(T\) benötigt \(N = T/\Delta t\) Schritte → globaler Fehler \(\mathcal{O}(\Delta t)\)

Konvergenzordnung: Ein Verfahren hat Ordnung \(p\), wenn der globale Fehler wie \(\mathcal{O}(\Delta t^p)\) skaliert. Lokaler Abbruchfehler: \(\mathcal{O}(\Delta t^{p+1})\). Euler: \(p = 1\).

Das Heun-Verfahren

Idee: durchschnittliche Steigung über das Intervall verwenden — Steigung am Ende muss zunächst geschätzt werden.

Prädiktor (Euler-Schritt): \[ \tilde{\mathbf{y}}(t + \Delta t) = \mathbf{y}(t) + \Delta t \, \mathbf{g}(\mathbf{y}(t), t) \]

Korrektor (Mittelwert der Steigungen): \[ \mathbf{y}(t + \Delta t) = \mathbf{y}(t) + \frac{\Delta t}{2} \left[ \mathbf{g}(\mathbf{y}(t), t) + \mathbf{g}(\tilde{\mathbf{y}}(t + \Delta t), t + \Delta t) \right] \]

• Konvergenzordnung 2: Halbierung von \(\Delta t\) → Fehler um Faktor 4 kleiner

Konvergenzordnung des Heun-Verfahrens

Taylor-Entwicklung der exakten Lösung (\(y' = g\), Kettenregel): \[ y(t + \Delta t) = y + \Delta t \, g + \frac{\Delta t^2}{2} \left( g \frac{\partial g}{\partial y} + \frac{\partial g}{\partial t} \right) + \mathcal{O}(\Delta t^3) \]

Entwicklung der Heun-Korrektorformel liefert exakt denselben Ausdruck: \[ y_{i+1} = y + \Delta t \, g + \frac{\Delta t^2}{2} \left( g \frac{\partial g}{\partial y} + \frac{\partial g}{\partial t} \right) + \mathcal{O}(\Delta t^3) \]

• Übereinstimmung bis einschließlich \(\mathcal{O}(\Delta t^2)\)
• Lokaler Abbruchfehler \(\mathcal{O}(\Delta t^3)\) → Verfahren zweiter Ordnung

Euler vs. Heun: harmonischer Oszillator

• Euler spiralt nach außen, Heun folgt dem Kreis deutlich besser

Heun vs. Leapfrog

Beide: Ordnung 2, Prädiktor-Korrektor-Struktur — aber:

• Runge-Kutta (Heun): universeller Integrator für beliebige \(\mathbf{y}' = \mathbf{g}(\mathbf{y}, t)\) — keine Annahmen über die Struktur
• Leapfrog: spezialisiert auf \(\ddot{x} = f(x)\) (konservative mechanische Systeme, Hamilton-Systeme)
• Leapfrog ist symplektisch → keine systematische Energiedrift
• Heun ist nicht symplektisch → spiralt im Phasenraum, nur langsamer als Euler
• Hamiltonsche Systeme: Leapfrog; allgemeine, nicht-konservative Probleme: Runge-Kutta

Klassisches Runge-Kutta-Verfahren (RK4)

Gewichteter Mittelwert von vier Steigungsschätzungen: \[ \mathbf{y}_{i+1} = \mathbf{y}_i + \frac{\Delta t}{6} \left[ \mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4 \right] \]

\[ \begin{aligned} \mathbf{k}_1 &= \mathbf{g}(\mathbf{y}_i, t_i) \\ \mathbf{k}_2 &= \mathbf{g}\left(\mathbf{y}_i + \tfrac{\Delta t}{2} \mathbf{k}_1, t_i + \tfrac{\Delta t}{2}\right) \\ \mathbf{k}_3 &= \mathbf{g}\left(\mathbf{y}_i + \tfrac{\Delta t}{2} \mathbf{k}_2, t_i + \tfrac{\Delta t}{2}\right) \\ \mathbf{k}_4 &= \mathbf{g}\left(\mathbf{y}_i + \Delta t \, \mathbf{k}_3, t_i + \Delta t\right) \end{aligned} \]

RK4: Interpretation der Stufen

• \(\mathbf{k}_1\): Steigung am Anfang (wie Euler)
• \(\mathbf{k}_2\): Steigung in der Mitte, Position über halben Euler-Schritt mit \(\mathbf{k}_1\)
• \(\mathbf{k}_3\): Steigung in der Mitte, mit verbesserter Schätzung aus \(\mathbf{k}_2\)
• \(\mathbf{k}_4\): Steigung am Ende, voller Schritt mit \(\mathbf{k}_3\)
• Gewichte \(1/6, 2/6, 2/6, 1/6\): minimieren den Fehler optimal
• Lokaler Abbruchfehler \(\mathcal{O}(\Delta t^5)\) → Verfahren vierter Ordnung

Verwandtschaft mit der Simpson-Regel

\[ \int_a^b f(x) \, \mathrm{d}x \approx \frac{b-a}{6} \left[ f(a) + 4 f\left(\frac{a+b}{2}\right) + f(b) \right] \]

• Mittelpunkt vierfach gewichtet → parabolische Interpolation
• Im RK4: Mittelwerte \(\mathbf{k}_2\), \(\mathbf{k}_3\) tragen zusammen Gewicht \(4/6\), Randwerte \(\mathbf{k}_1\), \(\mathbf{k}_4\) je \(1/6\)
• Zwei Mittelpunktschätzungen statt einer → weitere Fehlerreduktion

Konvergenz und Effizienz im Vergleich

• RK4: Halbierung von \(\Delta t\) → Fehler um Faktor 16 kleiner
• Trotz 4 Auswertungen pro Schritt: bei gleichem Aufwand um Größenordnungen genauer

Die Funktion solve_ivp

scipy.integrate.solve_ivp: robuste Implementierung von Runge-Kutta-Verfahren mit automatischer Schrittweitenkontrolle

Methode	Ordnung	Bemerkung
RK23	3(2)	für weniger anspruchsvolle Probleme
RK45	5(4)	Standard (Dormand-Prince)
DOP853	8(5,3)	höchste Genauigkeit

• Notation „5(4)“: propagierte Lösung der Ordnung 5 mit eingebettetem Fehlerschätzer der Ordnung 4 → Grundlage der adaptiven Zeitschrittsteuerung

Struktur eines solve_ivp-Aufrufs

from scipy.integrate import solve_ivp

def dgl(t, y, param1, param2):
    # WICHTIG: t ist das erste Argument, y das zweite!
    dydt = ...
    return dydt

solution = solve_ivp(
    dgl,                      # Die DGL-Funktion
    t_span=(t_start, t_end),  # Zeitintervall
    y0=[y1_0, y2_0],          # Anfangsbedingungen
    args=(param1, param2),    # Zusätzliche Parameter
    method='RK45',            # Integrationsverfahren
    rtol=1e-3, atol=1e-6      # Toleranzen
)
t = solution.t   # Zeitpunkte
y = solution.y   # Lösungswerte

• Achtung: solve_ivp erwartet dgl(t, y) — Zeit vor Zustand!

Beispiel: gedämpfter harmonischer Oszillator

\[ m \ddot{y} + c \dot{y} + k y = 0 \qquad\Rightarrow\qquad \dot{y} = v, \quad \dot{v} = -\frac{c}{m} v - \frac{k}{m} y \]

Verhalten abhängig von der kritischen Dämpfung \(c_\mathrm{krit} = 2\sqrt{km}\):

• \(c < c_\mathrm{krit}\): gedämpfte Schwingung
• \(c = c_\mathrm{krit}\): aperiodischer Grenzfall
• \(c > c_\mathrm{krit}\): überdämpft, kriecht zum Gleichgewicht

Gedämpfter Oszillator mit solve_ivp

Anwendung: Lotka-Volterra-Gleichungen

Idealisiertes Räuber-Beute-System (\(x\): Beute, \(y\): Räuber): \[ \frac{\mathrm{d}x}{\mathrm{d}t} = \alpha x - \beta xy, \qquad \frac{\mathrm{d}y}{\mathrm{d}t} = \delta xy - \gamma y \]

• \(\alpha\): Wachstumsrate der Beute, \(\beta\): Fressrate
• \(\delta\): Umwandlungseffizienz, \(\gamma\): Sterberate der Räuber
• Klassisches Beispiel eines nichtlinearen Systems

Lotka-Volterra: numerische Lösung

• Periodische Oszillationen: Räubermaxima folgen den Beutemaxima zeitversetzt
• Phasenraum: geschlossene Kurve um den Gleichgewichtspunkt

Zusammenfassung

• Runge-Kutta-Verfahren: mehrere Funktionsauswertungen pro Zeitschritt → höhere Ordnung als Euler (= RK-Verfahren 1. Ordnung)
• Heun (RK2): Prädiktor-Korrektor mit 2 Auswertungen, quadratische Konvergenz
• RK4: vier Steigungen, gewichtet nach der Simpson-Regel, lokaler Fehler \(\mathcal{O}(\Delta t^5)\) — dramatisch genauer bei moderatem Mehraufwand
• solve_ivp mit Standardmethode RK45: robuste Wahl für die meisten praktischen Anwendungen (adaptive Schrittweite)
• Höhere Ordnung lohnt sich: bei gleichem Rechenaufwand um Größenordnungen kleinerer Fehler