Thread Subject: Script for series expansion of sine

Hello folks,

To begin with, I browsed through the archives looking for relevant posts to my question and found some related topics. However, the posts were ill-structured and, appropriately, weren't answered. Hopefully I can be a little more specific.

I'm trying to write a script to calculate the series expansion of sine with the inputs x (argument) and n. Here's what I've thrown together so far:

% SINESERIES: computes sin(x) from series expansion.
% A script to evaluate the series expansion of sine with the formula:
% sin(x) = x-(x^3/3!)+(x^5/5!)-...
x = input('Enter the argument x: ');
n = input('Enter the interval n: ');
N = 1:1:n; % Creates row vector with n elements, 1 at a time.
k = 2*N-1; % For convenience with prod() function
j = N-1; % For cleanup.
sinseries = ((-1).^(j)).*((x.^(k))./(prod(k)))
sum(sinseries)

My calculations aren't turning out at all correctly, short of x = 1, n = 1. I've been toying with the equation to see if I translated it wrong, but I think the problem lies in how I'm expressing the series. If anybody has some advice, I'd greatly appreciate it. Thanks.

"Sean " <el_sean@yahoo.com> wrote in message <gsohrh$27i$1@fred.mathworks.com>...
> Hello folks,
>
> To begin with, I browsed through the archives looking for relevant posts to my question and found some related topics. However, the posts were ill-structured and, appropriately, weren't answered. Hopefully I can be a little more specific.
>
> I'm trying to write a script to calculate the series expansion of sine with the inputs x (argument) and n. Here's what I've thrown together so far:
>
> % SINESERIES: computes sin(x) from series expansion.
> % A script to evaluate the series expansion of sine with the formula:
> % sin(x) = x-(x^3/3!)+(x^5/5!)-...
> x = input('Enter the argument x: ');
> n = input('Enter the interval n: ');
> N = 1:1:n; % Creates row vector with n elements, 1 at a time.
> k = 2*N-1; % For convenience with prod() function
> j = N-1; % For cleanup.
> sinseries = ((-1).^(j)).*((x.^(k))./(prod(k)))
> sum(sinseries)
>
> My calculations aren't turning out at all correctly, short of x = 1, n = 1. I've been toying with the equation to see if I translated it wrong, but I think the problem lies in how I'm expressing the series. If anybody has some advice, I'd greatly appreciate it. Thanks.

  At first glance I would say the culprit here is prod(k). You are expecting its results to be a vector consisting of all the products from 1 up to a number that increases by 2 each time. Instead you are getting a single scalar answer consisting all the products in k. It would only work for n - 1. What you need is cumprod instead.

  You can see this if you follow each calculation step by step and compare it with how you would do it by hand.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsok1b$ruu$1@fred.mathworks.com>...
> "Sean " <el_sean@yahoo.com> wrote in message <gsohrh$27i$1@fred.mathworks.com>...
> > I'm trying to write a script to calculate the series expansion of sine with the inputs x (argument) and n. Here's what I've thrown together so far:
> >
> > % SINESERIES: computes sin(x) from series expansion.
> > % A script to evaluate the series expansion of sine with the formula:
> > % sin(x) = x-(x^3/3!)+(x^5/5!)-...
> > x = input('Enter the argument x: ');
> > n = input('Enter the interval n: ');
> > N = 1:1:n; % Creates row vector with n elements, 1 at a time.
> > k = 2*N-1; % For convenience with prod() function
> > j = N-1; % For cleanup.
> > sinseries = ((-1).^(j)).*((x.^(k))./(prod(k)))
> > sum(sinseries)

  I should have pointed out that there are more efficient ways of carrying out that computation. Suppose n = 5. You could do the computation this way:

 s = x*(1-x^2/(2*3)*(1-x^2/(4*5)*(1-x^2/(6*7)*(1-x^2/(8*9)))));

Moreover you can easily devise a for-loop that accomplishes this calculation in the same sequence of operations for a general n. You would have to work from the high end and go backwards. It obviously involves far fewer multiplications and divisions and is more accurate numerically in the bargain. In spite of the for-loop overhead, this method is bound to be faster for sufficiently large n, since it is an order n algorithm whereas your original one is order n^2.

  Another observation: Don't try this with values of x substantially larger than 1 unless you have set n very high. High values of x will cause the series to converge very slowly at first and only when the ratio x/k grows substantially smaller than 1 does it begin to converge at a decent rate. In the actual computation of the sine function, things are done in a very different manner than this.

Roger Stafford

"Sean " <el_sean@yahoo.com> wrote in message <gsohrh$27i$1@fred.mathworks.com>...
> Hello folks,
>
> To begin with, I browsed through the archives looking for relevant posts to my question and found some related topics. However, the posts were ill-structured and, appropriately, weren't answered. Hopefully I can be a little more specific.
>
> I'm trying to write a script to calculate the series expansion of sine with the inputs x (argument) and n. Here's what I've thrown together so far:
>
> % SINESERIES: computes sin(x) from series expansion.
> % A script to evaluate the series expansion of sine with the formula:
> % sin(x) = x-(x^3/3!)+(x^5/5!)-...
> x = input('Enter the argument x: ');
> n = input('Enter the interval n: ');
> N = 1:1:n; % Creates row vector with n elements, 1 at a time.
> k = 2*N-1; % For convenience with prod() function
> j = N-1; % For cleanup.
> sinseries = ((-1).^(j)).*((x.^(k))./(prod(k)))
> sum(sinseries)
>
> My calculations aren't turning out at all correctly, short of x = 1, n = 1. I've been toying with the equation to see if I translated it wrong, but I think the problem lies in how I'm expressing the series. If anybody has some advice, I'd greatly appreciate it. Thanks.

Sean,

The notes here http://www.derekroconnor.net/NA/LE/LE-2006-2.pdf explain why your code gives the wrong answers.

I hope you read and digest all 14 pages of this homework solution.

Regards,

Derek O'Connor.

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsok1b$ruu$1@fred.mathworks.com>...
> At first glance I would say the culprit here is prod(k). You are expecting its results to be a vector consisting of all the products from 1 up to a number that increases by 2 each time. Instead you are getting a single scalar answer consisting all the products in k. It would only work for n - 1. What you need is cumprod instead.

  Sean, I erred in telling you that just cumprod(k) would correct the error in your code. The prod(k) divisor is still certainly wrong, as I mentioned before. However, what you would need to rescue that code is this:

N = 1:n;
k = 2*N-1;
j = N-1;
p = cumprod(1:2*n-1);
sinex = sum((-1).^j.*x.^k./p(k));

  To make up for that lapse I'll spell out in detail the for-loop method I mentioned.

x2 = x^2/2;
s = 1;
for k = n-1:-1:1
 s = 1-x2/((2*k+1)*k)*s;
end
s = x*s;

Roger Stafford

Derek & Roger,

Thanks for your informative replies! I'm going to play around with this a bit, and get back to you.

-Sean