3.029 Spring 2022
Lecture 23 - 04/27/2022

Fun with Fourier Series

Fourier Series and Fourier Transforms are mathematical tools which allow us to express functions or data in an dual domain based on their frequencies

This is often a more compact notation, which allows us to take advantage of the dominant processes in a particular phenomenon or dataset

The importance of the Fourier transform in science and engineering cannot be overstated

Here are some examples of its use In Materials Science and Engineering:

Reciprocal space and diffraction (3.010)

Bloch waves in crystals (3.033)

Spectral methods for PDEs (3.030, 3.044)

You’ve already seen the formal definitions of Fourier Series and Fourier Transforms last semester in 3.019

We will quickly review the definitions, but will try and offer an alternative interactive exposition

Historical Note

Jean-Baptiste Joseph Fourier introduced the concept of a Fourier Series in 1807, while working on solving the heat equation

So let’s start our journey there

The 1D transient heat equation is given by

It is a linear partial differential equation

Let’s solve the equation for an “insulated” bar

i.e. the first spatial derivatives at the edges of the bar are zero

for various initial conditions

heatEquation=D[T[x,t],t]==D[T[x,t],{x,2}];
heatBCs={Derivative[1,0][T][0,t]==0,Derivative[1,0][T][1,t]==0};
heatIC[ic_]=T[x,0]==ic;

E.g. a cosine initial condition

heatSolution["cos(2πx)"][x_,t_]=DSolveValue[{heatEquation,heatBCs,heatIC[Cos[2π x]]},T[x,t],{x,t}]

Plot3D[heatSolution["cos(2πx)"][x,t],{x,0,1},{t,0,0.25},
PlotRange->All,MeshFunctions->{#2&},PlotStyle->None,Mesh->25,PlotPoints->25,
MeshStyle->Directive[Red,Thick],BoundaryStyle->Directive[Red,Thick],ImageSize->500,
AxesLabel->{"Position, x","Time, t","T(x,t)"},BaseStyle->Directive[16,Black],LabelStyle->Directive[16,Black]]

And a different cosine initial condition using a half the period

heatSolution["cos(4πx)"][x_,t_]=DSolveValue[{heatEquation,heatBCs,heatIC[Cos[4π x]]},T[x,t],{x,t}]

Plot3D[heatSolution["cos(4πx)"][x,t],{x,0,1},{t,0,0.25},
PlotRange->All,MeshFunctions->{#2&},PlotStyle->None,Mesh->25,PlotPoints->25,
MeshStyle->Directive[Blue,Thick],BoundaryStyle->Directive[Blue,Thick],ImageSize->500,
AxesLabel->{"Position, x","Time, t","T(x,t)"},BaseStyle->Directive[16,Black],LabelStyle->Directive[16,Black]]

What if we wanted to superimpose those initial conditions?

Well, we can always solve it directly

Plot3D[heatSolution["cos(2πx)+cos(4πx)"][x,t],{x,0,1},{t,0,0.25},
PlotRange->All,MeshFunctions->{#2&},PlotStyle->None,Mesh->25,PlotPoints->25,
MeshStyle->Directive[RGBColor[0,1/3,0],Thick],BoundaryStyle->Directive[RGBColor[0,1/3,0],Thick],ImageSize->500,
AxesLabel->{"Position, x","Time, t","T(x,t)"},BaseStyle->Directive[16,Black],LabelStyle->Directive[16,Black]]

Do you notice anything special about the symbolic solution?

The solution is in-fact identical to had we superimposed our two solutions independently

Plot3D[heatSolution["cos(2πx)"][x,t]+heatSolution["cos(4πx)"][x,t],{x,0,1},{t,0,0.25},
PlotRange->All,MeshFunctions->{#2&},PlotStyle->None,Mesh->25,PlotPoints->25,
MeshStyle->Directive[RGBColor[0,1/3,0],Thick],BoundaryStyle->Directive[RGBColor[0,1/3,0],Thick],ImageSize->500,
AxesLabel->{"Position, x","Time, t","T(x,t)"},BaseStyle->Directive[16,Black],LabelStyle->Directive[16,Black]]

This is a consequence of the linearity of the heat equation and it’s immensely useful

it means we can write the down the time-evolution of all* possible initial conditions as linear superpositions of our “eigen-solutions”

*there’s only one caveat - we need to be able to specify the initial condition as a superposition of simple sinusoidal functions as-well

this was Fourier’s motivation for developing Fourier Series

Consider the following initial condition

E.g. bringing together two bars with different temperatures

squareWavePlot=Plot[UnitStep[x(1/2-x)],{x,0,1},PlotStyle->Directive[Black,Thick],Frame->True,LabelStyle->Directive[Black,16],ImageSize->500,FrameLabel->{"Position, x","Initial Temperature, T(x,0)"},ExclusionsStyle->Dashed]

We can indeed model this using an appropriate Fourier Series

Manipulate[
DynamicModule[{fs,fcs,fss},
fs=FourierSeries[UnitStep[x(1/2-x)],x,n,FourierParameters->{1,2Pi}];
fcs=FourierCosSeries[UnitStep[x(1/2-x)],x,n,FourierParameters->{1,2Pi}];
fss=FourierSinSeries[UnitStep[x(1/2-x)],x,n,FourierParameters->{1,2Pi}];
Show[squareWavePlot,Plot[{fs,fcs,fss},{x,0,1},PlotStyle->{Red,Directive[RGBColor[0,1/3,0],Dashed],Directive[Blue,Dashed]}],PlotRange->5/4{-1,1}]],{{n,1,"Fourier Terms"},1,20,2},Paneled->False]

Notes:

Mathematica’s default convention for “Fourier Parameters” (which are just conventions for the formulas used) are slightly different than the ones we’re used to, so we explicitly specify them

?FourierParameters

We used three different Fourier Series expansions above

One using only Cosine terms `FourierCosSeries`, in Green

One using only Sine terms `FourierSinSeries`, in Blue

One using complex notation, `FourierSeries`, in Red

The accuracy of the approximation increases with more terms in the expansion

And indeed, if we let Mathematica turn its gears, we obtain a symbolic solution involving an infinite series of sine and cosine terms

heatSolution["unit-step"][m_][x_,t_]=DSolveValue[{heatEquation,heatBCs,heatIC[UnitStep[x(1/2-x)]]},T[x,t],{x,t}]/.∞->m

With[{f=Activate[heatSolution["unit-step"][100][x,t]]},
Plot3D[f,{x,0,1},{t,0,0.25},
PlotRange->All,MeshFunctions->{#2&},PlotStyle->None,Mesh->25,PlotPoints->25,
MeshStyle->Directive[Black,Thick],BoundaryStyle->Directive[Black,Thick],ImageSize->500,
AxesLabel->{"Position, x","Time, t","T(x,t)"},BaseStyle->Directive[16,Black],LabelStyle->Directive[16,Black]]]

‘Fourier’ Epicycles

Let’s see how far we can push this idea of “approximating periodic curves” using sine and cosine terms

First, note that it’s more convenient to work in complex notation, using Euler’s formula

Where the cosine and sine terms are given by:

In addition to the sine-cosine relation to exponentials, the above figure further demonstrates the “polar” notation of complex numbers, which can be interpreted as counter-clockwise notations in the complex plane

To investigate this “polar” mapping, let’s imagine the curve traced by a particle on a wheel with unit radius

Manipulate[
ParametricPlot[{Cos[t],Sin[t]},{t,-(10^(-6)),tmax},Frame->True,PlotStyle->Red,PlotRange->2{{-1,1},{-1,1}},LabelStyle->Directive[16,Black],ImageSize->250],{{tmax,0,"Extent of Trace"},0,2π},Paneled->False,SaveDefinitions->True]

And then add another wheel at the tip of this wheel, with half the radius, and rotating 7 times as fast

Manipulate[
ParametricPlot[{{Cos[t],Sin[t]},{Cos[tmax],Sin[tmax]}+1/2{Cos[7t],Sin[7t]},
{Cos[t],Sin[t]}+1/2{Cos[7t],Sin[7t]}},{t,-(10^(-6)),tmax},Frame->True,PlotStyle->{Directive[Red,Dashed],Directive[Blue,Dashed],Blue},PlotRange->2{{-1,1},{-1,1}},LabelStyle->Directive[16,Black],ImageSize->250],{{tmax,0,"Extent of Trace"},0,2π},Paneled->False,SaveDefinitions->True]

And then add yet another one at the tip of the second wheel with a third of the radius, rotating 17 times as fast, in the opposite direction

Manipulate[
ParametricPlot[{
{Cos[t],Sin[t]},
{Cos[tmax],Sin[tmax]}+1/2{Cos[7t],Sin[7t]},
{Cos[tmax],Sin[tmax]}+1/2{Cos[7tmax],Sin[7tmax]}+1/3{Sin[17t],Cos[17t]},
{Cos[t],Sin[t]}+1/2{Cos[7t],Sin[7t]}+1/3{Sin[17t],Cos[17t]}},{t,-(10^(-6)),tmax},Frame->True,PlotStyle->{Directive[Red,Dashed],Directive[Blue,Dashed],Directive[RGBColor[0,1/3,0],Dashed],RGBColor[0,1/3,0]},PlotRange->2{{-1,1},{-1,1}},LabelStyle->Directive[16,Black],ImageSize->250],{{tmax,0,"Extent of Trace"},0,2π},Paneled->False,SaveDefinitions->True]

And indeed, given enough little “wheels”, with varying radii, rotating at different rates, we can reproduce (almost) any closed periodic curve!

To see this, let’s express the final expression above in the complex plane

TrigToExp[#1+I#2&@@{Cos[t]+(1/2) Cos[7 t]+(1/3) Sin[17 t],(1/3) Cos[17 t]+Sin[t]+(1/2) Sin[7 t]}]

and recognize we can generalize this as

where  is the radius of each wheel

rotating with a phase

at an angular velocity

We now wish to extract these coefficients (Subscript[c, i]) given a closed periodic list of complex points

E.g. this elephant shape

ListLinePlot[,Frame->True,PlotStyle->Red,ImageSize->350,AspectRatio->1,FrameTicks->None]

The operation to extract the coefficients (Subscript[c, j]) given a list of complex points is very similar to extracting the coefficients of a Discrete Fourier Transform (DFT)

but ever so slightly different, so let’s quickly write our own

fourierCoefficients[data_,frequencyRange_]:=With[{n=Length[data]},
1/n*Table[Sum[data[[k]]*Exp[-I i k 2 Pi/n],{k,1,n}],{i,-frequencyRange,frequencyRange}]]

We first need to map our (x,y) points to the complex plane

complexElephantPoints=#1+#2I &@@@;

and then extract our coefficients for the first 10 frequencies as

fourierCoefficients[complexElephantPoints,10]

We can then implement the formula above as

epicyclesRepresentation[coefficients_,frequencyRange_][t_]:=Sum[coefficients[[j]]Exp[I(j-frequencyRange-1)t],{j,1,2 frequencyRange+1}]

epicyclesRepresentation[fourierCoefficients[complexElephantPoints,10],10][t]

and finally convert back to (x,y) and plot

parametricEquationsEpicycles[data_,frequencyRange_]:=Block[{coefficients,representation},
coefficients=fourierCoefficients[data,frequencyRange];
representation=epicyclesRepresentation[coefficients,frequencyRange][t];
Function[t,Evaluate[ReIm[ComplexExpand[representation]]]]
]

Manipulate[
With[{f=parametricEquationsEpicycles[complexElephantPoints,n]},
ParametricPlot[f[t],{t,0,2π},Axes->False,PlotStyle->Red,ImageSize->350,AspectRatio->1]],{{n,10,"Numbe of Terms"},5,100,5},Paneled->False]

We can also visualize the individual epicycles

epicyclesVisualization[data_,frequencyRange_]:=Block[{coefficients,radii,phase,ordering,edgePoints,circles},
coefficients=fourierCoefficients[data,frequencyRange];
radii=Abs/@coefficients;
phase=Arg/@coefficients;
ordering={frequencyRange+1}~Join~Riffle[Range[frequencyRange+2,2 frequencyRange+1],Reverse[Range[1,frequencyRange]]];

edgePoints[t_]=
ReIm[Accumulate[Table[coefficients[[j]]Exp[I(j-frequencyRange-1)t],{j,ordering}]]];
circles[t_]=Table[Circle[edgePoints[t][[j-1]],radii[[ordering[[j]]]]],{j,2,2frequencyRange+1}];

Function[t,Evaluate[{circles[t],Line[edgePoints[t]],Point[edgePoints[t]]}]]

]

Manipulate[
With[{
f=parametricEquationsEpicycles[complexElephantPoints,n],
gr=epicyclesVisualization[complexElephantPoints,n]},
Animate[
ParametricPlot[f[s],{s,-(10^(-6)),t},PlotStyle->Red,AspectRatio->Automatic,Axes->False,Epilog->gr[t],PlotRange->{{-400,400},{-400,200}},ImageSize->500],{{t,0,"Time"},0,2π},AnimationRunning->False,AnimationRate->.1,Paneled->False]],{{n,10,"Numbe of Terms"},5,100,5},Paneled->False]

The “person-curves” in Mathematica are built using a similar way

EntityList[(EntityClass["person curves", EntityClass["PopularCurve","Person"], "EntityClass[\"PopularCurve\", \"Person\"]", "popular curves"])]