Experiments with Fourier Series

Part 1. Synthesis of Periodic Functions

Enter the next block of commands to display the graph of the constant function.  Since we will be plotting graphs with similar names, we suppress the spelling warnings.

Off[General::spell1];

f[x_] :=  1;
graph1a  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
graph1b  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph1a, 
        DisplayFunction -> $DisplayFunction];

Add a first-order sine term to f and replot.  Each plot will display the previous graph and the new one.

f[x_] :=  1  +  2.7*Sin[x];
graph2a  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
graph2b  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph1b, graph2a, 
        DisplayFunction -> $DisplayFunction];

Now an order-two term:

f[x_] :=  1  +  2.7*Sin[x]  +  3.5*Cos[2*x];
graph3a  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
graph3b  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph2b, graph3a, 
        DisplayFunction -> $DisplayFunction];

Next an order-three term:

f[x_] :=  
        1  +  2.7*Sin[x]  +  3.5*Cos[2*x] - 2.4*Sin[3*x];
graph4a  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
graph4b  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph3b, graph4a, 
        DisplayFunction -> $DisplayFunction];

Now an order-four term:

f[x_] :=  
        1  +  2.7*Sin[x]  +  3.5*Cos[2*x] - 2.4*Sin[3*x] + 
          3*Cos[4*x];
graph5a  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
graph5b  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[graph4b, graph5a, 
        DisplayFunction -> $DisplayFunction];

Choose your own order-5 term, and plot the old and new graphs.  You need to replace the word "term" in the line below with your own choice.

f[x_] :=  
        1  +  2.7*Sin[x]  +  3.5*Cos[2*x] - 2.4*Sin[3*x] + 
          3*Cos[4*x]  +  term;
graph6  =  
        Plot[f[x], {x, -4*Pi, 4*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph6, graph5b, 
        DisplayFunction -> $DisplayFunction];

Answer to question in step 3:

Part 2. Fourier Approximations to a Triangular Wave

Enter the definition of f, and display the graph.

f[x_] :=  
        Which[x >= -Pi  &&  x < Pi, Pi - Abs[x], 
          x >= -3*Pi  &&  x < -Pi, Pi - Abs[x + 2*Pi], 
          x >= Pi  &&  x <= 3*Pi, Pi - Abs[x - 2*Pi]];
graph7  =  
        Plot[f[x], {x, -3*Pi, 3*Pi}, 
          PlotStyle -> RGBColor[0, 0, 1]];

Calculate the coefficients b[k].  Note that Mathematica will not symbolically integrate piecewise functions defined previously.  One approach is to use use NIntegrate.  However, with this approach integrals that are actually equal to zero may generate loss of precision warnings as below.

Clear[k];
b[k_] :=  (1/Pi)*
          NIntegrate[f[x]*Sin[k*x], {x, -Pi, Pi}];
b[1]

One way to persuade Mathematica to calculate the Fourier coefficients symbolically is to give it the explicit integrand over our range of integration.  Another approach is to express the function in terms of built-in mathematical functions that are valid over the whole domain.  For example, piecewise functions may be re-written using the Heaviside function (UnitStep in Mathematica) so that Mathematica does not have to figure out before-hand which branch to use.  Convince yourself that the function defined below is identical to f[x].

f1[
          x_] := (π - Abs[-2 π + x])*
            UnitStep[-(x - π) (x - 3 π)] + (π - 
                Abs[2 π + x])*
            UnitStep[-(x + 3 π) (x + π)] + (π - 
                Abs[x])*
            UnitStep[-(x - π) (x + π)];

We will use this version of f in our subsequent calculations.

f[x_] :=  f1[x];

We can then now symbolically determine b[k].

(1/Pi)*
        Integrate[(Pi - Abs[x])*Sin[k*x], {x, -Pi, Pi}, 
          Assumptions -> {Im[k] == 0}];
b[k_] = %

Of course, we could have anticipated this result since the product of an even and and odd function is odd, and we are integrating over a symmetric interval about the origin.

Calculate a[0], and give a geometrical explanation of this coefficient in the text cell below.

a[0]  =  (2/Pi)*Integrate[f[x], {x, 0, Pi}]

Answer:

Calculate the coefficients a[k] for k > 0.

(2/Pi)*
        Integrate[f[x]*Cos[k*x], {x, 0, Pi}, 
          Assumptions -> {k > 0  &&  Im[k] == 0}];
a[k_] = %

Compare the graph of the approximation (red) with the graph of the function (blue).

n  =  1;
s[x_] :=  a[0]/2  +  Sum[a[k]*Cos[k*x], {k, n}];
graph8  =  
        Plot[s[x], {x, -3*Pi, 3*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph7,  graph8, 
        DisplayFunction -> $DisplayFunction];

Answers to the questions in step 5:

Part 3. Fourier Approximations to a Saw-Tooth Wave

Graph the function f.  Note that Mathematica will connect the graph across vertical asymptotes.

f[x_] :=  
        Which[x ≥ -Pi  &&  x < Pi, x, 
          x ≥ -3*Pi  &&  x < -Pi, x + 2*Pi, 
          x ≥ Pi  &&  x ≤ 3*Pi, 
          x - 2*Pi];
graph9  =  
        Plot[f[x], {x, -3*Pi, 3*Pi}, 
          PlotStyle -> RGBColor[0, 0, 1]];

Here is f in terms of UnitStep.  Convince yourself this is correct.

f[x_] :=  
        x*UnitStep[-(x - Pi) (x + Pi)] + (x + 2*Pi)*
            UnitStep[-(x + 3*Pi) (x + Pi)] + (x - 2*Pi)*
            UnitStep[-(x - Pi) ( x - 3*Pi)];

Explain why a[0] must be 0.

Answer:

Calculate the coefficients a[k].

(1/Pi)*Integrate[f[x]*Cos[k*x], {x, -Pi, Pi}] // 
        Simplify;
a[k_]  =  %

Answer to the question in step 2:

Calculate the coefficients b[k] for a range of k values, and conjecture the general formula for b[k].

(1/Pi)*Integrate[x*Sin[k*x], {x, -Pi, Pi}] // 
        Simplify;
b[k_]  =  %

Compare the approximations to the function.

n  =  1;
s[x_] :=  Sum[b[k]*Sin[k*x], {k, n}]; 
graph10  =  
        Plot[s[x], {x, -3*Pi, 3*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[graph9, graph10, 
        DisplayFunction -> $DisplayFunction];

Answer to the question in step 5:

Part 4: Fourier Approximations and the Series (∑_(k = 1)^∞) (1/(k^2)).

Graph the function f.

f[x_] :=  
        Which[x ≥ -Pi  &&  x < Pi, Pi^2 - x^2, 
          x ≥ -3*Pi  &&  x < -Pi, 
          Pi^2 - (x + 2*Pi)^2, 
          x ≥ Pi  &&  x ≤ 3*Pi, 
          Pi^2 - (x - 2*Pi)^2];
graph11  =  
        Plot[f[x], {x, -3*Pi, 3*Pi}, 
          PlotStyle -> RGBColor[0, 0, 1]];

Here is f in terms of UnitStep.  Convince yourself this is correct.

f[
          x_] :=  (Pi^2 - x^2)*
            UnitStep[-(x - Pi) (x + 
                    Pi)] + (Pi^2 - (x + 2*Pi)^2)*
            UnitStep[-(x + 3*Pi) (x + 
                    Pi)] + (Pi^2 - (x - 2*Pi)^2)*
            UnitStep[-(x - Pi) ( x - 3*Pi)];

Calculate a[0].

a[0]  =  (2/Pi)*Integrate[Pi^2 - x^2, {x, 0, Pi}]

Calculate the coefficients a[k].

(2/Pi)*Integrate[(Pi^2 - x^2)*Cos[k*x], {x, 0, Pi}] // 
        Simplify;
a[k_] =  %

Calculate the coefficients b[k].

Clear[k, b];
(1/Pi)*Integrate[(Pi^2 - x^2)*Sin[k*x], {x, -Pi, Pi}]  // 
        Simplify;
b[k]  =  %

Compare the approximations to the function.

Clear[n];
n  =  1;
s[x_] :=  a[0]/2  +  Sum[a[k]*Cos[k*x], {k, n}];
graph12  =  
        Plot[s[x], {x, -3*Pi, 3*Pi}, PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
\(\(Show[graph11, graph12, 
        DisplayFunction -> $DisplayFunction];\)
\)

Part 5. Summary

Answers to the questions in Part 5: