How do I use parameter epsabs in scipy.integrate.quad in Python?

Question

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?

Answer 1

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.