Forelesning 6: Numerisk integrasjon

Det bestemte integralet

def f(x):
    return x**2

a = 0  # Nedre integrasjonsgrense
b = 2  # Øvre integrasjonsgrende
n = 10000 # Antall rektangler
dx = (b - a)/n # Bredden av hvert rektangel

A = 0 # Arealet av rektanglene
x = a

for i in range(n):
    # Summere opp arealet av hvert rektangel
    A = A + f(x) * dx
    x = x + dx

print(A)
    

2.666266679999545

def rektangelmetoden_venstre(f, a, b, n = 100000):
    dx = (b - a)/n # Bredden av hvert rektangel
    A = 0 # Arealet av rektanglene
    x = a
    for i in range(n):
        # Summere opp arealet av hvert rektangel
        A = A + f(x) * dx
        x = x + dx
    return A

def rektangelmetoden_høyre(f, a, b, n = 100000):
    dx = (b - a)/n # Bredden av hvert rektangel
    A = 0 # Arealet av rektanglene
    x = a + dx
    for i in range(n):
        # Summere opp arealet av hvert rektangel
        A = A + f(x) * dx
        x = x + dx
    return A

def rektangelmetoden_midtpunkt(f, a, b, n = 100000):
    dx = (b - a)/n # Bredden av hvert rektangel
    A = 0 # Arealet av rektanglene
    x = a + dx/2
    for i in range(n):
        # Summere opp arealet av hvert rektangel
        A = A + f(x) * dx
        x = x + dx
    return A

print(rektangelmetoden_venstre(f, 1, 2))
print(rektangelmetoden_høyre(f, 1, 2))
print(rektangelmetoden_midtpunkt(f, 1, 2))

2.333318333360897
2.3333483333608975
2.3333333333359163

def trapesmetoden(f, a, b, n = 100000):
    dx = (b - a)/n # Bredden av hvert rektangel
    A = (f(a) + f(b))/2*dx # Arealet av rektanglene
    x = a
    for i in range(1, n):
        # Summere opp arealet av hvert rektangel
        A = A + f(x) * dx
        x = x + dx
    return A

trapesmetoden(f, 1, 2)

2.3333033337608957

Bruke biblioteksfunksjoner til å integrere

Benytter scipy-biblioteket.

from scipy.integrate import trapezoid, simpson, quad

x1 = np.linspace(0,1,10000)
y1 = f(x)

trapes = trapezoid(y = y1, x = x1)
trapes

0.3333333350003334

simps = simpson(y = y1, x = x1)
simps

0.3333333333333333

kvadratur = quad(f, 0, 1)
kvadratur

(0.3333333333333333, 3.700743415417188e-15)