Romberg-integrasjon

Romberg-integrason er en metode innen numerisk analyse til å beregne verdien av bestemte integral. Den ble utviklet i 1955 av den norsk-tyske matematiker Werner Romberg. Med utgangspunkt i den enkleste trapesintegrasjon gjøres denne nøyaktigere ved suksessive anvendelser av Richardson-ekstrapolasjon.

Matematisk formulering

De fleste integrasjoner lar seg ikke gjøre eksakt. Man må derfor ofte benytte seg av metoder som bare gir et tilnærmet riktig svar. Dette gjelder også bestemte integral. Men da disse for gitte grenser har en viss numerisk verdi, kan numerisk analyse benyttes. Når en slik integrasjonsmetode i tillegg gir et svar som kan gjøres så nøyaktig som man måtte ønske det, vil det i praksis si at integrasjonen er gjennomført.

Integralet av en funksjon y = f(x) mellom grensene a og b skrives som

$I=\int _{a}^{b}f(x)\,dx$

Den enkleste, numeriske metode for utregning av dette integralet er trapesintegrasjon. Hele intervallet b - a diskretiseres ved å dele det opp i N mindre intervall, hvert med utstrekning h = (b - a)/N. Integralet over hvert slikt delintervall er da tilnærmet lik med arealet av et lite trapes. Legges alle slike delintegrasjoner sammen, får man den tilnærmete verdien

$T={\frac {h}{2}}(f(x_{0})+2f(x_{1})+2f(x_{2})+2f(x_{3})+\dotsb +2f(x_{N-1})+f(x_{N}))$

for hele integralet I  hvor x₀ = a og x_N = b. Nøyaktigheten til svaret ved bruk av denne metoden, forventes å bli bedre desto finere diskretiseringen av integrasjonen gjøres. Det vil si å gjøre N større som tilsvarer å gjøre h mindre.

Iterasjon

Gjøres antall intervall n større på en tilfeldig måte, må man forvente å regne ut funksjonsverdiene for helt nye verdier av argumentet. Dette vil vanligvis ta tid og gjøre beregningen lite praktisk. Første del av Romberg-integrasjon består i øke dette antallet nøyaktig med en faktor 2. Siden intervallene skal fortsatt ha samme lengde, betyr det at bare halvparten av alle funksjonsverdiene må regnes ut. De andre kjennes fra den forrige beregningen.

Da antall intervaller på den måten fordobles ved hver iterasjon, vil man ha i alt ha N = 2ⁿ intervall etter n iterasjoner hvis man starter ut med det ene intervallet h₀ = b - a for n = 0. Resultatet T(n) fra trapesapproksimasjonen kan da skrives på den iterative formen som

$T(n)={\frac {1}{2}}T(n-1)+h_{n}\sum _{k=1}^{2^{n-1}}f(a+(2k-1)h_{n})$

hvor h_n = (b - a)/2ⁿ er intervallengden etter n slike iterasjoner. I utgangspunktet n = 0 har man bare ett intervall. Hele integralet er da gitt ved arealet av et trapes som tilsvarer verdien

$T(0)={\frac {1}{2}}(b-a)(f(a)+f(b))$

I de fleste tilfeller er dette en meget dårlig tilnærming. Den blir raskt bedre etter noen få iterasjoner.

Ekstrapolasjon

Ved bruk av Euler-Maclaurin-formelen kan man vise generelt at nøyaktigheten E(n) = I - T(n)  av integralet etter n iterasjoner kan skrives som ekspansjonen

$E(n)=E_{1}(n)h^{2}+E_{2}(n)h^{4}+\cdots +E_{k}(n)h^{2k}+\cdots$

hvor koeffisientene E_k(n) er gitt ved differenser mellom høyere deriverte av funksjonen f(x) i integralets øvre og nedre grenser. Andre del av Romberg-integrasjon består i å benytte Richardson-ekstrapolasjon til å øke nøyaktigheten av hvert iterativt resultat T(n) = R(n,0) som er det laveste Romberg-resultat. Da nøyaktigheten av dette er av orden h², vil man etter en slik ekstrapolasjon få resultatet

${\begin{aligned}R(n,1)&={\frac {1}{2^{2}-1}}(2^{2}R(n,0)-R(n-1,0))={4 \over 3}T(n)-{1 \over 3}T(n-1)\\&={h_{n} \over 3}{\Big [}f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\cdots +4f(x_{2^{n}-1})+f(x_{2^{n}}){\Big ]}\end{aligned}}$

som nå vil være nøyaktig til orden h⁴. Dette er i overensstemmelse med hva man finner ved Simpson-integrasjon hvor man approksimerer funksjon y = f(x) med et diskret set med parabelbuer.

Da unøyaktigheten ved integrasjonen bare inneholder ledd med like potenser av h, kan man øke den til orden h⁶ ved å foreta en ny Richardson-ekstrapolasjon. Slik kan man fortsette til nøyaktigheten er blitt av orden h^2m  hvor 1≤ m ≤ n. I hvert stepp er derfor delresultatene fra Romberg-integrasjonen forbundet ved den rekursive formelen

$R(n,m)={\tfrac {1}{2^{2m}-1}}{\Big [}2^{2m}R(n,m-1)-R(n-1,m-1){\Big ]}$

med verdiene R(n,0) = T(n)  som utgangspunkt. Slik kan nøyaktigheten systematisk økes og er bare begrenset av avrundingsfeil i regneprosessen.

Delresultatene som inngår i denne integrasjonen, kan fremstilles på en oversiktlig måte i en såkalt «Romberg-tabell». Hver horisontal linje viser ser de midlertidige verdiene som fremkommer ved Richardson-ekstrapolasjon for en oppdeling av hele integrasjonsintervallet i N = 2ⁿ mindre intervaller. Da må man også benytte et delresulatet fra linjen over. På den måten griper denne Romberg-integrasjonen helt tilbake til verdien R(0,0). De tre første linjene i tabellen inneholder derfor elementene

${\begin{array}{l | c | c | r}n&O(h^{2})&O(h^{4})&O(h^{6})\\\hline 0&R(0,0)&&\\1&R(1,0)&R(1,1)&\\2&R(2,0)&R(2,1)&R(2,2)\\\end{array}}$

som hvert avhenger av nærmeste elementet til venstre for seg samt elementet over dette.

Eksempel

Som et enkelt eksempel kan man integrere funksjonen y = x⁵ mellom grensene a = 0 og b = 1. Det eksakte resultatet er da trivielt E = 1/6 = 0.1666667. Laveste approksimasjon er da R(0,0) = (1 - 0)/2 = 0.5, mens to trapeser gir R(1,0) = 17/64 = 0.265625. Ved en Richardson-ekstrapolasjon finnes herav resultatet R(1,1) = 3/16 = 0.187500 som er Simpson-resultatet. Ved enda en fordobling av delintervaller til N = 4 har man dermed de tre første linjene i tabellen:

${\begin{array}{l | c | c | r}n&O(h^{2})&O(h^{4})&O(h^{6})\\\hline 0&0.500000&&\\1&0.265625&0.187500&\\2&0.192383&0.167969&0.166667\\\end{array}}$

Bare etter to iterasjoner og to tilsvarende ekstrapolasjoner er resultatet fra Romberg-integrasjonen nøyaktig med fem riktige desimaler.