Ряд Тейлора Разница между exp(-x) и exp(+x)
Я пытаюсь написать программу, которая вычисляет ряд Тейлора 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
Как видите, он работает просто для положительного расчета, но не для отрицательного. У кого-нибудь есть идеи о том, почему ошибки округления могут пойти не так, просто добавляя отрицание на каждый второй термин? Кроме того, есть ли что-нибудь, что я мог бы реализовать, чтобы решить эту проблему?
Заранее спасибо!!
3 ответа
Я не думаю, что это на самом деле имеет какое-либо отношение к ошибке аппроксимации с плавающей точкой, я думаю, что есть еще один более значительный источник ошибок.
Как вы говорите сами, ваш метод действительно прост. Вы принимаете приближение ряда Тейлора в 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 на два и вызывая функцию рекурсивно, затем возводя в квадрат результат. Или как то так.
Для достижения наилучших результатов вам, вероятно, следует использовать установленный метод.
Редактировать:
Я решил написать код, чтобы увидеть, насколько хорошо работает моя идея. Похоже, что он идеально соответствует реализации библиотеки 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
Редактировать:
Для кого это интересно:
В комментариях был задан вопрос о том, насколько значительны ошибки с плавающей запятой в исходной программе. Мое первоначальное предположение состояло в том, что они были незначительны - я сделал тест, чтобы проверить, правда это или нет. Оказывается, это не так, и происходит значительная ошибка с плавающей запятой, но даже без ошибок с плавающей запятой существуют существенные ошибки, вносимые только одной серией Тейлора. Истинное значение ряда Тейлора на 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);
}
Использование полиномов Тейлора, вероятно, не очень хорошая идея для этого; см. http://www.embeddedrelated.com/showarticle/152.php для большой статьи об использовании многочленов Чебышева для приближения функций.
Рикандросс указал источник ошибки в этом случае, а именно, что расширение Тейлора для 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
Одной из наиболее распространенных операций, приводящих к ошибке округления, является вычитание почти равных чисел. При вычислении exp(-100) члены высшего порядка почти равны, поэтому ошибка округления может быстро накапливаться. Это усугубляется тем фактом, что более поздние итерации работают вне диапазона типа long double. По данным Microsoft, диапазон длинных двойных составляет от 1,7E +/- 308 (15 цифр), тогда как 200! 7,89 х 10^374.
Если вы просто ищете реализацию exp, у библиотеки math.h есть реализация.