Tilbake
7.4

7.4 Stabilitetsfunksjonen $R(z)$ og stabilitetsintervall

Bruk metoden på testlikningen $y'=\lambda y$, finn $R(z)$ med $z=\lambda h$, og bestem det stabile steglengdeintervallet — også for systemer via egenverdiene.

50 min
11 oppgaver
Stabilitetsfunksjonen $R(z)$stabilitetsintervall
Din fremgang i kapitlet
0 / 11 oppgaver
Forkunnskaper: kap. 7.1 — metodeformlene og eksemplet i oppgave 7 der eksplisitt Euler brøt sammen på y=10yy'=-10y. Dette kapitlet forklarer nøyaktig hva som skjedde der.

Du trenger dessuten å finne egenverdiene til en 2×22\times2-matrise: løs det(AλI)=0\det\left(A-\lambda I\right)=0, altså andregradslikningen

λ2(trA)λ+detA=0.\lambda^{2}-\left(\operatorname{tr}A\right)\lambda+\det A=0.

Kapitlet er ikke forutsetning for noe senere kapittel, men det samme resonnementet dukker opp igjen i Del 8, der stabilitetskravet for det eksplisitte varmeskjemaet blir r12r\le\tfrac12.

Når metoden lyver

I kap. 7.1 løste vi y=10yy'=-10y med eksplisitt Euler og h=0,25h=0{,}25, og fikk

y1=1,5.y_1=-1{,}5.

Den eksakte løsningen e10te^{-10t} er positiv og går raskt mot null. Metoden ga et negativt tall med større tallverdi enn startverdien — og fortsetter du, vokser yn\left|y_n\right| uten grense.

Det er ikke unøyaktighet. Det er kollaps. En unøyaktig metode bommer litt; en ustabil metode gir svar som ikke har noe med løsningen å gjøre.

Og legg merke til hva som ikke hjalp: en bedre metode. RK4 med samme hh på samme problem bryter også sammen, bare litt senere. Problemet er ikke ordenen — det er steglengden.

Spørsmålet dette kapitlet svarer på er: hvor stor kan hh være før metoden begynner å lyve?

Grepet er å teste metoden på det enkleste problemet som finnes. Likningen

y=λyy'=\lambda y

har løsningen y=y0eλty=y_0e^{\lambda t}, som for λ<0\lambda<0 går mot null. Bruker vi en ett-skritts-metode på den, får vi alltid

yn+1=R(z)yn,z=λh,y_{n+1}=R(z)\,y_n,\qquad z=\lambda h,

altså en ren multiplikasjon med et tall R(z)R(z). Da er kravet åpenbart: skal yny_n gå mot null, må R(z)1\left|R(z)\right|\le 1.

Det er hele kapitlet. Resten er å regne ut RR for hver metode og se hvilke zz som slipper gjennom.

Flashcard- og repetisjonsstoff — hopp trygt over ved førstegangslesing. Definisjonsboksene er samtidig kortene i flashcard-bunken. Tidsanslaget på 50 minutter gjelder kjernestoffet.

Løkke 1 — Testlikningen og R(z)R(z) (~13 min)

Testlikningen
Den enkleste differensiallikningen som fanger det vesentlige:

y=λy,y(0)=y0,y'=\lambda y,\qquad y(0)=y_0,

med eksakt løsning y(t)=y0eλty(t)=y_0e^{\lambda t}.

Hvorfor akkurat denne. For λ<0\lambda<0 dør løsningen ut, og det er nettopp den oppførselen en metode kan svikte på. For λ>0\lambda>0 vokser løsningen, og da vokser en numerisk løsning også — det er ikke interessant.

Hvorfor det er nok å se på én likning. Et lineært system y=Ay\mathbf y'=A\mathbf y kan diagonaliseres, og da faller det fra hverandre i uavhengige testlikninger med λ\lambda lik hver egenverdi. Det er hele koblingen til systemtilfellet i løkke 3.

λ\lambda kan være kompleks. Da beskriver testlikningen en svingning som dempes, og stabilitetsområdet blir et område i planet i stedet for et intervall.

Variabelen z=λhz=\lambda h

Produktet av likningens λ\lambda og metodens steglengde.

Hvorfor akkurat dette produktet dukker opp: alle metodene ganger hh med f(t,y)=λyf(t,y)=\lambda y, så hh og λ\lambda opptrer alltid sammen. Ingen metode kan skille dem.

Konsekvensen er praktisk viktig: stabilitet er ikke en egenskap ved metoden alene, og ikke ved problemet alene. Det er hλh\lambda som avgjør. En metode som er ustabil for λ=10\lambda=-10 og h=0,25h=0{,}25, er helt fin for λ=1\lambda=-1 og samme hh.

Merk at zz er dimensjonsløs. Størrelsen 1/λ1/|\lambda| er en tidsskala i problemet, og zz måler hvor mange slike tidsskalaer ett skritt spenner over.

Stabilitetsfunksjonen R(z)R(z)
Faktoren løsningen ganges med per skritt når metoden brukes på testlikningen:

yn+1=R(z)yn,z=λh.y_{n+1}=R(z)\,y_n,\qquad z=\lambda h.

Slik finner du den: sett f(t,y)=λyf(t,y)=\lambda y inn i metodeformelen og faktoriser ut yny_n. Det som står igjen, er R(z)R(z).

Etter nn skritt er

yn=R(z)ny0,y_n=R(z)^{\,n}\,y_0,

så oppførselen over tid avgjøres helt av om R(z)\left|R(z)\right| er større eller mindre enn 1.

Sammenlikn med det eksakte: den eksakte løsningen ganges med eλh=eze^{\lambda h}=e^{z} per skritt. R(z)R(z) er metodens tilnærming til eze^{z}, og for en metode av orden pp stemmer de til og med leddet zpz^{p} i rekkeutviklingen.

Det finnes ingen ferdig R(z)R(z) på det utdelte formelarket — utledningen må kunnes.

Absolutt stabilitet
Kravet

R(z)1.\left|R(z)\right|\le 1.

Innholdet: løsningen vokser ikke fra skritt til skritt. Er kravet oppfylt, er yn=R(z)ny0\left|y_n\right|=\left|R(z)\right|^{n}\left|y_0\right| ikke-voksende, og for R(z)<1\left|R(z)\right|<1 går den mot null — slik den eksakte løsningen gjør når λ<0\lambda<0.

Er kravet brutt, vokser den numeriske løsningen eksponentielt selv om den eksakte dør ut. Det er nøyaktig det som skjedde med y=10yy'=-10y og h=0,25h=0{,}25: der er z=2,5z=-2{,}5 og R(z)=1+z=1,5R(z)=1+z=-1{,}5, med R=1,5>1\left|R\right|=1{,}5>1.

Merk ordet «absolutt». Det skiller dette fra andre stabilitetsbegreper og betyr at vi ser på selve tallverdien til vekstfaktoren, uavhengig av hvor nøyaktig metoden er.

Kravet må kunnes.

Stabilitetsintervallet på den reelle aksen
Mengden av reelle, negative zz som oppfyller R(z)1\left|R(z)\right|\le 1. Den er alltid på formen

zz0.z^{*}\le z\le 0.

MetodeR(z)R(z)Intervall
Eksplisitt Euler1+z1+z[2, 0][-2,\ 0]
Heun (forbedret Euler)1+z+z221+z+\tfrac{z^{2}}{2}[2, 0][-2,\ 0]
Kutta 31+z+z22+z361+z+\tfrac{z^{2}}{2}+\tfrac{z^{3}}{6}[2,51, 0][-2{,}51,\ 0]
RK41+z+z22+z36+z4241+z+\tfrac{z^{2}}{2}+\tfrac{z^{3}}{6}+\tfrac{z^{4}}{24}[2,785, 0][-2{,}785,\ 0]
Bakover-Euler11z\dfrac{1}{1-z}hele z0z\le 0

Legg merke til at Euler og Heun har samme intervall, til tross for ulik orden. Høyere orden gir ikke automatisk større stabilitetsområde — det er to uavhengige egenskaper.
Og legg merke til hvor lite RK4 vinner: fra 22 til 2,7852{,}785, altså 39 prosent. Alle eksplisitte metoder har et endelig stabilitetsintervall, og det er den grunnleggende begrensningen som gjør implisitte metoder nødvendige for stive problemer.
✏️Utled $R(z)$ for Euler, Heun og bakover-Euler

Bruk hver av de tre metodene på testlikningen y=λyy'=\lambda y og finn stabilitetsfunksjonen.

Eksplisitt Euler. Formelen er yn+1=yn+hf(tn,yn)y_{n+1}=y_n+h\,f\left(t_n,y_n\right) med f(t,y)=λyf(t,y)=\lambda y:

yn+1=yn+hλyn=(1+hλ)yn=(1+z)yn.y_{n+1}=y_n+h\lambda y_n=\left(1+h\lambda\right)y_n=(1+z)\,y_n.

R(z)=1+z\boxed{R(z)=1+z}

Heun. Stigningstallene blir

k1=λyn,k2=λ(yn+hk1)=λ(yn+hλyn)=λyn(1+z).k_1=\lambda y_n,\qquad k_2=\lambda\left(y_n+hk_1\right)=\lambda\left(y_n+h\lambda y_n\right)=\lambda y_n\left(1+z\right).

yn+1=yn+h2(k1+k2)=yn+hλ2(yn+yn(1+z))=yn+z2yn(2+z).y_{n+1}=y_n+\frac{h}{2}\left(k_1+k_2\right)=y_n+\frac{h\lambda}{2}\left(y_n+y_n(1+z)\right)=y_n+\frac{z}{2}\,y_n\left(2+z\right).

yn+1=yn(1+z+z22).y_{n+1}=y_n\left(1+z+\frac{z^{2}}{2}\right).

R(z)=1+z+z22\boxed{R(z)=1+z+\frac{z^{2}}{2}}

Bakover-Euler. Formelen er implisitt:

yn+1=yn+hλyn+1.y_{n+1}=y_n+h\lambda y_{n+1}.

Samle yn+1y_{n+1} på venstre side:

yn+1(1hλ)=yn  yn+1=yn1z.y_{n+1}\left(1-h\lambda\right)=y_n\ \Longrightarrow\ y_{n+1}=\frac{y_n}{1-z}.

R(z)=11z\boxed{R(z)=\frac{1}{1-z}}

Sammenlikn alle tre med eze^{z}. Rekkeutviklingen er

ez=1+z+z22+z36+e^{z}=1+z+\frac{z^{2}}{2}+\frac{z^{3}}{6}+\dots

- Euler treffer til og med z1z^{1}orden 1
- Heun treffer til og med z2z^{2}orden 2
- Bakover-Euler: 11z=1+z+z2+z3+\dfrac{1}{1-z}=1+z+z^{2}+z^{3}+\dots, som treffer til og med z1z^{1}orden 1

Dette er en gratis ordenskontroll. Utvikler du R(z)R(z) i rekke og sammenlikner med eze^{z}, får du metodens orden uten å røre ordensbetingelsene fra kap. 7.2.

(Forbehold: kontrollen virker fordi testlikningen er så enkel. En metode kan treffe eze^{z} til høy orden og likevel ha lavere orden for generelle ff — akkurat som i oppgave 8 i kap. 7.2. Bruk den som indikasjon, ikke som bevis.)

Merk hvor mye pen struktur det er i disse tre svarene. For eksplisitte metoder blir RR et polynom — den avkuttede rekka for eze^{z}. For implisitte blir den en rasjonal funksjon, og det er nettopp brøkformen som gir den ubegrensede stabiliteten.

📝Oppgave 1

(Innstegsoppgave — ren innsetting.) Eksplisitt Euler brukes på y=4yy'=-4y.

a) Hva er R(z)R(z)?
b) Regn ut RR for h=0,1h=0{,}1, h=0,5h=0{,}5 og h=0,6h=0{,}6.
c) For hvilke av disse er metoden stabil?

📝Oppgave 2
O

Vis at stabilitetsintervallet for eksplisitt Euler på den reelle aksen er 2z0-2\le z\le 0, og finn maks steglengde for y=8yy'=-8y.

Løkke 2 — Stabilitetsintervallene og figuren (~13 min)

R(z)R(z) for eksplisitte Runge–Kutta-metoder
For en eksplisitt metode av orden pp med s=ps=p steg er

R(z)=1+z+z22!++zpp!,R(z)=1+z+\frac{z^{2}}{2!}+\dots+\frac{z^{p}}{p!},

altså Taylor-polynomet til eze^{z} av grad pp.

Hvorfor. Metoden er eksakt til orden pp, så R(z)R(z) må stemme med eze^{z} til og med leddet zpz^{p}. Og siden hvert stigningstall bare kan bidra med én potens av zz, kan RR ikke ha høyere grad enn ss.

Konsekvensen: RR er alltid et polynom for en eksplisitt metode, og et polynom går mot uendelig når z\left|z\right|\to\infty. Derfor har enhver eksplisitt metode et endelig stabilitetsintervall — det finnes alltid en hh som er for stor.

Det er den grunnleggende asymmetrien mellom eksplisitte og implisitte metoder, og hele grunnen til at stive problemer krever implisitte løsere.

A-stabilitet
En metode er A-stabil når stabilitetsområdet inneholder hele venstre halvplan:

Re(z)0  R(z)1.\operatorname{Re}(z)\le 0\ \Longrightarrow\ \left|R(z)\right|\le 1.

Innholdet: metoden er stabil for enhver steglengde, uansett hvor negativ λ\lambda er.

Bakover-Euler er A-stabil. Med R(z)=11zR(z)=\dfrac{1}{1-z} og Re(z)0\operatorname{Re}(z)\le 0 er 1z1\left|1-z\right|\ge 1, så

R(z)=11z1.\left|R(z)\right|=\frac{1}{\left|1-z\right|}\le 1.\qquad ✔

(For reell z0z\le 0 ser du det direkte: 1z11-z\ge 1, så brøken er høyst 1.)

Ingen eksplisitt metode kan være A-stabil, siden RR da er et polynom og vokser uten grense.

Prisen for A-stabilitet er at hvert skritt krever en likningsløsning. Gevinsten er at du kan velge hh etter nøyaktighet i stedet for etter stabilitet — og for et stivt problem er det forskjellen mellom tusen og en million skritt.

Maks steglengde hz/λh\le\left|z^{*}\right|/\left|\lambda\right|
Oversettelsen fra stabilitetsintervall til steglengde. Er intervallet [z,0]\left[z^{*},0\right] og λ<0\lambda<0 reell, gir z=λhzz=\lambda h\ge z^{*} at

hzλ.h\le\frac{\left|z^{*}\right|}{\left|\lambda\right|}.

Metodez\left|z^{*}\right|Maks hh for y=5yy'=-5y
Euler20,40{,}4
Heun20,40{,}4
RK42,7852{,}7850,5570{,}557
Bakover-Euler\inftyingen grense

Merk at RK4 bare vinner 39 prosent på steglengden, til tross for fire ganger så mange funksjonsevalueringer per skritt. For et stivt problem er RK4 altså dårligere enn Euler per evaluering — det er en av de mest kontraintuitive og viktigste innsiktene i numerisk ODE-løsning.
Formelen må kunnes.
✏️Stabilitetsintervallet for Heun — og hvorfor det er $[-2,0]$

Vis at Heuns metode har stabilitetsintervallet 2z0-2\le z\le 0 på den reelle aksen, og forklar hvorfor det er det samme som for Euler.

Steg 1 — funksjonen. Fra eksempel 1: R(z)=1+z+z22R(z)=1+z+\dfrac{z^{2}}{2}.

Steg 2 — den øvre grensen R1R\le 1.

1+z+z221  z+z220  z(1+z2)0.1+z+\frac{z^{2}}{2}\le 1\ \Longleftrightarrow\ z+\frac{z^{2}}{2}\le 0\ \Longleftrightarrow\ z\left(1+\frac z2\right)\le 0.

Produktet av to faktorer er negativt eller null mellom nullpunktene z=0z=0 og z=2z=-2:

2z0.-2\le z\le 0.

Steg 3 — den nedre grensen R1R\ge -1.

1+z+z221  2+z+z220  z2+2z+40.1+z+\frac{z^{2}}{2}\ge -1\ \Longleftrightarrow\ 2+z+\frac{z^{2}}{2}\ge 0\ \Longleftrightarrow\ z^{2}+2z+4\ge 0.

Diskriminanten er 416=12<04-16=-12<0, så andregradsuttrykket har ingen reelle nullpunkter og er alltid positivt (koeffisienten foran z2z^{2} er positiv).

Denne betingelsen er altså oppfylt for alle zz og gir ingen ny begrensning.

Steg 4 — konklusjon.

2z0\boxed{-2\le z\le 0}

Kontroll i endepunktet: R(2)=12+2=1R(-2)=1-2+2=1, altså R=1\left|R\right|=1 ✔ — nøyaktig på grensen. Og R(2,1)=12,1+2,205=1,105>1R(-2{,}1)=1-2{,}1+2{,}205=1{,}105>1

Hvorfor det er samme intervall som for Euler. Svaret ligger i steg 3: for Euler er R=1+zR=1+z, og der er det den nedre grensen 1+z11+z\ge -1 som gir z2z\ge -2. For Heun er den nedre grensen aldri aktiv, og det er i stedet den øvre grensen R1R\le 1 som gir z2z\ge -2.

Samme tall, helt ulik årsak. Det er ren tilfeldighet at grensene faller sammen — og det illustrerer at orden og stabilitetsintervall er uavhengige egenskaper.

Praktisk konsekvens: Heun koster dobbelt så mange funksjonsevalueringer som Euler og gir ingen som helst gevinst i tillatt steglengde. For et stivt problem er Heun altså rent tap, selv om den er dobbelt så nøyaktig for ikke-stive problemer.

Merk hvordan tallverdien håndteres. Ulikheten R1\left|R\right|\le 1 må alltid deles i to: R1R\le 1 og R1R\ge -1. Å bare sjekke den ene er den vanligste feilen i utledningen — og for Heun ville du da fått riktig svar av gal grunn.

📝Oppgave 3
O
Utled R(z)R(z) for midtpunktsmetoden

k1=f(tn,yn),k2=f(tn+h2, yn+h2k1),yn+1=yn+hk2,k_1=f\left(t_n,y_n\right),\quad k_2=f\left(t_n+\tfrac h2,\ y_n+\tfrac h2 k_1\right),\quad y_{n+1}=y_n+h\,k_2,

og finn stabilitetsintervallet på den reelle aksen.

📝Oppgave 4
O

Problemet y=25yy'=-25y skal løses fram til t=1t=1.

a) Finn maks hh for eksplisitt Euler, og hvor mange skritt det gir.
b) Samme for RK4, der z=2,785\left|z^{*}\right|=2{,}785.
c) Sammenlikn antall ff-evalueringer, og kommenter.

— naturlig pausepunkt (~26 min brukt) —

Du kan utlede R(z)R(z) og finne intervallet for én likning. De to siste løkkene handler om systemer, som er den formen sjangeren faktisk kommer i på settene.

Løkke 3 — Systemer og egenverdier (~13 min)

En enkelt skalarlikning er sjelden interessant i praksis. Det interessante er systemer — og der er stabilitetsanalysen nesten like enkel, hvis du vet hvor du skal se.

Stabilitet for lineære systemer
For y=Ay\mathbf y'=A\mathbf y med en diagonaliserbar matrise AA:

Skriv om i egenvektorbasis. Da faller systemet fra hverandre i uavhengige skalarlikninger

wj=λjwj,w_j'=\lambda_j w_j,

én for hver egenverdi λj\lambda_j til AA.

Metoden virker uavhengig på hver komponent, siden den er lineær. Derfor er kravet

R(λjh)1for ALLE j.\left|R\left(\lambda_j h\right)\right|\le 1\qquad\text{for ALLE } j.

Alle egenverdiene må inn i stabilitetsområdet samtidig. Det er én steglengde hh, og den må passe for hver av dem.

Framgangsmåten på eksamen, i tre steg:

1. Finn egenverdiene til AA (løs λ2(trA)λ+detA=0\lambda^{2}-\left(\operatorname{tr}A\right)\lambda+\det A=0 for en 2×22\times2-matrise).
2. Finn kravet på hh for hver egenverdi.
3. Ta det strengeste kravet.

Den mest negative egenverdien styrer
For reelle negative egenverdier og et intervall [z,0]\left[z^{*},0\right] er kravet

hzλjfor hver j,h\le\frac{\left|z^{*}\right|}{\left|\lambda_j\right|}\qquad\text{for hver } j,

og det strengeste kravet kommer fra den største λj\left|\lambda_j\right| — altså den mest negative egenverdien:

hzmaxjλj.h\le\frac{\left|z^{*}\right|}{\max_j\left|\lambda_j\right|}.

⚠ Dette er fella i sjangeren. Det er lett å gripe den egenverdien som står først, eller den som ser «viktigst» ut. Det er alltid den med størst tallverdi — den raskeste komponenten — som binder.

Intuisjonen: den raskt dempede komponenten er kanskje helt uinteressant for svaret, men den er der, og en eksplisitt metode må ta små nok skritt til å håndtere den. Det er nettopp definisjonen på et stivt problem.

Merk at det er tallverdien som teller, ikke fortegnet eller rekkefølgen. Med λ=1\lambda=-1 og λ=6\lambda=-6 er det 6-6 som styrer.

Stivhetsforholdet
Forholdet mellom den raskeste og den langsomste komponenten:

S=maxjλjminjλj.S=\frac{\max_j\left|\lambda_j\right|}{\min_j\left|\lambda_j\right|}.

Er SS stor (typisk over 1000), er problemet stivt: steglengden bestemmes av en komponent som forsvinner nesten umiddelbart, mens du må integrere lenge for å følge den langsomme.

Eksempel: med λ1=1\lambda_1=-1 og λ2=1000\lambda_2=-1000 er S=1000S=1000. Den raske komponenten er borte etter t0,005t\approx 0{,}005, men eksplisitt Euler må likevel bruke h0,002h\le 0{,}002 hele veien — også når det bare er den langsomme komponenten igjen.

Løsningen er en implisitt metode. Bakover-Euler er A-stabil og kan bruke hh etter nøyaktighet, ikke etter stabilitet.

Merk at et lite stivhetsforhold ikke betyr at problemet er lett — det betyr bare at stabilitet ikke er flaskehalsen.

✏️Eksamensnivå: stabilitet for et system
Systemet y=Ay\mathbf y'=A\mathbf y har

A=(5222).A=\begin{pmatrix}-5 & 2\\ 2 & -2\end{pmatrix}.

a) Finn egenverdiene til AA.
b) Finn maks steglengde for eksplisitt Euler.
c) Finn maks steglengde for RK4, med z=2,785\left|z^{*}\right|=2{,}785.
d) Hva ville bakover-Euler krevd?

a) Egenverdiene. For en 2×22\times2-matrise er karakteristisk likning

λ2(trA)λ+detA=0.\lambda^{2}-\left(\operatorname{tr}A\right)\lambda+\det A=0.

trA=5+(2)=7,detA=(5)(2)22=104=6.\operatorname{tr}A=-5+(-2)=-7,\qquad \det A=(-5)(-2)-2\cdot 2=10-4=6.

λ2+7λ+6=0.\lambda^{2}+7\lambda+6=0.

λ=7±49242=7±52.\lambda=\frac{-7\pm\sqrt{49-24}}{2}=\frac{-7\pm 5}{2}.

λ1=1,λ2=6\boxed{\lambda_1=-1,\qquad \lambda_2=-6}

Kontroll: summen 1+(6)=7=trA-1+(-6)=-7=\operatorname{tr}A ✔ og produktet (1)(6)=6=detA(-1)(-6)=6=\det A

b) Eksplisitt Euler. Kravet er 1+λjh1\left|1+\lambda_j h\right|\le 1 for begge egenverdiene, altså 2λjh0-2\le\lambda_j h\le 0.

λ1=1:h21=2,\lambda_1=-1:\quad h\le\frac{2}{1}=2,
λ2=6:h26=130,3333.\lambda_2=-6:\quad h\le\frac{2}{6}=\frac13\approx 0{,}3333.

Det strengeste kravet gjelder:

h13\boxed{h\le\frac13}

Det er λ2=6\lambda_2=-6 — den mest negative — som binder, akkurat som regelen sier.

Kontroll ved h=13h=\tfrac13: z1=13z_1=-\tfrac13 gir R=23R=\tfrac23 ✔, og z2=2z_2=-2 gir R=1R=-1, altså R=1\left|R\right|=1 ✔ — akkurat på grensen.

Kontroll ved h=0,4h=0{,}4: z2=2,4z_2=-2{,}4 gir R=1,4R=-1{,}4, altså R=1,4>1\left|R\right|=1{,}4>1 ✘ — ustabil, selv om λ1\lambda_1 er helt fornøyd.

c) RK4.

λ1=1:h2,785,\lambda_1=-1:\quad h\le 2{,}785,
λ2=6:h2,78560,4642.\lambda_2=-6:\quad h\le\frac{2{,}785}{6}\approx 0{,}4642.

h0,464\boxed{h\le 0{,}464}

Igjen er det λ2\lambda_2 som binder, og gevinsten mot Euler er de samme 39 prosentene.

d) Bakover-Euler. Metoden er A-stabil: for enhver λj\lambda_j med Reλj0\operatorname{Re}\lambda_j\le 0 og enhver h>0h>0 er

R(λjh)=11λjh1,\left|R\left(\lambda_j h\right)\right|=\frac{1}{\left|1-\lambda_j h\right|}\le 1,

siden 1λjh11-\lambda_j h\ge 1 når λj0\lambda_j\le 0 og h>0h>0.

Ingen stabilitetsgrense på hh. Steglengden kan velges etter nøyaktighetskravet alene.

Er problemet stivt? Stivhetsforholdet er

S=61=6,S=\frac{6}{1}=6,

som er lite. Dette er altså ikke et stivt problem, og eksplisitte metoder er helt greie her — kravet h13h\le\tfrac13 er ikke plagsomt.

Hadde egenverdiene vært 1-1 og 6000-6000, ville S=6000S=6000, og eksplisitt Euler måtte brukt h3,3104h\le 3{,}3\cdot 10^{-4} mens den interessante komponenten lever i flere sekunder. Da hadde bakover-Euler vært eneste fornuftige valg.

Tidsbruk på eksamen: a) 4 min, b) 4 min, c) 2 min, d) 3 min — omtrent 13 minutter.

📝Oppgave 5
O
Systemet y=Ay\mathbf y'=A\mathbf y har

A=(3214).A=\begin{pmatrix}-3 & 2\\ 1 & -4\end{pmatrix}.

a) Finn egenverdiene.
b) Finn maks steglengde for eksplisitt Euler.
c) Regn ut stivhetsforholdet og vurder om problemet er stivt.

Løkke 4 — Stabilitet mot nøyaktighet (~11 min)

Til slutt en distinksjon som er lett å blande sammen, og som skiller besvarelsene.

Stabilitet er ikke nøyaktighet

To helt ulike krav, som begge begrenser hh:

StabilitetNøyaktighet
KreverR(λh)1\left|R(\lambda h)\right|\le 1lokal feil \le Tol
Avhenger avλh\lambda hhp+1h^{p+1} og problemets glatthet
Bryter man detløsningen eksplodererløsningen blir litt gal
Bedre orden hjelpernesten ikkemye

Det viktigste å forstå: en stabil løsning kan være svært unøyaktig, og en unøyaktig løsning kan være helt stabil. De to er uavhengige.
Bakover-Euler på y=10yy'=-10y med h=0,5h=0{,}5: R=11+5=16\displaystyle R=\frac{1}{1+5}=\frac16, altså stabil ✔. Men y1=160,167\displaystyle y_1=\frac16\approx 0{,}167 mens eksakt er e50,0067e^{-5}\approx 0{,}006725 ganger for stort. Stabil, og likevel ubrukelig.
Konklusjonen for praksis: velg hh etter det strengeste av de to kravene. For ikke-stive problemer er det nøyaktighet; for stive problemer med en eksplisitt metode er det stabilitet.

R(z)R(z) direkte fra Butcher-tabellen
For en eksplisitt metode kan RR leses ut av tabellen uten å regne stigningstall:

R(z)=1+zb ⁣(IzA)11,R(z)=1+z\,\mathbf b^{\!\top}\left(I-zA\right)^{-1}\mathbf 1,

der 1\mathbf 1 er vektoren med bare ettall.

For eksplisitte metoder forenkles det, siden AA er nilpotent: rekka (IzA)1=I+zA+z2A2+\left(I-zA\right)^{-1}=I+zA+z^{2}A^{2}+\dots har bare endelig mange ledd, og

R(z)=1+zibi+z2i,jbiaij+z3i,j,kbiaijajk+R(z)=1+z\sum_i b_i+z^{2}\sum_{i,j}b_ia_{ij}+z^{3}\sum_{i,j,k}b_ia_{ij}a_{jk}+\dots

Kjenner du igjen summene? Det er nesten ordensbetingelsene fra kap. 7.2, bare med cjc_j byttet ut med 1.

På eksamen er det som regel raskere å sette inn f=λyf=\lambda y direkte, slik vi gjorde i eksempel 1. Formelen er verdt å kjenne fordi den viser hvorfor RR blir et polynom av grad høyst ss for eksplisitte metoder.

Dette er «kjenne»-stoff.

Kontrollrutine for en O-oppgave

Seks sjekker:

1. Er f=λyf=\lambda y satt inn i ALLE stigningstallene? Et glemt kik_i gir feil grad på RR.
2. Er yny_n faktorisert helt ut? Det som står igjen, skal være en funksjon av zz alene.
3. Stemmer RR med eze^{z} til orden pp? Det er en gratis kontroll på utledningen.
4. Er R1\left|R\right|\le 1 delt i BEGGE ulikhetene, R1R\le 1 og R1R\ge -1?
5. For systemer: er egenverdiene kontrollert mot sporet og determinanten?
6. Er den MEST negative egenverdien brukt?

Punkt 3 er den beste kontrollen. Har du utledet R(z)=1+z+z22R(z)=1+z+\tfrac{z^{2}}{2} for en andreordens metode, stemmer det med eze^{z} til og med z2z^{2} ✔. Får du 1+z+z21+z+z^{2}, har du regnefeil.

✏️Eksamensnivå: skalar og system i samme oppgave
a) Utled stabilitetsfunksjonen til Heuns metode og finn stabilitetsintervallet.
b) Finn den største steglengden som gir stabil løsning av y=5yy'=-5y.
c) Systemet y=Ay\mathbf y'=A\mathbf y har egenverdier 2-2 og 10-10. Finn maks hh for Heun, og for RK4 med z=2,785\left|z^{*}\right|=2{,}785.
d) Vurder om systemet i c) er stivt, og anbefal en metode.
a) Utledningen. Sett f(t,y)=λyf(t,y)=\lambda y i Heuns metode:

k1=λyn,k_1=\lambda y_n,
k2=λ(yn+hk1)=λyn(1+z),z=λh.k_2=\lambda\left(y_n+hk_1\right)=\lambda y_n\left(1+z\right),\qquad z=\lambda h.

yn+1=yn+h2(k1+k2)=yn+hλ2yn(1+(1+z))=yn+z2yn(2+z).y_{n+1}=y_n+\frac h2\left(k_1+k_2\right)=y_n+\frac{h\lambda}{2}\,y_n\left(1+(1+z)\right)=y_n+\frac z2\,y_n\left(2+z\right).

yn+1=yn(1+z+z22)  R(z)=1+z+z22.y_{n+1}=y_n\left(1+z+\frac{z^{2}}{2}\right)\ \Longrightarrow\ R(z)=1+z+\frac{z^{2}}{2}.

Kontroll mot ez=1+z+z22+z36+e^{z}=1+z+\tfrac{z^{2}}{2}+\tfrac{z^{3}}{6}+\dots: de stemmer til og med z2z^{2}, som svarer til orden 2 ✔

Stabilitetsintervallet. Del R1\left|R\right|\le 1 i to:

Øvre: 1+z+z221z(1+z2)02z01+z+\tfrac{z^{2}}{2}\le 1\Leftrightarrow z\left(1+\tfrac z2\right)\le 0\Leftrightarrow -2\le z\le 0.

Nedre: 1+z+z221z2+2z+401+z+\tfrac{z^{2}}{2}\ge -1\Leftrightarrow z^{2}+2z+4\ge 0, som har diskriminant 416<04-16<0 og derfor alltid er oppfylt.

2z0\boxed{-2\le z\le 0}

b) Maks hh for y=5yy'=-5y. Med λ=5\lambda=-5 er z=5hz=-5h, og kravet z2z\ge -2 gir

5h2  h25=0,4.-5h\ge -2\ \Longrightarrow\ h\le\frac{2}{5}=0{,}4.

h0,4\boxed{h\le 0{,}4}

Kontroll ved h=0,4h=0{,}4: z=2z=-2, R=12+2=1R=1-2+2=1, altså R=1\left|R\right|=1 ✔ — grensen. Ved h=0,5h=0{,}5: z=2,5z=-2{,}5, R=12,5+3,125=1,625>1R=1-2{,}5+3{,}125=1{,}625>1

c) Systemet. Den mest negative egenverdien er λ=10\lambda=-10.

Heun (z=2\left|z^{*}\right|=2):

h210=0,2.h\le\frac{2}{10}=0{,}2.

RK4 (z=2,785\left|z^{*}\right|=2{,}785):

h2,78510=0,2785.h\le\frac{2{,}785}{10}=0{,}2785.

(Fra λ=2\lambda=-2 alene ville kravene vært h1h\le 1 og h1,39h\le 1{,}39 — begge mye slappere.)

d) Er systemet stivt?

S=102=5.S=\frac{10}{2}=5.

Nei. Et stivhetsforhold på 5 er svært moderat. De to komponentene lever på tidsskalaer som skiller med en faktor 5 — den raske dør ut på t0,5t\approx 0{,}5, den langsomme på t2,5t\approx 2{,}5.

Anbefaling: RK4 med eksplisitt steglengde. Begrunnelsen har to deler:

1. Stabiliteten er ikke plagsom. Kravet h0,28h\le 0{,}28 er langt fra dramatisk, og du trenger uansett en hh i den størrelsesorden for å følge den raske komponenten nøyaktig.
2. Nøyaktigheten er det som teller her, og der er RK4 langt overlegen — vi så i kap. 7.1 at forspranget er flere størrelsesordener.

Når ville jeg snudd? Hadde egenverdiene vært 2-2 og 2000-2000 (S=1000S=1000), ville RK4 måtte bruke h1,4103h\le 1{,}4\cdot 10^{-3} gjennom hele integrasjonen — kanskje tusenvis av skritt for å følge en komponent som er borte etter et øyeblikk. Da er bakover-Euler eller en annen A-stabil metode det riktige valget, selv om ordenen er lavere.

Tidsbruk på eksamen: a) 6 min, b) 3 min, c) 4 min, d) 4 min — omtrent 17 minutter for en oppgave på 10 poeng.

📝Oppgave 6
O

Bakover-Euler brukes på y=λyy'=\lambda y med λ<0\lambda<0.

a) Utled R(z)R(z).
b) Vis at R(z)1\left|R(z)\right|\le 1 for alle reelle z0z\le 0.
c) Regn ut RR for λ=10\lambda=-10 og h=0,5h=0{,}5, og sammenlikn med den eksakte faktoren eze^{z}.

📝Oppgave 7
O
En metode har stabilitetsfunksjonen

R(z)=1+z+z22+z36.R(z)=1+z+\frac{z^{2}}{2}+\frac{z^{3}}{6}.

a) Hvilken orden har metoden, og hvor mange steg må den minst ha?
b) Vis at z=2z=-2 ligger inne i stabilitetsintervallet, mens z=3z=-3 ikke gjør det.
c) Bestem den nedre enden av stabilitetsintervallet med to desimalers nøyaktighet.

Repetisjonsoppgaver
Din fremgang
0 / 4 oppgaver
Symbol- og formelliste

Dette kapitlet er skrevet av Anthropics toppmodeller (Claude Opus og Claude Fable) og er foreløpig ikke manuelt gjennomgått — kvalitetskontrollen gjøres av uavhengige KI-agenter, og innmeldte feil rettes fortløpende. Funnet en feil? Meld fra, så retter vi den. Les mer om hvordan innholdet lages.

Skolesaga er en uavhengig læringsressurs og er ikke tilknyttet eller godkjent av Norges teknisk-naturvitenskapelige universitet. Dette er ikke offisielt studiemateriell. Les mer.