SciPy Integration

SciPy Integration#

The SciPy library provides a lot of methods for computing integrals:

from scipy import integrate
import numpy as np

For a function where you know \(f(x)\), quad is a good general purpose integrator.

Let’s compute

\[I = \int_0^{2\pi} \sin^2(x) dx\]

I, err = integrate.quad(lambda x: np.sin(x)**2, 0.0, 2.0 * np.pi)

print(I)

3.141592653589793

print(err)

2.3058791671639882e-09

Notice that it gives us the answer as well as an estimate of the error.

Sometimes our integrand takes arguments. Let’s integrate

\[I = \int_{-1}^1 A e^{-(x/\sigma)^2} dx\]

def integrand(x, A, sigma):
    return A * np.exp(-x**2 / sigma**2)

I, err = integrate.quad(integrand, -1.0, 1.0, args=(1.0, 2.0))
print(I, err)

1.8451240256511698 2.0484991765669867e-14

NumPy defines np.inf which can be used for integrating to infinity. So to compute:

\[I = \int_{-\infty}^{\infty} A e^{-(x/\sigma)^2} dx\]

we do

I, err = integrate.quad(integrand, -np.inf, np.inf,
                        args=(1.0, 1.0))
print(I, err)

1.7724538509055159 1.4202636780944923e-08

For a function where we only have the value at some samples \(\{ f_i\}\), we can use Simpson’s rule

Let’s compute

\[I = \int_0^{2\pi} \sin^2(x) dx\]

with \((x_i, f_i)\) defined at \(N\) points.

N = 17
x = np.linspace(0.0, 2.0*np.pi, N, endpoint=True)
y = np.sin(x)**2

I = integrate.simpson(y, x=x)
print(I)

3.141592653589793