Taylor series approximation of e^x at x =-20

I'm trying to evaluate the Taylor polynomials for the function e^x at x = -20. My results do not look right and I don't know what's wrong with my for loop. Also, I can't seem to plot my data correctly with one being the approximate and the actual one on the same graph. Here's my code:

    n = 1;
    x = -20;
    %terms = 0;
    terms = zeros(1, n+1);
    for j = 0:1:n
        terms(j+1)=(x.^j)/factorial(j);
        %display(terms)
    end;

    %display(terms)
    termsSum = sum(terms);

    display(terms)

    display('The approximate estimate at x = -20')
    display(termsSum)

    %true value at x=-20
    display('TRUE VALUE')
    display(exp(x))

    %absolute error between the approximated and true value
    display('The error: ')

    t = exp(x)
    y = termsSum
    h = abs(t-y)
    display(h)

    %plot the approximation
    %plot(-22:0.001:-18,termsSum,'o')

    plot(terms,h, '*')

Answer 1

[Matt J]

You're better off approximating exp(+20) with the Taylor series, and then taking the reciprocal of that.

    n = 50;
    x = -20;
    %terms = 0;
    terms = zeros(1, n+1);
    for j = 0:n
        terms(j+1)=(abs(x).^j)/factorial(j);

    end;

    termsSum = sum(terms);

    t = exp(x);
    y = termsSum.^sign(x);
    h = abs(t-y);

This avoids the numerical instability of summing over large terms with opposite signs. As for plotting, it's not clear to me what you're plotting things as a function of. n? x?

[cee878]

@Matt J Are you referring to the Taylor series function? Also, I was trying to plot n vs error, but I don't think I did that correctly.

[Matt J]

Something like this, then?

termsSum = cumsum(terms);

    t = exp(x);
    y = termsSum.^sign(x);
    h = abs(t-y);

    semilogy(0:n,h);

[cee878]

I wanted to use a for loop to really understand how to program the Taylor series rather than using the function. Looking at my code, which part is the one that's messing up with my iterations?

[Matt J]

Sorry, I don't follow you. Here is my complete recommendation for how to fix your code, combining everything we've talked about so far. It still has the for-loop, as you can see, although note my side comment as to how you could eliminate the loop with vectorized commands. The key is to use abs(x).^j/factorial(j), instead of just x.^j/factorial(j), and then compensate for the missing negative sign later.

    n = 50;
    x = -20;

    terms = zeros(1, n+1);
    for j = 0:n
        terms(j+1)=(abs(x).^j)/factorial(j);
    end;
    %loop could be replaced by, terms=abs(x).^(0:n)./factorial(0:n)

    termsSum = cumsum(terms);

    t = exp(x);
    y = termsSum.^sign(x);
    h = abs(t-y);
    semilogy(0:n,h);

[cee878]

What is cumsum?

[Matt J]

It is a MATLAB comand, see

http://www.mathworks.com/help/matlab/ref/cumsum.html

which computes partial sums of a vector (or along rows,columns, etc... of an array or matrix).

[cee878]

Yeah, I read about it. And I edited my code as the following:

clear all
clc
%%format bank
%Using Taylor polynomials to evaluate f(x)=e^x at x=-20
n = 100;
x = -20;
%terms = 0;
terms = zeros(1, n+1);
for j = 0:1:n
  terms(j+1)=(x.^j)/factorial(j);
end;

%display(terms)
termsSum = sum(terms); %holds the estimate

display(terms)

display('The approximate estimate at x = -20')
display(termsSum)

%true value at x=-20
display('TRUE VALUE')
display(exp(x))

%absolute error between the approximated and true value
display('The error: ')

display('THE ERROR')
display(abs(exp(x)-termsSum))

%plot the approximation and exact graph
semilogy(0:n,termsSum,'o')
semilogy(x,exp(x),'*')

Now, I want to plot all of the points in my saved vector to compare them to the actual value of e^-20, but it seems to not be working. I know I'm suppose to see a polynomial with degree 100. I only see a point when I'm plotting. Could you tell me the problem with this?

[Matt J]

As you've computed them x, exp(x), and termSum are all just scalars. To make a plot, you need a sequence of points.

[cee878]

Since terms is a vector, won't that hold a sequence of points for me?

[Matt J]

But you are not plotting your "terms" vector. I see only x, exp(x), and termSum in your calls to semilogy.

Also, I'm not sure what the benefit would be of plotting the individual terms of the Taylor series. It is the partial sums of the Taylor series - not the sequence of terms themselves - that is supposed to converge to exp(x).

Finally, I'm not sure why you're sticking with your original implementation of the for-loop. I have told you that it is numerically unstable, and so has John D'Errico in your duplicate thread. Did you try the alternative I suggested?

[cee878]

I guess I'm not quite sure what you mean by unstable, like how every term is changing signs in the Taylor polynomial changing signs? I wanted to graph the absolute error vs the "n" order of my polynomial so that I could see how at the approximated value of using Taylor polynomials at x = -20 will ever be close e^-20 regardless of how I choose n. And sorry for the trouble @Matt J, I'm kind of new to Matlab.

[Matt J]

Yes, but did you try the code I gave you several comments back? When I tested it, it seemed to do exactly what you describe.

[John D'Errico]

Matt points out the problem with computing either exp(-20) or exp(20). Both will be problematic in a series, due to the magnitude of those terms. However, there are nice ways to greatly improve the convergence of a series.

For example, you can use a simple identity or two to fix things. Consider this one:

exp(x) = (exp(x/k))^k

For example, if k = 2, this reduces to

exp(x) = (exp(x/2))^2 = exp(x/2)*exp(x/2)

In the case of x=-20, we can thus compute exp(-20) using k = 32.

exp(-20) = (exp(-20/32))^32

Lets see how this can be used to compute exp(-20) in efficient fashion.

xhat = -20/32;
n = 15; % 16 total terms in the series
ii = (0:n).';
t = xhat.^ii./factorial(ii)
t =
          1
     -0.625
    0.19531
   -0.04069
  0.0063578
-0.00079473
 8.2784e-05
-7.3914e-06
 5.7746e-07
-4.0101e-08
 2.5063e-09
 -1.424e-10
 7.4169e-12
-3.5658e-13
 1.5919e-14
-6.6329e-16

sum(t)
ans =
      0.53526

format long g
sum(t)
ans =
        0.53526142851899

So the above is a very good approximation to exp(-20/32). Raise it to the 32nd power...

sum(t)^32
ans =
    2.06115362243854e-09

exp(-20)
ans =
    2.06115362243856e-09

and we got almost full precision, using only 16 total terms. Had I used a larger value of k, the convergence would have been faster.

There are lots of ways of accelerating convergence of such a series. Generally, if we can make x smaller (closer to zero), then the series will converge with fewer terms.

[Matt J]

Note that my proposed method is a special case of John's method with k=-1. Running that case, you should see that it works, but it took me large n (~100) to get good precision.