I have a code which uses as part of the process the scipy.integrate.quad package. This should produce values (the blue squares) which follow the red line in the figure below. Some of the blue squares do not follow the trend which can only be derivative of error in the scipy.integrate.quad procedure. Does anybody know why the quad package might fail for only some of the responses? Might it have something to do with floating point arithmetic or some other underlying issue?
Note: I've seen this problem before in other functions and I'm sure the source of error is the quad package and not the rest of my code.
A likely reason for quad
returning a value near 0 is that the function is sharply localized relative to the size of the interval of integration, and the algorithm simply misses it. An example:
from scipy.integrate import quad
import numpy as np
np.around([quad(lambda x: np.exp(-x**2), -100, 100*n)[0] for n in range(1, 10)], 3)
returns [1.772, 1.772, 1.772, 0., 0., 1.772, 0., 1.772, 0.])
where zeros are incorrect. Here the Gaussian is supported near 0, and when the interval of integration is [-100, 300] or such, quad
never samples its value near 0. (For some other intervals it does, by chance.)
A remedy: use the parameter points
to indicate where the function is located, approximately. With
quad(lambda x: np.exp(-x**2), -100, 100*n, points=[-10, 10])
the above code returns 1.772 for every n.