4

考虑以下函数：

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

1 回答 1

8

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 在新问题中击败了我：Tricking numpy/python into representation very large and very small numbers。）

最后，正如 Simon Byrne 在Tricking numpy/python into representation 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)
于 2015-08-31T00:04:26.523 回答