pythonscipynumerical-integrationquad

Why does scipy.integrate.quad return 0 instead of the correct value of the integral?


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.

Result showingintegration error


Solution

  • 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.