Forelesning 6: Numerisk integrasjon

Det bestemte integralet

def f(x):
    return x    # Integrert = 0.5x**2

a = 0 # Startpunkt 
b = 2 # Sluttpunkt
A = 0 # Areal
N = 100000000 # Antall rektangler
dx = (b - a)/N # Bredden på rektanglene

x = a
for i in range(N):
    A = A + f(x)*dx # A = høyde x bredde
    x = x + dx

print(A)

1.9999999818801717

def f(x):
    return x

def integrate(f, a, b, N):
    dx = (b - a) / N  # Bredden på rektanglene
    A = 0  # Areal
    x = a
    
    for i in range(N):
        A += f(x) * dx  # A = høyde * bredde
        x += dx
    
    return A
    
resultat = integrate(f, 0, 2, 1000)
print(resultat)

def venstretilnærming(f, a, b, N = 100000):
    dx = (b - a)/N 
    x = a
    A = 0
    for i in range(N):
        A = A + f(x)*dx
        x = x + dx
    return A

def høyretilnærming(f, a, b, N = 100000):
    dx = (b - a)/N 
    x = a + dx
    A = 0
    for i in range(N):
        A = A + f(x)*dx
        x = x + dx
    return A

høyretilnærming(f, 0, 2)

2.0000199999987003

def midtpunktstilnærming(f, a, b, N = 100000):
    dx = (b - a)/N 
    x = a + dx/2
    A = 0
    for i in range(N):
        A = A + f(x)
        x = x + dx
    return A*dx

midtpunktstilnærming(f, 0, 2)

1.9999999999987004

Andre tilnærminger

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

trapesmetoden(f, 0, 2)

1.9999999999987024

Bruke biblioteksfunksjoner til å integrere

Benytter scipy-biblioteket.

from scipy.integrate import trapezoid, simpson, quad
import numpy as np

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

trapesmetoden = trapezoid(y,x)
trapesmetoden

2.0

simpsons = simpson(y,x)
simpsons

/var/folders/z_/zd2_19g1205dvcvdhgk10p680000gp/T/ipykernel_88515/2046446227.py:1: DeprecationWarning: You are passing x=[0.00000e+00 2.00002e-05 4.00004e-05 ... 1.99996e+00 1.99998e+00
 2.00000e+00] as a positional argument. Please change your invocation to use keyword arguments. From SciPy 1.14, passing these as positional arguments will result in an error.
  simpsons = simpson(y,x)

2.0000000000000004

quad(f, a, b)

(2.0, 2.220446049250313e-14)