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(-σ2x2)
One can plot this function for various values of σ using matplotlib.
$ ipython3 In : import matplotlib.pyplot as plt In : import numpy as np In : x=np.linspace(-3,3,200) In : for s in [1,2,4,8]: ...: plt.plot(x,s*np.exp(-x*x*s*s),label=("s=%f"%(s))) ...: In : 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
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?