3.4 Numerisk integrasjon

Integralet som ikke lot seg løse

Med substitusjon, delvis integrasjon og delbrøkoppspalting i verktøykassen føles det kanskje som om alle integraler nå kan knekkes. Så møter du dette:

$$\int_0^2 e^{-x^2}\,dx$$

Du prøver substitusjon — ingen kjerne-derivert. Delvis integrasjon — det blir bare verre. Sannheten er nedslående og fascinerende på én gang: $e^{-x^2}$ har ingen antiderivert som kan skrives med elementære funksjoner. Det er bevist. Og dette er ikke et obskurt unntak — det er normalsangen for funksjonene i statistikk ($e^{-x^2}$ ligger bak normalfordelingen), signalbehandling ($\displaystyle \frac{\sin x}{x}$) og geometri ($\sqrt{1+x^3}$). I tillegg kommer alle situasjoner der funksjonen bare finnes som måledata, uten formel i det hele tatt.

Men integralet eksisterer jo — arealet under kurven er der, like virkelig som alltid. Vi kan bare ikke nå det med antideriverte. Løsningen er å vende tilbake til Riemanns rektangler fra kapittel 2.4, men nå med et praktisk mål: ikke en grenseverdi på papiret, men et tall med kontrollert nøyaktighet. Det er numerisk integrasjon — og kappløpet står mellom rektangler, trapeser og parabler.

Rektangelmetoden: enkel, men grov

Utgangspunktet kjenner du: del $[a,b]$ i $n$ delintervaller med bredde $\displaystyle h = \frac{b-a}{n}$, og tilnærm arealet med rektangler. Bruker vi venstre endepunkt i hvert delintervall, får vi venstre rektangelmetode:

$$\int_a^b f(x)\,dx \approx h\left[f(x_0) + f(x_1) + \cdots + f(x_{n-1})\right]$$

Høyre endepunkt gir høyremetoden, og midtpunktet gir midtpunktmetoden — samme idé, smartere valg.

Testbenken vår blir $\displaystyle \int_0^2 x^2\,dx = \frac{8}{3} \approx 2{,}667$, der fasiten er kjent. Venstre metode med $n = 4$: $h = 0{,}5$, og

$$0{,}5\,[f(0) + f(0{,}5) + f(1) + f(1{,}5)] = 0{,}5 \cdot 3{,}5 = 1{,}75$$

En feil på hele $34\,\%$ — rektanglene henger etter den voksende grafen. Midtpunktmetoden med samme $n$ bruker $x = 0{,}25, 0{,}75, 1{,}25, 1{,}75$ og gir

$$0{,}5\,[0{,}0625 + 0{,}5625 + 1{,}5625 + 3{,}0625] = 2{,}625$$

— bare $1{,}6\,\%$ feil! Hemmeligheten er at midtpunktet balanserer: rektangelet stikker over grafen på den ene halvdelen og under på den andre, og feilene kansellerer nesten. Samme antall regneoperasjoner, dramatisk bedre svar. Lærdommen som preger hele kapittelet: hvor du måler, betyr mer enn hvor mye du regner.

Trapesmetoden: skrå tak

Rektanglenes svakhet er de flate takene. Hvorfor ikke la taket følge kurven — i det minste som en rett linje mellom nabopunktene? Da blir hver stripe et trapes med areal $\displaystyle \frac{h}{2}[f(x_i) + f(x_{i+1})]$, og summen over alle stripene gir trapesmetoden:

$$T_n = \frac{h}{2}\left[f(x_0) + 2f(x_1) + 2f(x_2) + \cdots + 2f(x_{n-1}) + f(x_n)\right]$$

Endepunktene teller én gang, de indre punktene to (de deles av to nabotrapesar).

På testbenken med $n = 4$: $T_4 = 0{,}25\,[0 + 2(0{,}25) + 2(1) + 2(2{,}25) + 4] = 0{,}25 \cdot 11 = 2{,}75$ — feil $3{,}1\,\%$. Bedre enn venstremetoden, men faktisk litt bak midtpunktmetoden her; for konvekse kurver skjærer trapesets tak konsekvent over grafen.

Hvor fort blir det bedre når $n$ øker? Feilestimatet sier $\displaystyle |E_T| \leq \frac{(b-a)^3}{12n^2}M_2$, der $M_2$ er maksimum av $|f''|$ — og det vesentlige er $n^2$ i nevneren: dobler du $n$, deles feilen på fire. For $\int_0^1 e^x\,dx$ ser man det svart på hvitt: feilen går fra $0{,}0356$ ved $n=2$ til $0{,}0089$ ved $n=4$ til $0{,}0022$ ved $n=8$ — faktor fire hver gang. Midtpunktmetoden deler samme orden, $O(1/n^2)$, med konstanten $24$ i stedet for $12$.

Trapesmetodens adelsmerke er måledata: har du bare tabellverdier fra et eksperiment — hastigheter, strømninger, sensoravlesninger — er trapesmetoden den naturlige måten å integrere dem på, uten noen formel for $f$.

Simpsons metode: parablene vinner

Rette linjer var bedre enn flate tak — hva med å la taket krumme? Simpsons metode legger en parabel gjennom tre og tre nabopunkter, og parabler følger kurvers krumning langt bedre enn korder. Formelen (som krever at $n$ er et partall):

$$S_n = \frac{h}{3}\left[f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + \cdots + 4f(x_{n-1}) + f(x_n)\right]$$

med det karakteristiske koeffisientmønsteret $1, 4, 2, 4, 2, \ldots, 4, 1$.

Testbenken, $n = 4$: $\displaystyle S_4 = \frac{0{,}5}{3}[0 + 4(0{,}25) + 2(1) + 4(2{,}25) + 4] = \frac{0{,}5}{3} \cdot 16 = \frac{8}{3}$ — eksakt! Ingen tilfeldighet: Simpsons metode er eksakt for alle polynomer til og med grad 3, ett hakk bedre enn parabel-opphavet skulle tilsi (en heldig symmetri visker ut tredjegradsfeilen).

For andre funksjoner er feilordenen $O(1/n^4)$: dobling av $n$ deler feilen på seksten. Sammenlikningen er nådeløs — for samme regnearbeid vinner Simpson nesten alltid.

Tilbake til den umulige utfordringen: $\int_0^1 \sqrt{1+x^3}\,dx$, uten elementær antiderivert. Med $n = 4$ gir midtpunkt $1{,}106$, trapes $1{,}115$ og Simpson $1{,}109$ — mot den «sanne» verdien $1{,}1114$ er Simpson nærmest. Og feilestimatene, som $\displaystyle |E_S| \leq \frac{(b-a)^5}{180n^4}M_4$, gjør tilnærmingen til vitenskap: vi kan garantere nøyaktigheten på forhånd, eller regne ut hvilken $n$ som trengs for ønsket presisjon. Tommelfingerreglene til slutt: trapes for tabelldata, Simpson for glatte funksjoner med valgfrie punkter, midtpunkt for raske overslag.

Oppsummering: kappløpet om arealet

Når antideriverten ikke finnes — som for $e^{-x^2}$, $\displaystyle \frac{\sin x}{x}$ og $\sqrt{1+x^3}$ — eller funksjonen bare finnes som måledata, overtar numerikken. Alle metodene deler grunnplanen fra Riemann: del $[a,b]$ i $n$ biter med bredde $\displaystyle h = \frac{b-a}{n}$, og summer bidrag. Forskjellen ligger i taket over hver stripe. Rektangelmetoden (flatt tak) er enklest men grovest, med venstre/høyre-variantene som systematisk bommer på monotone funksjoner; midtpunktvarianten redder mye ved å la feilene kansellere. Trapesmetoden (skrått tak, mønster $1,2,2,\ldots,2,1$) deler feilorden $O(1/n^2)$ med midtpunktmetoden — dobling av $n$ firedeler feilen — og er førstevalget for tabelldata. Simpsons metode (parabeltak, mønster $1,4,2,4,\ldots,4,1$, partall $n$) er eksakt for polynomer opp til grad $3$ og konvergerer som $O(1/n^4)$: dobling av $n$ deler feilen på seksten.

Feilestimatene — $\displaystyle \frac{(b-a)^3}{12n^2}M_2$ for trapes, $\displaystyle \frac{(b-a)^5}{180n^4}M_4$ for Simpson — løfter det hele fra gjetning til garanti: nøyaktigheten kan planlegges før regningen starter. Men hvem orker å summere hundre funksjonsverdier for hånd? Ingen — og ingen trenger det. I neste kapittel skriver vi metodene som Python-kode og lar maskinen løpe kappløpet for oss.