评估收敛发散项的乘积

如何解决评估收敛发散项的乘积

我想在 Python 中以数字方式计算以下表达式

exp((x^2))*Erfc(x).

其中 Erfc 表示互补误差函数 (https://en.wikipedia.org/wiki/Error_function#Complementary_error_function)。 Erfc(x) 可从 Scipy Special 获得。即使对于大 x，$exp(x^2)$ 发散，乘积也具有有限值。准确地说

Lim x->\infinity exp((x^2))*Erfc(x) = 1/(sqrt(Pi) * x).

不幸的是，Python 不计算这个表达式，而是在 x 的大值时返回 NaN。有没有办法实现这个函数，使得对于大 x，Python 仍然产生正确的答案？ （我的代码和输出如下。）

import numpy as np
import scipy 
from scipy import special as spl

def my_fun(x):

    return np.exp(x**2)*spl.erfc(x);


print (my_fun(10))

print (my_fun(50))

print (my_fun(100))

0.05614099274382259

    Warning (from warnings module):
      File "/Users/spjoy/Working/Projects/Wigner Crystal/Python/num_eval.py",line 7
        return np.exp(x**2)*spl.erfc(x);
    RuntimeWarning: overflow encountered in exp
    
    Warning (from warnings module):
      File "/Users/spjoy/Working/Projects/Wigner Crystal/Python/num_eval.py",line 7
        return np.exp(x**2)*spl.erfc(x);
    RuntimeWarning: invalid value encountered in double_scalars
    nan
    nan

解决方法

SciPy 具有该函数的实现：scipy.special.erfcx。

例如

In [8]: import numpy as np

In [9]: from scipy.special import erfcx,erfc

一个适中的值，其中 exp(x**2)*erfc(x) 不受上溢和下溢的影响，以验证我们使用 erfcx(x) 获得相同的结果：

In [10]: x = 2.5

In [11]: np.exp(x**2) * erfc(x)
Out[11]: 0.21080636406114353

In [12]: erfcx(x)
Out[12]: 0.2108063640611436

x 的大值，其中 exp(x**2) 溢出，但 erfcx(x) 给出结果：

In [13]: x = 50.0

In [14]: np.exp(x**2) * erfc(x)
<ipython-input-14-bb75732dec2d>:1: RuntimeWarning: overflow encountered in exp
  np.exp(x**2) * erfc(x)
<ipython-input-14-bb75732dec2d>:1: RuntimeWarning: invalid value encountered in double_scalars
  np.exp(x**2) * erfc(x)
Out[14]: nan

In [15]: erfcx(x)
Out[15]: 0.011281536265323772

In [16]: 1/(np.sqrt(np.pi)*x)  # Verify the asymptotic approximation.
Out[16]: 0.011283791670955126