Skip to main content

Uke 2 – fikspunktiterasjon

Mange likninger kan vi ikke løse med en formel, og slike likninger dukker opp overalt i ingeniørfag: en likevekt i en simulering, tidspunktet der to bevegelser møtes, eller en parameter som får en modell til å passe med målinger. Ta \(x^3+x-1=0\): Den har én løsning mellom 0 og 1, men ingen av formlene du kjenner, gir den direkte. Likevel finner en datamaskin løsningen med 16 riktige sifre på et øyeblikk. Trikset er å gjette og så forbedre gjetningen steg for steg, med samme regel hver gang. Det kaller vi en iterativ metode. Hver bruk av regelen er en iterasjon, og vi får en følge av anslag \(x_0,x_1,x_2,\ldots\). Denne uken lærer du å lage slike regler, å se om de virker, og å vite når du kan stoppe.

ViktigUkens spørsmål

Kan vi finne løsningen av \(x^3+x-1=0\) ved å gjenta en enkel regneregel – og hvordan vet vi når svaret er nøyaktig nok?

Vi skal prøve tre regler med nøyaktig samme løsning. De oppfører seg helt forskjellig, og resten av uken handler om hvorfor.

Vi bytter altså ut spørsmålet «kan vi skrive ned løsningen direkte?» med «kan vi lage en regel som gir stadig mindre feil?» Metoden fortsetter til et stoppkriterium er oppfylt, eller til den viser tegn på å mislykkes.

Modellproblemet vårt er \(f(x)=0\) med

\[f(x)=x^3+x-1.\]

En løsning kalles også et nullpunkt for \(f\). Før vi lager en iterasjonsregel, finner vi ut hvor løsningen ligger og om den er entydig. Cellen tegner \(f\) på \([0,1]\) og skriver ut fortegnene i endepunktene.

Grafen og fortegnene viser at en løsning ligger mellom \(0\) og \(1\). Siden

\[ f'(x)=3x^2+1>0, \]

er funksjonen strengt voksende og kan derfor ikke ha mer enn én løsning. Likningen har dermed nøyaktig én løsning i \([0,1]\).

Vi kaller den eksakte løsningen \(r\), altså \(r^3+r-1=0\). Vi kjenner ikke alle desimalene i \(r\). Som referanse bruker vi

\[r_{\mathrm{ref}}=0.6823278038280193.\]

Den er riktig i alle 16 desimaler, og algoritmene får aldri se den; vi bruker den bare til å måle feilen \(|x_n-r_{\mathrm{ref}}|\) når vi sammenligner metodene.

Gå i dybden: hvor nøyaktig er referanseverdien?

\(r_{\mathrm{ref}}\) er et binary64-tall. Sammenligner vi det med roten beregnet med 70 desimalsifres presisjon, har det lagrede binary64-tallet en absolutt referansefeil

\[ |r_{\mathrm{ref}}-r|<6\cdot10^{-17}. \]

Feil som er mindre enn referansefeilen, kan vi ikke måle pålitelig med denne referansen. Når en feilkurve senere kommer ned mot \(10^{-16}\), har vi derfor nådd grensen for hva referansen kan avgjøre.

For å iterere trenger vi en regel som tar inn \(x_n\) og gir \(x_{n+1}\). Vi kan isolere \(x\) på én side av likningen og få en ekvivalent likning

\[ x=g(x) \]

og deretter bruke høyresiden som oppdateringsregel:

\[ x_{n+1}=g(x_n). \]

Omskrivingen bestemmer dermed algoritmen. For modellproblemet kan vi for eksempel bruke

\[x=1-x^3,\qquad x=\sqrt[3]{1-x},\qquad x=\frac{2x^3+1}{3x^2+1}.\]

Alle tre er ekvivalente med modellikningen: Vi har bare flyttet ledd, tatt kubikkrot, eller lagt til \(2x^3+1\) på begge sider og delt med \(3x^2+1\), som alltid er positivt. Derfor har alle tre reglene \(r\) som fikspunkt.

Gå i dybden: hvorfor har de tre omskrivingene samme løsning?

Hvorfor beskriver alle tre den samme løsningen? Vi kan kontrollere det ved å starte med modellikningen og regne om:

\[ \begin{aligned} x^3+x-1=0 &\iff x=1-x^3,\\[2mm] x^3+x-1=0 &\iff x^3=1-x \iff x=\sqrt[3]{1-x},\\[2mm] x^3+x-1=0 &\iff 3x^3+x=2x^3+1 \iff x(3x^2+1)=2x^3+1 \iff x=\frac{2x^3+1}{3x^2+1}. \end{aligned} \]

I den siste linjen la vi først til \(2x^3+1\) på begge sider og satte \(x\) utenfor parentes. Deretter kan vi dele på \(3x^2+1\) fordi dette tallet alltid er positivt. Tegnet \(\iff\) betyr «hvis og bare hvis»: Regningen kan gå begge veier uten at vi legger til eller mister løsninger. Derfor har modellikningen og hver av de tre omskrevne likningene nøyaktig de samme løsningene.

Vi kaller nå de tre høyresidene \(g_1\), \(g_2\) og \(g_3\). Siden \(r\) er løsningen av modellikningen, gir omskrivingene

\[ r=g_1(r),\qquad r=g_2(r),\qquad r=g_3(r). \]

Et tall som ikke endrer seg når vi setter det inn i en regel \(g\), kalles et fikspunkt for regelen. Dermed er \(r\) det samme fikspunktet for alle tre reglene. Men for en \(x_n\) som ennå ikke er lik \(r\), gir de tre reglene vanligvis tre forskjellige neste verdier.

På papiret er reglene like gode. Men datamaskinen bruker ikke likningen – den bruker regelen, ett steg om gangen. Spørsmålet er om de tre følgene faktisk havner i \(r\). Funksjonen iterate, som gjentar regelen, ligger ferdig i oppsettet; vi skriver den selv i 2.1. Målet er å se om tre ekvivalente omskrivinger også gir tre like gode algoritmer.

Forsøk: én likning, tre regler

  • Kjør cellen og se hvilke følger som legger seg på den stiplede linjen.
  • Velg en annen startverdi i \([0,1]\), for eksempel \(0.2\) eller \(0.9\).
  • Endre antall steg fra 30 til 100, og undersøk om \(g_2\) når helt fram til løsningen.

Cellen definerer \(g_1\), \(g_2\) og \(g_3\), kjører 30 steg fra samme startverdi \(x_0=0.5\), skriver ut \(x_0,\ldots,x_8\) og den siste verdien, og tegner \(x_n\) mot \(n\). Den stiplede linjen er løsningen \(r\).

Bruk plottet til å svare på spørsmålene:

  • Hvilke av de tre følgene kommer nærmere den stiplede linjen?
  • Hvilken følge får minst absolutt feil etter åtte iterasjoner?
  • Hvilken følge ligger vekselvis over og under løsningen?
  • Blir avstanden til løsningen mindre når verdiene veksler side?

Skriv svarene kort i dine egne notater.

Hva viser forsøket? \(g_1\) hopper ut mot \(0\) og \(1\) og fortsetter å veksle mellom dem. \(g_2\) nærmer seg \(r\) fra begge sider, men sakte. \(g_3\) ligger på linjen etter tre–fire steg; etter fire steg er feilen bare omtrent \(3\cdot10^{-13}\). Både \(g_1\) og \(g_2\) veksler side, men bare for \(g_2\) krymper avstanden til løsningen; for \(g_1\) blir den stort sett større. Samme løsning, tre helt forskjellige algoritmer. Resten av uken forklarer hvorfor.

ViktigHovedidé

Å omskrive en likning er ikke bare algebra. Valget av \(x=g(x)\) bestemmer den numeriske algoritmen. Ekvivalente formler kan gi helt forskjellig oppførsel.

Etter denne uken skal du kunne

  • forklare hva et fikspunkt er, og hvordan en fikspunktiterasjon lager en følge av tilnærminger,
  • bruke iterasjonsverdier til å avgjøre om en metode nærmer seg en løsning, beveger seg bort, veksler mellom sider eller går i en syklus,
  • lese et spindelvevplott og bruke det til å forklare hvordan iterasjonen beveger seg,
  • bruke helningen \(g'(r)\) til å avgjøre om et fikspunkt trekker til seg nærliggende iterater eller skyver dem bort,
  • forklare forskjellen på konvergens nær et fikspunkt og konvergens fra et større område av startverdier,
  • sammenligne hvor raskt metoder nærmer seg løsningen ved hjelp av feilreduksjon og antall iterasjoner,
  • beregne og tolke absolutt feil, relativ feil, residual og skrittlengde,
  • programmere en fikspunktiterasjon som stopper på en kontrollert måte og oppgir hvorfor den stoppet.

I del 2.0 skrev vi om likningen \(f(x)=0\) til en ekvivalent likning \(x=g(x)\). Denne omskrivingen gjør mer enn å flytte ledd: Den gir en regel som kan ta et anslag \(x_n\) og lage det neste anslaget \(x_{n+1}=g(x_n)\). Målet er å finne en verdi som regelen ikke lenger flytter på.

NoteDefinisjon

Et tall \(r\) er et fikspunkt for \(g\) når

\[g(r)=r.\]

Navnet beskriver akkurat det som skjer: Setter vi \(r\) inn i \(g\), får vi det samme tallet tilbake. Verdien endres ikke når vi bruker \(g\).

For å finne fikspunktene til en gitt funksjon løser vi altså likningen \(g(x)=x\). For \(g(x)=x^2\) gir det \(x(x-1)=0\), altså fikspunktene \(0\) og \(1\).

Gå i dybden: fikspunktene til \(g(x)=x^2\) for hånd

Ta \(g(x)=x^2\) som et enkelt eksempel. Da må et fikspunkt oppfylle

\[ x^2=x. \]

Vi samler alt på én side og faktoriserer:

\[ x^2-x=x(x-1)=0. \]

Dermed er \(x=0\) og \(x=1\) de eneste fikspunktene. Vi kan kontrollere dem direkte: \(g(0)=0^2=0\) og \(g(1)=1^2=1\). Et tall som for eksempel \(2\) er ikke et fikspunkt, fordi \(g(2)=4\neq2\).

Når vi har valgt både en funksjon \(g\) og en startverdi \(x_0\), lager vi resten av følgen ved å bruke samme regel om igjen:

\[ x_{n+1}=g(x_n),\qquad n=0,1,2,\ldots \]

Dette kalles fikspunktiterasjon. Først beregner vi \(x_1=g(x_0)\), deretter \(x_2=g(x_1)\), og slik fortsetter vi. Treffer et iterat et fikspunkt eksakt, blir alle senere iterater lik dette fikspunktet: Hvis \(x_n=r\), får vi \(x_{n+1}=g(r)=r\).

Vanligvis treffer vi ikke fikspunktet eksakt. Da er spørsmålet om følgen kommer nærmere og nærmere et fikspunkt når vi fortsetter å iterere.

Ett iterat kan tilfeldigvis ligge nær løsningen. Det er derfor hele følgen som avgjør om metoden virker.

Vi sier at \(x_n\) konvergerer mot \(r\), og skriver \(x_n\to r\), når avstanden \(|x_n-r|\) går mot null. Mer presist: Uansett hvor liten positiv grense du velger, er avstanden under grensen fra et visst steg og utover. Iteratene trenger ikke bevege seg i én retning; det er avstanden til \(r\) som må bli mindre.

Følgene nedenfor viser fire typiske forløp, og ordene vi bruker om dem:

  • \(1,1/2,1/4,1/8,\ldots\) nærmer seg \(0\) fra én side. Avstanden til \(0\) halveres for hvert ledd: Følgen konvergerer.
  • \(1,-1/2,1/4,-1/8,\ldots\) ligger vekselvis på hver side av \(0\), men avstanden blir stadig mindre: Følgen alternerer, men konvergerer likevel mot \(0\).
  • \(1,-1,1,-1,\ldots\) gjentar de samme to verdiene: Følgen går i en syklus. Avstanden til én bestemt verdi går ikke mot null, så følgen konvergerer ikke.
  • \(1,2,4,8,\ldots\) beveger seg stadig lenger bort og vokser uten grense: Følgen divergerer.

Generelt sier vi at en følge divergerer når den ikke konvergerer, altså når den ikke nærmer seg én bestemt verdi. Også syklusen over divergerer i denne betydningen.

Alternering beskriver altså retningen på bevegelsen, ikke om metoden lykkes.

Grensen er et fikspunkt. Hvis følgen \(x_{n+1}=g(x_n)\) nærmer seg et tall \(r\), og \(g\) er kontinuerlig, er \(g(r)=r\). I spindelvevplottet i 2.2 ser du det: Trappen kan bare stoppe der kurven krysser diagonalen. Men dette resultatet lover ikke at følgen konvergerer – bare at grensen er et fikspunkt hvis den finnes.

Gå i dybden: hvorfor må grensen være et fikspunkt?

Anta nå at iteratene \(x_n\) nærmer seg et tall \(r\). Da nærmer også de neste iteratene \(x_{n+1}\) seg \(r\); vi har bare flyttet oss ett steg fram i den samme følgen. Samtidig er hvert neste iterat gitt av \(x_{n+1}=g(x_n)\).

Når \(g\) er kontinuerlig ved \(r\), vil verdiene \(g(x_n)\) nærme seg \(g(r)\) når \(x_n\) nærmer seg \(r\). Men \(g(x_n)\) og \(x_{n+1}\) er de samme tallene, og vi har nettopp sagt at disse nærmer seg \(r\). En følge kan ikke nærme seg to forskjellige tall. Derfor må

\[ r=g(r). \]

Dermed er \(r\) et fikspunkt. Argumentet forklarer hvorfor en konvergent fikspunktiterasjon gir en løsning av \(x=g(x)\). Det lover fortsatt ikke at enhver startverdi eller enhver omskriving vil få følgen til å konvergere.

Vi gjør formelen \(x_{n+1}=g(x_n)\) til et program: Vi lagrer den gjeldende verdien, bruker \(g\) på den én gang og gjentar. Koden nedenfor beregner et bestemt antall steg og beholder hele historikken, slik at vi kan undersøke forløpet etterpå. Vi prøver den på \(g_2\) fra \(x_0=0.5\).

Listen begynner med \(x_0\). Hver runde i løkken bruker det siste leddet til å beregne ett nytt ledd, så seks steg gir verdiene \(x_0,x_1,\ldots,x_6\). Historikken lar oss se om verdiene nærmer seg en grense, alternerer, gjentar seg eller beveger seg bort.

Øv på dette: følg tre steg av en fikspunktiterasjon.

Men når kan vi stoppe, og hvordan kontrollerer vi svaret? Det svarer vi på i 2.4 og 2.5.

Gå i dybden: et første blikk på residual og stopp

Hvis den siste beregnede verdien \(x\) ser ut til å være nær et fikspunkt, kontrollerer vi den med fikspunktresidualens størrelse

\[ |g(x)-x|. \]

For den opprinnelige likningen \(f(x)=0\) er den tilsvarende residualstørrelsen \(|f(x)|\). Absoluttverdien gjør begge til ikke-negative mål på hvor godt likningen er oppfylt. En liten residual er likevel ikke alene et bevis på at iterasjonen har konvergert eller at den absolutte feilen \(|x-r|\) er liten. Dette skillet blir presisert i del 2.4.

En datamaskin kan ikke beregne uendelig mange ledd. En praktisk metode trenger derfor en regel for når den skal stoppe, og et maksimalt antall iterasjoner som sikkerhetsgrense. Senere denne uken bygger vi inn disse kontrollene og lar programmet oppgi om det stoppet fordi endringen var liten, eller fordi det brukte opp alle tillatte iterasjoner.

Følg de tre første stegene

Her skal du følge den samme oppdateringsregelen steg for steg. I en fikspunktiterasjon beregnes hvert nytt ledd fra leddet rett før:

\[ x_{n+1}=g(x_n). \]

Startverdien \(x_0\) er utgangspunktet og teller ikke som en utført iterasjon. Første steg gir \(x_1=g(x_0)\), andre steg gir \(x_2=g(x_1)\), og tredje steg gir \(x_3=g(x_2)\).

I denne oppgaven er

\[ g(x)=\sqrt[3]{1-x} \]

Når uttrykket under rottegnet er negativt, skal den reelle kubikkroten brukes. Behold flere desimaler i mellomregningene, slik at tidlig avrunding ikke forplanter seg gjennom de neste stegene. Avrund bare \(x_3\) til tre desimaler.

En liste med iterater skjuler sammenhengen mellom hvert input og den tilhørende outputen. Et spindelvevplott viser den samme beregningen geometrisk, med to grafer:

  • Grafen \(y=g(x)\): For inputen \(x_n\) har punktet \((x_n,g(x_n))=(x_n,x_{n+1})\) høyden til neste iterat.
  • Diagonalen \(y=x\): Der er de to koordinatene like, så punktet \((x_{n+1},x_{n+1})\) gjør outputen til ny input.

Et fikspunkt ligger der grafene krysser hverandre, for der er \(g(r)=r\). Vi starter på diagonalen i \((x_0,x_0)\) og går

  1. loddrett til kurven for å lese av \(x_1=g(x_0)\),
  2. vannrett til diagonalen, slik at \(x_1\) blir ny \(x\)-koordinat,
  3. loddrett til kurven igjen for å få \(x_2=g(x_1)\), og så videre.

Slik dannes en trapp eller spiral. Trekker den seg inn mot et skjæringspunkt, nærmer iteratene seg et fikspunkt. Hjelperen cobweb_row tegner ett slikt plott per regel, side om side: kurven, diagonalen og den røde trappen. (For ett enkelt plott finnes cobweb_plot, som vi bruker i 2.5.) Tegnekoden ligger i oppsettet; det er trappen vi skal lese.

I 2.0 så vi på tallene; nå skal vi se dem tegnet. Vi tegner de tre reglene \(g_1\), \(g_2\) og \(g_3\) som spindelvev, med samme start \(x_0=0.5\). Målet er å se forskjellen på en trapp som trekker seg inn, en som går sakte i spiral innover, og en som skyves ut – og å finne ut hva ved kurven som avgjør det.

Forsøk: tre spindelvev

  • Kjør cellen og følg den røde stien i hvert panel.
  • Tell omtrent hvor mange trinn hver trapp trenger før den ser ut til å stå stille.
  • Bytt startverdien til \(0.9\) og undersøk om bildet for \(g_1\) endrer seg.
  • Prøv flere steg før du konkluderer: langsom bevegelse kan se ut som stillstand i et kort plott.

Cellen tegner de tre spindelvevene side om side: 20 steg for \(g_1\) og \(g_2\), og 8 steg for \(g_3\).

Hva viser forsøket? For \(g_1\) går trappen utover til en firkant mellom \(0\) og \(1\). For \(g_2\) går den i en sakte spiral innover. For \(g_3\) går den nesten rett inn. Se på helningen der kurven krysser diagonalen: \(g_1\) er bratt, \(g_2\) er slakere, og \(g_3\) er flat. Det er denne sammenhengen helningsregelen i neste fane handler om.

Tips

Et plott kan avdekke mønstre, foreslå hypoteser og finne moteksempler. Det viser bare hva som skjer for de valgte startverdiene og det valgte antallet iterasjoner, og er derfor ikke et bevis.

Spindelvevet viser hva som skjer for noen valgte startverdier. Nå vil vi forstå hvorfor: Hvor mye nærmere svaret bringer hvert steg oss?

La \(r\) være fikspunktet og \(x_n\) anslaget etter \(n\) iterasjoner. Feilen er

\[ e_n=x_n-r. \]

Hvis \(e_n>0\), ligger anslaget over \(r\); hvis \(e_n<0\), ligger det under. For eksempel betyr \(e_n=0.02\) at anslaget ligger \(0.02\) over fikspunktet. Det sentrale spørsmålet er: Hvordan avhenger neste feil \(e_{n+1}\) av den vi allerede har? Et fikspunkt er attraktivt når nærliggende iterater trekkes mot det, og frastøtende når de skyves bort.

Feilen etter neste steg er \(x_{n+1}-r\). Siden \(x_{n+1}=g(x_n)\) og fikspunktet oppfyller \(r=g(r)\), får vi den eksakte sammenhengen

\[ e_{n+1}=x_{n+1}-r=g(x_n)-g(r). \]

For å sammenligne ny og gammel feil kan vi, så lenge \(e_n\neq0\), dele den eksakte sammenhengen på \(e_n=x_n-r\):

\[ \frac{e_{n+1}}{e_n} =\frac{g(x_n)-g(r)}{x_n-r}. \]

Høyresiden er stigningstallet til sekanten gjennom punktene \((r,g(r))\) og \((x_n,g(x_n))\). Når \(x_n\) ligger nær \(r\), ligger punktene nær hverandre, og denne sekanthelningen er omtrent lik tangenthelningen \(g'(r)\). Derfor er

\[ \frac{e_{n+1}}{e_n}\approx g'(r). \]

Ganger vi med \(e_n\), får vi koblingen mellom feilene i to påfølgende iterater:

\[e_{n+1}\approx g'(r)e_n.\]

Tar vi absoluttverdier, får vi en enkel regel for størrelsen på neste feil:

\[|e_{n+1}|\approx |g'(r)|\,|e_n|.\]

ViktigLokal helningsregel
  • \(|g'(r)|<1\): attraktivt fikspunkt.
  • \(|g'(r)|>1\): frastøtende fikspunkt.
  • \(|g'(r)|=1\): regelen avgjør ikke.

Regelen gjelder lokalt, altså når \(x_n\) allerede er nær \(r\). Den sier ikke hva som skjer for en vilkårlig startverdi langt unna.

Når \(|g'(r)|=1\), avgjør ikke helningen alene; da er det leddene av høyere orden som avgjør, se Gå i dybden.

Gå i dybden: når helningen er akkurat 1

Tilfellet \(|g'(r)|=1\) viser hvorfor helningsregelen ikke alltid kan avgjøre saken alene. Se på to regler med samme fikspunkt \(r=0\):

\[ g_-(x)=x-x^3,\qquad g_+(x)=x+x^3. \]

Begge har

\[ g_-(0)=g_+(0)=0,\qquad g_{-}'(0)=g_{+}'(0)=1. \]

Helningsregelen gir derfor den samme første vurderingen for begge. Plottene under bruker små positive startfeil, men ulikt antall iterasjoner og ulike vinduer for å gjøre den langsomme bevegelsen synlig.

Her er \(r=0\), så feilen er ganske enkelt \(e_n=x_n\). Setter vi dette inn i reglene, får vi den nøyaktige feilutviklingen

\[ \begin{aligned} g_-:&\qquad e_{n+1}=e_n-e_n^3=e_n(1-e_n^2),\\ g_+:&\qquad e_{n+1}=e_n+e_n^3=e_n(1+e_n^2). \end{aligned} \]

For \(g_-\) og \(0<e_n<1\) ligger faktoren \(1-e_n^2\) mellom 0 og 1. Feilen blir derfor litt mindre i hvert steg, og den røde trappen beveger seg langsomt mot 0. For \(g_+\) er faktoren \(1+e_n^2\) større enn 1. Feilen blir litt større i hvert steg, og trappen beveger seg bort fra 0. Det er tredjegradsleddet, som helningen i selve fikspunktet ikke fanger opp, som skiller de to metodene.

Størrelsen \(|g'(r)|\) forteller om feilen krymper eller vokser. Fortegnet forteller hvilken side neste iterat havner på: En positiv helning bevarer vanligvis fortegnet til feilen, mens en negativ helning bytter fortegnet til feilen og gir bevegelse fram og tilbake over fikspunktet.

Hvor bratte er så de tre reglene våre i \(r\)? Det regner vi ut for hånd.

NoteHåndeksempel: helningene til \(g_1\), \(g_2\) og \(g_3\) i \(r\)

Utgangspunkt. \(r\approx0.6823\) oppfyller \(r^3+r-1=0\). Det gir to nyttige fakta: \(1-r=r^3\) og \(r^2\approx0.4656\).

Prøv selv. Deriver \(g_1\) og \(g_2\), og sett inn \(x=r\). Hvilken helning er størst i absoluttverdi?

Du trenger to derivasjonsregler. Kjerneregelen: Er \(g(x)=h(u(x))\), så er \(g'(x)=h'(u(x))\cdot u'(x)\). Brøkregelen: Er \(g(x)=\dfrac{u(x)}{v(x)}\), så er \(g'(x)=\dfrac{u'(x)v(x)-u(x)v'(x)}{v(x)^2}\).

Regnegangen.

  • \(g_1(x)=1-x^3\) gir \(g_1'(x)=-3x^2\). Altså er \(g_1'(r)=-3r^2\approx-3\cdot0.4656\approx-1.40\).
  • \(g_2(x)=(1-x)^{1/3}\) gir med kjerneregelen \(g_2'(x)=\tfrac13(1-x)^{-2/3}\cdot(-1)=-\dfrac{1}{3(1-x)^{2/3}}\). I \(r\) er \(1-r=r^3\), så \((1-r)^{2/3}=(r^3)^{2/3}=r^2\). Dermed er \(g_2'(r)=-\dfrac{1}{3r^2}\approx-\dfrac{1}{1.40}\approx-0.72\).
  • \(g_3(x)=\dfrac{2x^3+1}{3x^2+1}\) gir med brøkregelen \[g_3'(x)=\frac{6x^2(3x^2+1)-(2x^3+1)\cdot6x}{(3x^2+1)^2}=\frac{18x^4+6x^2-12x^4-6x}{(3x^2+1)^2}=\frac{6x\,(x^3+x-1)}{(3x^2+1)^2}.\] Telleren inneholder selve likningen \(x^3+x-1\), som er null i \(r\). Altså er \(g_3'(r)=0\).

Dette tar vi med oss. \(|g_1'(r)|\approx1.40\), \(|g_2'(r)|\approx0.72\) og \(g_3'(r)=0\). De to første er hverandres omvendte tall. Det er ingen tilfeldighet: \(g_2\) er den omvendte funksjonen til \(g_1\), så grafene er speilbilder om linjen \(y=x\). Det samme fikspunktet frastøter for den ene regelen og tiltrekker for den andre. \(g_3\) er flat i \(r\) fordi telleren i \(g_3'\) inneholder likningen vi løser.

For \(g_1\) vokser en liten feil altså med faktor omtrent \(1.40\) per steg, for \(g_2\) blir omtrent \(72\,\%\) av feilen igjen, og for \(g_3\) forsvinner den delen av feilen som er proporsjonal med \(e_n\): Nær fikspunktet bestemmes feilen da av mindre ledd av høyere orden.

Øv på dette: forutsi konvergens med helningen, og velg en parameter som gjør helningen null.

Forutsi konvergens med helningen

Likningen

\[ x^2-x-1=0 \]

har den positive løsningen \(r=\frac{1+\sqrt5}{2}\approx1.618\). Den kan skrives om til fikspunktform på flere måter. Vi ser på to av dem:

\[ x=1+\frac1x \qquad\text{og}\qquad x=x^2-1. \]

Begge har \(r\) som fikspunkt (kontroller det gjerne: begge følger av \(x^2=x+1\)). Kall høyresidene \(g_a\) og \(g_b\). Bruk helningsregelen fra 2.2 til å avgjøre hvilken regel som trekker nærliggende iterater inn mot \(r\). Sammenhengen \(r^2=r+1\) gjør regningen kortere.

Design en raskere iterasjon

Samme likning kan gi mange fikspunktiterasjoner, og parameteren i oppdateringsregelen kan gjøre forskjellen mellom langsom og rask konvergens. Her skal du velge parameteren med helningsregelen, uten å måtte prøve mange verdier numerisk.

En likning \(f(x)=0\) kan gjøres om til familien av fikspunktiterasjoner

\[ g_\lambda(x)=x-\lambda f(x), \]

der \(\lambda>0\) bestemmer hvor stor korreksjon vi tar i hvert steg. Før vi undersøker farten, kontrollerer vi at omskrivingen beholder løsningene. Hvis \(r\) er et nullpunkt for \(f\), er \(f(r)=0\), og innsetting gir

\[ g_\lambda(r) =r-\lambda f(r) =r-\lambda\cdot0 =r. \]

Dermed er hvert nullpunkt for \(f\) et fikspunkt for \(g_\lambda\), uansett hvilken positiv \(\lambda\) vi velger. Også motsatt vei gjelder: Hvis \(g_\lambda(r)=r\), får vi \(r-\lambda f(r)=r\), og siden \(\lambda>0\), må \(f(r)=0\). Nullpunktene og fikspunktene er altså de samme.

Parameteren flytter derfor ikke løsningen, men den endrer helningen og dermed den lokale konvergensfaktoren

\[ |g_\lambda'(r)|. \]

For \(f(x)=x^2-2\) er

\[ g_\lambda'(x)=1-2\lambda x. \]

En parameter som gjør \(g_\lambda'(r)=0\), fjerner den lineære delen av feilutviklingen og gir lokal faktor null. I denne oppgaven bruker vi det positive nullpunktet av \(f\).

Helningsregelen fra 2.2 beskriver hva som skjer når et iterat allerede er nær fikspunktet. Den sier ikke om en startverdi langt unna noen gang kommer dit.

  • En metode har lokal konvergens mot \(r\) dersom alle startverdier tilstrekkelig nær \(r\) gir konvergens mot \(r\). Det er dette \(|g'(r)|<1\) lover – uten å si hvor stort nærområdet er.
  • En metode har global konvergens på et område, for eksempel et intervall \([a,b]\), dersom den konvergerer for alle startverdier i området. «Global» må alltid kobles til et oppgitt område.

Mengden av startverdier som faktisk gir konvergens mot et bestemt fikspunkt, kalles fikspunktets tiltrekningsområde. Vi kan få et første bilde av det ved å prøve mange startverdier numerisk.

I eksperimentet nedenfor velger vi 400 startverdier mellom \(-1\) og \(1.5\). For hver startverdi utfører vi nøyaktig 100 iterasjoner med \(g_2\) og tegner ett punkt:

\[ (x_0,x_{100}). \]

Den vannrette koordinaten viser altså hvor iterasjonen startet, mens den loddrette koordinaten viser hvor den befant seg etter 100 steg. Den stiplede linjen markerer fikspunktet \(r\). Målet er å finne ut hvor langt unna \(r\) vi kan starte og fortsatt havne i \(r\) – og hva som skjer med startverdiene utenfor.

Forsøk: hvor kan vi starte?

  • Kjør cellen og finn området der punktene ligger på den stiplede linjen.
  • Regn ut final_value(g2, -0.01) og final_value(g2, 0.01) selv.
  • Endre antall steg fra 100 til 101, og se hva som skjer med de to ytre nivåene (forklaringen står i Gå i dybden under figuren).

Funksjonen final_value gjentar regelen et fast antall ganger og gir bare den siste verdien. Løkken bruker den på alle 400 startverdiene.

Hva viser forsøket? Den midterste, vannrette delen av figuren ligger oppå den stiplede linjen. Det betyr at startverdier mellom omtrent \(0\) og \(1\) har gitt \(x_{100}\approx r\). Eksperimentet tyder derfor på at intervallet \((0,1)\) ligger i tiltrekningsområdet til fikspunktet.

De to ytre delene ligger nær \(-0.536\) og \(1.154\), ikke nær fikspunktet. Disse to verdiene hører sammen: \(g_2(-0.536)\approx1.154\) og \(g_2(1.154)\approx-0.536\). Iterasjonen nærmer seg dermed en 2-syklus og alternerer mellom verdiene.

Gå i dybden: hvorfor bytter de ytre nivåene plass, og hvorfor åpent intervall?

Fordi figuren viser resultatet etter et partall steg, vises den negative delen av syklusen for startverdier til venstre og den positive delen for startverdier til høyre. Med 101 steg ville de to ytre nivåene ha byttet plass. Endepunktene \(0\) og \(1\) danner dessuten 2-syklusen \(0\mapsto1\mapsto0\); intervallet vi observerer konvergens fra, er derfor det åpne intervallet \((0,1)\).

Legg merke til hva som skjer rundt \(x_0=0\): \(0.01\) havner i \(r\), mens \(-0.01\) havner i 2-syklusen. En ørliten endring i startverdien kan gi et helt annet resultat. Denne følsomheten for startverdien skal du undersøke i det komplekse planet i prosjektet.

«Verdien etter 100 steg» er altså ikke automatisk en løsning. Og forsøket er ikke et bevis: Vi har prøvd 400 startverdier og 100 steg, så en startverdi som konvergerer langsomt, kan bli feilklassifisert, og smale områder mellom punktene kan bli oversett. Et bevis for global konvergens må kontrollere alle startverdier i det oppgitte området.

I 2.2 fant vi den lokale feilmodellen for fikspunktiterasjon

\[ |e_{n+1}|\approx |g'(r)|\,|e_n|, \qquad e_n=x_n-r. \]

Én iterasjon multipliserer altså feilen med omtrent samme faktor når \(x_n\) er nær \(r\). Vi kaller faktoren \(S\):

\[|e_{n+1}|\approx S|e_n|,\qquad 0<S<1.\]

Da sier vi at metoden har lineær konvergens, og \(S\) kalles konvergensfaktoren: andelen av feilen som blir igjen etter ett steg. For fikspunktiterasjon er \(S\approx|g'(r)|\), men samme mønster dukker opp i mange andre metoder, og da kan vi anslå \(S\) fra forholdet \(|e_{n+1}|/|e_n|\). Over \(k\) steg gir modellen

\[ |e_{n+k}|\approx S^k|e_n|. \]

Jo mindre \(S\), desto raskere går det.

«Fart» kan også bety antall iterasjoner til en ønsket nøyaktighet, eller faktisk regnetid. Her måler vi bare feilreduksjonen per iterasjon, altså forholdet \(|e_{n+1}|/|e_n|\).

Gå i dybden: tre mål på fart

Det finnes flere meningsfulle måter å sammenligne «fart» på:

  • Antall iterasjoner til ønsket nøyaktighet: Hvor mange steg trengs før \(|x_n-r|\) er mindre enn en valgt grense?
  • Feilreduksjon per iterasjon: Hvor stor andel av feilen blir igjen fra ett steg til det neste?
  • Faktisk regnetid eller arbeidsmengde: Hvor lang tid tar beregningen, eller hvor mange operasjoner krever den?

Disse målene kan rangere metoder forskjellig. En metode kan bruke få iterasjoner, men ha dyre iterasjoner. I resten av denne delen undersøker vi bare den matematiske feilreduksjonen per iterasjon. Vi måler absolutt feil \(|e_n|=|x_n-r|\) og sammenligner forholdet \(|e_{n+1}|/|e_n|\). Vi måler ikke faktisk regnetid her.

Gå i dybden: hvordan leser vi tallet \(S\)?
  • \(S\approx0.9\): omtrent \(90\%\) av feilen er igjen.
  • \(S\approx0.5\): feilen halveres omtrent.
  • \(S\approx0.1\): vi vinner omtrent ett nytt korrekt desimalsiffer per steg.
Note

«Lineær» betyr ikke at \(x_n\) danner en rett linje mot \(n\). Det betyr at hver iterasjon fjerner omtrent samme andel av feilen.

En annen feilmodell er

\[|e_{n+1}|\approx C|e_n|^2.\]

Dette kalles kvadratisk konvergens. For små feil gjør kvadreringen at feilen kan avta langt mer per iterasjon enn ved en fast lineær faktor. Når modellen gjelder og konstanten \(C\) ikke motvirker reduksjonen, dobles antall korrekte sifre omtrent i hvert steg etter hvert som feilen blir liten. Med \(C=1\) blir feilene for eksempel \(10^{-2}\to10^{-4}\to10^{-8}\to10^{-16}\): antall riktige desimaler dobles i hvert steg.

Tenk at en iterasjon beregner en temperatur som skal sendes videre til en regulator. Programmet kan skrive ut mange desimaler, men det betyr ikke at alle er pålitelige. Hvor mange desimaler kan vi stole på?

Med \(p\) desimaler ligger de mulige avrundede verdiene med avstand \(10^{-p}\) fra hverandre; med tre desimaler er for eksempel avstanden mellom \(2.718\) og \(2.719\) lik \(10^{-3}\). Rundt hver verdi ligger et intervall av tall som avrundes til nettopp den, og grensen mellom to nabointervaller ligger halvveis mellom de to verdiene.

Derfor er halvparten av avstanden, \(\frac12 10^{-p}\), den naturlige feilgrensen. Dersom tilnærmingen \(x_n\) ligger mindre enn denne avstanden fra den eksakte løsningen \(r\), sier vi at den har minst \(p\) korrekte desimaler:

\[ |x_n-r|<\frac12\,10^{-p}. \]

«Minst tre korrekte desimaler» betyr altså \(|x_n-r|<0.0005\). Vi teller ikke sifre som ser like ut på skjermen; vi kontrollerer den faktiske avstanden til løsningen. Denne forskjellen er viktig nær et avrundingsskille, der to svært nære tall kan få ulik siste desimal etter avrunding.

Gå i dybden: les korrekte desimaler fra en figur

Nå bruker vi en figur til å gjøre den samme kvalitetskontrollen for tre beregnede resultater. For et resultat \(x_n\) definerer vi feilen som

\[ e=x_n-r, \]

der \(r\) er den eksakte løsningen. Fortegnet forteller hvilken side av løsningen vi ligger på. Her er vi bare interessert i avstanden til løsningen og plotter derfor den absolutte feilen \(|e|=|x_n-r|\).

Tenk at tre beregninger, A, B og C, gir hver sin tilnærming til det samme tallet \(r\):

\[ \begin{aligned} x_A&=r+3.2\cdot10^{-2}, &e_A=x_A-r&= 3.2\cdot10^{-2},\\ x_B&=r-4.0\cdot10^{-4}, &e_B=x_B-r&=-4.0\cdot10^{-4},\\ x_C&=r+2.0\cdot10^{-7}, &e_C=x_C-r&= 2.0\cdot10^{-7}. \end{aligned} \]

B ligger litt under den eksakte løsningen og får derfor negativ feil, mens A og C ligger over. Når vi skal vurdere nøyaktigheten, er retningen uviktig. Vi bruker avstandene \(|e_A|\), \(|e_B|\) og \(|e_C|\).

Som i vanlige feilplott bruker vi en logaritmisk \(y\)-akse: store feil ligger øverst, og mindre feil ligger lenger ned. Hver vannrette kravlinje viser \(|e|=\frac12 10^{-p}\). Et punkt under linjen for \(p=3\) har liten nok feil til at vi kan stole på minst tre desimaler. Vi leser svaret ved å finne den nederste kravlinjen punktet har passert.

Følg hvert punkt nedover fra toppen. Det blå punktet har passert kravet til én korrekt desimal, men ikke kravet til to. Det oransje har passert tre-desimallinjen, men ikke fire-desimallinjen. Det grønne har passert seks-desimallinjen, men ikke sju-desimallinjen. Vi kan derfor stole på henholdsvis 1, 3 og 6 desimaler i de tre resultatene.

Denne figuren stopper ved \(10^{-8}\) fordi alle tre eksempelfeilene ligger i dette området. I et konvergensplott for binary64 fortsetter vi vanligvis aksen ned mot omtrent \(10^{-16}\) for tall av størrelse 1. Der flater feilkurven ofte ut fordi maskinen ikke kan lagre tettere tall på denne skalaen.

Nå kan vi bruke kriteriet på modellproblemet. Feilene til \(g_2\) og \(g_3\) spenner fra tideler til nesten \(10^{-16}\), og på en vanlig loddrett akse blir de minste verdiene presset sammen nær null. Et semilogaritmisk plott, ofte kalt et semilogy-plott i Python, løser dette ved å bruke

  • en vanlig, lineær akse for iterasjonsnummeret \(n\), og
  • en logaritmisk akse for den absolutte feilen \(|e_n|\).

På den logaritmiske aksen betyr like store loddrette avstander like store faktorer, ikke like store differanser. Ett hovedtrinn ned fra \(10^{-3}\) til \(10^{-4}\) betyr at feilen er delt på \(10\). En rett, synkende linje betyr derfor at feilen reduseres med omtrent samme faktor i hvert steg. En kurve som blir stadig brattere, kan tyde på raskere konvergens, for eksempel kvadratisk konvergens.

Som i 2.0 måler vi feilen mot referanseverdien \(r_{\mathrm{ref}}\). Målet er å se forskjellen på lineær og kvadratisk konvergens i ett bilde, og å måle konvergensfaktoren \(S\) for \(g_2\).

Forsøk: hvor fort går det?

  • Kjør cellen og sammenlign formen på de to kurvene.
  • Les av omtrent hvor mange steg \(g_2\) trenger for å vinne én ny desimal.
  • Bytt startverdien til \(0.9\) og undersøk om formen på kurvene endrer seg.

Cellen tegner den absolutte feilen \(|x_n-r_{\mathrm{ref}}|\) for \(g_2\) og \(g_3\) i 20 steg fra \(x_0=0.5\). Punkter der feilen er nøyaktig null, kan ikke tegnes på en logaritmisk akse og blir utelatt.

På den logaritmiske aksen går \(g_2\) nesten rett nedover. En rett linje betyr at feilen ganges med omtrent det samme tallet i hvert steg. Kurven for \(g_3\) blir stadig brattere: Etter fire steg er feilen omtrent \(3\cdot10^{-13}\), og etter det er den nede ved grensen for hva referansen kan måle. Den neste cellen regner ut denne faktoren for \(g_2\), altså forholdet \(|e_n|/|e_{n-1}|\). Målet er å sjekke om det nærmer seg \(|g_2'(r)|\approx0.72\) fra 2.2.

Hva viser forsøket? Den siste kolonnen er andelen av forrige feil som er igjen etter ett nytt steg. I starten veksler forholdet mellom omtrent \(0.61\) og \(0.82\). \(g_2'(r)\) er negativ, så iteratene havner annenhver gang over og under \(r\), og sekanthelningen fra 2.2 er litt ulik på de to sidene. Etter hvert legger forholdet seg nær \(|g_2'(r)|\approx0.716\). Konvergensfaktoren \(S\) kan altså både regnes ut med derivasjon og måles fra en tabell. Med \(S\approx0.72\) trengs omtrent sju steg per ny desimal, fordi \(0.716^7\approx0.1\).

Øv på dette: les korrekte desimaler fra et semilogy-plott og mål konvergensfaktoren fra en feiltabell.

Les nøyaktighet fra plott

Husk fra 2.3: på en logaritmisk akse betyr like store avstander like store faktorer.

Vi sammenligner tre fikspunktiterasjoner som alle har fikspunktet \(r=0\):

Kurve Oppdateringsregel Startverdi Forventet oppførsel nær \(0\)
Blå \(x_{n+1}=0.8x_n\) \(x_0=0.1\) langsom lineær konvergens
Oransje \(x_{n+1}=0.1x_n\) \(x_0=0.1\) rask lineær konvergens
Grønn \(x_{n+1}=x_n^2\) \(x_0=0.3\) kvadratisk konvergens

For de lineære reglene blir feilen multiplisert med henholdsvis \(0.8\) og \(0.1\) i hvert steg. For den grønne kurven er \(e_{n+1}=e_n^2\), så antallet korrekte desimaler vokser mye raskere når feilen først er liten.

Cellen regner ut de tre feilhistoriene og tegner dem på en logaritmisk akse.

Bruk kriteriet

\[ |e_n|<\frac12 10^{-p} \]

og finn det største heltallet \(p\) som oppfyller ulikheten. Dette er antallet korrekte desimaler ved iterasjon \(n\).

«Nye desimaler fra steg 4 til 5» betyr

\[ p_5-p_4, \]

altså antallet korrekte desimaler etter steg 5 minus antallet etter steg 4. Les først omtrentlige feil fra plottet, og bruk tierpotensene på den logaritmiske aksen til å velge riktig grense.

Den interaktive grafen nedenfor viser den raske lineære iterasjonen \(x_{n+1}=0.1x_n\) med \(x_0=0.1\). Det er den samme oransje kurven som i semilogy-plottet, forlenget med noen steg.

Den loddrette koordinaten er \(\log_{10}(|e_n|)\). For eksempel betyr høyden \(-4\) at feilen er \(10^{-4}\). Klikk på datapunktet som har sju korrekte desimaler etter kriteriet \(|e_n|<\frac12 10^{-7}\). Det valgte punktet markeres med rødt.

Den grønne iterasjonen \(x_{n+1}=x_n^2\) med \(x_0=0.3\) konvergerer kvadratisk mot fikspunktet \(r=0\). Her finnes det ikke nødvendigvis et iterat med akkurat sju korrekte desimaler: Ved kvadratisk konvergens kan antallet korrekte desimaler øke med flere i ett steg.

Grafen viser den samme grønne kurven som i semilogy-plottet. Klikk på det første punktet som har minst sju korrekte desimaler. Sammenlign feilen med grensen \(\frac12 10^{-7}=5\cdot10^{-8}\). Det valgte punktet markeres med rødt.

Hvor lenge må vi iterere?

En fikspunktiterasjon gir disse absolutte feilene:

\(n\) \(0\) \(1\) \(2\) \(3\)
\(\lvert e_n\rvert\) \(4.0\cdot10^{-2}\) \(2.4\cdot10^{-2}\) \(1.44\cdot10^{-2}\) \(8.64\cdot10^{-3}\)

Anta at modellen \(|e_{n+1}|\approx S|e_n|\) fra 2.3 fortsetter å gjelde med samme \(S\). Kriteriet for fire korrekte desimaler er \(|e_n|<\frac12 10^{-4}\).

I uke 1, del 1.0 målte vi absolutt feil \(|x-r|\) og relativ feil \(|x-r|/|r|\) (når \(r\neq0\)), der \(r\) er den eksakte løsningen og \(x\) en tilnærming. Begge svarer på spørsmålet: Hvor langt ligger svaret fra den riktige løsningen? Det kaller vi foroverfeil – vi går «forover» gjennom beregningen, fra problemet vi ga programmet, fram til feilen i resultatet.

AdvarselNår den absolutte feilen ikke kan beregnes

I en virkelig beregning kjenner vi vanligvis ikke \(r\). Dermed kan vi heller ikke beregne \(|x-r|\). Hvis vi kjente løsningen, trengte vi ofte ikke metoden.

Ved testing kan vi i stedet måle feilen mot en referanseverdi, \(|x-r_{\mathrm{ref}}|\). Da må referansefeilen være oppgitt og mindre enn feilene vi ønsker å undersøke.

Siden vi sjelden kjenner løsningen, trenger vi mål som kan beregnes uten den. For likningen \(f(x)=0\) setter vi tilnærmingen \(x\) inn i venstresiden. Tallet vi får,

\[ R(x)=f(x), \]

kalles residualen. Vi bruker vanligvis residualstørrelsen \(|R(x)|=|f(x)|\). Den måler hvor langt den innsatte verdien er fra å gi null i selve likningen. Residualen kan beregnes uten å kjenne løsningen \(r\), men den er ikke det samme som absolutt feil. For \(x=g(x)\) er fikspunktresidualen \(g(x)-x\).

Residualen leder til en annen måte å vurdere resultatet på: bakoverfeil. I stedet for å spørre hvor feil svaret er, spør vi

Hvor mye måtte vi endre det opprinnelige problemet for at den beregnede verdien \(x\) skulle vært et eksakt svar?

For en likning er svaret enkelt: Flytter vi grafen loddrett med residualen \(f(x)\), blir \(x\) et eksakt nullpunkt. Residualstørrelsen \(|f(x)|\) er altså en bakoverfeil.

Gå i dybden: bakoverfeil i mer generelle problemer

For likningen \(f(t)=0\) har den beregnede verdien \(x\) residual \(R(x)=f(x)\). Hvis vi endrer likningen til

\[ \widetilde f(t)=f(t)-f(x)=0, \]

blir \(x\) en eksakt løsning, fordi \(\widetilde f(x)=0\). Residualstørrelsen \(|f(x)|\) er derfor størrelsen på en konstant, loddrett endring av grafen som gjør det beregnede svaret eksakt. Med denne bestemte måten å endre problemet på kan residualen tolkes som en additiv bakoverfeil.

I mer generelle problemer betyr bakoverfeil den minste tillatte endringen i dataene, koeffisientene eller modellen som gjør det beregnede svaret eksakt. Da er den rå residualen ikke nødvendigvis hele bakoverfeilen; den må ofte skaleres og knyttes til hvilke deler av problemet som faktisk kan endres.

Geometrisk måler absolutt foroverfeil og residual avstander i to forskjellige retninger. Den absolutte feilen \(|x-r|\) er den vannrette avstanden fra tilnærmingen \(x\) til roten \(r\). Residualstørrelsen \(|f(x)|\) er den loddrette avstanden fra punktet \((x,f(x))\) på grafen til \(x\)-aksen.

Figuren nedenfor bruker den flate funksjonen

\[ f(x)=0.02(x-3), \]

som har roten \(r=3\). Dra det blå punktet langs \(x\)-aksen. Det oransje linjestykket viser absolutt feil, mens det lilla viser residualstørrelsen. Ved \(x=8\) er den absolutte feilen \(|8-3|=5\), men residualen er bare \(|f(8)|=0.1\).

Ved \(x=8\) kan vi også lese figuren som forover- og bakoveranalyse:

  • Forover: Det beregnede svaret ligger 5 enheter fra den egentlige roten.
  • Bakover: Hvis vi flytter funksjonsgrafen ned med \(0.1\), altså bruker \(\widetilde f(t)=f(t)-0.1\), blir \(x=8\) en eksakt rot i det endrede problemet.

Et resultat kan dermed ha liten bakoverfeil, men stor foroverfeil. Det skjer her fordi den flate grafen gjør at en liten loddrett endring flytter roten langt vannrett. En bratt graf kan gi det motsatte forholdet. Residualen forteller hvor godt likningen er oppfylt på funksjonens skala; den er ikke automatisk et direkte mål på avstanden til roten.

Endringen \(|x_{n+1}-x_n|\) kalles skrittlengde. Den svarer på spørsmålet: Hvor mye flyttet algoritmen seg i siste steg?

For fikspunktiterasjon er

\[|x_{n+1}-x_n|=|g(x_n)-x_n|.\]

Skrittlengden er derfor også fikspunktresidualen i \(x_n\), men tolkningene er ulike. Et lite skritt kan bety at vi er nær løsningen, at konvergensen er langsom, eller at flyttallsavrunding hindrer videre bevegelse.

Størrelse Formel Spørsmål
Absolutt feil (foroverfeil) \(|x-r|\) Hvor langt er vi fra løsningen?
Relativ feil (relativ foroverfeil) \(|x-r|/|r|\) Hvor stor er feilen relativt til svaret?
Residual (grunnlag for bakoverfeil) \(|f(x)|\) Hvor godt er likningen oppfylt?
Fikspunktresidual \(|g(x)-x|\) Hvor godt er fikspunktlikningen oppfylt?
Skrittlengde \(|x_{n+1}-x_n|\) Hvor mye flytter algoritmen seg?

Øv på dette: skrittlengde og feil i to tallformater, og tre kjøringer med ulik stoppårsak i 2.5.

En toleranse er ikke en feil. Den er en tallfestet grense vi velger for et bestemt mål, for eksempel absolutt feil, residual eller skrittlengde. Vi må alltid oppgi både målet og toleransen.

En absolutt toleranse atol setter en grense i samme skala som tallet:

\[|x_{n+1}-x_n|\leq\mathrm{atol}.\]

En relativ toleranse rtol skalerer grensen med størrelsen på svaret:

\[|x_{n+1}-x_n|\leq\mathrm{rtol}|x_{n+1}|.\]

Relativ toleranse fungerer dårlig nær null; absolutt toleranse alene tar ikke hensyn til svært store svar. Derfor kombineres de ofte:

\[|x_{n+1}-x_n|\leq\mathrm{atol}+\mathrm{rtol}|x_{n+1}|.\]

Cellen regner ut den samlede toleransen for tre svar av svært ulik størrelse. Se hvilken del som bestemmer grensen ved \(x=0\), og hvilken som bestemmer den ved \(x=10^6\).

Note
  • atol bestemmer hvilke absolutte endringer vi bryr oss om, særlig nær null.
  • rtol tilpasser kravet til størrelsen på svaret.
  • Toleransen er et brukerkrav, ikke en garanti for like liten absolutt feil.

Toleranse og tallformat

Fra uke 1, del 1.3 vet du at nabotallene rundt \(1\) ligger omtrent \(2.2\cdot10^{-16}\) fra hverandre i binary64. Ber vi om en skrittlengde under \(10^{-20}\), gir det ingen flere riktige sifre. Neste lagrede skritt er enten null eller omtrent \(10^{-16}\), og i verste fall veksler iterasjonen mellom to nabotall til det maksimale antallet iterasjoner er brukt opp. Velg derfor toleransen etter tallformatet.

Gå i dybden: binary64 og toleranser

Python float bruker vanligvis IEEE 754 binary64. Rundt \(1\) er avstanden til neste representerbare tall omtrent

\[2^{-52}\approx2.22\cdot10^{-16}.\]

Du kjenner maskinepsilon fra uke 1, del 1.3. Cellen minner om hva den betyr: \(1+\varepsilon\) kan lagres, men \(1+\varepsilon/2\) blir avrundet tilbake til \(1\).

ViktigIkke en absolutt grense på \(10^{-16}\)

Binary64 har omtrent 15–16 signifikante desimalsifre. Rundt \(1\) gir dette en oppløsning på omtrent \(10^{-16}\). Den absolutte avstanden mellom nabotall vokser for store tall og minker for små tall. Derfor passer relativ toleranse naturlig med flyttallsrepresentasjonen.

Som vi så over, gir flere iterasjoner ikke automatisk flere pålitelige sifre. Toleransen må velges i samsvar med tallformatet og nøyaktigheten som problemet faktisk krever.

Når tallformatet setter en grense

Vi bruker iterasjonen for \(\sqrt{2}\) som vi bygger i 2.5,

\[ x_{n+1}=\frac12\left(x_n+\frac{2}{x_n}\right), \qquad x_0=1, \]

og utfører nøyaktig de samme regneoperasjonene i to tallformater:

  • binary32, representert med numpy.float32, har omtrent 7 signifikante desimalsifre,
  • binary64, representert med vanlig Python float, har omtrent 15–16 signifikante desimalsifre.

Referanseverdien beregnes med 80 desimalsifres presisjon og brukes bare til å måle feil. Hvert iterat i de to kjøringene blir derimot lagret i det angitte tallformatet.

Begge kurvene faller raskt i starten fordi iterasjonen konvergerer kvadratisk nær løsningen (det forklarer vi i 2.5). Etter noen steg blir hver kurve vannrett. Dette feilplatået betyr at tallformatet ikke kan lagre en bedre tilnærming med denne iterasjonen. Flere steg gjentar bare samme lagrede verdi.

Legg merke til forskjellen mellom skrittlengde og feil mot referanseverdien: Når to påfølgende lagrede verdier er like, er skrittlengden \(|x_{n+1}-x_n|=0\). Feilen mot den mer presise referansen er likevel ikke null. Et stopp på grunn av null skrittlengde betyr derfor «ingen videre endring i dette tallformatet», ikke «eksakt løsning».

Les følgende fra semilogy-plottet:

  1. Starten på platået er den første iterasjonen der feilen har nådd nivået den beholder i alle senere steg.

  2. På platået finner du antall korrekte desimaler med kriteriet

    \[ |x_n-r|<\frac12 10^{-p}. \]

  3. Første nullskritt er iterasjonen \(n\) der den nye lagrede verdien for første gang er lik den forrige, altså \(x_n=x_{n-1}\).

En formel alene er ikke en ferdig algoritme. Programmet må også vite hvor det skal starte, når det skal stoppe, og hva det skal gjøre dersom iterasjonen ikke oppfører seg som forventet.

Funksjonen under gjentar fire handlinger: Den beregner et nytt iterat, kontrollerer at beregningene gir endelige tall, lagrer skrittlengde og residual, og sammenligner skrittlengden med toleransen. En liten hjelpefunksjon samler resultatene slik at vi slipper å gjenta den samme return-blokken. Linjen residual = lambda value: … lager en kort funksjon på én linje: Hvis likningen f er gitt som argument, måler den \(|f(x)|\); ellers måler den fikspunktresidualen \(|g(x)-x|\).

Resultatet er en ordbok. x er siste tilnærming, mens reason og iterations forteller hvorfor og når programmet stoppet. Listene history, steps og residuals lar oss kontrollere utviklingen i etterkant. Ved regnefeil eller ikke-endelige tall returnerer funksjonen siste godkjente iterat; iterations teller fullførte steg. Et mislykket steg lagres ikke. Hvis residualen ikke kan beregnes allerede i startpunktet, markeres den med nan.

Øv på dette: les tre kjøringer som en kontrollrapport.

Stoppmeldingen «lite skritt» sier bare at verdien sluttet å endre seg. Derfor returnerer fixed_point både skrittlengden og residualen, og vi bør se på begge. I 2.4 så vi at de svarer på forskjellige spørsmål. maxiter hindrer en divergent eller syklisk metode i å kjøre for alltid, og kontrollen av nan og inf hindrer at ugyldige tall blir rapportert som en løsning. Figurene i Gå i dybden-boksen under og i oppgaven «Kan resultatet sendes videre?» viser begge størrelsene på logaritmisk akse – samme type figur som du lager for Newton i prosjekt 2.

Gå i dybden: skrittlengde og residual for \(g_2\)

Vi bruker først \(g_2\) på modellproblemet. Stoppmeldingen alene sier bare at skrittlengden ble liten. Cellen kjører derfor fixed_point, skriver ut resultatet med residualen og plotter både skrittlengde og residual gjennom iterasjonen.

Skrittlengde og residual avtar ikke nødvendigvis i samme takt. Det er derfor nyttig å returnere begge, selv om stoppregelen i denne funksjonen bare bruker skrittlengden.

Til slutt ser vi på en regel der helningen i fikspunktet er null, som for \(g_3\). Målet er å se hvor fort det da går.

For å beregne \(\sqrt a\) kan vi løse \(x^2-a=0\). Omskrivingen \(x=a/x\) gir

\[x_{n+1}=\frac{a}{x_n}.\]

La oss prøve regelen med \(a=2\) og \(x_0=1\).

Følgen går i syklusen \(1,2,1,2,\ldots\), så skrittlengden blir aldri mindre enn \(1\). For å få en oppdatering som demper denne vekslingen, tar vi gjennomsnittet av \(x\) og \(a/x\):

\[x_{n+1}=\frac12\left(x_n+\frac{a}{x_n}\right).\]

Hvorfor gjennomsnittet? Produktet av de to tallene er alltid \(x\cdot\frac ax=a\). Er \(x\) for liten, er \(\frac ax\) for stor, og omvendt. Den riktige roten \(\sqrt a\) ligger derfor alltid mellom \(x\) og \(\frac ax\), og midtpunktet er et naturlig nytt anslag.

Cellen lager regelen for en vilkårlig \(a\), kjører fixed_point med \(a=2\) og skriver ut både iteratet og residualen \(|x^2-2|\) i hvert steg.

Hvorfor går det så fort? Regelen er \(g(x)=\tfrac12\bigl(x+\tfrac ax\bigr)\), så \(g'(x)=\tfrac12\bigl(1-\tfrac{a}{x^2}\bigr)\). I \(x=\sqrt a\) er \(a/x^2=1\), og dermed \(g'(\sqrt a)=0\) – akkurat som for \(g_3\) i 2.2. Den lineære delen av feilen forsvinner (helningsregelen fra 2.2), og vi får kvadratisk konvergens (2.3). I tabellen over faller residualen \(|x_n^2-2|\) fra \(6.9\cdot10^{-3}\) til \(6.0\cdot10^{-6}\) og \(4.5\cdot10^{-12}\): eksponenten dobles omtrent i hvert steg, og det gjør antall riktige sifre også.

Gå i dybden: hvorfor kvadreres residualen?

For regelen \(x_{n+1}=\dfrac12\left(x_n+\dfrac{a}{x_n}\right)=\dfrac{x_n^2+a}{2x_n}\) kan vi regne ut neste residual eksakt:

\[ \begin{aligned} x_{n+1}^2-a &=\frac{(x_n^2+a)^2}{4x_n^2}-a =\frac{x_n^4+2ax_n^2+a^2-4ax_n^2}{4x_n^2}\\[1mm] &=\frac{x_n^4-2ax_n^2+a^2}{4x_n^2} =\frac{(x_n^2-a)^2}{4x_n^2}. \end{aligned} \]

Residualen blir altså kvadrert og delt på \(4x_n^2\). Med \(a=2\) og \(x_n\approx\sqrt2\) er \(4x_n^2\approx8\):

\[ \frac{(6.9\cdot10^{-3})^2}{8}\approx6.0\cdot10^{-6}, \qquad \frac{(6.0\cdot10^{-6})^2}{8}\approx4.5\cdot10^{-12}, \]

akkurat tallene i tabellen. Dette er kvadratisk konvergens fra 2.3, nå for residualen.

Regelen er et spesialtilfelle av Newtons metode, som du møter i prosjekt 2.

Spindelvevet viser hvorfor det går så fort: Der kurven krysser diagonalen, er den vannrett, så trappen får nesten ingen trinn å gå.

Til slutt blir forbedringen mindre enn avstanden til neste representerbare binary64-tall. Da kan den nye verdien bli avrundet til det samme flyttallet. Skrittlengden blir da null, og flere iterasjoner gir nøyaktig samme lagrede verdi.

Svar på ukens spørsmål. Ja – men bare hvis regelen er valgt med omhu. \(g_1\) skyver oss bort fordi \(|g_1'(r)|\approx1.40>1\). \(g_2\) kommer fram, men sakte, fordi omtrent \(72\,\%\) av feilen blir igjen i hvert steg. \(g_3\) og regelen for \(\sqrt2\) er raske fordi \(g'(r)=0\). Små skritt og liten residual er nyttige kontroller, men gir ikke alene en garanti for liten løsningsfeil. For en slik garanti trenger vi også en feilgrense eller kunnskap om hvordan problemet forsterker residualen. Tallformatet begrenser dessuten nøyaktigheten vi kan oppnå.

ViktigHovedlærdom

En numerisk algoritme er en formel sammen med en startverdi, et konvergensområde, et mål på framgang, et stoppkriterium og en måte å oppdage feil på. Python lar oss undersøke valgene; analyse forklarer observasjonene.

Gå i dybden: ukens oppsummering i punkter
  • En fikspunktiterasjon trenger en startverdi og en oppdateringsregel.
  • Valget av oppdateringsregel bestemmer om og hvor raskt iterasjonen konvergerer.
  • Feil, residual og skrittlengde måler forskjellige sider av beregningen.
  • Et stoppkriterium må alltid kobles til et oppgitt mål og en toleranse.
  • Historikk og stoppårsak gjør det mulig å kontrollere resultatet, også når metoden ikke konvergerer.

I prosjektet bruker du Newtons metode i det komplekse planet og tegner tiltrekningsområdene.

Kan resultatet sendes videre?

En beregningsmodul skal finne den positive løsningen av

\[ f(x)=x^2-2=0. \]

Resultatet skal brukes videre i et kontrollsystem. Før det skjer, må en ingeniør kontrollere at programmet faktisk har funnet en løsning. Systemkravet er

\[ |f(x)|<10^{-8}. \]

Programmet har prøvd tre oppdateringsregler med samme startverdi \(x_0=1\) og de samme innstillingene atol=1e-10, rtol=1e-10 og maxiter=10:

\[ \begin{aligned} \text{A:}\quad x_{n+1}&=\frac12\left(x_n+\frac{2}{x_n}\right),\\ \text{B:}\quad x_{n+1}&=\frac{2}{x_n},\\ \text{C:}\quad x_{n+1}&=x_n-10^{-12}(x_n^2-2). \end{aligned} \]

Alle tre reglene har \(\sqrt2\) som fikspunkt, men de beveger seg svært forskjellig på vei dit. Regel C tar med vilje svært små skritt. Det gjør den nyttig for å undersøke om «lite skritt» alene er et trygt stoppkriterium.

Koden under kjører de tre reglene med fixed_point fra fanen «Implementasjon». For hver kjøring viser plottet

  • skrittlengden \(|x_n-x_{n-1}|\), som måler hvor mye den lagrede verdien endret seg,
  • residualen \(|f(x_n)|=|x_n^2-2|\), som måler hvor godt den siste verdien oppfyller likningen.

En kjøring kan bare godkjennes dersom residualen er under systemkravet. En stoppmelding er ikke nok i seg selv.

Les plottet som en kontrollrapport:

  1. Se først på stoppårsaken over hvert panel.
  2. Kontroller deretter om residualen ender under den stiplede kravlinjen.
  3. Se til slutt på utviklingen. Verdier som gjentar seg på to nivåer, tyder på en 2-syklus. Et lite skritt sammen med stor residual tyder på at programmet har stoppet før det nådde løsningen.

Bruk kodene A = 1, B = 2 og C = 3.

Stoppårsaken sier hva programmet observerte, ikke nødvendigvis at svaret er riktig. Bruk disse kodene:

  1. lite skritt
  2. maksimalt antall iterasjoner
  3. ikke-endelig verdi

Nå skal du velge hva en ingeniør bør gjøre med hver kjøring:

  1. godkjenn resultatet,
  2. bytt oppdateringsregel,
  3. krev både liten skrittlengde og liten residual før resultatet godkjennes.

For B hjelper det ikke bare å øke maxiter: regelen fortsetter å sende \(1\) til \(2\) og \(2\) tilbake til \(1\). For C er oppdateringsregelen på vei mot riktig fikspunkt, men så sakte at skrittlengden blir liten lenge før residualen er god nok.

Forbedre stoppkontrollen

Funksjonen fixed_point godkjenner nå et resultat når skrittlengden alene er liten. I en beregning med et eget residualkrav bør begge kravene være oppfylt. I Python kobler vi to krav med and: uttrykket a and b er sant bare når både a og b er sanne.

Bruk den redigerbare cellen under som en liten testbenk. Endre bare linjen med return, slik at alle tre kontrollene skriver bestått.

Når alle testene består, skiller stoppkontrollen mellom «verdien endrer seg nesten ikke» og «verdien oppfyller likningen godt nok». Det er denne forskjellen som gjør at A kan godkjennes, mens C må undersøkes videre.