# 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 `and`(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. # #