Drei Methoden, drei Bereiche
Für die Testgleichung dy/dx = λy haben drei explizite ODE-Methoden die folgenden Stabilitätsbereiche in der komplexen hλ-Ebene:
Eulersche Methode (erster Ordnung): Der Stabilitätsbereich ist die Scheibe |1 + hλ| ≤ 1, ein Kreis mit Radius 1, zentriert bei (-1, 0). Das reelle negative hλ muss in [-2, 0] liegen.
Runge-Kutta 2 (Mittelpunkt-Methode) (zweiter Ordnung): Der Stabilitätsbereich ist |1 + hλ + (hλ)²/2| ≤ 1. Größer als die Scheibe von Euler, aber immer noch beschränkt.
Runge-Kutta 4 (vierte Ordnung): Der Stabilitätsbereich erfüllt |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Das reelle negative hλ reicht bis etwa -2,785. Der Bereich ist erheblich größer als der von Euler.
Rückwärts-Euler (implizit): Der Stabilitätsbereich ist die gesamte komplexe Ebene außer der Scheibe |1 - hλ|⁻¹ > 1, äquivalent |1/(1-hλ)| ≤ 1. Für λ in der linken Halbebene (Re(λ) < 0) ist dies unbedingt stabil — keine Beschränkung auf h aus Stabilitätsgründen.
Die Verstärkungsfunktion
Für jede Runge-Kutta-Methode ist der Verstärkungsfaktor pro Schritt R(hλ) eine Polynomapproximation von e^(hλ):
- Euler: R(z) = 1 + z (bei Grad 1 gekürzt)
- RK2: R(z) = 1 + z + z²/2 (bei Grad 2 gekürzt)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (bei Grad 4 gekürzt)
Der Stabilitätsbereich ist {z ∈ ℂ : |R(z)| ≤ 1}. Die Verstärkung der wahren Lösung: |e^z| = e^(Re(z)). Für Re(z) < 0 (stabile ODE) nimmt die wahre Lösung ab. Die numerische Methode ist stabil, wenn |R(z)| ≤ 1 — was dem Abnahmeverhalten entspricht.
Rein imaginäre Eigenwerte: Oszillierende Systeme
Viele physikalische Systeme haben rein imaginäre Eigenwerte: λ = iω (Schwingungen ohne Dämpfung). Das Feder-Masse-System, Bahnmechanik, Pendeldynamik.
Für λ = iω: hλ = ihω liegt auf der imaginären Achse.
Eulers Stabilität auf der imaginären Achse: |1 + ihω|² = 1 + (hω)² > 1 für alle h > 0. Euler ist instabil für jede Schrittweite bei rein imaginären Eigenwerten. Die berechnete 'Schwingung' wächst ohne Grenzen.
RK4s Stabilität auf der imaginären Achse: Der Stabilitätsbereich erstreckt sich auf der imaginären Achse auf ungefähr |hω| ≤ 2,83. Für hinreichend kleine h bewältigt RK4 ungedämpfte Schwingungen. Euler kann das nicht.
Dieses geometrische Merkmal erklärt, warum Euler bei konservativen Systemen (Feder-Masse, Bahnen, Wellengleichungen) auch mit kleinem h fehlschlägt, während RK4 sie gut bewältigt.
Die Geometrie steifer Probleme
Ein steifes ODE-System hat Eigenwerte mit sehr unterschiedlichen Größen. Die Steifheitsquote: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Warum Steifheit teuer für explizite Lösungsmethoden ist:
Die Stabilität erfordert h·max|λᵢ| ≤ C (wobei C von der Methode abhängt). Der negativste Eigenwert legt die Schranke fest.
Die Genauigkeit für die langsame Dynamik erfordert h·min|λᵢ| ≥ ε (die langsamste Mode angemessen auflösen).
Wenn κ groß ist, zwingen diese beiden Anforderungen ein winziges h: klein genug für die Stabilität der schnellen Mode, groß genug, um die langsame Mode zu erfassen. Die Anzahl der Schritte skaliert mit κ.
Geometrisches Bild im Eigenwertspektrum: Die Eigenwerte der Jacobi-Matrix ∂f/∂y bilden eine Menge von Punkten in der komplexen Ebene. Der Stabilitätsbereich eines expliziten Lösungsmittels muss alle Punkte h·λᵢ enthalten. Wenn Eigenwerte von -1 bis -1000 reichen, muss der Stabilitätsbereich einen Bereich von 1000 entlang der realen Achse abdecken — es ist ein sehr kleines h erforderlich.
Implizite Lösungsmethoden: Der Stabilitätsbereich von Rückwärts-Euler deckt die gesamte linke Halbebene ab. Alle Eigenwerte mit Re(λ) < 0 befinden sich automatisch im Stabilitätsbereich, unabhängig von h. Die Beschränkung von h ergibt sich nur aus Genauigkeit, nicht aus Stabilität.
Steifheitsquote & Kosten
Betrachte ein chemisches Reaktionsnetzwerk mit schnellen Reaktionen (Zeitskala 10⁻⁶ s) und langsamen Reaktionen (Zeitskala 1 s).
Steifheitsquote: κ = 10⁶ / 1 = 10⁶.
Mit RK4 (Stabilitätsgrenze h·|λ_max| ≤ 2,785): h_max = 2,785 / 10⁶ ≈ 2,8 × 10⁻⁶ s.
Um über 10 s Reaktionszeit zu integrieren: steps = 10 / (2,8 × 10⁻⁶) ≈ 3,6 × 10⁶.
Mit Rückwärts-Euler (unbedingt stabil): h kann für die Genauigkeit der langsamen Reaktionen gewählt werden. h = 10⁻² s (100 Proben über die 1-s-Skala). Steps = 10 / 10⁻² = 1000.
Kostenverhältnis: explizite 3,6 Millionen Schritte gegen implizite 1000 Schritte — ein Faktor von 3600. Jeder implizite Schritt erfordert das Lösen eines linearen Systems (Kosten pro Schritt sind höher), aber die Gesamtkosten sind für sehr steife Probleme viel niedriger.
Warum n-dimensionale Röhren nicht das sind, was man denkt
In 2D ist eine 'Röhre' mit Radius ε um eine Kurve C die Menge der Punkte innerhalb der Entfernung ε von C. Der Querschnitt ist ein Kreis mit Radius ε. Das Volumen der Röhre wächst proportional mit ihrer Länge.
In n Dimensionen ändert sich die Röhrengeometrie grundlegend aufgrund des Phänomens aus Kapitel 9:
Das n-dimensionale Eckenparadoxon: In n-dimensionalem Raum liegt fast das gesamte Volumen eines n-dimensionalen Hyperwürfels in den Ecken — nicht in der zentralen Region. Mit zunehmendem n geht der Volumenanteil in der Entfernung ε vom Mittelpunkt gegen null für jedes feste ε.
Angewendet auf ODE-Lösungsröhren:
In 2D: Wenn die wahre Lösung durch den Mittelpunkt einer Röhre verläuft, befinden sich die meisten nahegelegenen Punkte nah bei der Kurve. Kleine Störungen halten dich nah bei der wahren Lösung.
In hohen Dimensionen: Die meisten Punkte innerhalb der Begrenzungsbox der Röhre sind tatsächlich weit weg von der echten Lösungskurve. Das 'Volumen' der Röhre wird von den Ecken dominiert — Regionen, die in mehreren Dimensionen gleichzeitig weit vom Mittelpunkt entfernt sind.
Folgerung für die Simulation: Mit 28 gekoppelten ODEs (Hamming's Navy-Abfangproblem) kann eine Störung der Größe ε in jeder Dimension eine Gesamtverschiebung von ε√28 ≈ 5,3ε von der wahren Lösung erzeugen. Die Röhre muss in Bezug auf die L2-Norm über alle Dimensionen verstanden werden, nicht nur als maximale Verschiebung in einer Dimension.
Stabilität in hohen Dimensionen: Ein System, bei dem jede Komponente unabhängig abklingt (jeder Eigenwert hat einen negativen Realteil), kann dennoch große kombinierte Verschiebungen zeigen, weil sich die Komponenten fehler in der L2-Norm addieren. Die 28-dimensionale Röhre ist nicht nur 28 unabhängige 1-dimensionale Röhren — die Geometrie koppelt sie.
Von Geometrie zum Design
Die geometrischen Erkenntnisse der Kapitel 18-20 kommen zusammen als eine Reihe von Designprinzipien für numerische Simulation:
Schrittweiteauswahl: h muss h·λ innerhalb des Stabilitätsbereichs für jeden Eigenwert platzieren. Bei steifen Systemen entfernen implizite Methoden die Stabilitätsbeschränkung und hinterlassen nur Genauigkeitsanforderungen.
Fehlerakkumulation in hohen Dimensionen: Der globale Fehler ist ein Vektor im n-dimensionalen Raum. Seine Norm wächst um das √n-fache des pro-Komponenten-Fehlers. Hochdimensionale Simulationen benötigen strengere Genauigkeitsanforderungen pro Schritt.
Rückkopplung als Stabilisator: Wenn die Simulation Rückkopplung einbezieht (die berechnete Ausgabe beeinflusst nachfolgende Eingaben, wie in einem Leitungssystem), dämpft konvergente Rückkopplung Fehler. Die Simulation kann ungenaue Eingaben für Größen in der Rückkopplungsschleife tolerieren.
Die Instabilität als Signal: Bei Problemen mit divergierenden Richtungsfeldern kann Instabilität ausgenutzt werden: Die Divergenzrichtung trägt Informationen über den Anfangsbedingungsfehler und ermöglicht korrigierende Anpassung.