Taylor Series for Exponential Function exp(-x)

Question

I have been working on a program for Taylor Series and used long double as the number format to allow calculation of large numbers. My program works pretty fine for positive exponents, yet it fails when it comes to negative exponents. The problem is that I am getting very large positive number when I calculate exp(-x) for some x. What may be the reason behind this? Thanks for your help in advance. You can see my code here:

#include <stdio.h>
#include <math.h>
//We need to write a factorial function beforehand, since we
//have factorial in the denominators.
//Remembering that factorials are defined for integers; it is
//possible to define factorials of non-integer numbers using
//Gamma Function but we will omit that.
//We first declare the factorial function as follows:
long double factorial (double);
//Long long integer format only allows numbers in the order of 10^18 so 
//we shall use the sign bit in order to increase our range.
//Now we define it,
long double
factorial(double n)
{
//Here s is the free parameter which is increased by one in each step and
//pro is the initial product and by setting pro to be 0 we also cover the
//case of zero factorial.
    int s = 1;
    long double pro = 1;
    if (n < 0)
        printf("Factorial is not defined for a negative number \n");
    else {
    while (n >= s) { 
    pro *= s;
    s++;
    }
    return pro;
    }
}

int main ()
{
    long double x[13] = { 1, 5, 10, 15, 20, 50, 100, -1, -5, -10, -20, -50, -100};
//Here an array named "calc" is defined to store 
//the values of x.

//int k;
////The upper index controls the accuracy of the Taylor Series, so
////it is suitable to make it an adjustable parameter. 
int p = 135;
long double series[13][p];
long double sum = 0;
int i, k;
for (i = 0; i <= 12;i++) {
for (k = 0; k <= p; k++){
    series[i][k] = pow(x[i], k)/( factorial(k));
    sum += series[i][k];
}
printf("Approximation for x = %Lf is %Lf \n", x[i], sum);
}
printf("%Lf \n", factorial(100));
}

Answer 1

You're not zeroing out your sum. You're just adding each new result to the previous.

Add sum = 0; as the first statement in the first for loop.

Moral: Always create a function for what should be a function. In this case, write a separate exp_taylor() function, rather than just writing it inline.

Answer 2

That's just the mathematical subject of numerical analysis for you. The MacLaurin series for e^x converges for all x, but let's see why it isn't useful for e^(-10).

e^x = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/720 + x^7/5040 + ... +x^n/n! + ...

e^(-10) = 1 - 10 + 100/2 - 1000/6 + 10000/24 -100000/120 + ...

What is the largest term in the series? 10^10/10!, which is approximately 2755.7319224. What is the true value of e^(-10) Approximately 0.00004539992. Adding up the series loses 9 digits of precision along the way, which you simply do not have.

Had you found e^(10) and taken the reciprocal, you would have been fairly safe. Had you computed e^(-10) directly by multiplying (1/e) 10 times, you would also be safe. But any series with alternating terms that can be far greater in magnitude than the true answer will cause these problems.

Even for functions with limited ranges, MacLaurin series are not used in practice. For example, one first takes the argument of a trig function and uses periodicity and trig identities to reduce the argument to the interval 0 < θ < π/4. Then, one often applies Chebychev approximation to reduce the error evenly. In other situations, continued fractions and Pade approximants are better than Trigonometric series. Bessel functions are best done by backwards recursion.

Look at a good numerical analysis book. Forman Acton's Numerical Methods that Usually Work is old-fashioned, but good.

Answer 3

You need to reset your sum to zero after each series is computed:

int main ()
{
  long double x[13] = { 1, 5, 10, 15, 20, 50, 100, -1, -5, -10, -20, -50, -100}; 
  const int p = 135;
  long double series[13][p];
  int i, k;
  for (i = 0; i <= 12;i++) {
    long double sum = 0;  // LOOK HERE
    for (k = 0; k <= p; k++){
      series[i][k] = pow(x[i], k)/( factorial(k));
      sum += series[i][k];
    }
    printf("Approximation for x = %Lf is %Lf \n", x[i], sum);
  }
}

Also note that you are expanding about x=0 and can expect significant error for x far from 0.