SciPy Integration

When a function cannot be integrated analytically, or is difficult to integrate analytically, people often resort to numerical integration methods. SciPy includes many routines for performing numerical integration. Most of them are found in the same scipy.integration library. The following table lists some commonly used functions.

Single Integral

The quad function is the workhorse of SciPy’s integral functions. Numerical integration is sometimes called orthogonal, hence the name “orthogonal.” It is often the default choice for performing a single integral of a function f ( x ) over a given fixed range from a to b.

int_{a}^{b} f(x)dx

The general form of quad is scipy.integration.quad(f, a, b) , where ‘f’ is the name of the function to be integrated, and ‘a’ and ‘b’ are the lower and upper limits, respectively. Let’s look at an example of a Gaussian function, integrated between 0 and 1.

We first need to define the function →f(x) = e^{-x^2}, which can be done with a lambda expression. We then call the quad method on this function.

import scipy.integrate
from numpy import exp
f= lambda x:exp(-x**2)
i = scipy.integrate.quad(f, 0, 1)
print i

The above program will produce the following output.

(0.7468241328124271, 8.291413475940725e-15)

The quad function returns two numbers: the first is the integral value, and the second is an estimate of the absolute error of the integral value.

Note — Because quad requires a function as its first argument, we cannot pass an exp directly as an argument. Quaternary functions accept positive and negative infinity as limits. The quad function can integrate standard predefined NumPy single-variable functions such as exp, sin, and cos.

Multiple Integrals

The mechanics of double and triple integrals are encapsulated in the functions dblquad, tplquad, and nquad. These functions integrate over four or six arguments, respectively. All limits of inner integrals need to be defined as functions.

Double Integral

dblquad has the general form scipy.integration.dblquad(func, a, b, gfun, hfun) . Here, func is the name of the function to be integrated, ‘a’ and ‘b’ are the lower and upper limits of the x-variable, respectively, and gfun and hfun are the names of the functions defining the lower and upper limits of the y-variable.

As an example, let’s perform the double integral method.

$int_{0}^{1/2}dyint_{0}^{sqrt{1-4y^2}}16xy16xy:dx

We use lambda expressions to define the functions f, g, and h. Note that even if g and h are constants, as they may be in many cases, they must be defined as functions, as we have done here for the lower bound.

import scipy.integrate
from numpy import exp
from math import sqrt
f = lambda x, y : 16*x*y
g = lambda x : 0
h = lambda y : sqrt(1-4*y**2)
i = scipy.integrate.dblquad(f, 0, 0.5, g, h)
print i

The above program will produce the following output.

(0.5, 1.7092350012594845e-14)

In addition to the routines described above, scipy.integration has several other integration routines, including nquad , which performs n-fold multiplication of integrals, and other routines that implement various integration algorithms. However, quad and dblquad will satisfy most of our numerical integration needs.