Drie Methoden, Drie Regio's
Voor de testvergelijking dy/dx = λy hebben drie expliciete ODE-methoden de volgende stabiliteitszones in het complexe hλ-vlak:
Eulermethode (eerste-orde): stabiliteitzone is de schijf |1 + hλ| ≤ 1, een cirkel met straal 1 gecentreerd op (-1, 0). Reële negatieve hλ moet in [-2, 0] liggen.
Runge-Kutta 2 (midpoint-methode) (tweede-orde): stabiliteitzone is |1 + hλ + (hλ)²/2| ≤ 1. Groter dan de schijf van Euler, maar nog steeds begrensd.
Runge-Kutta 4 (vierde-orde): stabiliteitzone voldoet aan |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Reële negatieve hλ strekt zich uit tot ongeveer -2.785. De zone is aanzienlijk groter dan die van Euler.
Achterwaartse Euler (impliciet): stabiliteitzone is het gehele complexe vlak behalve de schijf |1 - hλ|⁻¹ > 1, equivalent |1/(1-hλ)| ≤ 1. Voor λ in het linkerhalfvlak (Re(λ) < 0) is dit onvoorwaardelijk stabiel — geen beperking op h vanuit stabiliteit.
De Versterkingsfunctie
Voor elke Runge-Kutta-methode is de versterkingsfactor per stap R(hλ) een polynoomapproximatie van e^(hλ):
- Euler: R(z) = 1 + z (afgekapt op graad 1)
- RK2: R(z) = 1 + z + z²/2 (afgekapt op graad 2)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (afgekapt op graad 4)
De stabiliteitzone is {z ∈ ℂ : |R(z)| ≤ 1}. De werkelijke oplossingsversterkingsfactor: |e^z| = e^(Re(z)). Voor Re(z) < 0 (stabiele ODE) daalt de werkelijke oplossing. De numerieke methode is stabiel als |R(z)| ≤ 1 — overeenkomend met het dalende gedrag.
Zuiver Denkbeeldige Eigenwaarden: Oscillerende Systemen
Veel fysieke systemen hebben zuiver denkbeeldige eigenwaarden: λ = iω (oscillaties zonder demping). Het veer-massasysteem, orbitaalmecanica, slingerynamica.
Voor λ = iω: hλ = ihω ligt op de denkbeeldige as.
Stabiliteit van Euler op de denkbeeldige as: |1 + ihω|² = 1 + (hω)² > 1 voor alle h > 0. Euler is onstabiel voor elke stapgrootte op zuiver denkbeeldige eigenwaarden. De berekende 'oscillatie' groeit zonder grens.
Stabiliteit van RK4 op de denkbeeldige as: de stabiliteitzone strekt zich uit tot ongeveer |hω| ≤ 2.83 op de denkbeeldige as. Voor voldoende kleine h handelt RK4 ongedempte oscillaties af. Euler kan dit niet.
Deze geometrie is waarom Euler faalt op conservatieve systemen (veer-massa, banen, golfvergelijkingen) zelfs met kleine h, terwijl RK4 ze goed afhandelt.
De Geometrie van Stijve Problemen
Een stijf ODE-systeem heeft eigenwaarden met zeer verschillende magnitudes. De stijfheidsverhouding: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Waarom stijfheid duur is voor expliciete oplosvers:
Stabiliteit vereist h·max|λᵢ| ≤ C (waarbij C afhangt van de methode). De meest negatieve eigenwaarde bepaalt de grens.
Nauwkeurigheid voor de trage dynamica vereist h·min|λᵢ| ≥ ε (los de traagste modus adequate op).
Als κ groot is, dwingen deze twee vereisten een zeer kleine h: klein genoeg voor stabiliteit van de snelle modus, groot genoeg om de trage modus te bemonsteren. Het aantal stappen schalen met κ.
Geometrisch beeld in het eigenwaardespectrum: de eigenwaarden van de Jacobiaan ∂f/∂y vormen een reeks punten in het complexe vlak. De stabiliteitzone van een expliciete oplosver moet alle punten h·λᵢ bevatten. Als eigenwaarden van -1 tot -1000 reiken, moet de stabiliteitzone een bereik van 1000 langs de reële as bestrijken — wat zeer kleine h vereist.
Impliciete oplosvers: de stabiliteitzone van achterwaartse Euler bestrijkt het gehele linkerhalfvlak. Alle eigenwaarden met Re(λ) < 0 zijn automatisch binnen de stabiliteitzone ongeacht h. De beperking op h komt alleen uit nauwkeurigheid, niet uit stabiliteit.
Stijfheidsverhouding & Kosten
Beschouw een chemisch reactienetwerk met snelle reacties (tijdschaal 10⁻⁶ s) en trage reacties (tijdschaal 1 s).
Stijfheidsverhouding: κ = 10⁶ / 1 = 10⁶.
Met RK4 (stabiliteitslimiet h·|λ_max| ≤ 2.785): h_max = 2.785 / 10⁶ ≈ 2.8 × 10⁻⁶ s.
Om over 10 s reactietijd te integreren: stappen = 10 / (2.8 × 10⁻⁶) ≈ 3.6 × 10⁶.
Met achterwaartse Euler (onvoorwaardelijk stabiel): h kan worden gekozen voor nauwkeurigheid van de trage reacties. h = 10⁻² s (100 monsters over 1 s schaal). Stappen = 10 / 10⁻² = 1000.
Kostenverhouding: expliciete 3.6 miljoen stappen vs impliciete 1000 stappen — een factor van 3600. Elke impliciete stap vereist het oplossen van een lineair systeem (kosten per stap zijn hoger), maar de totale kosten zijn veel lager voor zeer stijve problemen.
Waarom n-Dimensionale Buizen Niet Zijn Wat Je Denkt
In 2D is een 'buis' met straal ε rond een curve C de reeks punten binnen afstand ε van C. De doorsnede is een cirkel met straal ε. Het volume van de buis groeit evenredig met zijn lengte.
In n dimensies verandert de buisgeometrie fundamenteel, vanwege het fenomeen uit Hoofdstuk 9:
De n-dimensionale hoekparadox: in n-dimensionale ruimte ligt vrijwel heel het volume van een n-dimensionale hyperkubus in de hoeken — niet in het centrale gebied. Naarmate n toeneemt, gaat het aandeel van het volume binnen afstand ε van het centrum naar nul voor elke vaste ε.
Toegepast op ODE-oplossingbuizen:
In 2D: als de werkelijke oplossing door het midden van een buis gaat, zijn de meeste nabijgelegen punten dicht bij de curve. Kleine verstoringen houden je dicht bij de werkelijke oplossing.
In hoge dimensies: de meeste punten in de begrenzingsdoos van de buis liggen eigenlijk ver van de werkelijke oplossingscurve. Het 'volume' van de buis wordt gedomineerd door de hoeken — regio's die in meerdere dimensies tegelijkertijd ver van het midden verwijderd zijn.
Gevolg voor simulatie: met 28 gekoppelde ODE's (Hammings Marine-onderscheppingsprobleem) kan een verstuwing van grootte ε in elke dimensie een totale verplaatsing van ε√28 ≈ 5.3ε van de werkelijke oplossing opleveren. De buis moet worden begrepen in termen van de L2-norm over alle dimensies, niet slechts de maximale verplaatsing in enige dimensie.
Stabiliteit in hoge dimensies: een systeem waarbij elke component onafhankelijk vervalt (elke eigenwaarde heeft negatief reëel deel) kan nog steeds grote gecombineerde verplaatsingen vertonen omdat de fouten van de componenten in L2-norm optellen. De 28-dimensionale buis is niet zomaar 28 onafhankelijke 1-dimensionale buizen — de geometrie koppelt ze.
Van Geometrie naar Ontwerp
De geometrische inzichten van Hoofdstukken 18-20 komen samen als een reeks ontwerpprincipes voor numerieke simulatie:
Stapgrootteselectie: h moet h·λ binnen de stabiliteitzone voor elke eigenwaarde plaatsen. Voor stijve systemen verwijderen impliciete methoden de stabiliteitsbeperking, waardoor alleen nauwkeurigheidsvereisten achterblijven.
Foutsamenstellng in hoge dimensies: de globale fout is een vector in n-dimensionale ruimte. De norm groeit als √n keer de per-componentfout. Hoog-dimensionale simulaties vereisen strenger per-stap nauwkeurigheidsvereisten.
Feedback als stabilizer: als de simulatie feedback verwerkt (de berekende output beïnvloedt volgende inputs, zoals in een geleideningssysteem), dempt convergente feedback fouten. De simulatie kan onnauwkeurige ingangen voor hoeveelheden binnen de feedbacklus verdragen.
De instabiliteit als signaal: voor problemen met divergente richtingsvelden kan instabiliteit worden benut: de divergentierichting bevat informatie over de beginvoorwaardefout, waardoor corrigerende aanpassingen mogelijk zijn.