Zeitschrittkontrolle

Differentialgleichungen – Kapitel 13

M. Ladecký, L. Pastewka

IMTEK, Universität Freiburg

Lernziele

Die Studierenden sollen…

• … das Konzept der automatischen Zeitschrittkontrolle beschreiben können.
• … das Prinzip eingebetteter Runge-Kutta-Verfahren zur Fehlerschätzung erklären können.
• … den Zusammenhang zwischen Toleranzparametern und Genauigkeit erklären können.
• … ein einfaches Runge-Kutta-Verfahren mit adaptiver Schrittweite implementieren können.
• … die Toleranzparameter von solve_ivp sinnvoll wählen können.

Motivation

• Bisher: Runge-Kutta-Verfahren mit konstantem Zeitschritt \(\Delta t\)
• Problem: Lösung ändert sich oft zeitweise schnell, zeitweise langsam → fester Zeitschritt muss für die schnellste Phase klein genug sein → in langsamen Phasen verschwenderisch
• Automatische Zeitschrittkontrolle: \(\Delta t\) wird dynamisch an die lokale Dynamik angepasst
• Zentrales Element: Fehlerschätzung — eingebettete Runge-Kutta-Verfahren liefern sie ohne zusätzliche Funktionsauswertungen

Eingebettete Runge-Kutta-Verfahren

Idee: Aus denselben Steigungen \(\mathbf{k}_1, \ldots, \mathbf{k}_s\) zwei Lösungen unterschiedlicher Ordnung berechnen: \[ \mathbf{y}_{n+1} = \mathbf{y}_n + \Delta t \sum_{i=1}^s b_i \mathbf{k}_i, \qquad \hat{\mathbf{y}}_{n+1} = \mathbf{y}_n + \Delta t \sum_{i=1}^s \hat{b}_i \mathbf{k}_i \]

Die Differenz dient als Schätzung des lokalen Fehlers: \[ \mathbf{e} = \mathbf{y}_{n+1} - \hat{\mathbf{y}}_{n+1} = \Delta t \sum_{i=1}^s (b_i - \hat{b}_i) \mathbf{k}_i \]

• Höherwertige Lösung wird für die Zeitpropagation verwendet
• Fehlerschätzung hat die niedrigere Ordnung → konservative Schätzung

Das RK23-Verfahren (Bogacki-Shampine)

Eingebettetes Verfahren der Ordnung 3(2) mit vier Stufen: \[ \begin{aligned} \mathbf{k}_1 &= \mathbf{g}(\mathbf{y}_n, t_n) \\ \mathbf{k}_2 &= \mathbf{g}\left(\mathbf{y}_n + \tfrac{\Delta t}{2} \mathbf{k}_1, t_n + \tfrac{\Delta t}{2}\right) \\ \mathbf{k}_3 &= \mathbf{g}\left(\mathbf{y}_n + \tfrac{3\Delta t}{4} \mathbf{k}_2, t_n + \tfrac{3\Delta t}{4}\right) \\ \mathbf{k}_4 &= \mathbf{g}(\mathbf{y}_{n+1}, t_{n+1}) \end{aligned} \]

• \(\mathbf{k}_4\) kann im nächsten Schritt wiederverwendet werden → effektiv nur drei neue Funktionsauswertungen pro Schritt (FSAL)

RK23: Lösung und Fehlerschätzer

Lösung dritter Ordnung (für die Zeitpropagation): \[ \mathbf{y}_{n+1} = \mathbf{y}_n + \frac{\Delta t}{9} \left( 2\mathbf{k}_1 + 3\mathbf{k}_2 + 4\mathbf{k}_3 \right) + \mathcal{O}(\Delta t^4) \]

Eingebettete Lösung zweiter Ordnung (für die Fehlerschätzung): \[ \hat{\mathbf{y}}_{n+1} = \mathbf{y}_n + \frac{\Delta t}{24} \left( 7\mathbf{k}_1 + 6\mathbf{k}_2 + 8\mathbf{k}_3 + 3\mathbf{k}_4 \right) + \mathcal{O}(\Delta t^3) \]

Fehlerschätzer aus der Differenz: \[ \mathbf{e} = \frac{\Delta t}{72} \left( -5\mathbf{k}_1 + 6\mathbf{k}_2 + 8\mathbf{k}_3 - 9\mathbf{k}_4 \right) + \mathcal{O}(\Delta t^4) \]

Schrittweitensteuerung: Toleranz

Vergleich des Fehlers \(\mathbf{e}\) mit einer benutzerdefinierten Toleranz: \[ \epsilon_\mathrm{tol} = \mathrm{atol} + \mathrm{rtol} \cdot |\mathbf{y}_{n+1}| \]

Skalierter Fehler (euklidische oder Maximumsnorm): \[ \epsilon = \left\| \frac{\mathbf{e}}{\epsilon_\mathrm{tol}} \right\| \]

• \(\epsilon \le 1\): Schritt akzeptiert
• \(\epsilon > 1\): Schritt verworfen, Wiederholung mit kleinerer Schrittweite
• In beiden Fällen: Schrittweite für den nächsten Schritt anpassen

Schrittweitensteuerung: Anpassung

Fehler eines Verfahrens der Ordnung \(p\) skaliert wie \(\Delta t^{p+1}\) → optimale neue Schrittweite: \[ \Delta t_\mathrm{neu} = \Delta t_\mathrm{alt} \cdot \epsilon^{-1/(p+1)} \]

In der Praxis mit Sicherheitsfaktor \(S < 1\) (typisch \(S = 0{,}9\)): \[ \Delta t_\mathrm{neu} = S \cdot \Delta t_\mathrm{alt} \cdot \epsilon^{-1/(p+1)} \]

• Sicherheitsfaktor vermeidet häufige Schrittablehnungen
• Zusätzlich: Begrenzung der Schrittweitenänderung pro Schritt

Testproblem: Van-der-Pol-Oszillator

Klassischer nichtlinearer Oszillator: \[ \ddot{x} - \mu(1 - x^2)\dot{x} + x = 0 \qquad\Rightarrow\qquad \dot{x} = v, \quad \dot{v} = \mu(1 - x^2)v - x \]

• \(|x| < 1\): negative Dämpfung — treibt das System an
• \(|x| > 1\): positive Dämpfung — bremst
• Ergebnis: stabiler Grenzzyklus mit Relaxationsoszillationen (langsame Phasen + schnelle Übergänge) — ideal für adaptive Verfahren

Adaptives RK23 in Aktion

• Schnelle Phasen: Schrittweite drastisch reduziert; langsame Phasen: große Schritte
• Skalierter Fehler bleibt stets unter der Toleranzgrenze 1

Vergleich mit solve_ivp

• solve_ivp verwendet standardmäßig RK45 (eingebettetes Paar 5(4))
• Prinzip der Schrittweitensteuerung identisch zum selbstimplementierten RK23
• Eigene RK23-Implementierung liefert vergleichbare Lösungen und ähnliche Schrittanzahl wie SciPy RK23
• RK45 benötigt weniger Schritte als RK23: höhere Ordnung erlaubt größere Zeitschritte
• Für die Praxis: professionell getestete Implementierung in SciPy verwenden

Einfluss der Toleranzwahl

• Kleinere Toleranz → mehr Schritte, aber genauere Ergebnisse
• Unterhalb von \(\sim 10^{-14}\): Rundungsfehler dominieren

Empfehlungen für die Toleranzwahl

• rtol steuert die relative Genauigkeit, atol die Genauigkeit für Werte nahe null
• Schnelle Visualisierungen, explorative Rechnungen: rtol=1e-3, atol=1e-6
• Produktionsrechnungen: rtol=1e-6, atol=1e-9
• Höchste Genauigkeit (deutlich höherer Aufwand): rtol=1e-10, atol=1e-13
• Grenze: Maschinengenauigkeit bei \(\sim 10^{-14}\)

Steife Differentialgleichungen

System mit sehr unterschiedlichen Zeitskalen, z. B.: \[ \dot{y}_1 = -1000\, y_1 + y_2, \qquad \dot{y}_2 = y_1 - y_2 \]

• Eigenwerte \(\lambda_1 \approx -1000\), \(\lambda_2 \approx -1\) → Zeitkonstanten \(\tau_1 = 0{,}001\) und \(\tau_2 = 1\) (Verhältnis 1000:1)
• Explizite Verfahren (RK23, RK45): sehr kleine Schritte nötig, um stabil zu bleiben — ineffizient
• Abhilfe: implizite Verfahren wie Radau oder BDF

Explizit vs. implizit beim steifen Problem

• Radau und BDF benötigen bei gleicher Genauigkeit deutlich weniger Funktionsauswertungen
• Faustregel: Wenn RK45 ungewöhnlich viele Schritte braucht → Problem ist vermutlich steif

Zusammenfassung

• Automatische Zeitschrittkontrolle: Schrittweite passt sich der lokalen Dynamik an — optimale Balance aus Genauigkeit und Aufwand
• Eingebettete Verfahren (RK23, RK45): zwei Lösungen unterschiedlicher Ordnung aus denselben Auswertungen; Differenz = lokaler Fehlerschätzer
• Schrittweitenanpassung: \(\Delta t_\mathrm{neu} = S \cdot \Delta t_\mathrm{alt} \cdot \epsilon^{-1/(p+1)}\)
• Toleranzen rtol/atol steuern Genauigkeit vs. Rechenaufwand
• Steife Probleme: implizite Verfahren (Radau, BDF) statt expliziter Runge-Kutta-Verfahren verwenden