Euler-Integration

Differentialgleichungen – Kapitel 11

M. Ladecký, L. Pastewka

IMTEK, Universität Freiburg

Lernziele

Die Studierenden sollen…

• … das explizite Euler-Verfahren als einfachstes numerisches Integrationsverfahren beschreiben können.
• … die Herleitung des Euler-Verfahrens aus der Annahme konstanter Steigung erklären können.
• … die Energiedrift des Euler-Verfahrens im Phasenraum erklären können.
• … das Euler-Verfahren auf skalare und vektorwertige Differentialgleichungen anwenden können.
• … das Leapfrog-Verfahren zur verbesserten Energieerhaltung anwenden können.

Motivation

• Nur ein kleiner Teil von Differentialgleichungen ist analytisch lösbar
• In der Praxis: numerischer Zugang erforderlich
• Einfachstes Verfahren: explizites Euler-Verfahren (Leonhard Euler)
• Zentrale Idee: Die DGL gibt die Steigung der Lösungskurve an → Lösung näherungsweise entlang der Tangente fortsetzen

Heuristische Herleitung

Anfangswertproblem: \[ \frac{\mathrm{d}y}{\mathrm{d}t} = g(y, t), \quad y(t_0) = y_0 \]

• Steigung der Lösungskurve am Punkt \((t_i, y_i)\): \(g(y_i, t_i)\)
• Annahme: Steigung bleibt über kleinen Zeitschritt \(\Delta t\) konstant
• Lineare Extrapolation zum nächsten Zeitpunkt \(t_{i+1} = t_i + \Delta t\):

\[ y_{i+1} = y_i + \Delta t \cdot g(y_i, t_i) \]

Das explizite Euler-Verfahren

\[ y_{i+1} = y_i + \Delta t \cdot g(y_i, t_i) \]

• “Explizit”: \(y_{i+1}\) folgt direkt aus den bekannten Größen \(y_i\), \(t_i\)
• Geometrische Interpretation: Schritt entlang der Tangente an die Lösungskurve
• Wahre Lösungskurve ist gekrümmt → Tangente weicht ab → Fehler

Das implizite Euler-Verfahren

Alternative: Steigung am Ende des Intervalls verwenden: \[ y_{i+1} = y_i + \Delta t \cdot g(y_{i+1}, t_{i+1}) \]

• Rechte Seite hängt von der unbekannten Größe \(y_{i+1}\) ab → in jedem Schritt eine (i. A. nichtlineare) Gleichung lösen
• Dafür: unbedingt stabil — auch bei großen Zeitschritten
• Für steife DGLs (sehr unterschiedliche Zeitskalen) oft die einzige praktikable Wahl

Skalares Beispiel: Relaxationsgleichung

\[ \frac{\mathrm{d}y}{\mathrm{d}t} = -\frac{y}{\tau}, \quad y(0) = y_0 \qquad\Rightarrow\qquad y(t) = y_0 \exp\left(-\frac{t}{\tau}\right) \]

Euler-Verfahren liefert die Rekursion: \[ y_{i+1} = \left(1 - \frac{\Delta t}{\tau}\right) y_i \quad\Rightarrow\quad y_n = \left(1 - \frac{\Delta t}{\tau}\right)^n y_0 \]

• Für \(\Delta t \ll \tau\): \((1 - \Delta t/\tau)^n \approx \exp(-t_n/\tau)\) — exakte Lösung
• Für größere \(\Delta t\): zunehmende Abweichung

Euler-Verfahren: Konvergenz

• Fehler skaliert linear mit \(\Delta t\) → Verfahren erster Ordnung

Vektorwertiges Problem: Feder-Masse-System

Bewegungsgleichung (2. Ordnung): \[ m \frac{\mathrm{d}^2 x}{\mathrm{d}t^2} = -k x \]

Umschreiben in System 1. Ordnung mit \(v = \mathrm{d}x/\mathrm{d}t\), \(\omega = \sqrt{k/m}\): \[ \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \mathbf{g}(\mathbf{y}, t), \quad \mathbf{y} = \begin{pmatrix} x \\ v \end{pmatrix}, \quad \mathbf{g}(\mathbf{y}, t) = \begin{pmatrix} v \\ -\omega^2 x \end{pmatrix} \]

Euler-Verfahren für Systeme

Gleiche Struktur wie im skalaren Fall: \[ \mathbf{y}_{i+1} = \mathbf{y}_i + \Delta t \cdot \mathbf{g}(\mathbf{y}_i, t_i) \]

Konkret für das Feder-Masse-System: \[ x_{i+1} = x_i + \Delta t \cdot v_i, \qquad v_{i+1} = v_i - \Delta t \cdot \omega^2 x_i \]

• Exakte Lösung: geschlossene Kreisbahn im Phasenraum \((x, v)\)
• Wie verhält sich das Euler-Verfahren?

Euler-Verfahren im Phasenraum

• Euler-Verfahren spiralt nach außen: dem System wird künstlich Energie zugeführt
• Effekt umso stärker, je größer \(\Delta t\) — numerische Energiedrift

Energiedrift: quantitative Analyse

Totale Energie: \[ E = \frac{1}{2} m v^2 + \frac{1}{2} k x^2 \]

Einsetzen der Euler-Formeln liefert nach einem Zeitschritt: \[ E_{i+1} = E_i \left( 1 + \omega^2 (\Delta t)^2 \right) \quad\Rightarrow\quad E_n = E_0 \left( 1 + \omega^2 (\Delta t)^2 \right)^n \]

• Für kleine \(\Delta t\): \(E_n \approx E_0 \exp(\omega^2 \Delta t \cdot t)\)
• Die Energie wächst exponentiell mit der Zeit!

Energiedrift: numerische Bestätigung

• Wachstumsrate skaliert linear mit \(\Delta t\) — wie theoretisch vorhergesagt

Das Leapfrog-Verfahren

• Auch Störmer-Verlet-Verfahren genannt
• Symplektisches Integrationsverfahren: deutlich bessere Energieerhaltung
• Prädiktor-Korrektor-Struktur: symmetrische Kombination zweier halber Euler-Schritte
• Drei Teilschritte:
  1. Prädiktor: halber Geschwindigkeitsschritt
  2. voller Positionsschritt
  3. Korrektor: halber Geschwindigkeitsschritt

Leapfrog: die drei Teilschritte

Prädiktor (halber Geschwindigkeitsschritt mit \(x_i\)): \[ v_{i+1/2} = v_i - \frac{\Delta t}{2} \omega^2 x_i \]

Positionsschritt (mit Zwischengeschwindigkeit): \[ x_{i+1} = x_i + \Delta t \cdot v_{i+1/2} \]

Korrektor (halber Geschwindigkeitsschritt mit \(x_{i+1}\)): \[ v_{i+1} = v_{i+1/2} - \frac{\Delta t}{2} \omega^2 x_{i+1} \]

Leapfrog: kompakte Darstellung

Elimination der Zwischengeschwindigkeit \(v_{i+1/2}\): \[ v_{i+1} = v_i - \frac{\Delta t}{2} \omega^2 (x_i + x_{i+1}) \]

• Mittelwert der Beschleunigungen am Anfang und Ende des Intervalls
• Entspricht der Trapezregel für numerische Integration
• Symmetrische Struktur → Zeitumkehrsymmetrie des Verfahrens

Symplektizität des Leapfrog-Verfahrens

Leapfrog-Schritt in Matrixform mit \(\lambda = (\omega \Delta t)^2\): \[ \begin{pmatrix} x_{i+1} \\ v_{i+1} \end{pmatrix} = \mathbf{M} \begin{pmatrix} x_i \\ v_i \end{pmatrix}, \quad \mathbf{M} = \begin{pmatrix} 1 - \frac{\lambda}{2} & \Delta t \\ -\omega^2 \Delta t \left(1 - \frac{\lambda}{4}\right) & 1 - \frac{\lambda}{2} \end{pmatrix} \]

\[ \det(\mathbf{M}) = \left(1 - \frac{\lambda}{2}\right)^2 + \lambda \left(1 - \frac{\lambda}{4}\right) = 1 \]

• \(\det(\mathbf{M}) = 1\) → flächenerhaltend im Phasenraum \((x, v)\)
• In 2D äquivalent zur Symplektizität → Langzeitstabilität
• Energie oszilliert um den wahren Wert statt systematisch zu driften

Euler vs. Leapfrog im Langzeitverhalten

• Euler: Energie nach \(t=100\) um mehr als das Hundertfache aufgebläht
• Leapfrog: Energie oszilliert leicht, bleibt aber im Mittel auf dem korrekten Niveau

Zusammenfassung

• Explizites Euler-Verfahren: \(y_{i+1} = y_i + \Delta t \cdot g(y_i, t_i)\) — einfachstes Verfahren, basiert auf konstanter Steigung über \(\Delta t\)
• Erste Ordnung: Fehler skaliert linear mit dem Zeitschritt; für dissipative Systeme bei kleinem \(\Delta t\) akzeptabel
• Energiedrift: Bei konservativen Systemen wächst die Energie exponentiell — Spirale nach außen im Phasenraum
• Leapfrog-Verfahren: symplektisch (\(\det \mathbf{M} = 1\)), erhält die Phasenraumstruktur, keine systematische Energiedrift
• Für Langzeitintegration konservativer Systeme: symplektische Verfahren statt Euler verwenden