Supplement for Lab 5: Fourier Series restart; with(plots): In this week's lab, you will learn to solve for the coefficients of the Fourier Series (all three versions) when fitting various functions. Many of the techniques we present here were based on a Maple Worksheet on Fourier Series developed by Anton Dzhamay from the University of Northern Colorado (http://www.unco.edu/adzham/) in 2004. Any credit for these techniques should be given to Dr. Dzhamay while any criticisms for typos should be directed at Juan Cabanela.
<Text-field style="Heading 1" layout="Heading 1"></Text-field> Implementing the Fourier Cosine Series In lecture we discussed the "full" Fourier Series which consists of both a cosine term and a sine term. In this supplement we will focus on just the Fourier cosine series for simplicity and let you worry about sine series in the lab itself. We start by setting up the variables n and m to always be treated as integers: assume(n,integer);assume(m,integer); For a function that is periodic over length L, we know that the cosine of 2\317\200n/L will keep cropping up over and over again, so rather than having to enter it over and over again, why not simply define a function c(x,n) that equals cos(n\317\200/L): c := (x,n) -> cos((2*Pi*n*x)/L); Let's be sure that stuck in memory by calling the new c function we defined: c(x,n); Notice that Maple is smart enough to apply the fact that m is an integer to simplify this function: c(L,m); We can compute the Fourier cosine coefficients for some function f(x) on the interval [0,L] using the Maple proc command to create a simple one line procedure for evaulating the given term in the Fourier cosine series. Notice that this looks a little different from what we did in lecture. I just dropped the intergral of the cosine-sine cross-term since over an entire period it always evaluates to zero and I divided by the remaining integral instead of coming up with seperate expressions for the special cases n=m=0 and so on: A:=proc(expr,var,n) simplify(int(expr*c(var,n),var=0..L)/int(c(var,n)*c(var,n),var=0..L)); end proc; Here I note that as we discussed in lecture, there are two possible "generic" expressions for the coefficients of the Fourier cosine series depending on whether LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJuRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEiMEYnRj5GPkYrRj4= or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZHLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIzAuRidGOS1GNjYtUSJ+RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSYwLjBlbUYnL0ZORldGUy1GLDYnUSVWZXJpRicvJSVib2xkR0YxRi8vRjNRLGJvbGQtaXRhbGljRicvJStmb250d2VpZ2h0R1ElYm9sZEYnLUYsNidRI2Z5RidGZm5GL0ZobkZqbi1GNjYwRlVGZm5GL0ZobkZqbkY7Rj5GQEZCRkRGRkZIRlZGWC1GLDYnUSV0aGF0RidGZm5GL0ZobkZqbkZgby1GLDYnUSR0aGVGJ0ZmbkYvRmhuRmpuRmBvLUYsNidRKWNvbXB1dGVkRidGZm5GL0ZobkZqbkZgby1GLDYnUSZ0ZXJtc0YnRmZuRi9GaG5Gam5GYG8tRiw2J1ElaGVyZUYnRmZuRi9GaG5Gam5GYG8tRiw2J1EmbWF0Y2hGJ0ZmbkYvRmhuRmpuRmBvLUYsNidRJXdoYXRGJ0ZmbkYvRmhuRmpuRmBvLUYsNidRJHlvdUYnRmZuRi9GaG5Gam5GYG8tRiw2J1Emd291bGRGJ0ZmbkYvRmhuRmpuRmBvLUYsNidRJ2V4cGVjdEYnRmZuRi9GaG5Gam5GYG8tRjY2MFElZnJvbUYnRmZuRi9GaG5Gam5GO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZkcUZgby1GLDYnUSV5b3VyRidGZm5GL0ZobkZqbkZgby1GLDYnUShsZWN0dXJlRidGZm5GL0ZobkZqbkZgby1GLDYnUSZub3Rlc0YnRmZuRi9GaG5Gam4tRjY2MFEiLkYnRmZuRi9GaG5Gam5GO0Y+RkBGQkZERkZGSEZWRlhGOQ== A(f(x),x,0); A(f(x),x,n); Now that we have a procedure for computing the nth term in the Fourier cosine series, let us write a procedure to compute the entire Fourier cosine series for a given function (note I am carefully dealing with the difference between the n=0 and n>0 solutions and using Sum instead of sum (as of Maple 12, the sum command actually tries to evaluate the sum instead of just producing a symbolic representation): cosineFS:=proc(expr,var,n) A(expr,var,0)+Sum(A(expr,var,m)*c(var,m),m=1..n); end proc; So what does the full series (to n equal infinity) look like? cosineFS(f(x),x,infinity); I'll also define a function to construct a polynomial of the first n terms in the cosine series by using the add command instead of sum: cosineFP:=proc(expr,var,n) A(expr,var,0)+add(A(expr,var,m)*c(var,m),m=1..n); end proc; The difference between sum and add is that sum tries to find a neat and tidy formula for the entire sum while add simply adds the terms up like they are numbers. For plotting the first several terms in the series, add is more useful because we don't expect the first several terms to add up to anything nice. So a polynomial consisting of the first five terms of the Fourier cosine series for f(x) can be found using: cosineFP(f(x),x,5);
<Text-field style="Heading 1" layout="Heading 1"></Text-field> Applying the Fourier Cosine Series procedure to functions. We have now written a procedure that should take any function f(x) and return the corresponding Fourier Cosine series to the nth order term. let's use it on a variety of functions and see how it works.
<Text-field style="Heading 2" layout="Heading 2"></Text-field> Fitting f(x) = 1 f:=x->1; Let's look at the nth coefficient of the Fourier cosine series for the sequence of the first 11 terms: seq(A[n]=A(f(x),x,n),n=0..10); Notice that all the coefficients except that for the n=0 term are zero. What does the A0 represent? Does this make sense? Let's actually look at the complete series cosineFS(f(x),x,infinity); Wow, that's pretty simple! Let's try something a little more involved.
<Text-field style="Heading 2" layout="Heading 2"></Text-field> Fitting g(x) = cos(2*pi*x/L)+0.5*cos(6*pi*x/L) It may not surprise you that you can fit a sum of two cosine waves with a Fourier cosine series! g:=x-> cos(2*Pi*x/L)+0.5*cos(6*Pi*x/L); L:=1; plot(g(x),x=0..L,y=-2..2); Let's again look at the first 5 coefficients of the Fourier cosine series and see if it makes sense: seq(A[n]=A(g(x),x,n),n=0..5); OK, notice we are recovering the cosine waves I added in! Let's actually look at the expression for the complete series: cosineFS(g(x),x,infinity); Wow! Maple seems to have oversimplified that! We will try constructing the series using add instead of sum: cosineFP(g(x),x,5); Here I have recovered exactly the function I put in because it is composed solely of cosine waves. This is one of the powerful aspects of Fourier series, they can pull out the distribution of "power" (amplitude) at different frequencies in any periodic function.
<Text-field style="Heading 2" layout="Heading 2"></Text-field> Fitting f(x) when it is a square wave We will call our square wave h(x) rather than f(x) to make sure Maple doesn't confuse it with the functions we have defined above. h:=x->piecewise(x<L/4,0.25,x>L/4 and x<3*L/4,0.75,x>3*L/4,0.25); I will have to set the range over which this function is periodic, so let me set L=1 and plot f(x) over that range: L := 1; plot(h(x),x=0..L,y=0..1); Let's look at the nth coefficient of the Fourier cosine series up to 10th order: seq(A[n]=A(h(x),x,n),n=0..10); OK, this is interesting, all the even terms above n=0 are zero. Why might that be? If you can understand why this is the case, you understand what is going on here! Hint: This square wave is symmetric about the center of this interval (x=0.5). Which cosine waves will show similar symmetry? Let's actually look at the expression for the complete series: cosineFS(h(x),x,infinity); And evaluate the first 5 terms: fit0 := cosineFP(h(x),x,0); fit1 := cosineFP(h(x),x,1);fit3 := cosineFP(h(x),x,3);fit5 := cosineFP(h(x),x,5); Let me try to plot this (remember I defined L = 1 earlier, so it is still assumed here). plot([h(x),fit0,fit1,fit3,fit5],x=0..L,legend=["h(x)","0th order fit","1st order fit","3rd order fit","5th order fit"]); By examining these first 5 terms, it should be clear that as we go to higher order Fourier cosine series, the fit of the entire series to this square wave is getting better and better. Let's see what happens as we keep increasing the number of terms from 0 to 20. I'm going to be sleazy and not use the animate command we taught you about previously because it doesn't allow for variable titles on plots. Instead we use the display which lets you run through a sequence of plots just like animate, but allows you to use the sprintf command to change the title to match the given order of the fit we are showing).: display( [plot(h(x),x=0..L,color=red,thickness=2,title=sprintf("Function f(x)=%a and ...",h(x))), seq(plot([h(x),cosineFP(h(x),x,n)],x=0..L,color=[red,blue],thickness=2,numpoints=1000, title=sprintf("... and its Fourier cosine series with %d terms, biggest n=%d and ...",n+1,n)),n=0..20) ],scaling=constrained,view=[0..L,0..1],insequence=true); What do you think will happen if I shift this boxcar wave so it is no longer centered in the interval 0 to L? Let's create a new function that is still periodic but not symmetric over the center of this interval by shifting this sqaure wave L/8 to the left. I'll plot the function as well, noting that we still have L = 1. hshifted:=x->piecewise(x<L/8,0.25,x>L/8 and x<5*L/8,0.75,x>5*L/8,0.25); plot(hshifted(x),x=0..L,y=0..1); Now lets try everything we tried above for the first square wave function with this new shifted one. We can start by looking at the first few coefficients of the Fourier cosine series: seq(A[n]=A(hshifted(x),x,n),n=0..10); cosineFS(hshifted(x),x,infinity); Well, the coefficients are different, let's check the quality of the fit... display( [plot(hshifted(x),x=0..L,color=red,thickness=2,title=sprintf("Function f(x)=%a and ...",hshifted(x))), seq(plot([hshifted(x),cosineFP(hshifted(x),x,n)],x=0..L,color=[red,blue],thickness=2,numpoints=1000, title=sprintf("... and its Fourier cosine series with %d terms, biggest n=%d and ...",n+1,n)),n=0..20) ],scaling=constrained,view=[0..L,0..1],insequence=true); Wow! That didn't seem to work at all! Why? If you can answer that question then you are well on your way to understanding Fourier Series.