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

当我从评估此函数的积分-14来-4使用scipy的quad功能,我得到以下结果:

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),并且erfcxat 的实现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时,你会得到灾难性的精度损失,而对于足够的负数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在新问题中击败了我:欺骗numpy/python代表非常大和非常小的数字.)

最后,正如Simon Byrne在Tricking numpy/python的回答中指出的那样,代表非常大和非常小的数字,要集成的函数可以表示为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)