restart; # An Example of a Fourier Series # In this notebook, you will study the Fourier series of the 2*Pi-periodic function which looks like # f(x) = x for 0 < x < 2*Pi. The graph of such a function is a sawtooth illustrated below. # The first question that comes to mind is: How do we represent this function in Maple? # The formula that we give uses the floor function, which is also called the greatest integer function. # For more on this subject, see Asmar, Section 2.1, Example 3 and Exercises 19--23. # Also, see the Mathematica Notebooks 2.2.periodic-func-1.nb and 2.2.periodic-func-2.nb. sawp:=(x,p)->x-p*floor(x/p); # sawtooth wave of period p, twp(x)=x on |x|

x; # So after reminding Maple that we intend n to be an integer, we get assume('n',integer);interface(showassumed=2);# display variable assumptions after eval a:=n->int(f(x)*cos(n*x),x=0..2*Pi)/Pi; a(0):=int(f(x),x=0..2*Pi)/(2*Pi); b:=n->int(f(x)*sin(n*x),x=0..2*Pi)/Pi; a(n); b(n); # The partial sums of the Fourier series are S:=(n,x)->a(0) + sum(a(k)*cos(k*x) + b(k)*sin(k*x),k=1..n); # For example, we compute the 4th partial sum. S(4,x); # No cosine terms appear, because the a(n)'s are 0 for n=1, 2, ... Our next step is to compare the # graph of the partial sum to the graph of the function itself. # Then we will see how well the partial sums are approximating the function f(x). # We plot the graphs of f(x) and the tenth partial sum over three periods. plot([f(x), S(10, x)],x=-3*Pi..3*Pi,color=[red,black]); # The approximation of the function by the partial sum of its Fourier series is good, # except near the points of discontinuity of f(x), at x = 0, -2*Pi, 2*Pi, -4*Pi, 4*Pi etc. # Near the points of discontinuity, the partial sums overshoot the graph of the function. # # The predictable over-shoot is called the Gibbs phenomenon. # # You can learn more about the Gibb's overshoot in Section 2.2, Exercises 21--24. # See also the Maple Worksheet 2.2.fourier-ser-1.ws. There are several other ways to # illustrate the graphs of the partial sums. We will show some now. # See the Mathematica Notebook 2.2.fourier-ser-1.nb for more methods of plotting. # The following method allows you to animate the graphics. # plots[animate]({f(x),S(N,x)},x=-3*Pi..3*Pi,N=1..15,view=-1..7,numpoints=100); # cpu killer X:=[seq(plot({f(x),S(N,x)},x=-3*Pi..3*Pi,view=-1..7,thickness=3),N=1..15)]: plots[display](X,insequence=true); # Graph animation shows the humps near the points of discontinuity moving toward the # discontinuties without diminishing in magnitude. This suggests that the humps will not # go away, or that the approximation will not improve near points of discontinuity, even if you # increase the terms in the partial sum. This illustrates the Gibbs phenomenon. # You can learn more about Gibb's over-shoot in the references given before. # To confirm this claim, let us compare the graph of f(x) with the graph of the 100th partial sum. p1:=plot(f(x),x=3*Pi..6*Pi,color=red,discont=true,thickness=3): p2:=plot(S(30,x),x=0..3*Pi,color=black,numpoints=400,thickness=3): p3:=plot(limit(f(x),x=2*Pi,left),x=0..6*Pi,color=green,thickness=2): p4:=plot(limit(f(x),x=2*Pi,right),x=0..6*Pi,color=green,thickness=4): plots[display]([p1,p2,p3,p4],title="Gibbs over-shoot is about 8.5%"); # Compute the height of the hump F:=unapply(S(100,x),x): plot(F(x),x=2*Pi-0.2..2*Pi); x0:=fsolve(diff(F(x),x)=0,x=127*Pi/64 .. 2*Pi); # Find the critical point y0:=evalf(F(x0)); # Height at the hump overshoot:=evalf(y0/limit(f(x),x=2*Pi,left)); printf("Over-shoot percentage = %a\n",(overshoot-1)*100);