SciPy Integration

SciPy Integration#

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

https://docs.scipy.org/doc/scipy/tutorial/integrate.html

from scipy import integrate
import numpy as np

Analytically defined function#

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

Pointwise defined function#

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