Linear Methods of Applied Mathematics

Calculations of some Fourier Series

(c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod.     All rights reserved.

Notes for the instructor.

This contains calculations and examples which correlate with chapter 4 of the WWW text by Harrell and Herod.

Students can be encouraged to cut and paste from this notebook to do homework.

Instructions

This notebook uses Mathematica  to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics.   The student needs only a basic knowledge of Mathematica  to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica  to work other problems of the student's own design.  

Calculations will be performed when the reader presses enter in a given calculation cell (bold print).  It is best to activate the cells in order, so that Mathematica  will be able to call on operators and functions defined in earlier cells.  Red color coding is used to warn the reader when a given calculation relies on an earlier one.

A useful substitution, which we shall often make, is:

TrigId = {Cos[Pi n_] -> (-1)^n, Sin[Pi n_] -> 0};

Some illustrative examples

A square pulse

Let's evaluate the Fourier series for the functions:
    f(x) = 1 for 0 <= x <= L/2, but 0 for L/2 <= x <= L.
and
     g(x) = x, 0 <= x <= L.  

We want to represent these functions in the form
a[0] + Sum[a[m] Cos[2 Pi m x/L, {m,1,N}] + Sum[b[n] Sin[2 Pi n x/L, {n,1,N}].  By the way, there is a procedure in Mathematica called Fourier, but it computes the discrete Fourier transform, which is a related but different operation.

The formulae for these coefficients are:

    a[0] := (1/L) Integrate[f[x], {x,0,L}]    (the average of f)
    a[m] := (2/L)  Integrate[f[x] Cos[2 Pi m x/L], {x,0,L}], m = 1, 2, ...
    b[n] := (2/L)  Integrate[f[x] Cos[2 Pi n x/L], {x,0,L}], n = 1, 2, ...   

We simplify matters by setting L = 1.

f[x_] := 1 /; 0 <= x <= 1/2
f[x_] := 0 /; 1/2 < x <= 1

Plot[f[x], {x, 0, 1}]

a[0] := Integrate[f[x], {x,0,1}]

a[0]

a[m_] := 2 Integrate[Cos[2 m Pi x], {x,0,1/2}]

a[m] /. TrigId

These coefficients are all zero.  Why?

b[n_] := 2 Integrate[Sin[2 n Pi x], {x,0,1/2}]

b[n] /. TrigId

Simplify[%]

If n is even, b[2 k] = 0. Thus only odd terms survive, in which case 
b[2 k + 1] = 2/(Pi (2 k + 1)). We can write the full series:

Clear[FullSeries]

FullSeries[x_,N_] := 1/2 +       Sum[2 Sin[2 Pi n x]/ (Pi n),{n,1,N,2}]

Explanation:  The list {n,1,N,2} tells Mathematica to increment n by 2 until it reaches N.  This takes care of the two forms of b[n] in an efficient way.

Plot[{FullSeries[x,1], f[x]}, {x,0,1}]

Plot[{FullSeries[x,3], f[x]}, {x,0,1}]

Plot[{FullSeries[x,7], f[x]}, {x,0,1}]

On the excitation of mechanical resonances

Suppose that experiments by K. Battle at the Wiener Staatsoper in Austria show that a crystal goblet will shatter if the intensity of the tone at frequency 1760 Hz (high A) exceeds .01.  (We use physical units in which the proportionality between the power and the square of the amplitude is 1, i.e., we define the intensity simply by a[k]^2 + b[k]^2.)  At Georgia 
Tech, not a particularly musical spot, we can generate a square pulse with amplitude A for 1/704 sec., then amplitude 0 for the same length of time, then A, etc., periodically with period 1/352.  What amplitude A will cause the glass to shatter in our laboratory?

The mapping from functions to their Fourier coefficients is linear.  Hence the coefficients can be obtained from the ones we just calculated. 

Step 1.  Scale the independent variable by replacing x with t and L with 1/330. 
Step 2.  Multiply all the coefficients by A.

The harmonic with frequency 1760 is the one with n=5, so the Fourier coefficients corresponding to this frequency are 
         a[5] = 0, and 
         b[5] = 2A/(5 Pi)
The intensity is 4 A^2/(25 Pi^2).  The intensity above which the glass will shatter is

A = Sqrt[25 Pi^2 / 400]

N[%]

A sawtooth function

For comparison, let us find another Fourier series, namely the one for the periodic extension of f(x) = x, 0 <= x <= L, sometmes designated x mod L.  Thus for x between L and 2L, the function is (x-L), for x between 2L and 3L it is (x-2L), etc.  Again, for simplicity we set L = 1.

The graph has a sawtooth shape:

Plot[Mod[x,1], {x,-2,3}]

a[0] :=  Integrate[x, {x,0,1}]

a[0]

Clear[a]

a[m_] := 2 Integrate[x Cos[2 m Pi x], {x,0,1}]

a[m] /. TrigId

% /. (-1)^(2 m) -> 1

These coefficients are all zero, again.  Why?

b[n_] := 2 Integrate[x Sin[2 n Pi x], {x,0,1}] /. TrigId /. (-1)^(2 n) -> 1

b[n]

Clear[FullSeries]

FullSeries[x_,N_] := 1/2 +       Sum[- (1/(n Pi)) Sin[2 Pi n x] ,{n,1,N}]

Let us see how the series converges:

Plot[{FullSeries[x,2], x}, {x,0,1}]

Plot[{FullSeries[x,6], x}, {x,0,1}]

Here is how it looks when we plot on a longer interval:

Plot[{FullSeries[x,4], Mod[x,1]}, {x,- 2,2}]

The Fourier series is converging nicely to the function except at the end-points of the interval, which are places where the full periodic saw-tooth function has jumps, and we saw something similar with the square pulse.  In both cases we see numerical evidence for the theorem that the Fourier series converges to f(x) where f(x) is continuous, and where it has a jump, the Fourier series converges to the average of the upper and the lower value at the jump.

A complicated function

Let's try out a different sort of example, where we need to integrate numerically: the function sin(Pi/x)

Clear[b]

b[n_] := NIntegrate[Sin[Pi/x] Sin[n Pi x], {x,-1,1}]

{b[1],b[2],b[3],b[4],b[5],b[6]}

We are clearly having convergence problems.  Let us put the series in by hand and look at 
it:

NewSeries[x_] := - 0.307133 Sin[Pi x] + 0.506251 Sin[2 Pi x] - 0.271449 Sin[3 Pi x] - 0.256129 Sin[4 Pi x] + 0.194751 Sin[5 Pi x] +  0.219291 Sin[6 Pi x]

Plot[NewSeries[x], {x, -1, 1}]
Plot[Sin[Pi/x], {x, -1, 1}]

A singular function

Here is another example, which is square-integrable but has an infinite jump.

Clear[b]
Clear[NewSeries]

b[n_] := NIntegrate[x^(-1/3) Sin[n Pi x], {x,-1,1}]

{b[1],b[2],b[3],b[4],b[5],b[6]}

NewSeries[x_] := 1.71069 Sin[Pi x] + 0.373604 Sin[2 Pi x] + 0.736895 Sin[3 Pi x] + 0.275216 Sin[4 Pi x] + 0.50107 Sin[5 Pi x] +  0.225203 Sin[6 Pi x]

Plot[NewSeries[x], {x, -1, 1}]

Automating the Fourier series

Here we set up some Mathematica  commands to compute the Fourier for a general function defined on an interval a < x < b.  Warning!  We do not recommend running these commands on a computer with a slow processor, since there are several nested definitions and integrations.  If you find yoursefl limited by the capabilities of your machine, it is recommended that you calculate the coefficients first, and then insert them into the series with a separate definition, as in the examples above.  We first give the formulae in case the integrals can be done in closed form:

Ave[f_, {a_,b_}] := (1/(b-a)) Integrate[(f /. x -> intvar1),                        {intvar1, a, b}]
A[f_,m_, {a_,b_}] := (2/(b-a)) Integrate[              (Cos[2 m Pi intvar2/(b-a)] f /. x -> intvar2),               {intvar2, a, b}]  /. TrigId
B[f_,n_, {a_,b_}] := (2/(b-a)) Integrate[              (Sin[2 n Pi intvar3/(b-a)] f /. x -> intvar3),                  {intvar3, a, b}]  /. TrigId

FullSeries[f_, infinity_, {x_,a_,b_}] := Ave[f,{a,b}] +       Sum[A[f,m,{a,b}] Cos[2 m Pi x/(b-a)], {m, 1, infinity}] +       Sum[B[f,n,{a,b}] Sin[2 n Pi x/(b-a)], {n, 1, infinity}]

An example:

FullSeries[x, 3, {x,0,1}]

Of course, it may be necessary to compute integrals numerically.  In that case, we use the commands:

NAve[f_, {a_,b_}] := (1/(b-a)) NIntegrate[(f /. x -> intvar1),                        {intvar1, a, b}]
NA[f_,m_, {a_,b_}] := (2/(b-a)) NIntegrate[              (Cos[2 m Pi intvar2/(b-a)] f /. x -> intvar2),               {intvar2, a, b}]  
NB[f_,n_, {a_,b_}] := (2/(b-a)) NIntegrate[              (Sin[2 n Pi intvar3/(b-a)] f /. x -> intvar3),                  {intvar3, a, b}]

NFullSeries[f_, infinity_, {x_,a_,b_}] := NAve[f,{a,b}] +       Sum[NA[f,m,{a,b}] Cos[2 m Pi x/(b-a)], {m, 1, infinity}] +       Sum[NB[f,n,{a,b}] Sin[2 n Pi x/(b-a)], {n, 1, infinity}]

An example:

NFullSeries[x^(2/3), 3, {x,0,1}]

Complex Fourier series

As in the previous sections, it will be convenient to program in some common trigonometric substitutions, but now we should include their exponential form:

ExpId = {E^(I Pi k) -> (-1)^k,E^(-I Pi k) -> (-1)^k}

In the last part of this notebook we will use the complex form of the Fourier series.  As you recall, the complex series is completely equivalent to the full real series, with both sines and cosines.  You may find it instructive to set the calculations up with the real functions for comparison, using the notebook for chapter 4 as a template.  The complex Fourier coefficients for a function f on an interval (a,b) are given by the formula:

FourierC[f_,k_, {a_,b_}] := (1/(b-a)) Integrate[              (Exp[- 2 k Pi I intvar2/(b-a)] f /. x -> intvar2),               {intvar2, a, b}]  /. TrigId /. ExpId

A simple example:  The constant function 1 for 0 < x < pi.

FourierC[1,k, {0,Pi}]

% /. k -> 2

%% /. k -> 3

FourierC[1,0, {0,Pi}]

Explanation of syntax

Recall that % refers to the previous calculation, %% to the one before that, etc.  The phrase /. k -> 3 tells Mathematica to evaluate the expression with the additional rule that k = 3.

So, we appear to be up and running; the only nonzero coefficient is the one for k=0, corresponding to the constant function.  This function is its own Fourier series.

A question to ponder:  What would have happened if I had written 
           %%% /. k -> 0 
in the line above?

Example:  The function x - x3 for -1 < x < 1.

FourierC[x - x^3,k, {-1,1}]

Simplify[%]

If we want a plot, it will be best to avoid complex numbers, and use the - completely equivalent - full Fourier series.

Fourier cosine series

Here we set up some Mathematica  commands to compute the Fourier cosine series.

Ave[f_, {a_,b_}] := (1/(b-a)) Integrate[(f /. x -> intvar1),                        {intvar1, a, b}]
AC[f_,m_, {a_,b_}] := (2/(b-a)) Integrate[              (Cos[m Pi intvar2/(b-a)] f /. x -> intvar2),               {intvar2, a, b}]  /. TrigId

CosSeries[f_, NN_, {x_,a_,b_}] := Ave[f,{a,b}] +       Sum[AC[f,m,{a,b}] Cos[m Pi x/(b-a)], {m, 1, NN}]

Examples:

CosSeries[Sin[Pi x], 5, {x,0,1}]

CosSeries[(x-1/2)^2, 5, {x,0,1}]

Fourier sine series

Here we set up some Mathematica  commands to compute the Fourier sine series.

BS[f_,n_, {a_,b_}] := (2/(b-a)) Integrate[              (Sin[n Pi intvar3/(b-a)] f /. x -> intvar3),                  {intvar3, a, b}]  /. TrigId

SinSeries[f_, NN_, {x_,a_,b_}] :=       Sum[BS[f,n,{a,b}] Sin[n Pi x/(b-a)], {n, 1, NN}]

Examples:

SinSeries[Sin[Pi x], 5, {x,0,1}]

SinSeries[(x-1/2)^2, 5, {x,0,1}]

Simplify[%]