Taylor Series Разница между exp (-x) и exp (+ x)

Question

Я пытаюсь написать программу, которая вычисляет серию Тейлора exp (-x) и exp (x) до 200 итераций, для больших x , (Ехр (х) = 1 + х + х^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

Я не думаю, что это действительно имеет отношение к ошибке приближения с плавающей запятой, я думаю, что есть еще один важный источник ошибок.

Как вы сами говорите, ваш метод очень прост. Вы принимаете приближение серии taylor на 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 на два и вызывать функцию рекурсивно, а затем возводить в квадрат результат. Или что-то вроде этого.

Для получения наилучших результатов вы, вероятно, должны использовать установленный метод.

Edit:

я решил написать код, чтобы увидеть, насколько хорошо на самом деле работает моя идея. Он, по-видимому, идеально сочетается с реализацией библиотеки C в диапазоне x = -10000 до x = 10000, по крайней мере, до такого количества бит точности, как я показываю. :)

Обратите также внимание на то, что мой метод является точным даже для значений 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 

Edit:

Ибо кто заинтересован:

Был некоторый вопрос, поднятый в комментариях о том, насколько значительным ошибки с плавающей точкой в ​​оригинале программа. Мое первоначальное предположение заключалось в том, что они были незначительными - я сделал тест, чтобы убедиться, что это правда или нет. Оказывается, это неверно, и происходит значительная ошибка с плавающей запятой, но даже без ошибок с плавающей запятой есть значительные ошибки, которые вводятся только по серии Тейлора. Истинное значение ряда Тейлора при 200 членах в x=-100, по-видимому, находится рядом с -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); 
} 

Answer 2

Одной из наиболее распространенных операций, вызывающих ошибку округления, является вычитание почти равных чисел. При вычислении exp (-100) члены более высокого порядка почти равны, поэтому ошибка округления может быстро накапливаться. Это усугубляется тем фактом, что более поздние итерации работают вне диапазона длинного двойного типа. Согласно microsoft, диапазон длинного двойника составляет от 1.7E +/- 308 (15 цифр), а 200! составляет 7,89 × 10 374.

Если вы просто ищете реализацию опыта, то math.h library имеет реализацию.

Answer 3

Использование полиномов Тейлора, вероятно, не является хорошей идеей для этого; см. http://www.embeddedrelated.com/showarticle/152.php за отличную статью об использовании полиномов Чебышева для аппроксимации функций.

rickandross указал на источник ошибки в этом случае, а именно, что разложение Тейлора для exp (-100) связано с различиями больших значений.

Существует простая модификация попытки Тейлора, которая дает разумные ответы на несколько тестовых примеров, которые я пробовал, а именно использование того факта, что exp (-x) = 1/exp (x). Эта программа:

#include <iostream> 
#include <cmath> 

double texp(double x) 
{ 
    double last=1.0; 
    double sum=1.0; 
    if(x<0) 
    return 1/texp(-x); 

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

    return sum; 
} 

void test_texp(double x) 
{ 
    double te=texp(x); 
    double e=std::exp(x); 
    double err=te-e; 
    double rerr=(te-e)/e; 
    std::cout << "x=" << x 
      << "\ttexp(x)=" << te 
      << "\texp(x)=" << e 
      << "\terr=" << err 
      << "\trerr=" << rerr 
      << "\n"; 
} 

int main() 
{ 
    test_texp(0); 
    test_texp(1); 
    test_texp(-1); 
    test_texp(100); 
    test_texp(-100); 
} 

дает этот вывод (обратите внимание, что двойная точность составляет около 2e-16):

x=0 texp(x)=1 exp(x)=1 err=0 rerr=0 
x=1 texp(x)=2.71828 exp(x)=2.71828 err=4.44089e-16 rerr=1.63371e-16 
x=-1 texp(x)=0.367879 exp(x)=0.367879 err=-5.55112e-17 rerr=-1.50895e-16 
x=100 texp(x)=2.68812e+43 exp(x)=2.68812e+43 err=1.48553e+28 rerr=5.52628e-16 
x=-100 texp(x)=3.72008e-44 exp(x)=3.72008e-44 err=-2.48921e-59 rerr=-6.69128e-16