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 å rekne ut bestemte integral numerisk. Sjølv om mange integral kan løysast analytisk, finst det utallige funksjonar som ikkje har ein lukka antiderivert. Då må vi ty til numeriske metodar.

Vi skal implementere tre klassiske metodar:
- Rektangelmetoden (den enklaste)
- Trapesmetoden (betre nøyaktigheit)
- Simpsons metode (svært god nøyaktigheit)

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

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

Rektangelmetoden

Den enklaste tilnærminga til numerisk integrasjon er å dele opp arealet under kurva i rektangel. Idéen er at summen av areala til rektangla tilnærmar det bestemte integralet.

Rektangelmetoden (venstrepunkt)

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

✏️Implementere rektangelmetoden i Python

Skriv ein Python-funksjon som reknar ut $\displaystyle\int_0^1 x^2\,dx$ ved hjelp av rektangelmetoden med $n = 1000$ delintervall. 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$ , noko som svarar til ein 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)

✏️Samanlikning av venstre- og midtpunktsmetoden

Samanlikn 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 meir nøyaktig!

📝Oppgave 3.5.1

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

Den eksakte verdien er $2$ .

a) Bruk $n = 100$ delintervall og finn tilnærma verdi.

b) Bruk $n = 1000$ delintervall og samanlikn feilen.

c) Kor mange delintervall treng du for at feilen skal vere mindre enn $0.0001$ ?

Oppgave 3.5.1

Medium

PythonAuto-lagret

Trapesmetoden

Trapesmetoden gir betre nøyaktigheit enn rektangelmetoden ved å tilnærme kurva med trapes i staden for rektangel. Eit trapes fangar opp meir av forma til kurva.

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 ein eller annan

\xi \in [a, b]

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

✏️Implementere trapesmetoden i Python

Skriv ein Python-funksjon for trapesmetoden og rekn ut $\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 blir redusert med faktor 100 når $n$ aukar med faktor 10, som stadfestar $O(1/n^2)$ -åtferda.

📝Oppgave 3.5.2

Samanlikn 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 metodane.

b) Rekn ut integralet med $n = 100$ for begge metodane.

c) Kor mange gonger meir nøyaktig er trapesmetoden?

Oppgave 3.5.2

Medium

PythonAuto-lagret

Simpsons metode

Simpsons metode tilnærmar kurva med parablar (andregradsfunksjonar) i staden for rette linjer. Dette gir endå betre nøyaktigheit enn trapesmetoden, særleg for glatte funksjonar.

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 ein eller annan

\xi \in [a, b]

Dette betyr at feilen avtek som $O(1/n^4)$ , altså mykje raskare enn trapesmetoden.

Bemerkelsesverdig: Simpsons metode er eksakt for polynom av grad 3 eller lågare!

✏️Implementere Simpsons metode i Python

Implementer Simpsons metode og rekn ut $\displaystyle\int_0^{\pi} \sin(x)\,dx$ med $n = 10$ . Samanlikn 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 berre $n = 10$ har Simpsons metode ein feil på ca. $7 \cdot 10^{-6}$ , medan trapesmetoden har feil på ca. $0.017$ . Simpsons metode er over 2000 gonger meir nøyaktig!

📝Oppgave 3.5.3

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

Dette integralet har ingen lukka antiderivert, så vi må bruke numeriske metodar!

a) Implementer Simpsons metode.

b) Finn ein tilnærma verdi med $n = 100$ .

c) Den "eksakte" verdien (rekna ut med mange desimalar) er ca. $0.746824$ . Kva er feilen?

Oppgave 3.5.3

Vanskelig

PythonAuto-lagret

Simpsons 3/8-regel

Det finst òg ein variant av Simpsons metode som bruker kubisk interpolasjon over tre delintervall.

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 treng vi ikkje alltid å implementere integrasjonsmetodar frå botnen av. Python-biblioteka NumPy og SciPy har ferdige, optimaliserte funksjonar for numerisk integrasjon.

SciPy integrate-modulen

n

✏️Integrasjon med SciPy

Bruk SciPy til å rekne ut $\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 handterer uendelege grenser automatisk og gir eit feilestimat!

✏️Trapesmetoden med NumPy på datapunkt

Gitt målepunkt $(x_i, y_i)$ frå eit eksperiment, rekn ut det tilnærma arealet under kurva.

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 staden for ein analytisk funksjon.

📝Oppgave 3.5.4

Bruk SciPy til å rekne ut følgjande integral:

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

Samanlikn med dei eksakte verdiane.

Oppgave 3.5.4

Medium

PythonAuto-lagret

Visualisering av numerisk integrasjon

Ein viktig del av å forstå numeriske metodar er å visualisere kva som skjer. Med Matplotlib kan vi teikne funksjonen og rektangla/trapesa som blir brukte til å tilnærme arealet.

✏️Visualisere rektangelmetoden

Lag ei visualisering som viser rektangla som blir brukte i rektangelmetoden for $f(x) = x^2$ på $[0, 2]$ med $n = 8$ delintervall.

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 ein figur som viser:
- Den blå kurva $f(x) = x^2$
- Lyseblå rektangel som representerer tilnærminga
- Tittel med både tilnærma og eksakt verdi

✏️Visualisere trapesmetoden

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

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 ei visualisering som samanliknar rektangelmetoden, trapesmetoden og Simpsons metode for $f(x) = x^3$ på $[0, 1]$ med $n = 4$ .

Vis tre subplot side om side som illustrerer kvar metode, og rekn ut feilen for kvar.

Oppgave 3.5.5

Vanskelig

PythonAuto-lagret

Monte Carlo-integrasjon

Ei heilt anna tilnærming til numerisk integrasjon er Monte Carlo-metoden, som bruker tilfeldige tal. Idéen er enkel: kast tilfeldige punkt i eit rektangel og tel kor mange som landar under kurva.

Monte Carlo-integrasjon

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

✏️Monte Carlo-integrasjon i Python

Bruk Monte Carlo-metoden til å estimere $\pi$ ved å rekne ut $\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 avtek som $1/\sqrt{N}$ : for å halvere feilen må vi firedoble talet på punkt.

📝Oppgave 3.5.6

Bruk Monte Carlo-integrasjon til å rekne ut volumet av ei einingssfære i 3D.

Volumet av ei 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 punkt i kuben $[-1, 1]^3$ og tel kor mange som er innanfor 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 metodar, er det viktig å forstå kor nøyaktig resultatet er og korleis feilen oppfører seg når vi aukar talet på delintervall.

📜Konvergensordenar for integrasjonsmetodar

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})$	Avheng 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

Rekn ut konvergensordenen til trapesmetoden for $\displaystyle\int_0^1 e^x\,dx$ ved å samanlikne 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 stadfestar at trapesmetoden har konvergensorden $O(1/n^2)$ .

📝Oppgave 3.5.7

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

a) Rekn ut feilen for $n = 10, 20, 40, 80, 160$ for alle tre metodane.

b) Lag eit log-log-plott av feil mot $n$ for alle tre metodane.

c) Stadfest konvergensordenane frå plottet (helninga skal vere -1, -2, -4 for dei respektive metodane).

Oppgave 3.5.7

Vanskelig

PythonAuto-lagret

Adaptiv integrasjon

Så langt har vi brukt jamt fordelte delintervall. Men nokre funksjonar varierer meir i nokre område enn andre. Adaptiv integrasjon justerer steglengda automatisk: små steg der funksjonen varierer mykje, store steg der han er jamn.

Adaptiv Simpson-integrasjon

S_1 = S(a, b)

✏️Implementere adaptiv Simpson

Implementer adaptiv Simpson-integrasjon for funksjonar som varierer mykje 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 utrekningane rundt toppen ved $x = 0.3$ .

📝Oppgave 3.5.8

Samanlikn talet på funksjonsevalueringar for vanleg Simpson og adaptiv Simpson for funksjonen

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

på intervallet $[0, 1]$ .

a) Implementer ein teljar for funksjonsevalueringar.

b) Finn kor mange evalueringar vanleg Simpson treng for å oppnå feil $< 10^{-8}$ .

c) Tel evalueringar 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øgare nøyaktigheit. Idéen er å bruke resultata frå trapesmetoden med forskjellige $n$ til å "ekstrapolere bort" feila.

Romberg-algoritmen

R_{i,j}

✏️Implementere Romberg-integrasjon

Implementer Romberg-integrasjon og rekn ut $\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 berre 5-6 iterasjonar har vi maskin-presisjon!

📝Oppgave 3.5.9

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

a) Samanlikn med scipy.integrate.quad.

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

Oppgave 3.5.9

Vanskelig

PythonAuto-lagret

Gauss-kvadratur

Gauss-kvadratur er ein familie av integrasjonsmetodar som oppnår maksimal nøyaktigheit for eit gitt tal på funksjonsevalueringar. I staden for å bruke jamt fordelte punkt, bruker Gauss-metodane optimalt plasserte punkt 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 å rekne ut $\displaystyle\int_0^1 e^{-x^2}\,dx$ med forskjellig tal på punkt.

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 berre 10 punkt oppnår vi maskin-presisjon!

📝Oppgave 3.5.10

Samanlikn Gauss-Legendre-kvadratur med Simpsons metode for følgjande integral:

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$ (oscillerande)

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

For kvar: finn talet på punkt/delintervall som gir feil $< 10^{-10}$ .

Oppgave 3.5.10

Vanskelig

PythonAuto-lagret

Dobbeltintegral numerisk

Numerisk integrasjon i to dimensjonar er ei naturleg utviding av metodane vi har sett. Vi kan rekne ut

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

ved å gjenta éin-dimensjonal integrasjon.

Dobbeltintegral med iterert integrasjon

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

✏️Rekne ut dobbeltintegral med SciPy

Rekn ut $\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 ikkje-rektangulært område

Rekn ut arealet av einingssirkelen ved å rekne ut $\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

Rekn ut 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$ punkt.

c) Visualiser paraboloiden og integrasjonsområdet.

Oppgave 3.5.11

Vanskelig

PythonAuto-lagret

Praktiske bruksområde

Numerisk integrasjon har utallige bruksområde i vitskap og ingeniørfag. Her ser vi på nokre eksempel.

✏️Rekne ut arbeid utført av ei variabel kraft

Ei fjør følgjer Hookes lov med $F(x) = kx$ for små deformasjonar, men for store deformasjonar blir ho ikkje-lineær: $F(x) = kx + \alpha x^3$ . Rekn ut arbeidet som krevst for å strekkje fjøra frå $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()

✏️Rekne ut buelengd

Rekn ut buelengda til kurva $y = \sin(x)$ frå $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

Fysikkoppgåve: Ein partikkel rører seg med farten

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

m/s.

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

b) Finn netto forflytning (endeleg 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ølgjande numeriske integrasjonsmetodar:

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 dei fleste 1D-oppgåver
2. Bruk scipy.integrate.dblquad for 2D-integral
3. Bruk Monte Carlo for dimensjon

\geq 3

4. Bruk Romberg eller Gauss for høg presisjon
5. Verifiser alltid med kjende integral eller ved å auke

n

6. Visualiser for å byggje intuisjon!

📝Oppgave 3.5.13

Prosjektoppgåve: Normalfordelinga

Sannsynstettleiken for normalfordelinga (Gauss-fordelinga) 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) Rekn ut sannsynet $P(-1 \leq X \leq 1)$ , altså arealet under kurva frå $-1$ til $1$ . Dette skal vere ca. $0.6827$ (68.27%).

c) Rekn ut $P(-2 \leq X \leq 2)$ og $P(-3 \leq X \leq 3)$ . Samanlikn med "68-95-99.7-regelen".

d) Visualiser normalfordelinga og skuggelegg området frå $-1$ til $1$ .

Oppgave 3.5.13

Vanskelig

PythonAuto-lagret

📝Oppgave 3.5.14

Bonusoppgåve: Lag ei interaktiv samanlikning

Lag eit program som:

a) Tek inn ein funksjon $f(x)$ , grenser $[a, b]$ , og tal på punkt $n$ .

b) Reknar ut integralet med alle metodane: rektangel, midtpunkt, trapes, Simpson, Romberg og Gauss-Legendre.

c) Viser ein tabell med resultata og feila (bruk scipy.quad som referanse).

d) Rangerer metodane etter nøyaktigheit for den gitte funksjonen.

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

Oppgave 3.5.14

Vanskelig

PythonAuto-lagret

Vanlege feil og fallgruver

1. Singularitetar:
Dersom integranden har singularitetar (uendelegheiter), vil standardmetodane feile. Bruk spesialiserte metodar eller del opp intervallet.

2. Oscillerande funksjonar:
Funksjonar som $\sin(1000x)$ krev svært mange delintervall. Bruk adaptiv integrasjon.

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

4. Gløyme partalskravet:
Simpsons metode krev at $n$ er eit partal.

5. Avrundingsfeil:
For svært små $h$ kan avrundingsfeil dominere. Bruk metodar av høgare orden i staden for å auke $n$ .

Oppsummering

I dette kapittelet har du lært:

- Rektangelmetoden: Tilnærmar integralet med $\sum f(x_i)\Delta x$ ; venstrepunkt, høgrepunkt eller midtpunkt kan brukast.
- 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 parablar gjennom tre og tre punkt, med feil av orden $\displaystyle \frac{1}{n^4}$ — eksakt for polynom av grad $\leq 3$ .
- Python-implementasjon: Alle metodane kan skrivast som korte løkker eller med NumPy/SciPy (scipy.integrate).
- Monte Carlo-integrasjon: Bruker tilfeldige punkt til å estimere integral, nyttig i høge dimensjonar.

Nøkkelomgrep

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

Viktige formlar

\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)

Repetisjonsoppgåver

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.