Numerisk integrasjon#

Læringsutbytte

Etter å ha arbeidet med dette temaet, skal du kunne:

  1. Forklare forskjellen på ulike tilnærminger til rektangelmetoden (venstre-, høyre- og midtpunktstilnærming).

  2. Forklare og utlede trapesmetoden.

  3. Implementere rektangelmetoden og trapesmetoden.

  4. Integrere funksjoner numerisk.

Integrasjon#

Du kjenner kanskje integrasjon som en metode å regne ut arealet under en graf eller volumet til legemer på. I tillegg kjenner du antakelig til integrasjon som den motsatte operasjonen av derivasjon. At derivasjon og integrasjon er motsatte operasjoner, er bevist gjennom analysens fundamentalteorem. Termen analyse brukes her om den greinen av matematikk som omhandler derivasjon og integrasjon (kalkulus). Siden datamaskinen kun kan operere med bestemte verdier, kan vi med numeriske algoritmer kun tilnærme det bestemte integralet. Vi antideriverer derfor ikke her.

Integraler i kjemi#

Det finnes mange anvendelser av integrasjon i kjemi. Integrasjon brukes til alt fra å beregne molekylgeometri i kvantekjemi til å beregne arealet til en topp i et NMR-spekter eller et kromatogram. Dessuten integrerer vi når vi løser differensiallikninger. I for eksempel NMR-spektre sier integralet under signalkurvene hva den relative forekomsten av H-atomer er i det bestemte miljøet. Numerisk integrasjon kan finne disse arealene for oss.

Her skal vi se på ulike metoder for å tilnærme det bestemte integralet.

Rektangelmetoden#

Vi kan enklest gjøre en tilnærming til det bestemte integralet ved å utnytte at det kan skrives som en grenseverdi av Riemann-summer. En Riemann-sum kan beskrives som en tilnærming til arealet under en graf ved hjelp av arealet til geometriske figurer. En vanlig tilnærming er å bruke rektangler:

Her benyttes \(N = 10\) rektangler for å tilnærme integralet av \(f(x) = \cos{(x)} +2\) for \(x\in[2,12]\). Bredden av rektanglene må være \((b-a)/N = (12-2)/10 = 1\). Vi ser også at høyden av hvert rektangel er \(f(x_n)\) der \(n\in[2,11]\), det vil si at vi lar venstresiden av rektangelet gå opp til grafen. Dersom vi regner ut arealet til hver av disse rektanglene, får vi 18.046675645664006, noe som ligger et lite stykke unna den analytiske verdien (\((\sin{(12)} + 2\cdot 12)-(\sin{(2)} + 2\cdot 2) \approx 18.554129655173885\). Dersom vi øker antall rektangler, her til 50, får vi naturlig nok et bedre estimat (her har vi inkludert resultatene i figuren):

Det er åpenbart, spesielt i den første figuren, at flere områder av rektanglene ligger utenfor grafen. De ligger både for langt over og under grafen flere steder. Men feilen blir ikke så stor som det kan se ut som fordi vi nettopp har områder både over og under grafen. Relativ feil er her ca. 2.7 % med 10 rektangler og 0.65 % med 50 rektangler.

La oss illustrere feilen vi får med en lineær funksjon, \(f(x) = 10x\):

Vi har i disse modellene målt rektangelhøyden på venstre ytterkant av rektangelet. Dette gir oss her en underestimering av integralet for hvert rektangel. Vi kan også måle høyden av rektanglene på høyre ytterkant. Da får vi en tilsvarende overestimering:

Vi ser her at vi får en stor underestimering med venstretilnærmingen og en stor overestimering med høyretilnærmingen, med en relativ feil på ca. 7.1 %. Det samme skjer for alle funksjoner i intervaller som er enten voksende eller synkende i hele intervallet. En måte å kompensere for dette, er å benytte høyden til rektangelet i midten istedenfor i endepunktet til venstre eller høyre:

Det aller beste resultatet får vi altså med denne midtpunktstilnærmingen. For lineære funksjoner vil vi få en eksakt verdi (selv med ett rektangel!) fordi rektanglene er like store over og under grafen. Men midtpunktstilnærmingen er generelt bedre også på for eksempel polynomer av høyere grad og trigonometriske funksjoner. La oss nå se på implementering av metodene.

Implementering av rektangelmetodene#

Rektangelmetoden (venstretilnærming)

Det bestemte integralet til en funksjon \(f(x)\) fra \(x = a\) til \(x = b\) kan tilnærmes ved arealet til \(n\) rektangler med bredden \(h = \frac{b-a}{n}\):

\[\int_a^b f(x) \ \mathrm{d}x \approx h \sum_{k=1}^{n} f(x_k)\]

Underveisoppgave

Programmet nedenfor gir en funksjon som bruker venstretilnærmingen av rektangelmetoden til å beregne det bestemte integralet av en funksjon f mellom a og b. Fyll inn det som mangler i metoden.

Løsningsforslag

def f(x):				#Definerer en funksjon som vi skal integrere.
    return x**3
    
def f_analytisk(x):     #Definerer analytisk verdi for sammenlikning.
    return (1/4)*x**4

def rektangelmetoden(f, a, b, n):
    A = 0.0							   
    h = (b-a)/n			#Bredden til rektanglene
    for k in range (n):	
        A = A + f(a + k*h)*h
    return A
    
print("Numerisk verdi:", rektangelmetoden(f, 0, 5, 1000))
print("Analytisk verdi:", f_analytisk(5)-f_analytisk(1))

Tilsvarende kan vi beskrive høyretilnærmingen slik:

Rektangelmetoden (høyretilnærming)

Det bestemte integralet til en funksjon \(f(x)\) fra \(x = a\) til \(x = b\) kan tilnærmes ved arealet til \(n\) rektangler med bredden \(h = \frac{b-a}{n}\):

\[\int_a^b f(x) \ \mathrm{d}x \approx h \sum_{k=1}^{n} f(x_{k+1})\]

Underveisoppgave

Implementer algoritmen for høyretilnærmingen som en Python-funksjon. Test og sammenlikn med venstretilnærmingen på integralet \(\int_2^8 f(x) = x^2 - 2x + 4 \ dx\).

Løsningsforslag

def rektangelmetoden_høyre(f, a, b, n):
    A = 0.0							   
    h = (b-a)/n
    for k in range(n):	
        A = A + f(a + (k+1)*h)*h
    return A

Den beste tilnærmingen med rektangler får vi med midtpunktstilnærmingen:

Rektangelmetoden (midtpunktstilnærming)

Det bestemte integralet til en funksjon \(f(x)\) fra \(x = a\) til \(x = b\) kan tilnærmes ved arealet til \(n\) rektangler med bredden \(h = \frac{b-a}{n}\):

\[\int_a^b f(x) \ \mathrm{d}x \approx h \sum_{k=1}^{n} f\left(\frac{1}{2}(x_k + x_{k+1})\right)\]

Underveisoppgave

Implementer algoritmen for høyretilnærmingen som en Python-funksjon. Test og sammenlikn med venstretilnærmingen på integralet \(\int_2^8 f(x) = x^2 - 2x + 4 \ dx\).

Løsningsforslag

def rektangelmetoden_midt(f, a, b, n):
    A = 0.0							   
    h = (b-a)/n
    for k in range (n):	
        A = A + f(a + (1+2*k)*(h/2))*h
    return A

Dersom vi har funksjoner med stor stor stigning eller minking (\(|f'(x) >> 0|\)), trenger vi mange rektangler for å få et godt resultat. Dette kan gi langsomme programmer dersom integrasjonen må gjentas flere ganger. Vi skal derfor nå se på noen forbedringer av rektangelmetoden.

Trapesmetoden#

Vi kan legge merke til at toppstykket i et rektangel er ei rett, horisontal linje. Ei slik linje kan representeres som et polynom av nullte grad, \(f(x) = ax^0 = a\), der \(a\) er et reellt tall. La oss se på mulighetene for å bytte ut dette toppstykket med et polynom av første grad, \(f(x) = ax^1 = ax\). Da får vi trapeser istedenfor rektangler. En algoritme for dette er litt mindre intuitiv og litt mer jobb å utlede, men vi spanderer på oss det. La oss ta utgangspunktet i trapesmetoden illustrert med ett trapes i intervallet \([a, b] = [2, 12]\) på \(f(x) = cos(x) + 2\):

La oss utlede en algoritme for trapesmetoden. Vi tar utgangspunkt i at arealet av et trapes er gitt ved følgende formel:

\[A_{trapes} = \frac{side 1+side 2}{2}\cdot h\]

De to sidene \(x_1\) og \(x_2\) er gitt ved henholdsvis \(f(a)\) og \(f(b)\). Høyden i trapeset blir stykket langs hele x-aksen, altså \(b-a\). Arealet blir derfor:

\[A = \frac{f(a) + f(b)}{2}\cdot(b-a)\]

La oss nå utvide til \(n\) trapeser. Da blir høyden av hvert trapes \(h = (x_1+x_2)/n\), noe som gir dette arealet for hvert \(i\)-te trapes:

\[A_i = \frac{f(x_i) + f(x_i)}{2}\cdot h\]

Summen blir da slik for \(n\) trapeser:

\[\frac{f(a)+f(a+h)}{2}h + \frac{f(a+h)+f(a+2h)}{2}h + \frac{f(a+2h)+f(a+3h)}{2}h + ... + \frac{f(a+ih)+f(b)}{2}h\]

Vi multipliserer så alle ledd med \(h\) og dividerer dem på 2. Dette setter vi utenfor uttrykket:

\[\frac{h}{2}\cdot (f(a)+f(a+h)+f(a+h)+f(a+2h)+f(a+2h)+ ... + f(a+ih) + f(b))\]

Trekker vi sammen like ledd, får vi:

\[\frac{h}{2}\cdot (f(a)+2f(a+h)+2f(a+2h)+2f(a+3h)+ ... + 2f(a+ih) + f(b))\]

Siden det bare er \(f(a)\) og \(f(b)\) som ikke er multiplisert med 2, kan vi forenkle:

\[h\left(\frac{f(a)+f(b)}{2} + (f(a+h)+f(a+2h)+f(a+3h)+ ... + f(a+ih))\right)\]

Den siste samlingen av ledd kan vi skrive som en sum. Da får vi trapesmetoden:

Trapesmetoden

Det bestemte integralet til en funksjon \(f(x)\) fra \(x = a\) til \(x = b\) kan tilnærmes ved arealet til \(n\) trapeser med bredden \(h = \frac{b-a}{n}\):

\[\int_a^b f(x) dx \approx h\left( \frac{f(a)+f(b)}{2} + \sum_{i=1}^{n-1} f(x_i) \right)\]

Underveisoppgave

Implementer trapesmetoden og test den på funksjonene \(f(x) = x^3 + 2x\) og \(f(x) = \sqrt(x)\) i intervallet \([2,4]\). Sammenlikn med resultatene du får fra rektangelmetoden (gjerne flere av dem!) med samme antall geometriske figurer.

Løsningsforslag

def trapesmetoden(f, a, b, n):
    h = (b-a)/n
    A = (f(a) + f(b))/2.0
    for k in range(1, n):
        A = A + f(a + k*h)
    return A*h

Dersom vi øker antallet trapeser, blir naturlig nok også tilnærmingen bedre, som vi kan se av denne figuren:

Simpsons metode#

Vi har nå sett på to tilnærminger til integralet som bruker henholdsvis nullte- og førstegradspolynomer som toppstykke på de geometriske figurene som vi beregner arealet av. For å få en enda bedre tilnærming, spesielt til oscillerende funksjoner, kan vi bruke et toppstykke av et polynom med høyere grad enn 1. Alle disse metodene kan utledes ved hjelp av interpolasjon, men vi skal ikke gjøre det her. Her nøyer vi oss med å vise og implementere algoritmen for tredjegradstilnærmingen. Denne algoritmen kalles Simpsons metode:

Simpsons metode

Det bestemte integralet til en funksjon \(f(x)\) fra \(x = a\) til \(x = b\) kan tilnærmes med \(n\) geometriske figurer, der \(n\) er et partall:

\[\int_a^b f(x) dx \approx \frac{h}{3} \left( f(a) + f(b) + 2\sum_{k = 1}^{\frac{n}{2} - 1} f(x_{2k}) + 4\sum_{k=1}^{\frac{n}{2}} f(x_{2k-1}) \right)\]

Metoden kan implementeres slik:

def simpsons(f, a, b, n):
    h = (b-a)/n
    total = f(a) + f(b)
    for k in range(1,n,2):
        A += 4 * f(a + k*h)
    for k in range(2,n,2):
        A += 2*f(a + k*h)
    return A*h/3

Oppgave

Studer koden ovenfor og prøv å finne igjen leddene i definisjonen av Simpsons metode.

Metodene vi har sett på, bygger på samme prinsipp, nemlig tilnærming av arealet under grafen ved hjelp av geometriske figurer med rektangelbase og et polynom av grad \(n\) som toppstykke. Siden prinsippet er det samme, kaller vi dem en familie av metoder (hyggelig, ikke sant?). Denne familien heter Newton-Cotes. Det vil si at for eksempel trapesregelen kalles en Newton-Cotes-metode av første grad, mens Simpsons metode er en Newton-Cotes-metode av andre grad. Det finnes mange andre metoder og familier innenfor numerisk integrasjon, men disse lar vi ligge foreløpig.

Bruk av biblioteker for å integrere#

Vi vender tilbake til scipy-biblioteket, som også inneholder metoder for numerisk integrasjon. Vi kan bruke biblioteket slik:

from scipy import integrate
import numpy as np

def f(x):
    return x**3 - 1

n = 1000
x = np.linspace(0,5,n)
y = f(x)

# Integrasjon
trapes = integrate.trapz(y,x) # Trenger arrayer som parameter
simpsons = integrate.simps(y,x) # Trenger arrayer som parameter
gauss_kvadratur = integrate.quad(f,0,5) # Trenger funksjon som parameter
print("Trapesmetoden:",trapes)
print("Simpsons metode:",simpsons)
print("Gauss kvadratur:",gauss_kvadratur) #Skriver ut svar og absolutt feil med Gauss kvadratur, en metode vi ikke har sett på

Trapesmetoden: 151.2501565629694
Simpsons metode: 151.25000015671972
Gauss kvadratur: (151.25, 1.6959623319719519e-12)

Underveisoppgave

Studer koden ovenfor. Hvordan fungerer de ulike funksjonene? Test gjerne ut funksjonene selv.

Multippel integrasjon#

Multippel integrasjon kan også gjennomføres med scipy-biblioteket. La oss prøve å løse følgende integraler med funksjonene dblquad og tplquad:

\[\int_0^{\frac{\pi}{2}} \int_{-1}^1 x \sin(y) - y e^x \ dx dy\]
\[\int_0^{3} \int_{0}^2 \int_{0}^1 xyz \ dx dy dz\]
from scipy import integrate
import numpy as np

def f(y,x):
    return x*np.sin(y) - y*np.exp(x)

def g(z, y, x):
    return x*y*z

numerisk_dobbel = integrate.dblquad(f, -1, 1, 0, np.pi/2)
numerisk_trippel = integrate.tplquad(g, 0, 1, 0, 2, 0, 3)

analytisk_dobbel = np.pi**2/8 * (1/np.exp(1)-np.exp(1))
analytisk_trippel = 9/2
print(f'Numerisk verdi av dobbeltintegralet: {numerisk_dobbel[0]}')
print(f'Analytisk verdi av dobbeltintegralet: {analytisk_dobbel}')

print(f'Numerisk verdi av trippelintegralet: {numerisk_trippel[0]}')
print(f'Analytisk verdi av trippelintegralet: {analytisk_trippel}')

Numerisk verdi av dobbeltintegralet: -2.899692718238082
Analytisk verdi av dobbeltintegralet: -2.8996927182380823
Numerisk verdi av trippelintegralet: 4.5
Analytisk verdi av trippelintegralet: 4.5

Oppgaver#

Oppgave 1

Løs puslespillet nedenfor. Programmet skal definere rektangelmetoden som en Python-funksjon.

Oppgave 2

Implementer algoritmen for rektangelmetoden som en Python-funksjon. Test metoden på integralet

\[\int_2^8 f(x) = x^2 - 2x + 4 \ dx\]

Oppgave 3

Implementer trapesmetoden og ulike tilnærmnger for rektangelmetoden som Python-funksjoner. Gjør en feilanalyse av metodene og sammenlikn svarene du får på integralet

\[\int_2^8 f(x) = x^2 - 2x + 4 \ dx\]

Oppgave 4

Utled rektangelmetoden ved hjelp av definisjonen av integralet og en figur.

Oppgave 5

Utled trapesformelen med utgangspunkt i definisjon av det bestemte integralet og formelen for arealet av et trapes. Bruk gjerne en figur.

Oppgave 6

En sannsynlighetsfordeling er en matematisk beskrivelse av hvor sannsynlig noe vil skje. Det finnes mange forskjellige sannsynlighetsfordelinger. I denne oppgava skal vi se på en fordeling \(f(x)\) som kalles en standard normalfordeling. Den er definert ved:

\[f(x) = \frac{1}{\sqrt{2\pi}}\cdot e^{-\frac{x^2}{2}}\]

Skriv et program som bruker trapes- eller rektangelmetoden for å integrere \(f(x)\) i området mellom \(a = -k\cdot\sigma = -k\) og \(b = k\cdot\sigma = k\) for \(k \in {1, 2, 3}\) og skriver ut resultatene. Du kan sette \(\sigma = 1\).

For å forsikre deg om at du bruker en grei verdi for \(n\), kan du se om resultatene blir omtrent lik 0.68, 0.95 og 0.997 for henholdsvis \(k \in {1, 2, 3}\).

Oppgave 7*

Noen funksjoner oppfører seg trøblete ved at de går mot minus eller pluss uendelig, eller ved at funksjonen varierer svært mye mot en bestemt verdi. La oss her ta denne funksjonen:

\[f(x) = \sin\left(\frac{1}{x}\right)\]

Funksjonen oppfører seg turbulent like ved 0. Plott funksjonen med (D_f = [-5, 0) \cup (0,5]) og studer grafen. Prøv å derivere og integrere med ulike metoder i punkter et stykke fra 0 og deretter punkter nær 0. Sammenlikn ulike metoder med den analytiske løsninga og kommenter resultatene.

Videoer#

Bruk videoene nedenfor som en innføring eller en reptisjon til numerisk integrasjon.