在 x=0x=0 附近计算 (ex−1)/x(ex−1)/x

计算科学 C++ C
2021-12-15 01:04:38

函数 f:x(ex1)/xx=0 附近有奇点。不过,这个奇点可以被解除:对于 x=1,应该有 f(x)=1,因为

ex=k=0xkk!
因此
(ex1)/x=k=1xk1k!
但是,形式 (ex1)/x 不仅没有在 x=0 处定义,而且在该点附近也是数值不稳定的;
为了在数值上计算非常小的 xf(x),可以使用泰勒展开式,即上述幂级数的截断。

:函数 f 有名字吗?换句话说,这是一个普遍的问题吗?

:是否有人知道可以很好地处理这种情况的 C/C++ 库,即使用接近 0 的适当度数的泰勒展开式和远离零的其他表示?

4个回答

可能可以从作为 C99 标准一部分的函数 expm1 开始,并在 x=0 附近准确计算 ex1

这是取消错误的一个实例。C 标准库(从 C99 开始)包含一个调用函数expm1来避免这个问题。如果您使用expm1(x) / x代替(exp(x) - 1.0) / x,您将不会遇到此问题(请参见下图)。 <code>fabs(expm1(x) / x - (exp(x) - 1.0) / x)</code>

这个特定问题的细节和解决方案在数值算法的准确性和稳定性的第 1.14.1 节中进行了详细讨论。W. Kahan 的论文第 19 页也解释了相同的解决方案,题为“浮点计算中舍入的无意识评估是多么无用?” . expm1GNU C 库中的实际实现与上面参考文献中描述的方法不同,并且在源代码中进行了详细记录

要回答您的第一个问题,不,该函数没有名称(至少不是广为人知的名称)。

正如其他人所提到的,计算函数的最佳方法是处理几种特殊情况。这就是任何库计算函数的方式。

  1. 情况 0:x = 0,返回 1。
  2. 案例1:$|x|<\delta$,返回$1+x/2$。根据您的数值精度和数据类型选择 $\delta$。对于,应该差不多对于浮动,大约是.|x|<δ, return 1+x/2. Choose δ empirically for your numerical precision and data type. For double2e-85e-4
  3. 其他情况:返回expm1(x)/x

使用截断的泰勒级数,您可以更复杂和更特殊的情况,但这可能不值得。实际上,情况1是否需要单独处理并不完全清楚,因为正如k20所指出的那样,取消是安全的。但是,分开处理会让我对它更有信心。

我记得这个问题早些时候在这个网站上被问过,令人惊讶的是答案是你只需要将特殊情况完全等于零。误差几乎为零。我没有链接。

是的,这个答案是完全错误的。我不知道为什么它被如此赞成,可能是因为它被如此权威地陈述。我找到了我想到的链接。它是在这里的数学堆栈交换上,而不是在 scicomp 堆栈交换上。expm1JM 在答案中给出了 -free 错误消除公式,并使用了转换u = exp(x)