Maple Division by Zero

Question

here is my work. i cannot figure out why when i make n=3 in the Fourier series, maple spits out a division by zero error.. i need some suggestions to get by this. The fourier series for n=1 and n=2 graph perfectly. This only starts to become a problem when n=3 and on. Im thinking it may be an issue with the way i define f(x), but ive tried numerous ways to define it, and i get the same result. New to maple, so go easy on me.

here is the link to some pictures of my work:

page1: http://tinypic.com/r/34iqmfa/6    
page2: http://tinypic.com/r/2rzruvm/6 

Answer 1

All the forget() calls in the code are just there to try and clear previously stored results, so that the timing comparisons might be more fair. You can remove all the calls to forget() and time(), if you wish.

I figured that you wanted your various ffsN procedures to add from 1 to N, rather than just add the Nth term (from N to N). So I made the last argument of ffs be i=1..N rather than i=N similar to what you had.

You can approach this in a purely exact symbolic way, or you can compute the coefficients numerically. And you can also try to re-use coefficients (exact symbolic, or floats) computed in previous calls that added fewer terms. This makes a measurable performance difference if you are going to compare the plots using "the sum of first 2 terms", "the sum of the first 5 terms", etc. One way to get such re-use is by using a recursive procedure with option remember, which is done below for both exact and float coefficient approaches.

I applied simplify to the exact symbolic coefficients, which collapses them to something more expected. (Think of orthogonality.) I also used value and inert Int * Sum rather than active int and sum or add, but that was just so I could also get float coefficients without doing any symbolic work.

This runs in Maple 15. Shout if it doesn't run, and then please state your version.

restart:
f:=x->sin(3*x)/x:
an:=1/Pi*Int(f(t)*cos(n*t),t=-Pi..Pi):
bn:=1/Pi*Int(f(t)*sin(n*t),t=-Pi..Pi):
ffs:=unapply(Sum(an*cos(n*x)+bn*sin(n*x),n=1..N),[x,N]):

a0:=1/(2*Pi)*int(f(t),t=-Pi..Pi):

a0+simplify(value(ffs(x,1)));
evalf(%);

seq(a0+simplify(value(ffs(x,i))),i=0..3);

seq(evalf(a0+ffs(x,i)),i=0..3);

forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
G5:=[seq(a0+simplify(value(ffs(x,i))),i=0..5)]:
plot(G5,x=-Pi..Pi);
time()-st;

forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
H5:=[seq(evalf(a0+ffs(x,i)),i=0..5)]:
plot(H5,x=-Pi..Pi);
time()-st;

ffsrecursive:=proc(x,N) option remember, system;
  if N<=0 then return evalf(a0);
  else
     return procname(x,N-1)+evalf(eval(an*cos(N*x)+bn*sin(N*x),n=N));
  end if;
end proc:

ffsrecursive(x,3);

ffsrecursiveexact:=proc(x,N) option remember, system;
  if N<=0 then return a0;
  else
     return procname(x,N-1)+simplify(value(eval(an*cos(N*x)+bn*sin(N*x),n=N)));
  end if;
end proc:

ffsrecursiveexact(x,3);

forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
S5exact:=[seq(ffsrecursiveexact(x,i),i=0..5)]:
plot(S5exact,x=-Pi..Pi);
time()-st;

forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
S5:=[seq(ffsrecursive(x,i),i=0..5)]:
plot(S5,x=-Pi..Pi);
time()-st;

# difference w.r.t. f(x)
forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
plot([seq(ffsrecursive(x,i)-f(x),i=15..19)],x=-Pi..Pi);
time()-st;

# difference w.r.t. f(x)
forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
plot([seq(a0+simplify(value(ffs(x,i)))-f(x),i=15..19)],x=-Pi..Pi);
time()-st;

# difference w.r.t. f(x)
forget(evalf): forget(int): forget(`evalf/int`):
forget(Si): forget(Ci): forget(ffsrecursive): forget(value):
st:=time():
H15_19diff:=[seq(a0+evalf(ffs(x,i))-f(x),i=15..19)]:
plot(H15_19diff,x=-Pi..Pi);
time()-st;

[addendum] The submitter asked for just one simple approach. Here is the exact, non-recursive approach (of the 4 above).

restart:

f:=x->sin(3*x)/x:
an:=1/Pi*Int(f(t)*cos(n*t),t=-Pi..Pi):
bn:=1/Pi*Int(f(t)*sin(n*t),t=-Pi..Pi):

ffs:=unapply(Sum(an*cos(n*x)+bn*sin(n*x),n=1..N),[x,N]):

a0:=1/(2*Pi)*int(f(t),t=-Pi..Pi):

G5:=seq(a0+simplify(value(ffs(x,i))),i=0..5);

plot([f(x), G5],x=-Pi..Pi,
     legend=[f(x),"n=0","n=1","n=2","n=3","n=4","n=5"],
     color=[black,gold,cyan,green,blue,magenta,red]);