3.5 Programmering av integrasjon

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(
    "Tilnærmet verdi:"
    f" {resultat:.6f}"
)
print(f"Eksakt verdi: {1/3:.6f}")
print(
    "Feil:"
    f" {abs(resultat - 1/3):.6f}"
)

Utskrift:

Tilnærmet verdi: 0.332833
Eksakt verdi: 0.333333
Feil: 0.000500

Feilen er ca. $0{,}0005$, noe som tilsvarer en relativ feil på ca. $0{,}15\,\%$.

Svar: Rektangelmetoden gir $\approx 0{,}332833$, mens den eksakte verdien er $\displaystyle \frac{1}{3} \approx 0{,}333333$. Feilen er ca. $0{,}0005$.

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"
        f"{abs(v - eksakt):.2e}\t"
        f"{abs(m - eksakt):.2e}"
    )

Utskrift:

n		Venstrepunkt	Midtpunkt
--------------------------------------------------
10		4.83e-02	8.33e-04
100		4.98e-03	8.33e-06
1000		5.00e-04	8.33e-08

Midtpunktsmetoden er dramatisk mer nøyaktig!

Svar: Med $n = 1000$ er feilen ca. $5{,}0 \cdot 10^{-4}$ med venstrepunktsmetoden og ca. $8{,}3 \cdot 10^{-8}$ med midtpunktsmetoden.

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}:"
        f" {resultat:.8f}, feil ="
        f" {feil:.2e}"
    )

Utskrift:

n =   10: 1.71971349, feil = 1.43e-03
n =  100: 1.71829615, feil = 1.43e-05
n = 1000: 1.71828197, feil = 1.43e-07

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

Svar: Trapesmetoden gir $1{,}71971349$ med $n = 10$ (feil $1{,}43 \cdot 10^{-3}$) og $1{,}71828197$ med $n = 1000$ (feil $1{,}43 \cdot 10^{-7}$).

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

Utskrift:

Eksakt verdi: 2
Trapes (n=10):  1.9835235375, feil = 1.65e-02
Simpson (n=10): 2.0001095173, feil = 1.10e-04

Med bare $n = 10$ har Simpsons metode en feil på ca. $1{,}1 \cdot 10^{-4}$, mens trapesmetoden har feil på ca. $0{,}017$. Simpsons metode er ca. 150 ganger mer nøyaktig.

Svar: Med $n = 10$ gir Simpsons metode $2{,}00011$ (feil ca. $1{,}1 \cdot 10^{-4}$) og trapesmetoden $1{,}98352$ (feil ca. $0{,}017$).

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(
    "SciPy quad:  "
    f" {resultat:.10f}"
)
print(
    f"Eksakt verdi: {eksakt:.10f}"
)
print(
    "Feilestimering fra quad:"
    f" {feil:.2e}"
)
print(
    "Faktisk feil:"
    f" {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!

Svar: quad gir $0{,}8862269255$, som er lik $\displaystyle \frac{\sqrt{\pi}}{2}$ med ti desimaler.

import numpy as np

# Simulerte måledata (f.eks.
# hastighet over tid)
# tid i sekunder
t = np.array([0, 1, 2, 3, 4, 5])
# hastighet i m/s
v = np.array([0, 5, 8, 10, 9, 6])

# Beregn tilbakelagt strekning
# (integral av hastighet)
strekning = np.trapezoid(v, t)

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

Utskrift:

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

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

Svar: Tilbakelagt strekning er ca. $35{,}0$ m.

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(
    'Rektangelmetoden: n ='
    f' {n}\nTilnærmet:'
    f' {integral_approx:.4f},'
    ' Eksakt:'
    f' {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

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(
    'Trapesmetoden: n ='
    f' {n}\nTilnærmet:'
    f' {integral_approx:.4f},'
    ' Eksakt: 2.0000'
)
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

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
# For reproduserbarhet
np.random.seed(42)

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

Utskrift:

N =     100: pi ≈ 3.094681, feil = 0.046912
N =    1000: pi ≈ 3.114703, feil = 0.026890
N =   10000: pi ≈ 3.153273, feil = 0.011680
N =  100000: pi ≈ 3.144075, feil = 0.002482
N = 1000000: pi ≈ 3.141744, feil = 0.000151

Feilen avtar i gjennomsnitt som $1/\sqrt{N}$: for å halvere feilen må vi firedoble antall punkter. Én enkelt kjøring følger ikke dette nøyaktig, fordi punktene er tilfeldige.

Svar: Med $N = 1\,000\,000$ punkter gir programmet estimatet $3{,}141744$ for $\pi$, med feil ca. $0{,}00015$.

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"
            f"{faktor:.2f}"
        )
    else:
        print(
            f"{n}\t\t{feil:.2e}\t"
            "-"
        )

    prev_feil = feil

Utskrift:

n		Feil		Faktor
----------------------------------------
10		1.43e-03	-
20		3.58e-04	4.00
40		8.95e-05	4.00
80		2.24e-05	4.00
160		5.59e-06	4.00
320		1.40e-06	4.00

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

Svar: Feilen blir 4 ganger mindre hver gang $n$ dobles, så trapesmetoden har konvergensorden 2.

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
        ):
            # Richardson-
            # ekstrapolering
            return (
                S_sum
                + (S_sum - S_ab)
                / 15
            )

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

Utskrift:

Adaptiv Simpson: 26.7794504459
SciPy quad:      26.7794504459
Forskjell:       1.82e-12

Den adaptive metoden konsentrerer beregningene rundt toppen ved $x = 0{,}3$.

Svar: Integralet er $\approx 26{,}7794504459$. Adaptiv Simpson og quad gir samme verdi (forskjell ca. $2 \cdot 10^{-12}$).

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(
    "Feil:    "
    f" {abs(resultat - eksakt):.2e}"
)

Utskrift:

Romberg-tabell:
----------------------------------------------------------------------
i=0: 1.8591409142
i=1: 1.7539310925  1.7188611519
i=2: 1.7272219046  1.7183188419  1.7182826879
i=3: 1.7205185922  1.7182841547  1.7182818422  1.7182818288
i=4: 1.7188411286  1.7182819741  1.7182818287  1.7182818285  1.7182818285
i=5: 1.7184216603  1.7182818376  1.7182818285  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!

Svar: Romberg gir $1{,}718281828459 = e - 1$ med feil ca. $4 \cdot 10^{-16}$.

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}:"
        f" {resultat:.12f}, feil"
        f" = {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]
    halvbredde = (b - a) / 2
    sentrum = (a + b) / 2
    x_trans = halvbredde * x + sentrum
    scale = halvbredde

    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}:"
        f" {resultat:.12f}, feil"
        " ="
        f" {abs(resultat - eksakt):.2e}"
    )

Utskrift:

Gauss-Legendre-kvadratur:
--------------------------------------------------
n =  2: 0.746594688283, feil = 2.29e-04
n =  3: 0.746814584191, feil = 9.55e-06
n =  4: 0.746824468131, feil = 3.35e-07
n =  5: 0.746824126766, feil = 6.05e-09
n =  6: 0.746824132890, feil = 7.77e-11
n =  8: 0.746824132812, feil = 6.44e-15
n = 10: 0.746824132812, feil = 1.11e-16

Manuell Gauss-Legendre med numpy.polynomial.legendre:
n = 3: 0.746814584191, feil = 9.55e-06
n = 5: 0.746824126766, feil = 6.05e-09
n = 7: 0.746824132812, feil = 7.89e-13

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

Svar: Gauss-Legendre med 10 punkter gir $0{,}746824132812$, med feil ca. $10^{-16}$.

import numpy as np
from scipy import integrate


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


# Beregn dobbeltintegralet
resultat, feil = (
    integrate.dblquad(
        f, 0, 1, 0, 1
    )
)

print(
    "Dobbeltintegral:"
    f" {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(
    "Iterert:        "
    f" {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(
    "Simpson 2D:     "
    f" {resultat_simpson:.10f}"
)

Utskrift:

Dobbeltintegral: 0.5577462854
Feilestimering:  8.29e-15
Iterert:         0.5577462854
Simpson 2D:      0.5577463616

Svar: Dobbeltintegralet er $\approx 0{,}5577462854$.

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(
    "Beregnet areal:"
    f" {areal:.10f}"
)
print(
    "Eksakt (pi):   "
    f" {np.pi:.10f}"
)
print(
    "Feil:          "
    f" {abs(areal - np.pi):.2e}"
)

# Alternativ: Monte Carlo
# For reproduserbarhet
np.random.seed(42)
N = 1000000
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
inside = np.sum(x**2 + y**2 <= 1)
# 4 = areal av kvadratet [-1,1] x
# [-1,1]
areal_mc = 4 * (inside / N)

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

Utskrift:

Beregnet areal: 3.1415926536
Eksakt (pi):    3.1415926536
Feil:           5.33e-15

Monte Carlo:    3.1399
Feil MC:        0.0017

Svar: dblquad gir arealet $3{,}1415926536 = \pi$. Monte Carlo med $10^6$ punkter gir $3{,}1399$ (feil ca. $0{,}0017$).

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

# Parametre
k = 100  # N/m (fjærstivhet)
# N/m^3 (ikke-lineær koeffisient)
alpha = 500
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(
    "  Arbeid ="
    f" {arbeid_linear:.4f} J"
)
forskjell = (
    (arbeid - arbeid_linear)
    / arbeid_linear
    * 100
)
print(f"\nForskjell: {forskjell:.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=(
        r'Ikke-lineær: $F = kx +'
        r' \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()

Svar: Arbeidet er $0{,}5125$ J for den ikke-lineære fjæra og $0{,}5000$ J for den lineære, en forskjell på $2{,}5\,\%$.

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"
    f" til 2*pi:"
)
print(f"L = {buelengde:.6f}")
print(
    "(Dette er ca."
    f" {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}"
)
forlengelse = (
    (buelengde - 2 * np.pi)
    / (2 * np.pi)
    * 100
)
print(f"Forlengelse: {forlengelse:.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(
    '$y = \sin(x)$, buelengde $L'
    f' = {buelengde:.4f}$'
)
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.show()

Svar: Buelengden er $L \approx 7{,}640396$.