python - 使用scipy在dblquad集成中获得的错误是什么意思？

Question

我想知道不同 Python 集成例程给出的错误的含义，例如dblquad。由于我们不知道确切的积分值，误差的估计是如何计算的？引用是什么？在我的一些计算中，我发现增加积分极限会使误差达到极高的值。由于这只是一个误差估计，是否值得依赖这样的结果？

Answer 1

简答

is it at all advisable to rely on such a result?

在大多数情况下，是的。但是，如果您认为积分例程表现异常并且您不信任它的输出，请尝试改变方法:例如，将积分区域划分为多个部分，分别对每个部分积分并查看结果是否相加。

说明

在数值积分中，人们通过使用两种方法计算积分(或具有两种步长的相同方法)并考虑结果之间的差异来估计误差。 Deep within Fortran source我们发现的 SciPy 的 quadpack 例程

abserr = dabs((resk-resg)*hlgth)

其中resg是10点高斯公式的结果，resk是21点克朗罗德公式的结果。请参阅维基百科文章 Gauss–Kronrod quadrature formula对于这些的数学意义。 (hlgth 是积分积分长度的一半；这里的长度是由于缩放。)

实际上，我引用的公式并不是最终的误差估计，它是一种非常粗略的第一种方法。两行之后我们看到

abserr = resasc*(0.2d+03*abserr/resasc)**1.5d+00)

这正是维基百科文章所说的:

The recommended error estimate is (200*|gauss - kronrod|)^1.5

这种对绝对误差的估计不能保证限制计算积分和实际积分(后者未知)之间的差异。 “推荐”的估计往往在实践中起作用，指数 1.5(该方法的收敛顺序)有一些数学理由，但我们永远不知道它是否真的涵盖了实际误差。

毕竟，该函数仅在其域内的有限多个点上进行评估。据我们所知，它可能恰好在那些点为 0，而在其他地方，在积分例程没有看到的地方，它可能恰好为 0。

例子

这是一个看起来很简单的函数的积分，但它被错误地计算了:

from scipy.integrate import quad
import numpy as np
quad(lambda x: np.exp(-x**2), -1e2, 1e3)

返回 (4.176612573788305e-60, 7.896357364711954e-60)。实际积分约为 1.77(圆周率的平方根)。错误估计 8e-60 是完全错误的，值 4e-60 也是如此。原因是这个函数局部化在0附近，积分的区间是[-100, 1000]，要大很多。 quad 算法并没有碰巧在具有相当大值的任何点对函数进行采样，因此它继续认为它在任何地方都几乎为零。