遇到溢出：exp() 乘以 erfc()答案

【问题标题】：overflow encountered: multipliation of exp() times erfc()遇到溢出：exp() 乘以 erfc()
【发布时间】：2020-04-08 15:49:53
【问题描述】：

我想评估以下公式：

c = exp{x}*erfc{y}

（参见下面代码中x 和y 的定义。）

问题是 x 和 y 变得非常大，我得到 exp{x} 的非常大的值和非常小的 erfc(y) 的值。

import numpy as np
import scipy as sci
k = 5.7e-3
D = 1.53e-8
R = 1.5e-5
r = 1e-6

t = np.linspace(0.0,12,10)

x = (r/R)  +  (D/(R*R) - k)*t
y = (r/(2*np.sqrt(D*t))) + np.sqrt(D*t)/R

exp_x = np.exp(x)
erfc_y = sci.special.erfc(y)
print("x = \n{} ".format(x))
print("y = \n{}".format(y))

print("exp(x) = \n{}".format(exp_x))
print("erfc(y) = \n{}".format(erfc_y))

print("exp(x) * erfc(y)= \n{}".format(exp_x*erfc_y))

我的想法是将评价改为

log{exp(x)*erfc(y)} = log{exp(x)} + log{erfc(y)} = x + log{erfc(y)}

之后，我可以计算

exp(x + log{erfc(y)})

但这里有问题：当我要计算时

log{erfc(y)} = log{1 - erfc(y)}

我遇到类似的问题，erfc 会接近 1，我会遇到精度问题。

有什么想法可以解决我的问题吗？

【问题讨论】：

作为一个快速的想法，您可以将log 近似应用于您的样本数据的一部分（假设您拥有erfc(y)==0 的唯一点是第一个）
是的，我也有同样的想法。我会测试这个，当我有时间的时候。但是，我能够用 Mathematica 解决我的问题。这不是我喜欢解决问题的方式，但它现在有效。
将log1p(x) 用于log(1+x)

标签： python numpy scipy

【解决方案1】：

这个问题的解决方案是Cadena给出的，在这篇文章中：

Cadena, F. (1989)。 使用个人计算机解决污染物传输模型的数值方法。教育中的计算机， 9(2), 34-6.

我无法获得文章，但解决方案也由

Lin, J. S. 和 Hildemann, L. M. (1995)。 用于预测垃圾填埋场挥发性有机化合物气体排放的非稳态分析模型。有害物质杂志, 40(3), 271-295.

Cadena 提出用另一个具有负幂项的指数函数来近似互补误差函数。它减少了导致溢出的原始指数函数所产生的高功率。

假设W是需要求值的乘法：

W=exp(y)*erfc(x)

W 可以近似为：

W = (a1*t+a2*t^2+a3*t^3+a4*t^4+a5*t^5)*exp(y-x^2), t=1/(1+p*x), if x >= 0

W = 2*exp(y)-(a1*t+a2*t^2+a3*t^3+a4*t^4+a5*t^5)*exp(y-|x|^2), t=1/(1+p*|x|), if x < 0

在哪里

a1 = 0.254829592
a2 = -0.284496736
a3 = 1.42141741
a4 = -1.453152027
a5 = 1.061405429
p = 0.3275911

【讨论】：

【解决方案2】：

使用scipy.special.erfcx() 得到exp(x**2)*erfc(x)。

然后你只需实现erfcx(x)*exp(x-x**2)，它再次给出erfc(x)*exp(x)。

如果您的论点x=y，这当然效果最好。

【讨论】：