Wolfram Alpha 和 scipy.integrate.quad 对同一个积分给出了不同的答案

Question

考虑以下函数：

import numpy as np
from scipy.special import erf

def my_func(x):
    return np.exp(x ** 2) * (1 + erf(x))

当我使用 scipy 的 quad 函数计算此函数从 -14 到 -4 的积分时，我得到以下结果：

In [3]: from scipy import integrate

In [4]: integrate.quad(my_func, -14, -4)
/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:289: UserWarning: 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)
Out[4]: (0.21896647054443383, 0.00014334175850538866)

也就是大约0.22.

然而，当我将这个积分提交给 Wolfram Alpha 时，我得到了一个非常不同的结果：

-5.29326 X 10 ^ 69.

怎么回事？我猜这与 scipy 给我的警告有关。在 python 中评估此积分的最佳方法是什么？

注意：增加 limit 会更改警告，但 scipy 结果不变：

In [5]: integrate.quad(my_func, -14, -4, limit=10000)
/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:289: UserWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  warnings.warn(msg)
Out[5]: (0.21894780966717864, 1.989164129832358e-05)

Answer 1

TL;DR：被积函数等价于 erfcx(-x)，scipy.special.erfcx 处 erfcx 的实现处理数值问题：

In [10]: from scipy.integrate import quad

In [11]: from scipy.special import erfcx

In [12]: quad(lambda x: erfcx(-x), -14, -4)
Out[12]: (0.6990732491815446, 1.4463494884581349e-13)

In [13]: quad(lambda x: erfcx(-x), -150, -50)
Out[13]: (0.6197754761443759, 4.165648376274775e-14)

您可以通过更改积分参数的符号和限制来避免 lambda 表达式：

In [14]: quad(erfcx, 4, 14)
Out[14]: (0.6990732491815446, 1.4463494884581349e-13)

问题是 1 + erf(x) 的负值 x 的数值评估。随着 x 减小，erf(x) 接近 -1。然后加 1，得到 catastrophic loss of precision，对于足够负的 x（特别是 x < -5.87），1 + erf(x) 在数值上为 0。

请注意，Wolfram Alpha 的默认行为也存在同样的问题。我必须点击 "More digits" 两次才能得到合理的答案。

解决方法是重新制定您的函数。您可以将 1+erf(x) 表示为 2*ndtr(x*sqrt(2))，其中 ndtr 是正态累积分布函数，可从 scipy.special.ndtr (see, for example, https://en.wikipedia.org/wiki/Error_function 获得）。这是您的函数的替代版本，以及将其与 scipy.integrate.quad:

集成的结果

In [133]: def func2(x):
   .....:     return np.exp(x**2) * 2 * ndtr(x * np.sqrt(2))
   .....: 

In [134]: my_func(-5)
Out[134]: 0.1107029852258767

In [135]: func2(-5)
Out[135]: 0.11070463773306743

In [136]: integrate.quad(func2, -14, -4)
Out[136]: (0.6990732491815298, 1.4469372263470424e-13)

点击 "More digits" 两次后 Wolfram Alpha 的答案是 0.6990732491815446...

这是使用数值稳定版本时函数图的样子：

为避免量级非常大的参数上溢或下溢，您可以在 log-space:

中进行部分计算

from scipy.special import log_ndtr

def func3(x):
    t = x**2 + np.log(2) + log_ndtr(x * np.sqrt(2))
    y = np.exp(t)
    return y

例如

In [20]: quad(func3, -150, -50)
Out[20]: (0.6197754761435517, 4.6850379059597266e-14)

（看起来@ali_m 在新问题中击败了我：。）

最后，正如 Simon Byrne 在的回答中指出的那样，要集成的函数可以表示为作为 erfcx(-x)，其中 erfcx 是缩放的互补误差函数。它可用 scipy.special.erfcx.

例如，

In [10]: from scipy.integrate import quad

In [11]: from scipy.special import erfcx

In [12]: quad(lambda x: erfcx(-x), -14, -4)
Out[12]: (0.6990732491815446, 1.4463494884581349e-13)

In [13]: quad(lambda x: erfcx(-x), -150, -50)
Out[13]: (0.6197754761443759, 4.165648376274775e-14)

Wolfram Alpha 和 scipy.integrate.quad 对同一个积分给出了不同的答案

Wolfram Alpha and scipy.integrate.quad give me different answers for the same integral

python

math

numpy

scipy

wolframalpha