Numerical Integration Goes Wrong

Sometimes one imagines that computers are so good at numerical integration, also known as quadrature, that one can forget all the details and leave them to it. This page shows that this is not quite so.

The scipy package provides several integration routines, and the integrand we consider here is the Gaussian function. We shall use the form

f(x) = σ exp(-σ²x²)

One can plot this function for various values of σ using matplotlib.

$ ipython3
In [1]: import matplotlib.pyplot as plt

In [2]: import numpy as np

In [3]: x=np.linspace(-3,3,200)

In [4]: for s in [1,2,4,8]:
   ...:     plt.plot(x,s*np.exp(-x*x*s*s),label=("s=%f"%(s)))
   ...:

In [5]: plt.show()

The integral of this from -∞ to +∞ is √π, and from 0 to +∞ is 0.5√π. The integral does not depend on the parameter σ, so we shall write some simple code to evaluate it for various values of σ.

#!/usr/bin/python3
from math import exp,sqrt,pi
from scipy import integrate

def gaussian(x):
    return(sigma*exp(-x*x*sigma*sigma))

print("0.5*sqrt(pi)=",0.5*sqrt(pi))

sigma=1
while(sigma<100):
    (y,err)=integrate.quad(gaussian,0,100)
    print("%5d %f %g"%(sigma,y,err))
    sigma=2*sigma

This prints the integral from 0 to 100, and the estimated error, for various values of σ. For the values of σ used, the difference between integrating to 100, and integrating to ∞, is negligible.

0.5*sqrt(pi)= 0.8862269254527579
    1 0.886227 9.88409e-11
    2 0.886227 9.88409e-11
    4 0.886227 9.88409e-11
    8 0.886227 9.88409e-11
   16 0.886227 9.88409e-11
   32 0.886227 9.88409e-11
   64 0.000000 3.99851e-20

So it worked fine up to σ=32, but was very wrong at σ=64. Very wrong in that the correct answer lies far outside the estimated error give in the final column. Further investigation shows a sharp change from σ=42 to σ=43.

As σ increases, the Gaussian function gets thinner, and taller. What may be surprising is that the integration routine does not evaluate the function at either limit. If it did try evaluating the function at x=0, it would return the answer σ, which would show that the integral is unlikely to be zero. If we assume σ=43, then at x=0.1 the function has the value 43 exp(-4.3*4.3), which is 4x10^-7, already quite close to zero, but certainly not zero. For all values of x<0.045 the function's value is greater than one.

What went wrong? The integration routine evaluated the function at merely 63 points between 0 and 100. It did bias its sampling somewhat around the lower limit, but the two smallest values of x it sampled were 0.1085 and 0.2171. How can we tell? Well, it is easy to add to the function gaussian(x) the line print(x) before the return statement.

If instead of integrating the above Gaussian, we integrate the same with the constant 0.000001 (1e-6) added, the integral over the range used should increase by 0.0001 (1e-4). Does it? Does it switch from the correct answer to forgetting the Gaussian peak at the same value of σ as before?

Differential Equations

The integration of differential equations can go astray due to rounding errors. The Physics section has an example of rounding errors in collisions.