General Taylor Polynomial code

Question

cee878 2016-2-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/267702-general-taylor-polynomial-code

评论： Walter Roberson 2016-2-14

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

2 个评论
显示无隐藏无

John D'Errico 2016-2-11

编辑：John D'Errico 2016-2-11

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 2016-2-12

General Zachary Taylor, Hero of the Mexican War?

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2016-2-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/267702-general-taylor-polynomial-code#answer_209488

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

John D'Errico 2016-2-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/267702-general-taylor-polynomial-code#answer_209494

编辑：John D'Errico 2016-2-11

在 MATLAB Online 中打开

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. :)

6 个评论
显示 4更早的评论隐藏 4更早的评论

John D'Errico 2016-2-12

编辑：John D'Errico 2016-2-12

在 MATLAB Online 中打开

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 2016-2-13

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 2016-2-14

在 MATLAB Online 中打开

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)), ...]

请先登录，再进行评论。