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

【问题标题】：What is the meaning of error obtained in dblquad integration using scipy?使用 scipy 在 dblquad 集成中获得的错误是什么意思？
【发布时间】：2017-05-30 11:18:37
【问题描述】：

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

【问题讨论】：

标签： python scipy integration

【解决方案1】：

简答

依赖这样的结果是否可取？

在大多数情况下，是的。但是，如果您认为集成例程的行为很奇怪并且您不信任它的输出，请尝试改变方法：例如，将集成区域划分为多个部分，分别对每个部分进行集成，然后查看结果是否相加。

说明

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

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

其中 resg 是 10 点高斯公式的结果，而 resk 是 21 点 Kronrod 公式的结果。有关这些的数学含义，请参阅 Wikipedia 文章 Gauss–Kronrod quadrature formula。（hlgth 是积分长度的一半；这里的长度是由于缩放。）

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

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

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

建议的误差估计为 (200*|gauss - kronrod|)^1.5

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

毕竟，该函数仅在其域内的有限多个点上进行评估。据我们所知，在这些点上它可能恰好为 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（pi 的平方根）。错误估计 8e-60 是完全错误的，值 4e-60 也是如此。原因是这个函数局部化在0附近，积分区间是[-100, 1000]，要大很多。 quad 算法没有碰巧在任何具有可观值的点对函数进行采样，因此它继续认为它几乎处处为零。

【讨论】：

感谢您解释和引用 abserr 公式。我的被积函数非常复杂，并且被评估为域 [0, 5e9] 中的 40k 点（0 - inf 的双积分）。其中 5e9 是无穷大的替代品。我得到的错误可能高达 1e6/1e11。这就是为什么我想知道它是否可以依赖。我没有收到任何溢出错误，但是我增加的限制越多（比如 5e11），错误就越高（1e11 from 1e6）。
写入应该非常高的整数值，例如 1e21，错误看起来很小 (1e6/1e11)。您认为在这种情况下错误可以吗？
非常大的集成限制可能是个问题；我怀疑你得到的数字是假的。我在回答中包含了一个发生这种事情的例子。据我所知，你的积分可能是发散的。如果没有，改变变量可能会有所帮助（也许极坐标会更好？）。或者，一种二元方法：在大小为 2^k 的区域上对不同的 k 值进行积分；比较结果。这不是我可以在 cmets 中解决的问题。
感谢您的示例和您的 cmets。我明白你的意思。我将尝试尽可能精确地指定集成域。积分已经在极坐标中。我没有使用二元方法，但我会检查一下。
如果它已经在极地，二元分解很方便：在 0