I am trying to compute the integrals more precise by specifying the parameter epsabs for scipy.integrate.quad, say we are integrating the function sin(x) / x^2 from 1e-16 to 1.0
from scipy.integrate import quad
import numpy
integrand = lambda x: numpy.sin(x) / x ** 2
integral = quad(integrand, 1e-16, 1.0)
This gives us
(36.760078801255595, 0.01091187908038005)
To make the results more precise, we specify the absolute error tolerance by epsabs
from scipy.integrate import quad
import numpy
integrand = lambda x: numpy.sin(x) / x ** 2
integral = quad(integrand, 1e-16, 1.0, epsabs = 1e-4)
The result is exactly the same and the error is still as large as 0.0109! Am I understanding the parameter epsabs wrong? What should I do differently to increase the precision of integral?
According to scipy manual quad function has limit argument to specify
An upper bound on the number of subintervals used in the adaptive algorithm.
By default the value of limit is 50.
You code return warning message
quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used.
warnings.warn(msg, IntegrationWarning)
You have to change limit argument, i.e.:
from scipy.integrate import quad
import numpy
integrand = lambda x: numpy.sin(x) / x ** 2
print(quad(integrand, 1e-16, 1.0, epsabs = 1e-4, limit=100))
Output:
(36.7600787611414, 3.635057215414274e-05)
There is no warning message in output. Number of subdivisions is under 100 and quad got required accuracy.