如何以数值稳定的方式计算 (exp(t) - 1)/t？答案

【问题标题】：How can I compute (exp(t) - 1)/t in a numerically stable way?如何以数值稳定的方式计算 (exp(t) - 1)/t？
【发布时间】：2021-06-09 21:17:54
【问题描述】：

表达式 (exp(t) - 1)/t 收敛到 1，而 t 趋于 0。但是，当进行数值计算时，我们会得到不同的结果：

In [19]: (exp(10**(-12)) - 1) * (10**12)                                        
Out[19]: 1.000088900582341

In [20]: (exp(10**(-13)) - 1) * (10**13)                                        
Out[20]: 0.9992007221626409

In [21]: (exp(10**(-14)) - 1) * (10**14)                                        
Out[21]: 0.9992007221626409

In [22]: (exp(10**(-15)) - 1) * (10**15)                                        
Out[22]: 1.1102230246251565

In [23]: (exp(10**(-16)) - 1) * (10**16)                                        
Out[23]: 0.0

有没有什么方法可以在不遇到这些问题的情况下计算这个表达式？我曾考虑过使用幂级数，但我对自己实现这一点持谨慎态度，因为我不确定要使用多少个术语等实现细节。

如果相关，我将 Python 与 scipy 和 numpy 一起使用。

【问题讨论】：

使用decimal 模块对你有用吗？
检查是否有函数expm1()，它以数字合理的方式计算exp(x)-1，也就是说，在原点附近没有灾难性的抵消。
你也可以犯同样的取消错误两次并计算(exp(t)-1)/((1+t)-1)。 // 是的，numpy 有 expm1。
@njuffa 这似乎有部分帮助（在 numpy 中有这样的函数），但是，在我的情况下，t 有时太小以至于1 / t 有时会返回inf，所以看来，我实际上需要一种计算“一次全部”表达的方法。
然后你需要用它的泰勒展开替换函数，例如1+t*t==1。如果你做一些更高级的事情，比如(exp(t)-1-t)/(t**2)，上面所有的技巧都会停止工作。

标签： numerical-methods numerical-stability

【解决方案1】：

cmets 中关于微小值的讨论似乎毫无意义。如果t 太小以至于导致下溢，则表达式为1“很久以来”。事实上，泰勒发展是

1 + t/2 + t²/6 + t³/24...

只要t < 1 ulp，浮点表示就是1。

除此之外，expm1(t)/t 会做得很好。

【讨论】：