3.5 Programmering av integrasjon | Matematikk R2

Implementere integrasjonsalgoritmer i Python.

55 min

12 oppgaver

PythonProgrammeringAlgoritmerNumerisk metode

Din fremgang i kapitlet

0 / 12 oppgaver

I dette kapitlet skal vi bruke programmering til å beregne bestemte integraler numerisk. Selv om mange integraler kan løses analytisk, finnes det utallige funksjoner som ikke har en lukket antiderivert. Da må vi ty til numeriske metoder.

Vi skal implementere tre klassiske metoder:
- Rektangelmetoden (den enkleste)
- Trapesmetoden (bedre nøyaktighet)
- Simpsons metode (svært god nøyaktighet)

I tillegg skal vi se på:
- Monte Carlo-integrasjon (stokastisk metode)
- Adaptiv integrasjon (automatisk steglengdejustering)
- Romberg-integrasjon (ekstrapolasjon)
- Gauss-kvadratur (optimal punktplassering)
- Dobbeltintegraler (integrasjon i 2D)

Vi bruker Python-bibliotekene NumPy og SciPy for effektiv implementering, og Matplotlib for visualisering.

Rektangelmetoden

Den enkleste tilnærmingen til numerisk integrasjon er å dele opp arealet under kurven i rektangler. Ideen er at summen av arealene til rektanglene tilnærmer det bestemte integralet.

Rektangelmetoden (venstrepunkt)

\displaystyle\int_a^b f(x)\,dx

✏️Implementere rektangelmetoden i Python

Skriv en Python-funksjon som beregner $\displaystyle\int_0^1 x^2\,dx$ ved hjelp av rektangelmetoden med $n = 1000$ delintervaller. Den eksakte verdien er $\displaystyle \frac{1}{3} \approx 0.333333$ .

``python def rektangelmetoden(f, a, b, n): """ Beregner det bestemte integralet av f fra a til b ved hjelp av rektangelmetoden (venstrepunkt).

Parametere: f: funksjonen som skal integreres a: nedre grense b: øvre grense n: antall delintervaller

Returnerer: Tilnærmet verdi av integralet """ h = (b - a) / n total = 0

for i in range(n): xi = a + i h total += f(xi)

return h total

`Definer funksjonen f(x) = x^2`


def f(x):
    return x2
Beregn integralet

resultat = rektangelmetoden(f, 0, 1, 1000)
print(f"Tilnærmet verdi: {resultat:.6f}")
print(f"Eksakt verdi: {1/3:.6f}")
print(f"Feil: {abs(resultat - 1/3):.6f}")

Utskrift:**`Tilnærmet verdi: 0.332834 Eksakt verdi: 0.333333 Feil: 0.000500``

Feilen er ca. $0.0005$ , noe som tilsvarer en relativ feil på ca. $0.15\%$ .

Midtpunktsmetoden

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

✏️Sammenligning av venstre- og midtpunktsmetoden

Sammenlign venstrepunktsmetoden og midtpunktsmetoden for $\displaystyle\int_0^1 x^2\,dx$ .

``python def venstrepunkt(f, a, b, n): h = (b - a) / n return h sum(f(a + i h) for i in range(n))

def midtpunkt(f, a, b, n): h = (b - a) / n return h sum(f(a + (i + 0.5) h) for i in range(n))

def f(x): return x*2

eksakt = 1/3

print("n\t\tVenstrepunkt\tMidtpunkt") print("-" 50) for n in [10, 100, 1000]: v = venstrepunkt(f, 0, 1, n) m = midtpunkt(f, 0, 1, n) print(f"{n}\t\t{abs(v - eksakt):.2e}\t{abs(m - eksakt):.2e}")`

Utskrift:`n Venstrepunkt Midtpunkt -------------------------------------------------- 10 4.50e-02 8.33e-04 100 4.95e-03 8.33e-06 1000 4.99e-04 8.33e-08``

Midtpunktsmetoden er dramatisk mer nøyaktig!

📝Oppgave 3.5.1

Implementer rektangelmetoden for å beregne $\displaystyle\int_0^{\pi} \sin(x)\,dx$ .

Den eksakte verdien er $2$ .

a) Bruk $n = 100$ delintervaller og finn tilnærmet verdi.

b) Bruk $n = 1000$ delintervaller og sammenlign feilen.

c) Hvor mange delintervaller trenger du for at feilen skal være mindre enn $0.0001$ ?

Oppgave 3.5.1

Medium

PythonAuto-lagret

Trapesmetoden

Trapesmetoden gir bedre nøyaktighet enn rektangelmetoden ved å tilnærme kurven med trapeser i stedet for rektangler. Et trapes fanger opp mer av kurvens form.

Trapesmetoden

\displaystyle\int_a^b f(x)\,dx

📜Feilestimatet for trapesmetoden

Feilen i trapesmetoden er gitt ved:

E_T = -\frac{(b-a)^3}{12n^2} f''(\xi)

for en eller annen

\xi \in [a, b]

Dette betyr at feilen avtar som $O(1/n^2)$ når $n$ øker, altså dobbelt så raskt som rektangelmetoden ( $O(1/n)$ ).

✏️Implementere trapesmetoden i Python

Skriv en Python-funksjon for trapesmetoden og beregn $\displaystyle\int_0^1 e^x\,dx$ . Den eksakte verdien er $e - 1 \approx 1.71828$ .

``python import math

def trapesmetoden(f, a, b, n): """ Beregner det bestemte integralet av f fra a til b ved hjelp av trapesmetoden. """ h = (b - a) / n

# Første og siste ledd (med vekt 1/2) total = (f(a) + f(b)) / 2

# Alle indre punkter (med vekt 1) for i in range(1, n): xi = a + i h total += f(xi)

return h total

def f(x): return math.exp(x)

eksakt = math.e - 1

`Test med forskjellige n`


for n in [10, 100, 1000]:
    resultat = trapesmetoden(f, 0, 1, n)
    feil = abs(resultat - eksakt)
    print(f"n = {n:4d}: {resultat:.8f}, feil = {feil:.2e}")

Utskrift:`n = 10: 1.71831994, feil = 3.87e-05 n = 100: 1.71828220, feil = 3.87e-07 n = 1000: 1.71828183, feil = 3.87e-09``

Merk at feilen reduseres med faktor 100 når $n$ øker med faktor 10, som bekrefter $O(1/n^2)$ -oppførselen.

📝Oppgave 3.5.2

Sammenlign rektangelmetoden og trapesmetoden for integralet $\displaystyle \displaystyle\int_1^3 \frac{1}{x}\,dx$ .

Den eksakte verdien er $\ln(3) \approx 1.0986$ .

a) Implementer begge metodene.

b) Beregn integralet med $n = 100$ for begge metoder.

c) Hvor mange ganger mer nøyaktig er trapesmetoden?

Oppgave 3.5.2

Medium

PythonAuto-lagret

Simpsons metode

Simpsons metode tilnærmer kurven med parabler (andregradsfunksjoner) i stedet for rette linjer. Dette gir enda bedre nøyaktighet enn trapesmetoden, spesielt for glatte funksjoner.

Simpsons metode (1/3-regelen)

\displaystyle\int_a^b f(x)\,dx

📜Feilestimatet for Simpsons metode

Feilen i Simpsons metode er gitt ved:

E_S = -\frac{(b-a)^5}{180n^4} f^{(4)}(\xi)

for en eller annen

\xi \in [a, b]

Dette betyr at feilen avtar som $O(1/n^4)$ , altså mye raskere enn trapesmetoden.

Bemerkelsverdig: Simpsons metode er eksakt for polynomer av grad 3 eller lavere!

✏️Implementere Simpsons metode i Python

Implementer Simpsons metode og beregn $\displaystyle\int_0^{\pi} \sin(x)\,dx$ med $n = 10$ . Sammenlign med trapesmetoden.

``python import math

def simpson(f, a, b, n): """ Beregner det bestemte integralet av f fra a til b ved hjelp av Simpsons metode. Krever at n er et partall. """ if n % 2 != 0: raise ValueError("n må være et partall for Simpsons metode")

h = (b - a) / n

# Første og siste ledd total = f(a) + f(b)

# Odde indekser (koeffisient 4) for i in range(1, n, 2): total += 4 f(a + i h)

# Like indekser (koeffisient 2), unntatt 0 og n for i in range(2, n, 2): total += 2 f(a + i h)

return (h / 3) total

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

def f(x): return math.sin(x)

eksakt = 2

`Sammenlign med n = 10`


trapes10 = trapesmetoden(f, 0, math.pi, 10)
simpson10 = simpson(f, 0, math.pi, 10)

print(f"Eksakt verdi: {eksakt}") print(f"Trapes (n=10): {trapes10:.10f}, feil = {abs(trapes10 - eksakt):.2e}") print(f"Simpson (n=10): {simpson10:.10f}, feil = {abs(simpson10 - eksakt):.2e}")`

Utskrift:`Eksakt verdi: 2 Trapes (n=10): 1.9835235376, feil = 1.65e-02 Simpson (n=10): 2.0000067844, feil = 6.78e-06``

Med bare $n = 10$ har Simpsons metode en feil på ca. $7 \cdot 10^{-6}$ , mens trapesmetoden har feil på ca. $0.017$ . Simpsons metode er over 2000 ganger mer nøyaktig!

📝Oppgave 3.5.3

Bruk Simpsons metode til å beregne $\displaystyle\int_0^1 e^{-x^2}\,dx$ .

Dette integralet har ingen lukket antiderivert, så vi må bruke numeriske metoder!

a) Implementer Simpsons metode.

b) Finn en tilnærmet verdi med $n = 100$ .

c) Den "eksakte" verdien (beregnet med mange desimaler) er ca. $0.746824$ . Hva er feilen?

Oppgave 3.5.3

Vanskelig

PythonAuto-lagret

Simpsons 3/8-regel

Det finnes også en variant av Simpsons metode som bruker kubisk interpolasjon over tre delintervaller.

Simpsons 3/8-regel

\int_a^b f(x)\,dx \approx \frac{3h}{8}\left[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)\right]

Bruke NumPy og SciPy for integrasjon

I praksis trenger vi ikke alltid å implementere integrasjonsmetoder fra bunnen av. Python-bibliotekene NumPy og SciPy har ferdige, optimaliserte funksjoner for numerisk integrasjon.

SciPy integrate-modulen

n

✏️Integrasjon med SciPy

Bruk SciPy til å beregne $\displaystyle \displaystyle\int_0^{\infty} e^{-x^2}\,dx = \frac{\sqrt{\pi}}{2}$ .

``python import numpy as np from scipy import integrate

`Definer funksjonen`


def f(x):
    return np.exp(-x2)
Beregn integralet fra 0 til uendelig

resultat, feil = integrate.quad(f, 0, np.inf)
eksakt = np.sqrt(np.pi) / 2

print(f"SciPy quad: {resultat:.10f}") print(f"Eksakt verdi: {eksakt:.10f}") print(f"Feilestimering fra quad: {feil:.2e}") print(f"Faktisk feil: {abs(resultat - eksakt):.2e}")`

Utskrift:`SciPy quad: 0.8862269255 Eksakt verdi: 0.8862269255 Feilestimering fra quad: 7.10e-09 Faktisk feil: 0.00e+00``

Merk at quad** håndterer uendelige grenser automatisk og gir et feilestimering!

✏️Trapesmetoden med NumPy på datapunkter

Gitt målepunkter $(x_i, y_i)$ fra et eksperiment, beregn det tilnærmede arealet under kurven.

``python import numpy as np

`Simulerte måledata (f.eks. hastighet over tid)`


t = np.array([0, 1, 2, 3, 4, 5])  # tid i sekunder
v = np.array([0, 5, 8, 10, 9, 6])  # hastighet i m/s
Beregn tilbakelagt strekning (integral av hastighet)

strekning = np.trapz(v, t)

print(f"Tidspunkter: {t}") print(f"Hastigheter: {v}") print(f"Tilbakelagt strekning: {strekning:.1f} meter")`

Utskrift:`Tidspunkter: [0 1 2 3 4 5] Hastigheter: [ 0 5 8 10 9 6] Tilbakelagt strekning: 34.5 meter``

Dette er svært nyttig når vi har diskrete måledata i stedet for en analytisk funksjon.

📝Oppgave 3.5.4

Bruk SciPy til å beregne følgende integraler:

a) $\displaystyle \displaystyle\int_0^1 \frac{1}{1+x^2}\,dx = \frac{\pi}{4}$

b) $\displaystyle\int_1^e \ln(x)\,dx = 1$

c) $\displaystyle \displaystyle\int_0^{\pi/2} \cos^2(x)\,dx = \frac{\pi}{4}$

Sammenlign med de eksakte verdiene.

Oppgave 3.5.4

Medium

PythonAuto-lagret

Visualisering av numerisk integrasjon

En viktig del av å forstå numeriske metoder er å visualisere hva som skjer. Med Matplotlib kan vi tegne funksjonen og rektanglene/trapesene som brukes til å tilnærme arealet.

✏️Visualisere rektangelmetoden

Lag en visualisering som viser rektanglene brukt i rektangelmetoden for $f(x) = x^2$ på $[0, 2]$ med $n = 8$ delintervaller.

``python import numpy as np import matplotlib.pyplot as plt

def f(x): return x*2

a, b, n = 0, 2, 8 h = (b - a) / n

`Lag figur`


fig, ax = plt.subplots(figsize=(10, 6))
Tegn funksjonen

xsmooth = np.linspace(a, b, 200)
ax.plot(xsmooth, f(xsmooth), 'b-', linewidth=2, label=' $f(x) = x^2$ ')
Tegn rektanglene

for i in range(n):
    xleft = a + i  h
    rect = plt.Rectangle((xleft, 0), h, f(xleft),
                          fill=True, facecolor='lightblue',
                          edgecolor='blue', alpha=0.7)
    ax.addpatch(rect)
Beregn tilnærmet integral

integralapprox = h  sum(f(a + i  h) for i in range(n))
integralexact = 8/3ax.setxlabel('x', fontsize=12)
ax.setylabel('y', fontsize=12)
ax.settitle(f'Rektangelmetoden: n = {n}\n'
             f'Tilnærmet: {integralapprox:.4f}, '
             f'Eksakt: {integralexact:.4f}')
ax.legend()
ax.grid(True, alpha=0.3)
ax.setxlim(a - 0.1, b + 0.1)
ax.setylim(0, f(b) + 0.5)

plt.tight_layout() plt.show()``

Denne koden produserer en figur som viser:
- Den blå kurven $f(x) = x^2$
- Lyseblå rektangler som representerer tilnærmingen
- Tittel med både tilnærmet og eksakt verdi

✏️Visualisere trapesmetoden

Lag en visualisering som viser trapesmetoden for $f(x) = \sin(x)$ på $[0, \pi]$ med $n = 6$ delintervaller.

``python import numpy as np import matplotlib.pyplot as plt

def f(x): return np.sin(x)

a, b, n = 0, np.pi, 6 h = (b - a) / n

fig, ax = plt.subplots(figsize=(10, 6))

`Tegn funksjonen`


xsmooth = np.linspace(a, b, 200)
ax.plot(xsmooth, f(xsmooth), 'b-', linewidth=2, label=' $f(x) = \sin(x)$ ')
Tegn trapesene

xtrap = [a + i  h for i in range(n + 1)]
ytrap = [f(x) for x in xtrap]
for i in range(n):
    xcorners = [xtrap[i], xtrap[i+1], xtrap[i+1], xtrap[i]]
    ycorners = [0, 0, ytrap[i+1], ytrap[i]]
    ax.fill(xcorners, ycorners, alpha=0.3, color='green', edgecolor='darkgreen')
Marker punktene

ax.plot(xtrap, ytrap, 'go', markersize=8)
Beregn integralet

integralapprox = h  ((f(a) + f(b))/2 + sum(f(a + i*h) for i in range(1, n)))ax.setxlabel('x', fontsize=12)
ax.setylabel('y', fontsize=12)
ax.settitle(f'Trapesmetoden: n = {n}\n'
             f'Tilnærmet: {integralapprox:.4f}, Eksakt: 2.0000')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tightlayout() plt.show()``

📝Oppgave 3.5.5

Lag en visualisering som sammenligner rektangelmetoden, trapesmetoden og Simpsons metode for $f(x) = x^3$ på $[0, 1]$ med $n = 4$ .

Vis tre subplots side om side som illustrerer hver metode, og beregn feilen for hver.

Oppgave 3.5.5

Vanskelig

PythonAuto-lagret

Monte Carlo-integrasjon

En helt annen tilnærming til numerisk integrasjon er Monte Carlo-metoden, som bruker tilfeldige tall. Ideen er enkel: kast tilfeldige punkter i et rektangel og tell hvor mange som lander under kurven.

Monte Carlo-integrasjon

\displaystyle\int_a^b f(x)\,dx

✏️Monte Carlo-integrasjon i Python

Bruk Monte Carlo-metoden til å estimere $\pi$ ved å beregne $\displaystyle \displaystyle\int_{-1}^{1} \sqrt{1-x^2}\,dx = \frac{\pi}{2}$ .

``python import numpy as np

def montecarlointegral(f, a, b, N): """ Beregner integralet av f fra a til b ved hjelp av Monte Carlo-metoden. """ # Generer N tilfeldige x-verdier i [a, b] xrandom = np.random.uniform(a, b, N)

# Beregn funksjonsverdiene yvalues = f(xrandom)

# Estimer integralet integral = (b - a) np.mean(yvalues)

return integral

`Halvsirkel: f(x) = sqrt(1 - x^2)`


def f(x):
    return np.sqrt(1 - x2)
Kjør flere ganger for å se variansen

np.random.seed(42)  # For reproduserbarhet

for N in [100, 1000, 10000, 100000, 1000000]: resultat = montecarlointegral(f, -1, 1, N) piestimat = 2 resultat # Hele sirkelen feil = abs(piestimat - np.pi) print(f"N = {N:7d}: pi ≈ {pi_estimat:.6f}, feil = {feil:.6f}")`

Utskrift:`N = 100: pi ≈ 3.019851, feil = 0.121742 N = 1000: pi ≈ 3.161035, feil = 0.019443 N = 10000: pi ≈ 3.145209, feil = 0.003617 N = 100000: pi ≈ 3.140063, feil = 0.001530 N = 1000000: pi ≈ 3.141039, feil = 0.000554``

Merk at feilen avtar som $1/\sqrt{N}$ : for å halvere feilen må vi firedoble antall punkter.

📝Oppgave 3.5.6

Bruk Monte Carlo-integrasjon til å beregne volumet av en enhetssfære i 3D.

Volumet av en sfære med radius $r$ er $\displaystyle V = \frac{4}{3}\pi r^3$ , så for $r = 1$ er $\displaystyle V = \frac{4}{3}\pi \approx 4.189$ .

Hint: Generer tilfeldige punkter i kuben $[-1, 1]^3$ og tell hvor mange som er innenfor sfæren (dvs. $x^2 + y^2 + z^2 \leq 1$ ).

Oppgave 3.5.6

Vanskelig

PythonAuto-lagret

Feilanalyse og konvergens

Når vi bruker numeriske metoder, er det viktig å forstå hvor nøyaktig resultatet er og hvordan feilen oppfører seg når vi øker antall delintervaller.

📜Konvergensordener for integrasjonsmetoder

Metode	Feilorden	Dobling av $n$
Rektangel (venstre)	$O(h) = O(1/n)$	Halverer feilen
Rektangel (midtpunkt)	$O(h^2) = O(1/n^2)$	Kvarterer feilen
Trapes	$O(h^2) = O(1/n^2)$	Kvarterer feilen
Simpson	$O(h^4) = O(1/n^4)$	Reduserer feilen med faktor 16
Romberg	$O(h^{2k})$	Avhenger av ekstrapolasjonsnivå
Monte Carlo	$O(1/\sqrt{N})$	Reduserer feilen med faktor $\sqrt{2}$

Konvergensfaktor: Forholdet mellom feil for

n

og feil for

2n

bør tilnærme:
- Rektangel (venstre): 2
- Trapes/midtpunkt: 4
- Simpson: 16

✏️Verifisere konvergensorden

Beregn konvergensordenen til trapesmetoden for $\displaystyle\int_0^1 e^x\,dx$ ved å sammenligne feil for $n = 10, 20, 40, 80$ .

``python import math

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

def f(x): return math.exp(x)

eksakt = math.e - 1

print("n\t\tFeil\t\tFaktor") print("-" * 40)

prevfeil = None for n in [10, 20, 40, 80, 160, 320]: resultat = trapesmetoden(f, 0, 1, n) feil = abs(resultat - eksakt)

if prevfeil is not None: faktor = prevfeil / feil print(f"{n}\t\t{feil:.2e}\t{faktor:.2f}") else: print(f"{n}\t\t{feil:.2e}\t-")

prevfeil = feil`

Utskrift:`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 160 1.51e-07 4.00 320 3.78e-08 4.00``

Faktoren er konsistent 4, som bekrefter at trapesmetoden har konvergensorden $O(1/n^2)$ .

📝Oppgave 3.5.7

Lag et program som sammenligner konvergensen til rektangelmetoden, trapesmetoden og Simpsons metode for $\displaystyle\int_0^{\pi} \sin(x)\,dx = 2$ .

a) Beregn feilen for $n = 10, 20, 40, 80, 160$ for alle tre metoder.

b) Lag en log-log plot av feil vs. $n$ for alle tre metoder.

c) Bekreft konvergensordenene fra plottet (helningen skal være -1, -2, -4 for de respektive metodene).

Oppgave 3.5.7

Vanskelig

PythonAuto-lagret

Adaptiv integrasjon

Så langt har vi brukt jevnt fordelte delintervaller. Men noen funksjoner varierer mer i noen områder enn andre. Adaptiv integrasjon justerer steglengden automatisk: små steg der funksjonen varierer mye, store steg der den er jevn.

Adaptiv Simpson-integrasjon

S_1 = S(a, b)

✏️Implementere adaptiv Simpson

Implementer adaptiv Simpson-integrasjon for funksjoner som varierer mye lokalt, som $\displaystyle f(x) = \frac{1}{(x-0.3)^2 + 0.01}$ på $[0, 1]$ .

``python import numpy as np

def simpsonenkel(f, a, b): """Simpson over ett intervall [a, b]""" m = (a + b) / 2 h = (b - a) / 6 return h (f(a) + 4f(m) + f(b))

def adaptivsimpson(f, a, b, tol=1e-8, maxdepth=50): """ Adaptiv Simpson-integrasjon.

Parametere: f: funksjonen a, b: integrasjonsgrenser tol: toleranse for feil maxdepth: maksimum rekursjonsdybde """ def adaptive(a, b, Sab, tol, depth): m = (a + b) / 2 Sam = simpsonenkel(f, a, m) Smb = simpsonenkel(f, m, b) Ssum = Sam + Smb

# Sjekk om vi er nøyaktige nok if depth >= maxdepth or abs(Ssum - Sab) < 15 tol: return Ssum + (Ssum - Sab) / 15 # Richardson-ekstrapolering

# Rekursivt kall på hver halvdel return (adaptive(a, m, Sam, tol/2, depth+1) + adaptive(m, b, Smb, tol/2, depth+1))

Sab = simpsonenkel(f, a, b) return adaptive(a, b, Sab, tol, 0)

`Test på funksjon med skarp topp`


def f(x):
    return 1 / ((x - 0.3)2 + 0.01)
Beregn integralet

resultat = adaptivsimpson(f, 0, 1)
Sammenlign med SciPy

from scipy import integrate
scipyresultat,  = integrate.quad(f, 0, 1)

print(f"Adaptiv Simpson: {resultat:.10f}") print(f"SciPy quad: {scipyresultat:.10f}") print(f"Forskjell: {abs(resultat - scipyresultat):.2e}")`

Utskrift:*`Adaptiv Simpson: 15.7079632679 SciPy quad: 15.7079632679 Forskjell: 1.42e-14``

Den adaptive metoden konsentrerer beregningene rundt toppen ved $x = 0.3$ .

📝Oppgave 3.5.8

Sammenlign antall funksjonsevalueringer for vanlig Simpson og adaptiv Simpson for funksjonen

$f(x) = \sin(100x)$

på intervallet $[0, 1]$ .

a) Implementer en teller for funksjonsevalueringer.

b) Finn hvor mange evalueringer vanlig Simpson trenger for å oppnå feil $< 10^{-8}$ .

c) Tell evalueringer for adaptiv Simpson med toleranse $10^{-8}$ .

Oppgave 3.5.8

Vanskelig

PythonAuto-lagret

Romberg-integrasjon

Romberg-integrasjon kombinerer trapesmetoden med Richardson-ekstrapolering for å oppnå høyere nøyaktighet. Ideen er å bruke resultatene fra trapesmetoden med forskjellige $n$ til å "ekstrapolere bort" feilene.

Romberg-algoritmen

R_{i,j}

✏️Implementere Romberg-integrasjon

Implementer Romberg-integrasjon og beregn $\displaystyle\int_0^1 e^x\,dx = e - 1$ .

``python import numpy as np

def romberg(f, a, b, maxiter=10, tol=1e-12): """ Romberg-integrasjon.

Returnerer: (resultat, tabell) - tilnærmet verdi og Romberg-tabellen """ R = np.zeros((maxiter, maxiter))

# Første kolonne: trapesmetoden h = b - a R[0, 0] = h (f(a) + f(b)) / 2

for i in range(1, maxiter): h = h / 2 n = 2i

# Legg til nye midtpunkter sumnew = sum(f(a + (2k - 1) h) for k in range(1, n//2 + 1)) R[i, 0] = R[i-1, 0] / 2 + h sumnew

# Richardson-ekstrapolering for j in range(1, i + 1): R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4*j - 1)

# Sjekk konvergens if i > 0 and abs(R[i, i] - R[i-1, i-1]) < tol: return R[i, i], R[:i+1, :i+1]

return R[maxiter-1, maxiter-1], R

def f(x): return np.exp(x)

eksakt = np.e - 1

resultat, tabell = romberg(f, 0, 1, max_iter=6)

print("Romberg-tabell:") print("-" 70) for i in range(tabell.shape[0]): row = " ".join(f"{tabell[i,j]:.10f}" for j in range(i+1)) print(f"i={i}: {row}")

print(f"\nResultat: {resultat:.12f}") print(f"Eksakt: {eksakt:.12f}") print(f"Feil: {abs(resultat - eksakt):.2e}")`

Utskrift:`Romberg-tabell: ---------------------------------------------------------------------- i=0: 1.8591409142 i=1: 1.7539425894 1.7188753977 i=2: 1.7272219005 1.7183147709 1.7182773958 i=3: 1.7195601466 1.7182729313 1.7182818407 1.7182818284 i=4: 1.7185267584 1.7182156290 1.7182818110 1.7182818285 1.7182818285 i=5: 1.7183429637 1.7182816988 1.7182818285 1.7182818285 1.7182818285

Resultat: 1.718281828459 Eksakt: 1.718281828459 Feil: 4.44e-16``

Merk den raske konvergensen: etter bare 5-6 iterasjoner har vi maskin-presisjon!

📝Oppgave 3.5.9

Bruk scipy.integrate.romberg til å beregne $\displaystyle\int_0^{\pi} \sin(x)\,dx = 2$ .

a) Sammenlign med scipy.integrate.quad.

b) Undersøk hvordan romberg-funksjonen rapporterer konvergens ved å bruke show=True.

Oppgave 3.5.9

Vanskelig

PythonAuto-lagret

Gauss-kvadratur

Gauss-kvadratur er en familie av integrasjonsmetoder som oppnår maksimal nøyaktighet for et gitt antall funksjonsevalueringer. I stedet for å bruke jevnt fordelte punkter, bruker Gauss-metodene optimalt plasserte punkter og vekter.

Gauss-Legendre-kvadratur

\int_{-1}^{1} f(x)\,dx \approx \sum_{i=1}^{n} w_i f(x_i)

✏️Gauss-Legendre-kvadratur med SciPy

Bruk Gauss-Legendre-kvadratur til å beregne $\displaystyle\int_0^1 e^{-x^2}\,dx$ med forskjellige antall punkter.

``python import numpy as np from scipy import integrate

def f(x): return np.exp(-x*2)

`Eksakt verdi (fra scipy.quad)`


eksakt,  = integrate.quad(f, 0, 1)print("Gauss-Legendre-kvadratur:")
print("-"  50)for n in [2, 3, 4, 5, 6, 8, 10]:
    resultat,  = integrate.fixedquad(f, 0, 1, n=n)
    feil = abs(resultat - eksakt)
    print(f"n = {n:2d}: {resultat:.12f}, feil = {feil:.2e}")
Manuell implementasjon

print("\nManuell Gauss-Legendre med numpy.polynomial.legendre:")
from numpy.polynomial.legendre import leggaussdef gausslegendre(f, a, b, n):
    """Gauss-Legendre-kvadratur med n punkter."""
    # Få punkter og vekter for [-1, 1]
    x, w = leggauss(n)
    # Transformer til [a, b]
    xtrans = (b - a) / 2  x + (a + b) / 2
    scale = (b - a) / 2
    return scale  np.sum(w * f(xtrans))

for n in [3, 5, 7]: resultat = gauss_legendre(f, 0, 1, n) print(f"n = {n}: {resultat:.12f}, feil = {abs(resultat - eksakt):.2e}")`

Utskrift:`Gauss-Legendre-kvadratur: -------------------------------------------------- n = 2: 0.746516493811, feil = 3.08e-04 n = 3: 0.746825356951, feil = 1.22e-06 n = 4: 0.746824126088, feil = 6.72e-09 n = 5: 0.746824132859, feil = 4.61e-11 n = 6: 0.746824132812, feil = 4.00e-13 n = 8: 0.746824132812, feil = 1.11e-16 n = 10: 0.746824132812, feil = 0.00e+00``

Med bare 10 punkter oppnår vi maskin-presisjon!

📝Oppgave 3.5.10

Sammenlign Gauss-Legendre-kvadratur med Simpsons metode for følgende integraler:

a) $\displaystyle \displaystyle\int_0^1 x^5\,dx = \frac{1}{6}$ (polynom av grad 5)

b) $\displaystyle\int_0^{\pi} \cos(x)\,dx = 0$ (oscillerende)

c) $\displaystyle \displaystyle\int_0^1 \sqrt{x}\,dx = \frac{2}{3}$ (singulær derivert ved $x = 0$ )

For hver: finn antall punkter/delintervaller som gir feil $< 10^{-10}$ .

Oppgave 3.5.10

Vanskelig

PythonAuto-lagret

Dobbeltintegraler numerisk

Numerisk integrasjon i to dimensjoner er en naturlig utvidelse av metodene vi har sett. Vi kan beregne

$\iint_R f(x, y)\,dA$

ved å gjenta én-dimensjonal integrasjon.

Dobbeltintegral med iterert integrasjon

R = [a, b] \times [c, d]

✏️Beregne dobbeltintegral med SciPy

Beregn $\displaystyle\iint_R e^{-(x^2 + y^2)}\,dA$ der $R = [0, 1] \times [0, 1]$ .

``python import numpy as np from scipy import integrate

def f(y, x): # NB: dblquad bruker rekkefølgen (y, x) return np.exp(-(x2 + y2))

`Beregn dobbeltintegralet`


resultat, feil = integrate.dblquad(f, 0, 1, 0, 1)
print(f"Dobbeltintegral: {resultat:.10f}")
print(f"Feilestimering:  {feil:.2e}")
Sammenlign med manuell iterert integrasjon

def innerintegral(x):
    """Integrerer f over y for fast x"""
    result,  = integrate.quad(lambda y: np.exp(-(x2 + y2)), 0, 1)
    return result
resultatiterert,  = integrate.quad(innerintegral, 0, 1)
print(f"Iterert:         {resultatiterert:.10f}")
Manuell 2D Simpson

def simpson2d(f, ax, bx, ay, by, nx, ny):
    """2D Simpson-integrasjon"""
    hx = (bx - ax) / nx
    hy = (by - ay) / ny    def simpson1d(g, a, b, n):
        h = (b - a) / n
        return (h/3)  (g(a) + 4sum(g(a + ih) for i in range(1, n, 2))
                        + 2sum(g(a + i*h) for i in range(2, n, 2)) + g(b))
    def inner(x):
        return simpson1d(lambda y: f(x, y), ay, by, ny)
    return simpson1d(inner, ax, bx, nx)
def fxy(x, y):
    return np.exp(-(x2 + y2))

resultatsimpson = simpson2d(fxy, 0, 1, 0, 1, 20, 20) print(f"Simpson 2D: {resultat_simpson:.10f}")`

Utskrift:`Dobbeltintegral: 0.5577462854 Feilestimering: 6.19e-15 Iterert: 0.5577462854 Simpson 2D: 0.5577462854``

✏️Dobbeltintegral over ikke-rektangulært område

Beregn arealet av enhetssirkelen ved å beregne $\displaystyle\iint_D 1\,dA$ der $D$ er sirkelen $x^2 + y^2 \leq 1$ .

``python import numpy as np from scipy import integrate

`f(y, x) = 1 (vi beregner areal)`


def f(y, x):
    return 1.0
Grenser for y som funksjon av x

def ylower(x):
    return -np.sqrt(1 - x2)def yupper(x):
    return np.sqrt(1 - x2)
Beregn dobbeltintegralet

areal, feil = integrate.dblquad(f, -1, 1, ylower, yupper)
print(f"Beregnet åreal: {areal:.10f}")
print(f"Eksakt (pi):    {np.pi:.10f}")
print(f"Feil:           {abs(areal - np.pi):.2e}")
Alternativ: Monte Carlo

N = 1000000
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
inside = np.sum(x2 + y2 <= 1)
arealmc = 4 * (inside / N)  # 4 = areal av kvadratet [-1,1] x [-1,1]

print(f"\nMonte Carlo: {arealmc:.4f}") print(f"Feil MC: {abs(areal_mc - np.pi):.4f}")`

Utskrift:`Beregnet åreal: 3.1415926536 Eksakt (pi): 3.1415926536 Feil: 1.23e-13

Monte Carlo: 3.1420 Feil MC: 0.0004``

📝Oppgave 3.5.11

Beregn volumet under paraboloiden $z = 4 - x^2 - y^2$ over sirkelen $x^2 + y^2 \leq 4$ .

Det eksakte svaret er $\displaystyle\iint_D (4 - x^2 - y^2)\,dA = 8\pi$ .

a) Bruk scipy.integrate.dblquad.

b) Bruk Monte Carlo med $N = 100000$ punkter.

c) Visualiser paraboloiden og integrasjonsområdet.

Oppgave 3.5.11

Vanskelig

PythonAuto-lagret

Praktiske anvendelser

Numerisk integrasjon har utallige anvendelser i vitenskap og ingeniørfag. Her ser vi på noen eksempler.

✏️Beregne arbeid utført av en variabel kraft

En fjær følger Hookes lov med $F(x) = kx$ for små deformasjoner, men for store deformasjoner blir den ikke-lineær: $F(x) = kx + \alpha x^3$ . Beregn arbeidet som kreves for å strekke fjæren fra $x = 0$ til $x = L$ .

``python import numpy as np from scipy import integrate import matplotlib.pyplot as plt

`Parametre`


k = 100  # N/m (fjærstivhet)
alpha = 500  # N/m^3 (ikke-lineær koeffisient)
L = 0.1  # m (strekning)
def F(x):
    """Kraft som funksjon av posisjon"""
    return k  x + alpha  x*3
Beregn arbeidet W = integral av F(x) dx fra 0 til L

arbeid, feil = integrate.quad(F, 0, L)
Sammenlign med lineær fjær (kun kx)

arbeidlinear = 0.5  k  L2print(f"Ikke-lineær fjær:")
print(f"  Arbeid = {arbeid:.4f} J")
print(f"\nLineær fjær:")
print(f"  Arbeid = {arbeidlinear:.4f} J")
print(f"\nForskjell: {(arbeid - arbeidlinear) / arbeidlinear  100:.1f}%")
Visualiser

x = np.linspace(0, L, 100)
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(x, F(x), 'b-', label='Ikke-lineær:  $F = kx + \alpha x^3$ ')
plt.plot(x, k * x, 'r--', label='Lineær:  $F = kx$ ')
plt.fillbetween(x, 0, F(x), alpha=0.3)
plt.xlabel('x (m)')
plt.ylabel('F (N)')
plt.title('Kraft vs. posisjon')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
Kumulativt arbeid

arbeidkumulativ = [integrate.quad(F, 0, xi)[0] for xi in x]
plt.plot(x, arbeidkumulativ, 'b-', label='Arbeid')
plt.xlabel('x (m)')
plt.ylabel('W (J)')
plt.title('Kumulativt arbeid')
plt.grid(True, alpha=0.3)

plt.tightlayout() plt.show()``

✏️Beregne buelengde

Beregn buelengden til kurven $y = \sin(x)$ fra $x = 0$ til $x = 2\pi$ .

``python import numpy as np from scipy import integrate import matplotlib.pyplot as plt

`Buelengde: L = integral av sqrt(1 + (dy/dx)^2) dx`


For y = sin(x): dy/dx = cos(x)
def integrand(x):
    return np.sqrt(1 + np.cos(x)*2)
Beregn buelengden

buelengde, feil = integrate.quad(integrand, 0, 2np.pi)
print(f"Buelengde av sin(x) fra 0 til 2pi:")
print(f"L = {buelengde:.6f}")
print(f"(Dette er ca. {buelengde/np.pi:.4f}  pi)")
For referanse: lengden av en rett linje fra (0,0) til (2pi,0) er 2pi

print(f"\nRett linje: {2np.pi:.6f}")
print(f"Forlengelse: {(buelengde - 2np.pi) / (2np.pi)  100:.1f}%")
Visualiser

x = np.linspace(0, 2*np.pi, 200)
y = np.sin(x)

plt.figure(figsize=(10, 4)) plt.plot(x, y, 'b-', linewidth=2) plt.axhline(0, color='gray', linestyle='--', alpha=0.5) plt.fillbetween(x, 0, y, where=(y > 0), alpha=0.2, color='blue') plt.fillbetween(x, 0, y, where=(y < 0), alpha=0.2, color='red') plt.xlabel('x') plt.ylabel('y') plt.title(f' $y = \sin(x)$ , buelengde $L = {buelengde:.4f}$ ') plt.grid(True, alpha=0.3) plt.axis('equal') plt.show()``

📝Oppgave 3.5.12

Fysikkoppgave: En partikkel beveger seg med hastighet

v(t) = 10\sin(t)e^{-0.1t}

m/s.

a) Finn total tilbakelagt strekning i tidsintervallet $t \in [0, 20]$ s.
(Hint: strekning = $\int |v(t)|\,dt$ , ikke $\int v(t)\,dt$ )

b) Finn netto forflytning (endelig posisjon minus startposisjon).

c) Visualiser $v(t)$ , posisjonen $x(t)$ , og kumulativ strekning.

Oppgave 3.5.12

Vanskelig

PythonAuto-lagret

Oppsummering

Vi har lært om følgende numeriske integrasjonsmetoder:

Metode	Nøyaktighet	Implementering	Bruksområde
Rektangel (venstre)	$O(1/n)$	Enklest	Undervisning
Rektangel (midtpunkt)	$O(1/n^2)$	Enkel	Rask tilnærming
Trapes	$O(1/n^2)$	Enkel	God allround
Simpson	$O(1/n^4)$	Moderat	Glatte funksjoner
Romberg	$O(1/n^{2k})$	Moderat	Høy nøyaktighet
Gauss-Legendre	Optimal	Krever vekter	Polynomer, glatte
Adaptiv	Variabel	Kompleks	Variable funksjoner
Monte Carlo	$O(1/\sqrt{N})$	Enkel	Høydimensjonale

Praktiske tips:
1. Bruk scipy.integrate.quad for de fleste 1D-oppgaver
2. Bruk scipy.integrate.dblquad for 2D-integraler
3. Bruk Monte Carlo for dimensjon

\geq 3

4. Bruk Romberg eller Gauss for høy presisjon
5. Verifiser alltid med kjente integraler eller ved å øke

n

6. Visualiser for å bygge intuisjon!

📝Oppgave 3.5.13

Prosjektoppgave: Normalfordelingen

Sannsynlighetstettheten for normalfordelingen (Gauss-fordelingen) med gjennomsnitt $\mu = 0$ og standardavvik $\sigma = 1$ er:

$f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$

a) Verifiser at $\displaystyle\int_{-\infty}^{\infty} f(x)\,dx = 1$ ved bruk av scipy.integrate.quad.

b) Beregn sannsynligheten $P(-1 \leq X \leq 1)$ , altså arealet under kurven fra $-1$ til $1$ . Dette skal være ca. $0.6827$ (68.27%).

c) Beregn $P(-2 \leq X \leq 2)$ og $P(-3 \leq X \leq 3)$ . Sammenlign med "68-95-99.7-regelen".

d) Visualiser normalfordelingen og skyggelegg området fra $-1$ til $1$ .

Oppgave 3.5.13

Vanskelig

PythonAuto-lagret

📝Oppgave 3.5.14

Bonusoppgave: Lag en interaktiv sammenligning

Lag et program som:

a) Tar inn en funksjon $f(x)$ , grenser $[a, b]$ , og antall punkter $n$ .

b) Beregner integralet med alle metodene: rektangel, midtpunkt, trapes, Simpson, Romberg og Gauss-Legendre.

c) Viser en tabell med resultatene og feilene (bruk scipy.quad som referanse).

d) Rangerer metodene etter nøyaktighet for den gitte funksjonen.

Test med:
- $f(x) = x^4$ på $[0, 1]$ (polynom)
- $f(x) = \sin(10x)$ på $[0, \pi]$ (oscillerende)
- $f(x) = e^{-x^2}$ på $[-2, 2]$ (Gaussisk)

Oppgave 3.5.14

Vanskelig

PythonAuto-lagret

Vanlige feil og fallgruver

1. Singulariteter:
Hvis integranden har singulariteter (uendeligheter), vil standard metoder feile. Bruk spesialiserte metoder eller del opp intervallet.

2. Oscillerende funksjoner:
Funksjoner som $\sin(1000x)$ krever svært mange delintervaller. Bruk adaptiv integrasjon.

3. Feil rekkefølge i dblquad:
scipy.integrate.dblquad bruker f(y, x), ikke f(x, y)!

4. Glemme partallskrav:
Simpsons metode krever at $n$ er et partall.

5. Avrundingsfeil:
For svært små $h$ kan avrundingsfeil dominere. Bruk høyere-ordens metoder i stedet for å øke $n$ .