In-class Demonstration: Fourier Approximations By Russell Blyth and Mike May, S.J. restart:with(plots): We start by using Maple to perform some tedious integrals for us. Verify that the functions used as basis vectors in the approximating subspace are orthogonal. Int(sin(i*Pi*x)*sin(j*Pi*x),x=-1..1) = int(sin(i*Pi*x)*sin(j*Pi*x),x=-1..1); Int(sin(i*Pi*x)*sin(i*Pi*x),x=-1..1) = int(sin(i*Pi*x)*sin(i*Pi*x),x=-1..1); Int(sin(i*Pi*x)*cos(j*Pi*x),x=-1..1) = int(sin(i*Pi*x)*cos(j*Pi*x),x=-1..1); Int(cos(i*Pi*x)*cos(j*Pi*x),x=-1..1) = int(cos(i*Pi*x)*cos(j*Pi*x),x=-1..1); Int(cos(i*Pi*x)*cos(i*Pi*x),x=-1..1) = int(cos(i*Pi*x)*cos(i*Pi*x),x=-1..1); That left us with some work still to do. Try again. This time we first tell Maple that i and j are integers (Maple assumes i and j are nonzero): assume(i,integer): assume(j,integer): interface(showassumed=0): Int(sin(i*Pi*x)*sin(j*Pi*x),x=-1..1) = int(sin(i*Pi*x)*sin(j*Pi*x),x=-1..1); Int(sin(i*Pi*x)*sin(i*Pi*x),x=-1..1) = int(sin(i*Pi*x)*sin(i*Pi*x),x=-1..1); Int(sin(i*Pi*x)*cos(j*Pi*x),x=-1..1) = int(sin(i*Pi*x)*cos(j*Pi*x),x=-1..1); Int(cos(i*Pi*x)*cos(j*Pi*x),x=-1..1) = int(cos(i*Pi*x)*cos(j*Pi*x),x=-1..1); Int(cos(i*Pi*x)*cos(i*Pi*x),x=-1..1) = int(cos(i*Pi*x)*cos(i*Pi*x),x=-1..1); We need to also integrate against 1: Int(sin(i*Pi*x)*1,x=-1..1) = int(sin(i*Pi*x)*1,x=-1..1); Int(cos(i*Pi*x)*1,x=-1..1) = int(cos(i*Pi*x)*1,x=-1..1); Int(1*1,x=-1..1) = int(1*1,x=-1..1); To compensate for the value of this final integral, we must remember to half the coefficient of this constant function later. Graph several of these integrands as a reality check - note that in each case it appears that the signed area from x=-1 to x=1 has the appropriate value (0 or 1): plot(sin(2*Pi*x)*sin(3*Pi*x),x=-1..1); plot(sin(3*Pi*x)*cos(4*Pi*x),x=-1..1); plot(cos(5*Pi*x)*cos(3*Pi*x),x=-1..1); plot(sin(3*Pi*x)*sin(3*Pi*x),x=-1..1); Example Find Fourier approximations for f(x) = x. We compute the Fourier coefficients for f(x) = x Int(x*sin(j*Pi*x),x=-1..1) = int(x*sin(j*Pi*x),x=-1..1); Int(x*cos(j*Pi*x),x=-1..1) = int(x*cos(j*Pi*x),x=-1..1); Int(x,x=-1..1) = int(x,x=-1..1); Note there are only sine terms (not surprisingly, since f is an odd function). Next use Maple to write a general expression for the order n approximation Fapprox := (x,n) -> sum(-2*(-1)^j/(j*Pi)*sin(j*Pi*x),j=1..n); We can now plot (say) a 10 term approximation against f(x). eval(Fapprox(x,10)); plot([Fapprox(x,10),x],x=-2..2,y = -2..2, title = " 10th Order Fourier approximation to y =x", titlefont = [HELVETICA,12], color = [green,blue]); We can also animate this plot. mindeg := 2: degsteps := 20: bydeg :=4: A := display(seq(plot([Fapprox(x,mindeg + bydeg*i),x], x=-2..2,y=-2..2, title="degree = "||(mindeg + bydeg*i), titlefont=[HELVETICA, 14]), i=0..degsteps), insequence=true): B := animate(func,x=-2..2,y=-2..2,frames=degsteps+1): print(" x vs its Fourier approximation"); display(A,B,view=[-2..2,-2..2]); It is also useful to plot the error between the function and the approximation. plot(Fapprox(x,10)-x,x=-1..1,y=-1..1, title = "error in 10th Order Fourier approximation to y = x", titlefont = [HELVETICA,12], color = red); Similarly, we can animate the plot of the error function. mindeg := 2: degsteps := 20: bydeg :=4: print("Error in nth order Fourier approximation to y=x"): display(seq( plot(Fapprox(x,mindeg + bydeg*i)-x, x=-1..1, y=-1..1, color=blue, title="order = "||(mindeg + bydeg*i), titlefont=[HELVETICA, 14]), i=0..degsteps), insequence=true); Higher order approximations match the function f even more closely!