负 x 的泰勒级数 e^x 的误差

Posted 2023-03-28

【中文标题】负 x 的泰勒级数 e^x 的误差

【英文标题】：Error of Taylor Series e^x for negative x

【发布时间】：2021-12-26 21:25:47

【问题描述】：

我在使用泰勒级数计算 e^x 时注意到，当我们为负 x 计算它时，绝对误差很大。是因为我们没有足够的精度来计算它吗？

（我知道为了防止它我们可以使用 e^(-x)=1/e^x)

#include <stdio.h>
#include <math.h>

double Exp(double x);

int main(void)

    double x;
    printf("x="); 
    scanf("%le", &x);

    printf("%le", Exp(x)); 
    return 0;


double Exp(double x)

    double h, eps = 1.e-16, Sum = 1.0; 
    int i = 2;

    h = x; 

    do
    
        Sum += h; 
        h *= x / i;
        i++;
     while (fabs(h) > eps);

    return Sum ; 

例如： x=-40 值是 4.24835e-18 但程序给我 3.116952e-01。绝对误差是 ~0.311

x=-50 值是 1.92875e-22 程序给我 2.041833e+03。绝对误差是 ~2041.833

【问题讨论】：

帮助您的读者看到您所看到的。您是如何使用您共享的代码注意到它的？

显示x 的值来说明情况，显示您获得的结果，并说明您观察到的绝对错误。

在每次迭代中为h 的值添加打印。它正在处理一些相当大的正数和负数，希望得到一个小的结果。 Subtraction is dangerous.

【参考方案1】：

该问题是由算法中间阶段的舍入误差引起的。 h 与 40/2 * 40/3 * 40 / 4 * ... 一样迅速增长，并在符号上振荡。连续迭代的i、h 和Sum 的值可以在下面找到（为简洁起见省略了一些数据点）：

x=-40
i=2 h=800 Sum=-39
i=3 h=-10666.7 Sum=761
i=4 h=106667 Sum=-9905.67
i=5 h=-853333 Sum=96761
i=6 h=5.68889e+06 Sum=-756572
...
i=37 h=-1.37241e+16 Sum=6.63949e+15
i=38 h=1.44464e+16 Sum=-7.08457e+15
i=39 h=-1.48168e+16 Sum=7.36181e+15
i=40 h=1.48168e+16 Sum=-7.45499e+15
i=41 h=-1.44554e+16 Sum=7.36181e+15
i=42 h=1.37671e+16 Sum=-7.09361e+15
i=43 h=-1.28066e+16 Sum=6.67346e+15
i=44 h=1.16423e+16 Sum=-6.13311e+15
i=45 h=-1.03487e+16 Sum=5.50923e+15
i=46 h=8.99891e+15 Sum=-4.83952e+15
...
i=97 h=-2610.22 Sum=1852.36
i=98 h=1065.4 Sum=-757.861
i=99 h=-430.463 Sum=307.534
...
i=138 h=1.75514e-16 Sum=0.311695
i=139 h=-5.05076e-17 Sum=0.311695
3.116952e-01

sum 的峰值幅度为7e15。这就是失去精度的地方。类型 double 可以以大约 1e-16 的精度表示。这给出了大约0.1 - 1 的预期绝对误差。 由于预期总和（exp(-40) 的值接近于零，最终的绝对误差接近于部分总和的最大绝对误差。

对于x=-50，总和的峰值是1.5e20，由于double 的有限表示，绝对误差大约为1e3 - 1e4，接近观察到的值。

如果不对算法进行重大更改以避免形成这些部分和，则无法解决太多问题。或者，将exp(-x) 计算为1/exp(x)。

【参考方案2】：

对于负 x，即使在 1.0 + x 的第一个总和中添加交替的 +/- 项也会产生计算问题，因为最终总和误差可能与 1.0 的最低有效位或大约 1 部分一样糟糕在 10¹⁶。这意味着x_min 中的Exp(x_min) == 1.0e-16 是最小有用计算值（例如x 约为-36）

一个简单的解决方案是形成一个好的Exp(positive_x) 和负值...

double Exp(double x) 
  if (x < 0) 
    return 1.0 / Exp(-x);
  
  ...

一个好的（和简单的）Exp(positive_x) 计算项直到term + 1.0 仍然是 1.0，因为额外的小项不会显着改变总和。适用于所有 x（非常小的错误），但在结果应该低于正常值时可以使用改进。

double my_exp(double x) 
  if (x < 0) 
    return 1.0 / my_exp(-x);
  
  double sum = 1.0;
  unsigned n = 1;
  double term = 1.0;
  do 
    term *= x / n++;
    sum += term;
    if (!isfinite(term)) 
      return term;
    
   while (1.0 != term + 1.0);
  return sum;