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

您可以通过更改积分参数和限制的符号来避免 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 的默认行为也存在同样的问题。我不得不点击“更多数字”两次才能得到合理的答案。

解决方法是重新表述您的函数。您可以将 1+erf(x) 表示为 2*ndtr(x*sqrt(2))，其中 ndtr 是正态累积分布函数，可从 scipy.special.ndtr 获得(参见，例如，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)

点击“更多数字”两次后 Wolfram Alpha 的答案是 0.6990732491815446...

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

---

为了避免具有非常大数量级的参数的溢出或下溢，您可以在对数空间中进行部分计算:

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 在新问题中打败了我:Tricking numpy/python into representing very large and very small numbers。)

---

最后，正如 Simon Byrne 在 Tricking numpy/python into representing very large and very small numbers 的回答中指出的那样，被积分的函数可以表示为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)