Estaba calculando e^x usando Taylor Series y noté que cuando lo calculamos para x negativo, el error absoluto es grande. ¿Es porque no tenemos suficiente precisión para calcularlo?
(Sé que para prevenirlo podemos usar 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 ; }Por ejemplo: x=-40 el valor es 4.24835e-18 pero el programa me da 3.116952e-01. El error absoluto es ~0.311
x=-50 el valor es 1.92875e-22 el programa me da 2.041833e+03. El error absoluto es ~2041.833
El problema es causado por errores de redondeo en la fase intermedia del algoritmo. La h crece rápidamente como 40/2 * 40/3 * 40 / 4 * ... y oscila en signo. Los valores para i , h y Sum para x=-40 para iteraciones consecutivas se pueden encontrar a continuación (algunos puntos de datos se omiten por brevedad):
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 La magnitud máxima de la suma es 7e15 . Aquí es donde se pierde la precisión. El tipo double se puede representar con una precisión de aproximadamente 1e-16 . Esto da un error absoluto esperado de alrededor de 0.1 - 1 . Como la suma esperada (valor de exp(-40) está cerca de cero, el error absoluto final está cerca del error absoluto máximo de las sumas parciales.
Para x=-50 , el valor máximo de la suma es 1.5e20 lo que da el error absoluto debido a la representación finita del double en alrededor de 1e3 - 1e4 lo que está cerca de uno observado.
No se puede arreglar mucho sin cambios significativos en el algoritmo para evitar formar esas sumas parciales. Alternativamente, calcule exp(-x) como 1/exp(x) .
Para x negativa, la suma de los términos +/- alternos crea problemas de cálculo incluso en la primera suma de 1.0 + x , ya que se puede esperar que el error de la suma final sea tan malo como el bit menos significativo de 1,0 o aproximadamente 1 parte en 10 16 . Esto implica que x_min como en Exp(x_min) == 1.0e-16 es el valor computacional útil mínimo (p. ej., x sobre -36)
Una solución simple es formar una buena Exp(positive_x) y para valores negativos...
double Exp(double x) { if (x < 0) { return 1.0 / Exp(-x); } ... Una Exp(positive_x) buena (y simple) calcula los términos hasta que un term + 1.0 sigue siendo 1.0, ya que los términos pequeños adicionales no cambian la suma de manera significativa. Funciona bien para todos los x (error muy pequeño), excepto que podría usar mejoras cuando el resultado debería ser inferior a lo normal.
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; }