I'm trying to write a code that uses a taylor series of sin(x) to estimate the value of sin(x), I'm using a while loop so that the loop runs as long as the absolute value of the previous estimate - the current estimate is greater than 1e-4.

x = input('Enter an angle in radians from 0 to 2pi: ');

run = 0;

est = 2;

previous_est = 1;

while abs(est-previous_est) > 1e-4

run = run + 1;

est(run+1) = (-1)^(run)*(x^(run+1))/factorial(run+1);

previous_est = est;

end

And then there is code to display the outputs, however for the estimates I keep getting large numbers in both the negative and positive.

Answer 1

[John D'Errico]

I think you don't understand how series work (converge/don't converge.) I think you also need to understand factorials, and double precision arithmetic. You need to understand that a series like this is the SUM of the terms, not just individual terms. Finally, you have the wrong formula for the individual terms. So...

1. What is the sine series? The n'th term S_n can be written as:

S_n = -1^(n-1)*x^(2*n-1)/factorial(2*n-1)

So, when n == 1, we have

-1^0*x^1/factorial(1) == x

etc.

You were using something closer to the formula for the terms in the series for exp, NOT sine. While they are related, the sine series uses essentially every other term, with an alternating sign on the terms.

2. You need to SUM the terms. Each individual term does not converge to sin(x). It is the SUM of the terms.

3. Next, you need to understand that for large values of x (x a bit larger than 2*pi), that these terms may become very large. So you might see massive subtractive cancellation happening, so that the series may never converge. Here as long as you are limited to the interval [0,2*pi], it will not get that bad.

x = 2*pi;
n = 1:20;

All the terms at once:

S_n = (-1).^(n-1).*x.^(2*n-1)./factorial(2*n-1)
S_n =
  Columns 1 through 13
       6.2832      -41.342       81.605      -76.706       42.059      -15.095         3.82     -0.71812      0.10423    -0.012032    0.0011309  -8.8235e-05   5.8057e-06
  Columns 14 through 20
  -3.2649e-07   1.5874e-08  -6.7384e-10   2.5191e-11  -8.3572e-13    2.477e-14  -6.5983e-16

And the cumulative sum...

cumsum(S_n)
ans =
  Columns 1 through 13
       6.2832      -35.059       46.547      -30.159         11.9      -3.1951      0.62488    -0.093246     0.010983   -0.0010482   8.2741e-05  -5.4944e-06   3.1127e-07
  Columns 14 through 20
  -1.5224e-08   6.4946e-10  -2.4376e-11    8.149e-13  -2.0824e-14   3.9452e-15   3.2854e-15

Of course, we know that sin(2*pi) is zero. But even MATLAb does not get that exactly true. Some tiny amount of slop is present.

sin(x)
ans =
  -2.4493e-16

So you can see that we need at least about 11 terms to compute an approximation to sin(2*pi), and that the first few terms are crazy in size compared to the final result. The net effect is called subtractive cancellation. Were x a bit larger, this effect would become so large that the sine series would simply never converge when computed in double precision arithmetic.

Now, I don't necessarily expect you to do this computation as I did without an explicit loop. In fact, you probably need to use a loop as one of the conditions of your assignment.

As a better example, I'll pick a number like 4.5 for x.

x = 4.5;

n = 1;
allterms = x;
tol = 1.e-4;
sineapprox = x;
while abs(allterms(end)) > tol
  n = n + 1;
  allterms(n) = -allterms(n-1)*x.^2/((2*n-2)*(2*n-1));
end
sineapprox = sum(allterms);

See that we can get the next term from the last by negating the previous term, multiplying by x^2, and then dividing by what it takes to go from factorial(2*n-3) to factorial(2*n-1).

format long g
sineapprox
sineapprox =
      -0.977531099355645

sin(4.5)
ans =
      -0.977530117665097

I was able to stop when the most recent term was less than tol, because after a point, all the terms will be rapidly decreasing, and changing in sign.

The loop needed 10 terms to get below the desired tolerance.

n
n =
    10

format short g
allterms
allterms =
        4.5      -15.188       15.377      -7.4141       2.0852     -0.38387     0.049829   -0.0048049   0.00035772  -2.1181e-05

Even for x as small as 4.5, see that the second and third term were both fairly large, pretty near each other in value, and of opposite sign. This series will never converge to full precision using doubles for x as large as 4.5. In fact, I would expect it to be correct to at most roughly 15 significant digits.