I believe the issue is due to np.exp(-x)
quickly becoming very small as x
increases, which results in evaluating as zero due to limited numerical precision. For example, even for x
as small as x=10**2*
, np.exp(-x)
evaluates to 3.72007597602e-44
, whereas x
values of order 10**3
or above result in 0
.
I do not know the implementation specifics of quad
, but it probably performs some kind of sampling of the function to be integrated over the given integration range. For a large upper integration limit, most of the samples of np.exp(-x)
evaluate to zero, hence the integral value is underestimated. (Note that in these cases the provided absolute error by quad
is of the same order as the integral value which is an indicator that the latter is unreliable.)
One approach to avoid this issue is to restrict the integration upper bound to a value above which the numerical function becomes very small (and, hence, contributes marginally to the integral value). From your code snipet, the value of 10**4
appears to be a good choice, however, a value of 10**2
also results in an accurate evaluation of the integral.
Another approach to avoid numerical precision issues is to use a module that performs computation in arbitrary precision arithmetic, such as mpmath
. For example, for x=10**5
, mpmath
evaluates exp(-x)
as follows (using the native mpmath
exponential function)
import mpmath as mp
print(mp.exp(-10**5))
3.56294956530937e-43430
Note how small this value is. With the standard hardware numerical precision (used by numpy
) this value becomes 0
.
mpmath
offers an integration function (mp.quad
), which can provide an accurate estimate of the integral for arbitrary values of the upper integral bound.
import mpmath as mp
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))
1.0
0.999999650469474
0.999999999996516
0.999999999999997
We can also obtain even more accurate estimates by increasing the precision to, say, 50
decimal points (from 15
which is the standard precision)
mp.mp.dps = 50;
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))
1.0
0.99999999999999999999999999999999999999999829880262
0.99999999999999999999999999999999999999999999997463
0.99999999999999999999999999999999999999999999999998
In general, the cost for obtaining this accuracy is an increased computation time.
P.S.: It goes without saying that if you are able to evaluate your integral analytically in the first place (e.g., with the help of Sympy
) you can forget all the above.