Tre metoder, tre regioner
För testekvationen dy/dx = λy har tre explicita ODE-metoder följande stabilitetsregioner i det komplexa hλ-planet:
Eulers metod (första ordningen): stabilitetsregionen är skivan |1 + hλ| ≤ 1, en cirkel med radie 1 centrerad vid (-1, 0). Verkligt negativt hλ måste ligga i [-2, 0].
Runge-Kutta 2 (mittpunktsmetod) (andra ordningen): stabilitetsregionen är |1 + hλ + (hλ)²/2| ≤ 1. Större än Eulers skiva, men fortfarande begränsad.
Runge-Kutta 4 (fjärde ordningen): stabilitetsregionen uppfyller |1 + hλ + (hλ)²/2 + (hλ)³/6 + (hλ)⁴/24| ≤ 1. Verkligt negativt hλ sträcker sig till ungefär -2,785. Regionen är väsentligt större än Eulers.
Bakåt Euler (implicit): stabilitetsregionen är hela det komplexa planet förutom skivan |1 - hλ|⁻¹ > 1, ekvivalent |1/(1-hλ)| ≤ 1. För λ i vänster halvplan (Re(λ) < 0) är detta ovillkorligt stabilt; ingen begränsning på h från stabilitet.
Amplifieringsfunktionen
För alla Runge-Kutta-metoder är amplifieringsfaktorn per steg R(hλ) en polynomisk approximation till e^(hλ):
- Euler: R(z) = 1 + z (trunkerad vid grad 1)
- RK2: R(z) = 1 + z + z²/2 (trunkerad vid grad 2)
- RK4: R(z) = 1 + z + z²/2 + z³/6 + z⁴/24 (trunkerad vid grad 4)
Stabilitetsregionen är {z ∈ ℂ : |R(z)| ≤ 1}. Den sanna lösningens amplifiering: |e^z| = e^(Re(z)). För Re(z) < 0 (stabil ODE) sönderfaller den sanna lösningen. Den numeriska metoden är stabil om |R(z)| ≤ 1; detta motsvarar sönderfallsbeteendet.
Rent imaginära egenvärden: Oscillatoriska system
Många fysiska system har rent imaginära egenvärden: λ = iω (oscillationer utan dämpning). Fjäder-massa-systemet, orbitmekanik, pendelynamik.
För λ = iω ligger hλ = ihω på den imaginära axeln.
Eulers stabilitet på den imaginära axeln: |1 + ihω|² = 1 + (hω)² > 1 för alla h > 0. Euler är instabil för vilken stegstorlek som helst på rent imaginära egenvärden. Den beräknade 'oscillationen' växer obegränsat.
RK4:s stabilitet på den imaginära axeln: stabilitetsregionen sträcker sig till ungefär |hω| ≤ 2,83 på den imaginära axeln. För tillräckligt litet h hanterar RK4 odämpade oscillationer. Euler kan inte.
Denna geometri är anledningen till att Euler misslyckas med konservativa system (fjäder-massa, omloppsbanor, vågekvationer) även med litet h, medan RK4 hanterar dem väl.
Den styva problemens geometri
Ett styvt ODE-system har egenvärden med mycket olika magnituder. Styvhetsförhållandet: κ = max|Re(λᵢ)| / min|Re(λᵢ)| >> 1.
Varför styvhet är dyrt för explicita lösare:
Stabilitet kräver h·max|λᵢ| ≤ C (där C beror på metoden). Det mest negativa egenvärdet sätter gränsen.
Noggrannhet för den långsamma dynamiken kräver h·min|λᵢ| ≥ ε (lösa den långsammaste moden på lämpligt sätt).
Om κ är stor tvingar dessa två krav ett litet h: litet nog för stabilitet för den snabba moden, tillräckligt stort för att sampla den långsamma moden. Antalet steg skaleras med κ.
Geometrisk bild i egenvärdespektrumet: egenvärden för Jacobian ∂f/∂y bildar en uppsättning punkter i det komplexa planet. En explicit lösares stabilitetsregion måste innehålla alla punkter h·λᵢ. Om egenvärden sträcker sig från -1 till -1000 måste stabilitetsregionen täcka ett intervall på 1000 längs den reella axeln; detta kräver mycket litet h.
Implicita lösare: bakåt Eulers stabilitetsregion täcker helt vänster halvplan. Alla egenvärden med Re(λ) < 0 är automatiskt inuti stabilitetsregionen oavsett h. Begränsningen på h kommer endast från noggrannhet, inte från stabilitet.
Styvhetsförhållande & kostnad
Betrakta ett kemiskt reaktionsnätverk med snabba reaktioner (tidsskala 10⁻⁶ s) och långsamma reaktioner (tidsskala 1 s).
Styvhetsförhållande: κ = 10⁶ / 1 = 10⁶.
Med RK4 (stabilitetsgräns h·|λ_max| ≤ 2,785): h_max = 2,785 / 10⁶ ≈ 2,8 × 10⁻⁶ s.
För att integrera över 10 s reaktionstid: steg = 10 / (2,8 × 10⁻⁶) ≈ 3,6 × 10⁶.
Med bakåt Euler (ovillkorligt stabil): h kan väljas för noggrannhet för de långsamma reaktionerna. h = 10⁻² s (100 sampel över 1 s skala). Steg = 10 / 10⁻² = 1000.
Kostnadförhållande: explicit 3,6 miljoner steg mot implicit 1000 steg; en faktor 3600. Varje implicit steg kräver att lösa ett linjärt system (kostnad per steg är högre), men den totala kostnaden är mycket lägre för mycket styva problem.
Varför n-dimensionella rör inte är vad du tror
I två dimensioner är ett 'rör' med radie ε omkring en kurva C uppsättningen av punkter inom avståndet ε från C. Tvärsnittet är en cirkel med radie ε. Rörets volym växer proportionellt med dess längd.
I n dimensioner förändras rörets geometri fundamentalt på grund av fenomenet från kapitel 9:
N-dimensionella hörnens paradox: i n-dimensionellt rum ligger nästan all volym i en n-dimensionell hyperkub i hörnen; inte i den centrala regionen. Då n ökar går fraktionen av volymen inom avståndet ε från mitten mot noll för all fixerad ε.
Tillämpat på ODE-lösningsrör:
I två dimensioner: om den sanna lösningen passerar genom rörets centrum ligger de flesta närliggande punkter nära kurvan. Små störningar håller dig nära den sanna lösningen.
I höga dimensioner: de flesta punkter inom rörets begränsningsbox ligger faktiskt långt från den sanna lösningskurvan. Rörets 'volym' domineras av hörnen; regioner långt från mitten i flera dimensioner samtidigt.
Konsekvens för simulering: med 28 kopplade ODEer (Hammings marinintercept-problem) kan en störning av storlek ε i varje dimension producera ett totalt förflyttning på ε√28 ≈ 5,3ε från den sanna lösningen. Röret måste förstås i termer av L2-normen över alla dimensioner, inte bara den största förflyttningen i någon dimension.
Stabilitet i höga dimensioner: ett system där varje komponent sönderfaller oberoende (varje egenvärde har negativ realdel) kan fortfarande visa stora kombinerade förflyttningar eftersom komponenternas fel adderas i L2-normen. Röret i 28 dimensioner är inte bara 28 oberoende 1-dimensionella rör; geometrin kopplar dem.
Från geometri till design
De geometriska insikterna från kapitlen 18-20 kommer tillsammans som en uppsättning designprinciper för numerisk simulering:
Stegstorlek-val: h måste placera h·λ inuti stabilitetsregionen för alla egenvärden. För styva system tar implicita metoder bort stabilitetsbegränsningen; det lämnar bara noggrannhetskrav.
Felackumulering i höga dimensioner: det globala felet är en vektor i n-dimensionellt rum. Dess norm växer som √n gånger per-komponent-felet. Höga dimensionella simuleringar behöver strängare per-steg-noggrannhetskrav.
Återkoppling som stabilisator: om simuleringen inkorporerar återkoppling (den beräknade utmatningen påverkar efterföljande ingångar, som i ett ledningssystem) dämpar konvergent återkoppling fel. Simuleringen kan tolerera oprecisa ingångar för kvantiteter inuti återkopplingsslingan.
Instabilitet som signal: för problem med divergerande riktningsfält kan instabilitet exploateras: divergensriktningen bär information om feltillståndsfelgränsen; detta möjliggör korrigerande justering.