9.4: Gaussian Quadratures

Previously, we have assumed in our analysis of numerical integration schemes that the discretization points are equally-spaced. This assumption is not strictly necessary; for instance, we can easily modify the mid-point rule and trapezium rule formulas to work with non-equally-spaced points.

However, if we are free to choose the discretization points for performing numerical integration, and these points need not be equally spaced, then it is possible to exploit this freedom to further improve the accuracy of the numerical integration. The idea is to optimize the placement of the discretization points, so as to minimize the resulting numerical error. This is the basic idea behind the method of integration by Gaussian quadratures.

We will not discuss the details of this numerical integration method. To use it, you can call the scipy.integrate.quad function. This function uses a low-level numerical library named QUADPACK, which performs quadrature integration with adaptive quadratures—i.e., it automatically figures out how many discretization points should be used, and where they should be located, in order to produce a result with the desired numericaly accuracy.

Because QUADPACK figures out the discretization points for itself, you have to pass quad a function representing the integrand, rather than an array of integrand values as with trapz or simps. The standard way to call the function is

t = quad(f, a, b)

which calculates the integral

\[t = \int_a^b f(x) dx.\]

The return value is a tuple of the form (t,err), where t is the value of the integral and err is an estimate of the numerical error. The quad function also accepts many other optional inputs, which can be used to specify additional inputs for passing to the integrand function, the error tolerance, the number of subintervals to use, etc. See the documentation for details.

The choice of whether to perform numerical integration using Simpson's rule (simps) or Gaussian quadratures (quad) is situational. If you already know the values of the integrands at a pre-determined set of discretization points (e.g., from the result of a finite-difference calculation), then use simps . If you can define a function that can quickly calculate the value of the integrand at any point, use quad.

Here is an example of using quad to compute the integral \(\int_{0}^\infty \frac{dx}{x^2 + 1}\):

from scipy import *
from scipy.integrate import quad

def f(x):
    return 1.0 / (x*x+1)

integ, _ = quad(f, 0, 1000)
print(integ)

(Note that quad returns two values; the first is the computed value of the integral, and the other is the absolute error, which we're not interested in, so we toss it in the "throwaway" variable named _. See the documentation for details.) Running the above program prints the result \(1.569796\dots,\) which differs from the exact result, \(\pi/2 = 1.570796\dots,\) by a relative error of \(0.06\%\).

Here is another example, where the integrand takes an additional parameter: \(\int_0^\infty x e^{-\lambda x}\, dx\):

from scipy import *
from scipy.integrate import quad

def f(x, lambd):
    return x * exp(-lambd * x)

integ, _ = quad(f, 0, 100, args=(0.5))
print(integ)

Running the program prints the result \(4.0\), which agrees with the exact result of \(1/\lambda^2\) for the chosen value of the parameter \(\lambda =0.5\).