0

I'm new here.

Description of scenario:

I have written MATLAB code to use the partial sum S_n(x)=Sum from j=0 to n of (-1)^j*x^(2j+1)/(2j+1)! to approximate sin x. I have written a computer program to evaluate this partial sum by methods

LS: computing and summing the terms from the largest term first to the smallest term last

and

SL: doing the calculation in the reverse order.

I have used my program to calculate S_n(x) for x = 0.1, 1, and 10 by both methods LS and SL, and compared the results against sin x. I have used n = 10, 100, and 1000.

Question: I am running into a problem in the last two lines of my code. It says "NaN" (Not a number), but also shows deviant numbers in the previous column. What am I doing wrong? Can anyone help me?

Here is my code:

The first function:

function ret = hw1_6_1(x,n)
ret=((-1)^n)*((x^((2*n)+1))/(factorial((2*n)+1)));
end

The second function:

function ret = partialsum(x,n,log)
%% Begin code
% j=n
% 1 to n+1
% LS 0 to n
% SL n to 0
%((-1)^n)*((x^((2*n)+1))/(factorial((2*n)+1)))
clear sum
ret = 0;

if log == 1
    for i=0:1:n
    ret= ret + hw1_6_1(x,i) ;
    i=i+1;
    end
elseif log == 0
    for i=n:-1:0
    ret= ret + hw1_6_1(x,i) ;
    i=i+1;
end
end
end

Finally, my main code:

%% Hw 1 Problem 6

% Running approximation of sin x
% LS and SL
% LS == log=1, SL == log=0
% x = 0.1, 1, 10
% n = 10,100,1000
clear all
clc
format long %displays more decimal points

%% For x = 0.1

sin1 = sin(0.1);

% LS 
% n = 10, 100, 1000 Generated in array format

LS1 = [partialsum(0.1,10,1);partialsum(0.1,100,1);partialsum(0.1,1000,1)];

% Now SL

SL1 = [partialsum(0.1,10,0);partialsum(0.1,100,0);partialsum(0.1,1000,0)];

%% For x = 1

sin2 = sin(1);

% LS
% n = 10, 100, 1000 Generated in array format

LS2 = [partialsum(1,10,1);partialsum(1,100,1);partialsum(1,1000,1)];

% Now SL

SL2 = [partialsum(1,10,0);partialsum(1,100,0);partialsum(1,1000,0)];

%% For x = 10

sin3 = sin(10);

% LS
% n = 10, 100, 1000 Generated in array format

LS3 = [partialsum(10,10,1);partialsum(10,100,1);partialsum(10,1000,1)];

% Now SL

SL3 = [partialsum(10,10,0);partialsum(10,100,0);partialsum(10,1000,0)];

%% Comparison stage

Sines = [sin1;sin2;sin3]

Approxs = [LS1 SL1 LS2 SL2 LS3 SL3]'  

My output:

    Sines =
   0.099833416646828
   0.841470984807897
  -0.544021110889370


Approxs =
   0.099833416646828   0.099833416646828   0.099833416646828
   0.099833416646828   0.099833416646828   0.099833416646828
   0.841470984807897   0.841470984807897   0.841470984807897
   0.841470984807897   0.841470984807897   0.841470984807897
   2.761090925979680  -0.544021110889270                 NaN
   2.761090925979687  -0.544021110889190                 NaN

Thanks in advance.

  • 4
    (k) * 10^1000 / 1000! = Inf / Inf = NaN – Aki Suihkonen Sep 30 '13 at 04:35
  • What about the 2.761090925979680 -0.544021110889270? –  Sep 30 '13 at 04:36
  • 3
    2.76 is because 10th degree polynomial diverges from sine at much earlier point -- try to plot sin(0:0.01:10) vs. LS(10,10,0:0.01:10); (May require vectorizing). – Aki Suihkonen Sep 30 '13 at 04:45
  • 2
    Just a comment: Try to avoid using `log` as a variable name, since it is actually the name of a [built-in function](http://www.mathworks.com.au/help/matlab/ref/log.html) – Colin T Bowers Sep 30 '13 at 05:29
  • @Anonymous - on the same note as [Colin](http://stackoverflow.com/questions/19086471/matlab-code-glitch-at-the-end#comment28216709_19086471), it is best [not to use `i` as a variable name in Matlab](http://stackoverflow.com/questions/14790740/using-i-and-j-as-variables-in-matlab) as well. – Shai Sep 30 '13 at 06:14

1 Answers1

2

1) sin_taylor(10) == 2.76 is because 20th degree polynomial approximation diverges from sine at much smaller value of x.

2) The sequence x^n / n! for |x|>1, is not monotonic, but has a maximum at about n=6; thus there will be some rounding differences compared to cases, where |x|<= 1 and one is summing up a rapidly vanishing sequence.

3) Both the nominator and the denominator (10^n, n!) will be represented as infinities, when n is large. Inf / Inf is defined as NaN, and NaN + anything == NaN. This is different to partial_sum(10, 100), where one sums up terms (x / inf == 0).

Aki Suihkonen
  • 19,144
  • 1
  • 36
  • 57