General Taylor Polynomial code

I'm trying to create a general code for evaluating Taylor polynomials for f(x)=x-sin(x)/x^3 and plot their relative error vs n iterations. I ran the code but I kept getting a N must be non-negative integers. Now, I'm confused. Also, I'm pretty new at Matlab. Some help would be great.

function[] = Taylor4(n,x)
n = input('Degree: ');
x = input('Value of x: ');
f=(1/x^2)-(sin(x)/x^3);
terms = zeros(1, n+1);
for j = 0:1:n
    terms(j+1)=(-1).^(j+1)*(x.^(2j-2))/factorial(2j+1);
end;

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

display('The iteration: ')
display(terms)

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

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

%absolute error between the approximated and true value
t = exp(x);
y = termsSum.^sign(x);
h = abs(t-y);
display('THE RELATIVE ERROR: ')
display(abs(f-termsSum))

%plot the Taylor approximation and exact value
semilogy(0:n,h); hold on %the Taylor polyonmial estimate when n = 10
semilogy(x,exp(x),'*') %exact value at x = -20

[John D'Errico]

Why would you create a function that has input arguments of n and x, and THEN use input on those same variables inside the function???????????

As far as your error goes, N never appears in that code. So we cannot guess what you did.

Show how you called the function, else how can we help you?

[Star Strider]

General Zachary Taylor, Hero of the Mexican War?

Answer 1

[Walter Roberson]

You have

      terms(j+1)=(-1).^(j+1)*(x.^(2j-2))/factorial(2j+1);

In MATLAB, the syntax 2j+1 means complex(1,2) -- that is, i or j immediately following a number indicates the number is to be multiplied by the imaginary constant to form a complex number.

MATLAB has no implicit multiplication of variables. If you want to multiply a variable by a value, you need to use either * or .* to indicate the multiplication.

Answer 2

[John D'Errico]

Other questions abound. Do you really want to approximate the function:

f(x)=x-sin(x)/x^3

That is what you wrote after all. Or do you want

f(x) = (1/x^2)-(sin(x)/x^3)

Which is what your code shows?

I'll assume the latter, since it is more interesting. Perhaps you just don't appreciate the need for parens in the first case.

In that event, then a quick check of your series approximation suggests you got that wrong.

syms x
f=(1/x^2)-(sin(x)/x^3);

taylor(f,'order',8)
ans =
- x^6/362880 + x^4/5040 - x^2/120 + 1/6

If we look at the j==0 term in your series approximation...

(-1).^(j+1)*(x.^(2*j-2))/factorial(2*j+1)

(after fixing the 2j issue) when j == 0, this reduces to

-1/factorial(1)

Clearly not 1/6. So you want to carefully re-derive that series approximation. Yes, I could do it for you, but that would take all the fun out of it for you. :)

[cee878]

My intention was to approximate f using Taylor polynomials and look at the exact value. When I derived the approximation on paper that's what I got for the for-loop when j=1 is the starting point.

[cee878]

Also, I fixed the j issue (thanks for the catch), and switched the starting point for-loop. However, after the first iteration I started to get * for my outputs. What am I coding incorrectly?

function[] = Taylor4(n,x)
n = input('Degree: ');
x = input('Value of x: ');
f=(1/x^2)-(sin(x)/x^3);
terms = zeros(1, n+1);
for j = 1:1:n
    terms(j+1)=(-1).^(j+1)*(x.^(2*j-2))/factorial(2*j+1);
end;

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

display('The iteration: ')
display(terms)

display('The approximate estimate: ')
display(termsSum1)

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

%absolute error between the approximated and true value
t = exp(x);
y = termsSum.^sign(x);
h = abs(t-y);
display('THE RELATIVE ERROR: ')
display(abs(f-termsSum))

%plot the Taylor approximation and exact value
semilogy(0:n,h); hold on %the Taylor polyonmial estimate when n = 10
semilogy(x,f,'*') %exact value at x = -20

[John D'Errico]

Ok, so lets first check your series. I'll let MATLAB do the work.

Note that since the loop now starts at 1, there is no need to offset your index into terms.

n = 4;
syms x
terms = zeros(1, n,'sym');
for j = 1:1:n
  terms(j)=(-1).^(j+1)*(x.^(2*j-2))/factorial(2*j+1);
end;

terms
terms =
[ 1/6, -x^2/120, x^4/5040, -x^6/362880]

Ok. So that worked very nicely, and it was consistent with the series that Taylor yields. So nice when mathematics works. :)

Next, you say that you are getting something strange for output. But I don't know what values you are using for x or n.

[John D'Errico]

So, looking further. It appears that you have used the code from your last assignment. BAD IDEA, at least if you are going to be sloppy about it. One thing I learned very early in my career was to be an absolute perfectionist in my code. (That does pay off when you are writing software.)

Did you really want to have this line in there?

t = exp(x);

That is a holdover from your last assignment.

Anyway, you will have the very same problem here as you did in the last assignment.

For example, if x = -20 again, then what happens with those terms?

terms
terms =
  0.16667  -3.3333   31.746  -176.37  641.33

Again, they explode. Not as fast as the exponential series, because you divided by x^2 in this problem, but you will need to go moderately far out before the terms start to approach zero.

So, go back to your code, and fix these problems. With x=-20, it appears as if you would need to go out to at least 33 or more terms to get anywhere close to double precision accuracy. So, after 33 terms, I get this when x = -20:

sum(terms)
ans =
       0.00238588184315855

f
f =
     0.00238588184365905

To test the result, lets see how well we did.

x = sym(-20);
vpa(x^-2 - sin(x)/x^3)
ans =
0.0023858818436590465432029875020193

f and the symbolic result are very close, with only the very last reported digit in error. But the series is off by quite a bit. Again the problem becomes subtractive cancellation error in the series. We can never achieve more than about 10 significant digits with that series for x=-20, at least in double precision.

(f - sum(terms))/f
ans =
      2.09775790119139e-10

[cee878]

Okay, I fixed it. However, conceptually I still don't really get what's the difference between "cumsum" and "sum". Could you explain that a little in my code?

[Walter Roberson]

v = [2 5 8];
sum(v)

would be 2+5+8 which would be the single value 15.

cumsum(v)

would be [(2), (2+5), (2+5+8)] = [2 7 15], a vector, in which each position is the cumulative total of the original vector, [sum(v(1)), sum(v(1:2)), sum(v(1:3)), ...]