Phys 238
Fourier Series Representation of Various Functions

Clear["Global`*"]

Fourier Expansion

Given f[x] periodic in [0, λ], the Fourier series representation (up to nmax terms) is given by

k[n_] := n 2 π / λ
fs[nmax_, x_] := a0/2 + Sum[a[[n]]Cos[k[n] x] + b[[n]] Sin[k[n] x], {n, 1, nmax}]

Define functions to compute the a and b coefficients.  We’ll calculate them all at once, and then build the series later.

af[n_] := (2/λ) NIntegrate[Cos[k[n] x] f[x], {x, 0, λ}]
bf[n_] := (2/λ) NIntegrate[Sin[k[n] x] f[x], {x, 0, λ}]

Different Functions to Consider

λ = 10; (* For plotting *)

Gaussian:

μ = λ/2;
σ = λ/10;
gauss[x_] := Exp[-(((Mod[x,λ]-μ))^2)/(2 (σ^2))]
Plot[gauss[x], {x, 0, 2 λ}, LabelStyle->Larger]

Triangle Wave:

triangle[x_] := Piecewise[{
{((2 Mod[x, λ])/λ),0 ≤ Mod[x, λ] ≤ λ/2},
{2 (1 - ((Mod[x, λ])/λ)), λ/2 < Mod[x, λ] ≤ λ}
}]

Plot[triangle[x], {x, 0, 2 λ}, LabelStyle->Larger]

Square Wave:

square[x_] := Piecewise[{
{1, 0 ≤ Mod[x, λ] ≤ λ/2},
{-1, λ/2 < Mod[x, λ] ≤ λ}
}]

Plot[square[x], {x, 0, 2 λ}, LabelStyle->Larger]

Sawtooth Wave:

saw[x_] := Mod[x, λ]
Plot[saw[x], {x, 0, 2 λ}, LabelStyle->Larger]

Pick a particular function:

Clear[f]
f[x_] :=square[x];
Plot[f[x], {x, 0, 2λ}]

Compute the Fourier Coefficients

For these integrals, Mathematica’s numerical integration will often raise various convergence errors.  You should read and understand the warnings before deciding whether to ignore them, but if you do wish to ignore them, you can do so with the Off[] function.

Off[NIntegrate::ncvb];Off[NIntegrate::slwcon]

nmax = 30;
Print["a0 = ", a0 = af[0]];
a = Table[af[n], {n, 1, nmax}];
b = Table[bf[n], {n, 1, nmax}];
Print[TableForm[Table[{n, a[[n]], b[[n]]}, {n, 1, 9}], TableHeadings->{None, {"n", "Cos Coeff.", "Sin Coeff."}}]];
ListPlot[{a, b}, Filling -> Axis, FillingStyle -> Black, PlotRange->All, PlotStyle->{Blue, Red}, PlotLegends -> {"Cosine", "Sine"}, LabelStyle->Larger]

c = Table[Sqrt[a[[n]]^2 + b[[n]]^2], {n, 1, nmax}];
ListPlot[c, Filling -> Axis, FillingStyle -> Black, PlotRange->All, AxesLabel->{"Mode n", "Amplitude"},  LabelStyle->Larger]

Reconstructing the original waveform

ymax = FindMaximum[f[x], {x, λ/2, 0, λ}][[1]];
ymax = Max[ymax, 1.0, Max[Abs[a]], Max[Abs[b]]];
ymin = FindMinimum[f[x], {x, λ/2, 0, λ}][[1]];
ymin = Min[ymin, 0.0, -Max[Abs[a]], -Max[Abs[b]]];Manipulate[Plot[Evaluate[{f[x], a[[n]] Cos[k[n]x],b[[n]] Sin[k[n] x], fs[n, x]}], {x, 0, 2λ},PlotLegends -> {"Original", "Cos term", "Sin term", "Series"},  PlotRange->{1.2ymin, 1.2ymax}, ImageSize -> Scaled[0.7]], {n, 0, nmax, 1, Appearance->"Open"}]

Pick a function built from random numbers

Clear[f]
randvalues = Table[{x, RandomReal[]}, {x, 0,λ, λ/20}];
f = Interpolation[randvalues, InterpolationOrder -> 2];
Plot[f[x], {x, 0, λ}, PlotRange->All]

Compute the Fourier Coefficients

nmax = 30;

Dynamic[n] (* Just to indicate progress *)
a0 = af[0];
a = Table[af[n], {n, 1, nmax}];
b = Table[bf[n], {n, 1, nmax}];
c = Table[Sqrt[a[[n]]^2 + b[[n]]^2], {n, 1, nmax}];
ListPlot[{a, b}, Filling -> Axis, FillingStyle -> Black, PlotRange->All, PlotStyle->{Blue, Red}, PlotLegends -> {"Cosine", "Sine"}, LabelStyle->Larger]
ListPlot[c, Filling -> Axis, FillingStyle -> Black, PlotRange->All, AxesLabel->{"Mode n", "Amplitude"},  LabelStyle->Larger]

Reconstructing the original waveform

ymax = FindMaximum[f[x], {x, λ/2, 0, λ}][[1]];
ymax = Max[ymax, 1.0, Max[Abs[a]], Max[Abs[b]]];
ymin = FindMinimum[f[x], {x, λ/2, 0, λ}][[1]];
ymin = Min[ymin, 0.0, -Max[Abs[a]], -Max[Abs[b]]];
Manipulate[Plot[Evaluate[{f[x], a[[n]] Cos[k[n]x],b[[n]] Sin[k[n] x], fs[n, x]}], {x, 0, λ}, PlotRange->{ymin, ymax}, ImageSize -> Scaled[0.8], AspectRatio->1/2], {n, 0, nmax, 1, Appearance->"Open"}]

Periodic Continuation

p1 = Plot[f[x], {x, 0, λ}];
p2 = Plot[f[x - λ], {x, λ, 3 λ/2}];
p3 = Plot[Evaluate[fs[nmax, x]], {x, λ/2, 3λ/2}, PlotStyle->Red];
Show[p3, p1, p2, PlotRange->All, ImageSize -> 600, AxesOrigin -> {0, 0}]