python - 如何在 Python 中使用 scipy.integrate.quad 中的参数 epsabs？

Question

我试图通过为 scipy.integrate.quad 指定参数 epsabs 来更精确地计算积分，假设我们正在积分函数 sin(x )/x^2 从 1e-16 到 1.0

from scipy.integrate import quad
import numpy

integrand = lambda x: numpy.sin(x) / x ** 2
integral = quad(integrand, 1e-16, 1.0)

这给了我们

(36.760078801255595, 0.01091187908038005)

为了使结果更加精确，我们通过 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)

结果一模一样，误差仍然大到0.0109!我对参数 epsabs 的理解是否错误？我应该采取哪些不同的措施来提高积分的精度？

Answer 1

根据scipy手册quad function有 limit 参数来指定

An upper bound on the number of subintervals used in the adaptive algorithm.

默认情况下，limit 的值为 50。您编写返回警告消息的代码

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)

您必须更改limit参数，即:

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))

输出:

(36.7600787611414, 3.635057215414274e-05)

输出中没有警告消息。分割数低于 100，quad 达到了所需的精度。