# FOURIER SERIES

#  In this notebook, we illustrate how Maple can be used to carry out
#  basic computations with periodic functions and Fourier series.
#  References in the notebook correspond to the book Partial
#  Differential Equations and Boundary Value Problems, by Nakhle Asmar,
#  Prentice Hall, 2000.

# Example 1 The Sawtooth function (Section 2.2, Example 1)

# Consider the function 
#   f(x) = (Pi-x)*(1/2)  if 0 < x, x < 2*Pi 
# and  
#   f(x+2*Pi) = f(x) otherwise.  
# To compute its Fourier series, all we need is the function's formula 
# on the interval [0, 2*Pi];. We will use the Euler formulas (5)--(6),
# Section 2.2.

f:= x->1/2*(Pi - x);

a:= n->1/Pi *int(f(x)*cos(n*x), x= 0.. 2*Pi);
b:= n->1/Pi *int(f(x)*sin(n*x), x= 0.. 2*Pi);
a(0) = 1/(2*Pi)*int(f(x), x= 0.. 2*Pi);

# To see the results, we will ask Maple to evaluate these coefficients.

a(0);
a(n);
b(n);

# Since Maple does not know that n; is an integer, we use the assume
# command to force the issue. The evaluations of the a(n) and b(n) will
# contain n~ because Maple is showing that there are assumptions on n.

assume('n',integer);
eval(a(0));
eval(a(n));
eval(b(n));

# As expected, a(n)=0 and b(n)=1/n.  We can now form the partial sums of
# the Fourier series:

S:=(N,x)->sum( sin(n*x)/n,n=1..N);

# We know that the partial sums of the Fourier series converge to the
# function. This important fact can be illustrated graphically as
# follows looking at the sum of only 20 terms.

plot(S(20,x),x=-2*Pi..2*Pi);

# The graph of the 20th partial sum, S(20, x);, is very close to the
# graph of the function on the interval [-2*Pi;, 2*Pi;], except near the
# points x = -2*Pi;, 0, and 2*Pi;.  What is happening near these points
# is a phenomenon, known as the Gibbs phenomenon, which will be
# discussed later.  At this point, let us consider a few other ways to
# visualize the partial sums.  Here is a method for plotting the partial
# sums, one at a time, S(1, x);, S(2, x);,..., S(14, x); . We use the
# numpoints option to increase the accuracy of the graph as we add more
# terms.

for k from 1 to 14 do
 plot(S(k,x),x=-2*Pi..2*Pi,numpoints=30*k) od;

# You can create an animation of these graphs in two ways. Both require
# first loading in the plots package of commands.

with(plots):

# The first way involves two commands. First a loop to create the plots,
# and then the display command to create the animation. Notice with the
# first command we end with : to avoid seeing the formal output which
# would be essentially a list of drawing commands.

for k from 1 to 14 do
 A[k]:=plot(S(k,x),x=-2*Pi..2*Pi,numpoints=30*k) od:
A:=[seq(A[k],k=1..14)]: 
display(A,insequence=true);

# We instead can use the animate command which is faster but
# occasionally can be tricky if Maple tries to fill the range of values
# with non-integer values.

animate( S(N,x),x=-2*Pi..2*Pi,N=1..14,numpoints=300);

# As you animate the graphs, look near the original function's
# discontinuities some of which are at x = 0;, x = -2*Pi; and x = 2*Pi;
# for example. Note how the humps on the graphs of the partial sums
# nearest the discontinuities move toward the discontintuity, but their
# amplitude does not seem to diminish.  You can try adding more terms of
# the partial sums; this effect will persist.  This is known as the
# Gibbs phenomenon.  (See Section 2.2, Exercises 21-24.) Two comments
# are in order here. 1--You can ask Mathematica to plot a subsequence of
# the partial sums by repeating the previous command line and specifying
# the step increase for n.  For example, if we want to see the partial
# sums S(2,x), s(4,x), s(6,x), up to S(12,x) we would do the following.

for k from 1 to 6 do
 A[k]:=plot(S(2*k,x),x=-2*Pi..2*Pi,numpoints=30*k) od:
A:=[seq(A[k],k=1..6)]: 
display(A,insequence=true);

# Our second comment has to do with the way Maple renders a graph. 
# Since we did not specify the range of the y-axis, Maple will choose
# what it thinks are the most appropriate scales for the x and y
# coordinates. These scales may vary with the picture.  In animating
# graphics, it is convenient to fix the range on the y-axis to better
# compare the graphs.  This can be done using "PlotRange", as follows.

for k from 1 to 6 do
 A[k]:=plot(S(2*k,x),x=-2*Pi..2*Pi,numpoints=30*k,view=-2..2) od:
A:=[seq(A[k],k=1..6)]: 
display(A,insequence=true);

# In the above graphs, we have specified the y-range as {-1.7,1.7}. Now,
# let us present one more useful way of plotting.  This one uses the
# versitile " seq command.  First, you ask Maple to make a list of the
# partial sums, then you ask Maple the list.  Here is an illustration.

psums := [seq(S(k,x),k=1..5)];
plot((psums), x=-2 *Pi..2* Pi);

# We can get Maple to graph a subsequence of partial sums by adjusting
# the last example. To get the odd partial sums, we could use

psums := [seq(S(2*k-1,x),k=1..5)];
plot((psums), x=-2 *Pi..2* Pi);

# The next method for plotting can be useful some times.  We will ask
# Maple to plot the graphs in an array, as follows.

B:=animate( S(N,x),x=-3*Pi..3*Pi,N=1..12,view=-2..2,axes=none):
display(B);

# Let us end this part with further comments about these graphs.
# 1-- As we include more sine terms in the partial sums, we see more
#     oscillations on the graphs.  
# 2-- The amplitude of the oscillations tend to deviate less and less
#     from the sawtooth function, and indeed, the partial sums seem to 
#     converge to the function, confirming the Fourier Series
#     Representation Theorem (Theorem 1, Section 2.2).
# 3-- Note that all partial sums are equal to 0 at the points of
#     discontintuity. This too confirms our expectation, since at these
#     points, 0 is the average of the limits of the function from the
#     left and the right at 0 and other discontinuities.
# 4-- Near the points of discontintuity, the graphs of the partial sums
#     overshoot the graph of the function causing the Gibbs phenomenon.
#