Numerical integration visualized

scipy.integrate.quad is a popular method of numerical integration with Python. Let’s see how it chooses the points at which to evaluate the function being integrated. Begin with a simple example, the exponential function.

exp1
exp(x) on [-1, 1]
The blue dots indicate the evaluation points, their y-coordinate being the order of evaluation. So, it begins with x=0, continues with something close to -1, etc. The function, drawn in red, is not to its true vertical scale.

The placement of dots gives away the method of integration: it is the Gauss-Kronrod quadrature with 10 Gauss nodes and 21 Kronrod nodes, abbreviated (G10, K21).  The Gauss nodes are included in the Kronrod nodes, and of course the function is not evaluated there again. The evaluation process is slightly out of order in that 0 is a Kronrod node but comes first, followed by 10 Gauss nodes, followed by 10 remaining Kronrod nodes. The process ends there, as the comparison of G10 and K21 results shows the necessary precision was reached.

The square root on [0, 1] is not such a nice function, so (G10, K21) does not reach the required precision at once. The interval is bisected again and again until it does.

sqrt
sqrt(x) on [0, 4]
  Surely the cube root is even worse. Or is it?

cbrt
cbrt(x) on [-1, 1]
The nodes are symmetric about the midpoint of the interval of integration. So for any odd function on an interval symmetric about 0 we get G10 = 0 and K21 = 0, and the process stops at once. To see the effect of the cube root singularity, one has to use a non-symmetric interval such as [-1, 2].

cbrt12
cbrt(x) on [-1, 2]
That’s still fewer subdivisions than for the square root: cancellation between the left and right neighborhoods of 0 still helps. Let’s look at smooth functions next.

sinxsq
sin(x^2) on [0, 7]
Rapid oscillations forces subdivisions here. There are also other reasons for subdivision, such as the loss of analyticity:

exp_1x
exp(-1/x) on [0, 1]
Although exp(-1/x) is infinitely differentiable on this interval, the fact that it is not analytic at 0 makes it a more difficult integrand that exp(x). Finally, an analytic function with no oscillation which still needs a bunch of subintervals:

runge
1/(1+x^2) on [-5, 5]
This is the standard example used to illustrate the Runge phenomenon. Although we are not interpolating here, numerical integration is also influenced by the function having a small radius of convergence of its Taylor series. After all, G10 can be thought of as degree-9 interpolation of the given function at the Gauss nodes (the zeros of a Legendre polynomial),  with the formula returning the integral of the interpolating polynomial.

The code used to plot these things:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
global eval_points
f = lambda x: np.exp(x)
def f_int(x):
    eval_points.append(x)
    return f(x)

eval_points = []
a, b = -1, 1
quad(f_int, a, b)
n = len(eval_points)
t = np.linspace(a, b, 1000)
y = f(t)
yy = n*(y - y.min())/(y.max() - y.min())
plt.plot(t, yy, 'r') 
plt.plot(eval_points, np.arange(0, n), '.')
plt.show()

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s