Forelesning 6: Numerisk integrasjon

Det bestemte integralet

def f(x):
    return 3*x**2

def analytisk_integrert(x):
    return x**3

N = 10000
a = 0
b = 3

A = 0
dx = (b - a)/N

for i in range(N):
    A = A + f(a + i*dx)*dx
    
print(A)
print(analytisk_integrert(b) - analytisk_integrert(a))

26.995950134999916
27

def rektangelmetoden_venstre(f, a, b, N=10000):
    A = 0
    dx = (b - a)/N
    for i in range(N):
        A = A + f(a + )(i*dx)
    return A*dx

rektangelmetoden_venstre(f, 0, 3)

26.99595013499999

def rektangelmetoden_høyre(f, a, b, N=10000):
    A = 0
    dx = (b - a)/N
    for i in range(N):
        A = A + f(a + (i+1)*dx)
    return A*dx

rektangelmetoden_høyre(f, 0, 3)

26.99595013499999

def rektangelmetoden_midtpunkt(f, a, b, N=10000):
    A = 0
    dx = (b - a)/N
    for i in range(N):
        A = A + f(a + (1+2*i)*(dx/2))*dx
    return A

rektangelmetoden_midtpunkt(f, a, b)

26.999999932499964

Andre tilnærminger

def trapesmetoden(f, a, b, N=10000):
    dx = (b - a)/N
    A = (f(a) + f(b))/2
    for i in range(1,N):
        A = A + f(a + i*dx)
    return A*dx

trapesmetoden(f, 0, 3)

27.000000134999997

Bruke biblioteksfunksjoner til å integrere

Benytter scipy-biblioteket.

from scipy import integrate # eventuelt from scipy.integrate import trapezoid
import numpy as np

def f(x):
    return 3*x**2

a = 0
b = 3
x = np.linspace(a, b, 100000)
y = f(x)

trapesmetoden = integrate.trapezoid(y, x)
simpsons = integrate.simpson(y, x)
gauss = integrate.quad(f, a, b)

print(f"trapesmetoden: {trapesmetoden}, simpsons metode: {simpsons}, gauss: {gauss}")

trapesmetoden: 27.000000001350028, simpsons metode: 27.00000000000001, gauss: (27.0, 2.9976021664879227e-13)

from sympy import *

x = symbols("x")
y = x**2 - cos(x/2) + 1/exp(2*x)

integrate(y,x)

\[\displaystyle \frac{x^{3}}{3} - 2 \sin{\left(\frac{x}{2} \right)} - \frac{e^{- 2 x}}{2}\]