c++ - exp(-x) 和 exp(+x) 之间的泰勒级数差异

Question

我正在尝试编写一个程序来计算 exp(-x) 和 exp(x) 的泰勒级数最多 200 次迭代，对于大 x。 (exp(x)=1+x+x^2/2+...).

我的程序非常简单，看起来应该可以完美运行。然而，它对 exp(-x) 发散，但对 exp(+x) 收敛得很好。到目前为止，这是我的代码:

long double x = 100.0, sum = 1.0, last = 1.0;

for(int i = 1; i < 200; i++) {
        last *= x / i;    //multiply the last term in the series from the previous term by x/n
        sum += last; //add this new last term to the sum
    }
cout << "exp(+x) = " << sum << endl;

x = -100.0; //redo but now letting x<0
sum = 1.0;
last = 1.0;

for(int i = 1; i < 200; i++) {
            last *= x / i;
            sum += last;
        }
    cout << "exp(-x) = " << sum << endl;

当我运行它时，我得到以下输出:

exp(+x) = 2.68811691354e+43 
exp(-x) = -8.42078025179e+24

当实际值为:

exp(+x) = 2.68811714182e+43 
exp(-x) = 3.72007597602e-44

如您所见，它对正计算工作得很好，但对负计算却不行。有没有人知道为什么舍入误差会因为仅在每个其他项上加上一个负数而变得如此错误？另外，有什么我可以实现的来解决这个问题吗？

提前致谢!

Answer 1

我认为这实际上与浮点逼近误差没有任何关系，我认为还有另一个更重要的误差源。

正如您自己所说，您的方法非常简单。您正在 x=0 对该函数进行泰勒级数逼近，然后在 x=-100 对其求值。

您实际期望这种方法有多准确，为什么？

在高层次上，您应该只期望您的方法在 x=0 附近的一个狭窄区域是准确的。泰勒近似定理告诉你，例如如果您采用 N 系列中的 x=0 项，您的近似值至少对 O(|x|)^(N+1) 是准确的。所以如果你有 200 个术语，你应该准确到例如在 10^(-60) 范围内或在 [-0.5, 0.5] 范围内。但是在 x=100 中，泰勒定理只会给你一个非常糟糕的界限。

从概念上讲，您知道随着 e^{-x} 趋于负无穷大，x 趋于零。但是你的近似函数是一个固定次数的多项式，任何非常数多项式都趋向于渐进地加无穷大或负无穷大。因此，如果您考虑 x 的整个可能值范围，则相对误差必须是无界的。

所以简而言之，我认为您应该重新考虑您的方法。您可能会考虑的一件事是，仅对满足 x 的 -0.5f <= x <= 0.5f 值使用泰勒级数方法。对于任何大于 x 的 0.5f，您尝试将 x 除以二并递归调用该函数，然后对结果进行平方。或者类似这样的东西。

为了获得最佳结果，您可能应该使用既定方法。

编辑:

我决定编写一些代码，看看我的想法实际效果如何。它似乎与 x = -10000 到 x = 10000 范围内的 C 库实现完美匹配，至少与我显示的精度一样多。 :)

另请注意，即使 x 的值大于 100，我的方法也是准确的，其中泰勒级数方法实际上也在正端失去了准确性。

#include <cmath>
#include <iostream>

long double taylor_series(long double x)
{
    long double sum = 1.0, last = 1.0;

    for(int i = 1; i < 200; i++) {
            last *= x / i;    //multiply the last term in the series from the previous term by x/n
            sum += last; //add this new last term to the sum
    }

    return sum;
}

long double hybrid(long double x)
{
    long double temp;
    if (-0.5 <= x && x <= 0.5) {
        return taylor_series(x);
    } else {
        temp = hybrid(x / 2);
        return (temp * temp);
    }
}

long double true_value(long double x) {
    return expl(x);
}

void output_samples(long double x) {
    std::cout << "x = " << x << std::endl;
    std::cout << "\ttaylor series = " << taylor_series(x) << std::endl;
    std::cout << "\thybrid method = " << hybrid(x) << std::endl;
    std::cout << "\tlibrary = " << true_value(x) << std::endl;
}

int main() {
    output_samples(-10000);
    output_samples(-1000);
    output_samples(-100);
    output_samples(-10);
    output_samples(-1);
    output_samples(-0.1);
    output_samples(0);
    output_samples(0.1);
    output_samples(1);
    output_samples(10);
    output_samples(100);
    output_samples(1000);
    output_samples(10000);
}

输出:

$ ./main 
x = -10000
    taylor series = -2.48647e+423
    hybrid method = 1.13548e-4343
    library = 1.13548e-4343
x = -1000
    taylor series = -2.11476e+224
    hybrid method = 5.07596e-435
    library = 5.07596e-435
x = -100
    taylor series = -8.49406e+24
    hybrid method = 3.72008e-44
    library = 3.72008e-44
x = -10
    taylor series = 4.53999e-05
    hybrid method = 4.53999e-05
    library = 4.53999e-05
x = -1
    taylor series = 0.367879
    hybrid method = 0.367879
    library = 0.367879
x = -0.1
    taylor series = 0.904837
    hybrid method = 0.904837
    library = 0.904837
x = 0
    taylor series = 1
    hybrid method = 1
    library = 1
x = 0.1
    taylor series = 1.10517
    hybrid method = 1.10517
    library = 1.10517
x = 1
    taylor series = 2.71828
    hybrid method = 2.71828
    library = 2.71828
x = 10
    taylor series = 22026.5
    hybrid method = 22026.5
    library = 22026.5
x = 100
    taylor series = 2.68812e+43
    hybrid method = 2.68812e+43
    library = 2.68812e+43
x = 1000
    taylor series = 3.16501e+224
    hybrid method = 1.97007e+434
    library = 1.97007e+434
x = 10000
    taylor series = 2.58744e+423
    hybrid method = 8.80682e+4342
    library = 8.80682e+4342

编辑:

对于感兴趣的人:

关于原始程序中浮点错误的严重程度，评论中提出了一些问题。我最初的假设是它们可以忽略不计——我做了一个测试，看看这是不是真的。事实证明这不是真的，并且存在重大的浮点错误，但即使没有浮点错误，泰勒级数本身也会引入重大错误。 x=-100 中 200 项的泰勒级数的真实值似乎接近 -10^{24}，而不是 10^{-44}。我使用 boost::multiprecision::cpp_rational 检查了这一点，它是建立在任意精度整数类型之上的任意精度有理数类型。

输出:

x = -100
    taylor series (double) = -8.49406e+24
                (rational) = -18893676108550916857809762858135399218622904499152741157985438973568808515840901824148153378967545615159911801257288730703818783811465589393308637433853828075746484162303774416145637877964256819225743503057927703756503421797985867950089388433370741907279634245166982027749118060939789786116368342096247737/2232616279628214542925453719111453368125414939204152540389632950466163724817295723266374721466940218188641069650613086131881282494641669993119717482562506576264729344137595063634080983904636687834775755173984034571100264999493261311453647876869630211032375288916556801211263293563
                           = -8.46257e+24
    library                = 3.72008e-44
x = 100
    taylor series (double) = 2.68812e+43
                (rational) = 36451035284924577938246208798747009164319474757880246359883694555113407009453436064573518999387789077985197279221655719227002367495061633272603038249747260895707250896595889294145309676586627989388740458641362406969609459453916777341749316070359589697827702813520519796940239276744754778199440304584107317957027129587503199/1356006206645357299077422810994072904566969809700681604285727988319939931024001696953196916719184549697395496290863162742676361760549235149195411231740418104602504325580502523311497039304043141691060121240640609954226541318710631103275528465092597490136227936213123455950399178299
                           = 2.68812e+43
    library                = 2.68812e+43

代码:

#include <cmath>
#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>

typedef unsigned int uint;
typedef boost::multiprecision::cpp_rational rational;

// Taylor series of exp

template <typename T>
T taylor_series(const T x) {
    T sum = 1, last = 1;

    for (uint i = 1; i < 200; i++) {
        last = last * (x / i);
        sum = sum + last;
    }
    return sum;
}

void sample(const int x) {
    std::cout << "x = " << x << std::endl;
    long double e1 = taylor_series(static_cast<long double>(x));
    std::cout << "\ttaylor series (double) = " << e1 << std::endl;
    rational e2 = taylor_series(static_cast<rational>(x));
    std::cout << "\t            (rational) = " << e2 << std::endl;
    std::cout << "\t                       = " << static_cast<long double>(e2) << std::endl;
    std::cout << "\tlibrary                = " << expl(static_cast<long double>(x)) << std::endl;
}

int main() {
    sample(-100);
    sample(100);
}