Differensialligninger

Hva slags ligninger snakker vi om?

Et initialverdiproblem har formen y' = f(t, y) med en startverdi y(t0) = y0. Vi vil finne y(t) videre. To strategier dominerer: analytisk når strukturen tillater det (lineære systemer, 2. ordens med konstante koeffisienter), og numerisk når den ikke gjør det.

Dette kapittelet binder sammen de to. Først bruker vi lineær algebra til å løse lineært system eksakt. Så bygger vi numeriske skjemaer som approksimerer y(t) i et rutenett av tidspunkter t0 < t1 < t2 < \ldots med {{term:h}} h = t{n+1} - t_n.

Diagonalisering — basis-skifte som forenkler alt

En matrise A er diagonaliserbar hvis A = PDP^{-1} der D er diagonal, eller ekvivalent: A har n lineært uavhengige egenvektorer. ƒdiagonalisering flytter regningen til den enkleste basen som finnes — den der A bare strekker hver akse med en faktor.

egenverdi / egenvektor λ og egenvektoren v oppfyller Av = λv, det vil si (A − λI)v = 0. For at det skal finnes en ikke-null v må determinanten være null, så egenverdiene er røttene av ƒkarakteristisk polynom. For 2x2 forenkles dette via spor og determinant til ƒspor og determinant (2x2) — du leser av spor og determinant og løser en andregradsligning.

Lineære systemer y' = A y

Når A er diagonaliserbar er løsningen rett fram: skift basis med w = P^{-1} y, så blir systemet w' = D w med n frakoblede ligninger wi' = λi w_i. Hver av dem løses med en eksponential. Tilbake i originale koordinater får vi ƒløsning av lineært system y' = a y.

Hvis du foretrekker matrise-formulering: ƒdiagonal eksponentiering sier at matrise-eksponentialen er like enkel som diagonalen. Dette er hele grunnen til at diagonalisering er så nyttig for dynamiske systemer — den fanger den lange-tids-oppførselen i n separate skalar-problemer.

Lineær 2. ordens

Ligningen ay'' + by' + cy = 0 har ƒkarakteristisk ligning (2. ordens) som karakteristisk ligning. Tre tilfeller:

- To distinkte reelle røtter r1, r2: y = C1 e^{r1 t} + C2 e^{r2 t}. - Dobbeltrot r: y = (C1 + C2 t) e^{rt}. Den ekstra t-faktoren kompenserer for at vi mangler en andre uavhengig løsning. - Komplekse røtter α ± iβ: y = e^{αt}(C1 cos βt + C2 sin βt). Realdelen styrer demping, imaginærdelen frekvens.

Hver kombinasjon dukker opp i fysikken: overdempet, kritisk dempet og underdempet svingning.

Et 2x2-eksempel: kobling og frakobling

La A = [[3, 1], [1, 3]]. Spor er 6, determinant 8, så det karakteristiske polynomet er λ² − 6λ + 8 = (λ − 2)(λ − 4). Egenverdiene er λ1 = 2 og λ2 = 4. Egenvektorer: for λ = 2 finner vi v1 = (1, −1)^T fra (A − 2I)v = 0; for λ = 4 finner vi v2 = (1, 1)^T. Da er P kolonnene (v1, v2) og løsningen av y' = Ay er y(t) = c1 e^{2t} (1, −1)^T + c2 e^{4t} (1, 1)^T. To frakoblede modi som vokser med ulik hastighet.

Hvis A i stedet har komplekse egenverdier α ± iβ får vi løsninger som spiraler — kombinasjoner av e^{αt} cos(βt) og e^{αt} sin(βt). Da er løsningen reell selv om egenverdiene er komplekse, fordi de komplekse leddene parres opp.

Numerikk: trinn for trinn

For ligninger uten lukket form bygger vi en sekvens u0, u1, u2, \ldots som approksimerer y(t0), y(t1), y(t2), \ldots. Den enkleste oppskriften er eulers metode: ƒeulers metode (eksplisitt). Geometrisk: vi sitter på (tn, un), regner stigningen f(tn, un), og går rett fram i den retningen et tidssteg.

Eksplisitt Euler har global feil O(h) — første orden. Halverer du h, halveres feilen. Det høres ok ut til du ser at RK4 gir O(h^4): samme halvering reduserer feilen med faktor 16. Forskjellen blir astronomisk over lange intervaller.

For å bli konkret: løs y' = y, y(0) = 1 (eksakt løsning e^t) på [0, 1] med Euler og RK4.

Med h = 0.1 gir Euler u_{10} ≈ 2.594 (eksakt 2.718, feil 0.124). RK4 med samme h gir 2.7183 — på fjerde desimal er løsningen rett. Det er O(h^4) i praksis.

Stabilitetsregion — hva tåler metoden?

Et nyttig analysert verktøy er stabilitetsregionen: settet av komplekse hλ-verdier der metoden anvendt på y' = λy ikke vokser ut av kontroll. Eksplisitt Euler har en liten sirkel rundt −1 i det komplekse planet. Tar du et h slik at hλ havner utenfor, eksploderer u_n. Implisitt Euler har derimot hele venstre halvplan som stabilitetsregion — den er A-stabil. Crank–Nicolson også. Disse to klassene er grunnen til at implisitte metoder dominerer for stive problemer: stabiliteten er ikke begrenset av hvor stort h du våger.

Implisitte metoder og stive problemer

Snur du tegningen og bruker stigningen i enden av steget får du implisitt euler: ƒimplisitt (bakover) euler. Nå står u_{n+1} på begge sider og må løses som ligning — for ikke-lineær f er det en Newton-iterasjon hvert steg. Hvorfor gidde? Fordi metoden er stabil for stive problemer: ligninger der løsningen har komponenter med svært ulik tidsskala.

Et klassisk eksempel: y' = -1000y, y(0) = 1. Eksplisitt Euler eksploderer med mindre h < 2/1000. Implisitt Euler holder seg pent uansett hvor stort h du velger. Prisen er én ligningløsning per steg; gevinsten er at h kan være tusen ganger større.

trapesmetoden / crank–nicolson er gjennomsnittet av eksplisitt og implisitt Euler: ƒtrapesmetoden / crank–nicolson. Den er implisitt, andre orden, og symmetrisk i tid — ofte førstevalget for diffusjons-typer.

Runge–Kutta og Butcher-tablå

runge–kutta-metode-familien tar ett steg med s mellom-evalueringer av f og kombinerer dem: ƒgenerell runge–kutta-metode. butcher-tablå-tablået oppsummerer hele metoden — knutene c, koeffisient-matrisen A og vektene b. En metode er eksplisitt metode hvis A er strengt nedre triangulær (hver Ki avhenger bare av tidligere Kj); ellers er den implisitt metode.

For at metoden i det hele tatt skal være konsistent — at lokal feil skal gå raskere mot null enn h — trengs ƒkonsistens av runge–kutta. Det er bare summen av vektene. Høyere orden krever flere algebraiske betingelser på c og a_{ij}.

Den mest brukte eksplisitte metoden er rk4: ƒklassisk rk4 med fire f-evalueringer som veies sammen. Global feil O(h^4). For ingeniørberegninger er det nesten alltid bra nok. midtpunktmetoden er den enklere kusinen — to evalueringer, orden 2: ƒmidtpunktmetoden (rk2).

Konvergens, konsistens og feilorden

To feilbegreper må holdes adskilt. lokal avbruddsfeil er bommen i ett enkelt steg gitt at starten av steget er eksakt — for en metode av orden p er den O(h^{p+1}). global avbruddsfeil er den akkumulerte feilen ved tid t etter mange steg. Sammenhengen er ƒsammenheng lokal/global feil: én h-faktor "spises" av at antall steg er proporsjonalt med 1/h.

konvergens betyr at un → y(tn) når h → 0. Det er ikke gratis — konsistens (lokal feil går mot null) må kombineres med nullstabilitet (gamle feil må ikke vokse ukontrollert) for å sikre det. For ettstegs-metoder følger nullstabilitet automatisk, så da er konsistens nok. For flerstegs-metoder må det sjekkes separat.

Hvilken metode velger man?

I praksis er det tre spørsmål: hvor stiv er ligningen, hvor mye nøyaktighet trenger jeg, og hvor dyr er én f-evaluering? For ikke-stive problemer med moderat krav til nøyaktighet er RK4 førstevalg — fire f-evalueringer per steg gir orden 4 og er enkel å implementere. For stive problemer går du til implisitt Euler eller Crank–Nicolson; ekstra ligningløsning per steg er en lav pris for å unngå tusenfold mindre h. For svært høyt presisjonskrav over lange intervaller bruker du adaptive RK-varianter (Dormand–Prince, Cash–Karp) som justerer h underveis.

Det er også verdt å huske at vi snakker om approksimasjoner, ikke svar. Lange simuleringer akkumulerer avrundingsfeil i tillegg til avbruddsfeil. Vil du ha y(100) regnet ut over en million steg, blir bit-røra etter hvert dominerende. Da må enten metode (mer presisjon per steg) eller aritmetikk (dobbel/kvad-presisjon) oppgraderes.