使用scipy在dblquad积分中获得的误差是什么意思？

Question

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

Answer 1

简答

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

In most cases, yes. But if you think the integration routine is behaving strangely and you don't trust its output, do try to change the approach: for example, divide the region of integration into parts, integrate over each part separately and see if the results add up.

Explanation

In numerical integration, one estimates the error by using two methods to compute the integral (or the same method with two step sizes), and considering the difference between the results. Deep within Fortran source of SciPy's quadpack routines we find

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

where resg is the result of the 10-point Gauss formula, and resk is the result of the 21-point Kronrod formula. See the Wikipedia article Gauss–Kronrod quadrature formula for the mathematical meaning of these. (The hlgth is half the length of the integral of integration; the length is here due to scaling.)

Actually, the formula I quoted is not the final error estimate, it's a very crude first approach to it. Two lines later we see

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

which is exactly what the Wikipedia article states:

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

This kind of estimate of the absolute error is not guaranteed to bound the difference between the computed and actual integral (the latter being unknown). The "recommended" estimate tends to work in practice, and there's some mathematical justification for the exponent 1.5 (the order of convergence of the method) but we never know if it really covers the actual error.

After all, the function was only evaluated at finitely many points within its domain. For all we know, it could happen to be 0 at those points and something huge elsewhere, at the points where the integration routine didn't look.

Example

Here is an integral of a simple-looking function which is evaluated incorrectly:

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算法没有碰巧在具有可观值的任何点对函数进行采样，因此它继续认为它在任何地方都几乎为零。