3.5 Programmering av integrasjon | Matematikk R2

Implementere integrasjonsalgoritmer i Python.

55 min

19 oppgaver

PythonProgrammeringAlgoritmerNumerisk metode

Du leser den tradisjonelle versjonen

Din fremgang i kapitlet

0 / 19 oppgaver

Kapitlets plass i kurset

Bygger på

3.4Numerisk integrasjon

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$ .

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):
        x_i = a + i * h
        total += f(x_i)

    return h * total

# Definer funksjonen f(x) = x^2
def f(x):
    return x**2

# 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$ .

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$ .

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):
        x_i = a + i * h
        total += f(x_i)

    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.

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
trapes_10 = trapesmetoden(f, 0, math.pi, 10)
simpson_10 = simpson(f, 0, math.pi, 10)

print(f"Eksakt verdi: {eksakt}")
print(f"Trapes (n=10):  {trapes_10:.10f}, feil = {abs(trapes_10 - eksakt):.2e}")
print(f"Simpson (n=10): {simpson_10:.10f}, feil = {abs(simpson_10 - 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}$ .

import numpy as np
from scipy import integrate

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

# 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.

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.

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
x_smooth = np.linspace(a, b, 200)
ax.plot(x_smooth, f(x_smooth), 'b-', linewidth=2, label='$f(x) = x^2$')

# Tegn rektanglene
for i in range(n):
    x_left = a + i * h
    rect = plt.Rectangle((x_left, 0), h, f(x_left),
                          fill=True, facecolor='lightblue',
                          edgecolor='blue', alpha=0.7)
    ax.add_patch(rect)

# Beregn tilnærmet integral
integral_approx = h * sum(f(a + i * h) for i in range(n))
integral_exact = 8/3

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title(f'Rektangelmetoden: n = {n}\n'
             f'Tilnærmet: {integral_approx:.4f}, '
             f'Eksakt: {integral_exact:.4f}')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(a - 0.1, b + 0.1)
ax.set_ylim(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.

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
x_smooth = np.linspace(a, b, 200)
ax.plot(x_smooth, f(x_smooth), 'b-', linewidth=2, label='$f(x) = \sin(x)$')

# Tegn trapesene
x_trap = [a + i * h for i in range(n + 1)]
y_trap = [f(x) for x in x_trap]

for i in range(n):
    x_corners = [x_trap[i], x_trap[i+1], x_trap[i+1], x_trap[i]]
    y_corners = [0, 0, y_trap[i+1], y_trap[i]]
    ax.fill(x_corners, y_corners, alpha=0.3, color='green', edgecolor='darkgreen')

# Marker punktene
ax.plot(x_trap, y_trap, 'go', markersize=8)

# Beregn integralet
integral_approx = h * ((f(a) + f(b))/2 + sum(f(a + i*h) for i in range(1, n)))

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title(f'Trapesmetoden: n = {n}\n'
             f'Tilnærmet: {integral_approx:.4f}, Eksakt: 2.0000')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
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}$ .

import numpy as np

def monte_carlo_integral(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]
    x_random = np.random.uniform(a, b, N)

    # Beregn funksjonsverdiene
    y_values = f(x_random)

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

    return integral

# Halvsirkel: f(x) = sqrt(1 - x^2)
def f(x):
    return np.sqrt(1 - x**2)

# Kjør flere ganger for å se variansen
np.random.seed(42)  # For reproduserbarhet

for N in [100, 1000, 10000, 100000, 1000000]:
    resultat = monte_carlo_integral(f, -1, 1, N)
    pi_estimat = 2 * resultat  # Hele sirkelen
    feil = abs(pi_estimat - 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$ .

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)

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

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

    prev_feil = 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]$ .

import numpy as np

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

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

    Parametere:
        f: funksjonen
        a, b: integrasjonsgrenser
        tol: toleranse for feil
        max_depth: maksimum rekursjonsdybde
    """
    def _adaptive(a, b, S_ab, tol, depth):
        m = (a + b) / 2
        S_am = simpson_enkel(f, a, m)
        S_mb = simpson_enkel(f, m, b)
        S_sum = S_am + S_mb

        # Sjekk om vi er nøyaktige nok
        if depth >= max_depth or abs(S_sum - S_ab) < 15 * tol:
            return S_sum + (S_sum - S_ab) / 15  # Richardson-ekstrapolering

        # Rekursivt kall på hver halvdel
        return (_adaptive(a, m, S_am, tol/2, depth+1) +
                _adaptive(m, b, S_mb, tol/2, depth+1))

    S_ab = simpson_enkel(f, a, b)
    return _adaptive(a, b, S_ab, tol, 0)

# Test på funksjon med skarp topp
def f(x):
    return 1 / ((x - 0.3)**2 + 0.01)

# Beregn integralet
resultat = adaptiv_simpson(f, 0, 1)

# Sammenlign med SciPy
from scipy import integrate
scipy_resultat, _ = integrate.quad(f, 0, 1)

print(f"Adaptiv Simpson: {resultat:.10f}")
print(f"SciPy quad:      {scipy_resultat:.10f}")
print(f"Forskjell:       {abs(resultat - scipy_resultat):.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$ .

import numpy as np

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

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

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

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

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

        # 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[max_iter-1, max_iter-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.

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.fixed_quad(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 leggauss

def gauss_legendre(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]
    x_trans = (b - a) / 2 * x + (a + b) / 2
    scale = (b - a) / 2

    return scale * np.sum(w * f(x_trans))

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]$ .

import numpy as np
from scipy import integrate

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

# 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 inner_integral(x):
    """Integrerer f over y for fast x"""
    result, _ = integrate.quad(lambda y: np.exp(-(x**2 + y**2)), 0, 1)
    return result

resultat_iterert, _ = integrate.quad(inner_integral, 0, 1)
print(f"Iterert:         {resultat_iterert:.10f}")

# Manuell 2D Simpson
def simpson_2d(f, ax, bx, ay, by, nx, ny):
    """2D Simpson-integrasjon"""
    hx = (bx - ax) / nx
    hy = (by - ay) / ny

    def simpson_1d(g, a, b, n):
        h = (b - a) / n
        return (h/3) * (g(a) + 4*sum(g(a + i*h) for i in range(1, n, 2))
                        + 2*sum(g(a + i*h) for i in range(2, n, 2)) + g(b))

    def inner(x):
        return simpson_1d(lambda y: f(x, y), ay, by, ny)

    return simpson_1d(inner, ax, bx, nx)

def f_xy(x, y):
    return np.exp(-(x**2 + y**2))

resultat_simpson = simpson_2d(f_xy, 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$ .

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 y_lower(x):
    return -np.sqrt(1 - x**2)

def y_upper(x):
    return np.sqrt(1 - x**2)

# Beregn dobbeltintegralet
areal, feil = integrate.dblquad(f, -1, 1, y_lower, y_upper)

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(x**2 + y**2 <= 1)
areal_mc = 4 * (inside / N)  # 4 = areal av kvadratet [-1,1] x [-1,1]

print(f"\nMonte Carlo:    {areal_mc:.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$ .

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)
arbeid_linear = 0.5 * k * L**2

print(f"Ikke-lineær fjær:")
print(f"  Arbeid = {arbeid:.4f} J")
print(f"\nLineær fjær:")
print(f"  Arbeid = {arbeid_linear:.4f} J")
print(f"\nForskjell: {(arbeid - arbeid_linear) / arbeid_linear * 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.fill_between(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
arbeid_kumulativ = [integrate.quad(F, 0, xi)[0] for xi in x]
plt.plot(x, arbeid_kumulativ, 'b-', label='Arbeid')
plt.xlabel('x (m)')
plt.ylabel('W (J)')
plt.title('Kumulativt arbeid')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

✏️Beregne buelengde

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

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, 2*np.pi)

print(f"Buelengde av sin(x) fra 0 til 2*pi:")
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 (2*pi,0) er 2*pi
print(f"\nRett linje: {2*np.pi:.6f}")
print(f"Forlengelse: {(buelengde - 2*np.pi) / (2*np.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.fill_between(x, 0, y, where=(y > 0), alpha=0.2, color='blue')
plt.fill_between(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$ .

Oppsummering

I dette kapittelet har du lært:

- Rektangelmetoden: Tilnærmer integralet med $\sum f(x_i)\Delta x$ ; venstrepunkt, høyrepunkt eller midtpunkt kan brukes.
- Trapesmetoden: $\displaystyle T = \frac{\Delta x}{2}\left(f(x_0) + 2f(x_1) + \cdots + 2f(x_{n-1}) + f(x_n)\right)$ , med feil av orden $\displaystyle \frac{1}{n^2}$ .
- Simpsons metode: Bruker parabler gjennom tre og tre punkter, med feil av orden $\displaystyle \frac{1}{n^4}$ — eksakt for polynomer av grad $\leq 3$ .
- Python-implementasjon: Alle metodene kan skrives som korte løkker eller med NumPy/SciPy (scipy.integrate).
- Monte Carlo-integrasjon: Bruker tilfeldige punkter til å estimere integraler, nyttig i høye dimensjoner.

Nøkkelbegreper

Begrep	Forklaring
Rektangelmetoden	Sum av rektangelarealer
Trapesmetoden	Sum av trapesarealer, feil $\sim 1/n^2$
Simpsons metode	Parabeltilnærming, feil $\sim 1/n^4$
Monte Carlo	Estimat basert på tilfeldige punkter

Viktige formler

\displaystyle \Delta x = \frac{b-a}{n}

- Trapes:

\displaystyle T = \frac{\Delta x}{2}\left(f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)\right)

- Simpson:

\displaystyle S = \frac{\Delta x}{3}\left(f(x_0) + 4f(x_1) + 2f(x_2) + \cdots + 4f(x_{n-1}) + f(x_n)\right)

Repetisjonsoppgaver

Din fremgang

0deloppgaver0 / 5 oppgaver

Dette kapitlet er skrevet av Anthropics toppmodeller (Claude Opus og Claude Fable) og er foreløpig ikke manuelt gjennomgått — kvalitetskontrollen gjøres av uavhengige KI-agenter, og innmeldte feil rettes fortløpende. Funnet en feil? Meld fra, så retter vi den. Les mer om hvordan innholdet lages.