3.5 Programmering av integrasjon

Maskinen overtar regnearbeidet

I forrige kapittel regnet du trapessummer for hånd og kjente på hvor fort det blir slitsomt — fire delintervaller er greit, fire hundre er utenkelig. Men metodene er jo bare oppskrifter: del opp, regn funksjonsverdier, summer med riktige vekter. Og oppskrifter er nøyaktig det datamaskiner er laget for å følge.

I dette kapittelet oversetter vi rektangel-, trapes- og Simpsons metode til Python — hver av dem blir en funksjon på under ti linjer — og slipper dem løs med tusenvis av delintervaller. Vi måler feilene og ser konvergensordenene fra teorien dukke opp i utskriftene: faktor 4 for trapes, faktor 16 for Simpson. Deretter åpner vi proffverktøyene i NumPy og SciPy, der én linje kode integrerer med maskinpresisjon, og avslutter med en metode som er like overraskende som den er nyttig: Monte Carlo-integrasjon, der tilfeldige tall beregner arealer — og estimerer $\pi$.

Metodene blir kode

La oss starte der forrige kapittel sluttet: venstre rektangelmetode. Som Python-funksjon:

def rektangelmetoden(f, a, b, n):
    h = (b - a) / n
    total = 0
    for i in range(n):
        total += f(a + i * h)
    return h * total

Mønsteret er metoden selv: steglengde $\displaystyle h = \frac{b-a}{n}$, en løkke som summerer funksjonsverdier, og multiplikasjon med $h$ til slutt. Med f = lambda x: x**2 og $n = 1000$ gir den $0{,}332834$ for $\displaystyle \int_0^1 x^2\,dx = \frac{1}{3}$ — feil $0{,}0005$, omtrent $0{,}15\,\%$.

Trapesmetoden er nesten lik; bare endepunktene behandles spesielt (vekt $\displaystyle \frac{1}{2}$):

def trapesmetoden(f, a, b, n):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    for i in range(1, n):
        total += f(a + i * h)
    return h * total

Og Simpson krever partall $n$ og to løkker — én for vekt $4$ (odde indekser), én for vekt $2$ (indre like indekser):

def simpson(f, a, b, n):
    if n % 2 != 0:
        raise ValueError("n må være et partall")
    h = (b - a) / n
    total = f(a) + f(b)
    for i in range(1, n, 2):
        total += 4 * f(a + i * h)
    for i in range(2, n, 2):
        total += 2 * f(a + i * h)
    return (h / 3) * total

Styrkeprøven: $\int_0^{\pi} \sin x\,dx = 2$. Med bare $n = 10$ gir trapesmetoden feil $1{,}65 \cdot 10^{-2}$, mens Simpson gir $6{,}78 \cdot 10^{-6}$ — over to tusen ganger mer nøyaktig for samme arbeid. Ti linjer kode gjenskaper hele forrige kapittels teori.

Konvergens du kan se — og proffverktøyene

Teorien lovet at trapesfeilen går som $O(1/n^2)$. Med en liten testløkke kan vi se det. For $\int_0^1 e^x\,dx = e - 1$ skriver vi ut feilen for $n = 10, 20, 40, 80, 160, 320$ og forholdet mellom påfølgende feil:

n      Feil        Faktor
10     3.87e-05    -
20     9.68e-06    4.00
40     2.42e-06    4.00
80     6.05e-07    4.00

Faktoren er nøyaktig $4$ hver gang $n$ dobles — konvergensorden $O(1/n^2)$, eksperimentelt bekreftet. For Simpson ville tilsvarende tabell vist faktor $16$. Slik feilanalyse er god praksis hver gang du bruker numerikk: kjør med økende $n$, og se at svaret stabiliserer seg.

I praksis slipper du ofte å skrive metodene selv. SciPy har ferdige, optimaliserte rutiner:

from scipy import integrate
resultat, feil = integrate.quad(f, 0, np.inf)

quad velger punkter adaptivt, håndterer til og med uendelige grenser, og returnerer både svaret og et feilestimat — for $\int_0^{\infty} e^{-x^2}\,dx$ treffer den fasiten $\displaystyle \frac{\sqrt{\pi}}{2}$ med feil under $10^{-8}$. For måledata finnes numpy.trapz: har du tidspunkter t = [0, 1, 2, 3, 4, 5] og målte hastigheter v = [0, 5, 8, 10, 9, 6], gir np.trapz(v, t) tilbakelagt strekning $34{,}5$ meter direkte fra tabellen — trapesmetoden på ekte data, én linje.

Arbeidsdelingen er typisk for beregningsorientert matematikk: skriv metodene selv én gang for å forstå dem, bruk biblioteket etterpå for å anvende dem.

Monte Carlo: tilfeldigheten som regner

Til slutt en metode fra en helt annen planet. Monte Carlo-integrasjon — oppkalt etter kasinobyen — bruker tilfeldige tall: trekk $N$ tilfeldige $x$-verdier uniformt fra $[a,b]$, regn gjennomsnittet av funksjonsverdiene, og gang med intervallbredden:

$$\int_a^b f(x)\,dx \approx (b-a) \cdot \frac{1}{N}\sum_{i=1}^{N} f(x_i)$$

Logikken: integralet er gjennomsnittsverdien av $f$ ganget med bredden, og et tilfeldig utvalg estimerer gjennomsnittet. I kode:

def monte_carlo(f, a, b, N):
    x = np.random.uniform(a, b, N)
    return (b - a) * np.mean(f(x))

Morsomste anvendelse: estimere $\pi$. Halvsirkelen $f(x) = \sqrt{1-x^2}$ over $[-1,1]$ har areal $\displaystyle \frac{\pi}{2}$, så $2 \cdot$ Monte Carlo-estimatet $\approx \pi$. Med $N = 100$ får vi kanskje $3{,}02$; med $N = 10^6$ rundt $3{,}1410$. (Beslektet er «dartskive-varianten»: kast tilfeldige punkter i kvadratet $[-1,1] \times [-1,1]$ og tell andelen som lander i enhetssirkelen — andelen ganger $4$ estimerer $\pi$.)

Konvergensen er ærlig talt treg: feilen går som $O(1/\sqrt{N})$, så fire ganger flere punkter gir bare halvert feil — langt bak Simpsons $O(1/n^4)$. Hvorfor bry seg, da? Fordi Monte Carlo har to superkrefter: feilordenen er uavhengig av dimensjonen (gull verdt for integraler i mange variabler, der rutenettmetoder kollapser), og metoden er robust og latterlig enkel å implementere. SciPy byr for øvrig også på dblquad for dobbeltintegraler — og verktøykassen vokser med behovet: adaptive metoder som forfiner seg selv der funksjonen varierer mest, Romberg-ekstrapolasjon som pumper opp trapesresultater, og Gauss-kvadratur med optimalt plasserte punkter. Prinsippet bak dem alle har du nå: summer smarte stikkprøver av funksjonen.

Oppsummering: ti linjer kode, tusen rektangler

Numerisk integrasjon gikk i dette kapittelet fra håndarbeid til maskinkraft. De tre klassiske metodene ble Python-funksjoner på under ti linjer hver: rektangelmetoden som en ren sum-løkke, trapesmetoden med halve vekter på endepunktene, og Simpson med vekslende vekter $4$ og $2$ (og partallskrav på $n$). Eksperimentene bekreftet teorien i sanntid: trapesfeilen falt med faktor $4$ per dobling av $n$ ($O(1/n^2)$), og Simpson slo trapes med en faktor over $2000$ på $\int_0^{\pi}\sin x\,dx$ allerede ved $n = 10$.

Proffverktøyene tar over derfra: scipy.integrate.quad integrerer adaptivt med feilestimat på kjøpet — også over uendelige intervaller — og np.trapz integrerer rene måledata. Og Monte Carlo viste numerikkens ville side: tilfeldige punkter som estimerer integraler (og $\pi$), med treg, men dimensjonsuavhengig konvergens $O(1/\sqrt{N})$ som gjør metoden uunnværlig i mange dimensjoner.

Arbeidsregelen du tar med videre: forstå metoden ved å skrive den selv, bruk biblioteket i produksjon, og verifiser alltid mot kjente integraler eller økende $n$. Med både analytiske og numeriske verktøy i kassen er du klar for integralregningens mest håndgripelige anvendelse: å regne ut volumet av ting — fra vinglass til planeter.