Fourier and Laplace transforms

Lecture given for third year Bachelor students in physics

Gerald Kneller

Code

SetOptions[$FrontEnd,InitializationCellWarning->False]

Complex numbers

Roots

rootArray = Function[{n, r, ϕ}, 
   Table[N[{{0, 0}, {r^(1/n)*Cos[(ϕ + 2*k*Pi)/n], 
       r^(1/n)*Sin[(ϕ + 2*k*Pi)/n]}}], {k, 0, n - 1}]]

Matrix form

Id=IdentityMatrix[2];

Imag={{0,-1},{1,0}};

Imag.Imag+Id//MatrixForm

Fourier transform and series

Π[t_,σ_]=If[Abs[t]<σ,1,0];

FTΠ[ω_,σ_]=FourierTransform[Π[t,σ], t, ω, FourierParameters -> {1, -1},  Assumptions -> {σ > 0}]

Δ1[t_,σ_]=Π[t,σ]/(2σ);

FTΔ1[ω_,σ_]=FourierTransform[Δ1[t,σ], t, ω, FourierParameters -> {1, -1},  Assumptions -> {σ > 0}]

Δ2[t_,σ_]=PDF[NormalDistribution[0,σ],t]

FTΔ2[ω_,σ_]=FullSimplify[FourierTransform[Δ2[t,σ], t, ω, FourierParameters -> {1, -1},  Assumptions -> {σ > 0}],  Assumptions -> {σ > 0}]

Δ3[t_,σ_]=PDF[CauchyDistribution[0,σ],t]

FTΔ3[ω_,σ_]=FullSimplify[FourierTransform[Δ3[t,σ], t, ω, FourierParameters -> {1, -1},  Assumptions -> {σ > 0}],  Assumptions -> {σ > 0}]

W[t_, T_, ϵ_] = T*Δ1[t,ϵ]

G[t_, T_, r_] = (1 - r^2)/(1 + r^2 - 2*r*Cos[2*Pi*(t/T)])

Brief review of complex numbers

Historical context and definition

As many concepts in mathematics, complex numbers have been introduced as an extension of existing concepts and definitions, namely to be able to solve algebraic equations for any set of coefficients. The standard example is the quadratic equation

(x^2)+px+q=0,

where p and q are arbitrary real numbers, . The equation  may be solved by writing it in the form

(((x+(p/2)))^2) +q-((p^2)/4)=0,

that shows that the equation has two solutions,

(Subscript[x, 1,2])=-(p/2)±√(((p^2)/4)-q).

Depending on the choice p and q, the square root can have a negative argument (p^2)/4-q<0, and no solution exists in the ensemble of real numbers. It seems that (1501– 1576) and  Rafael Bombelli (1526 - 1572) were the first to conceive complex numbers to be able to solve algebraic equations as () for any set of coefficients p and q and the name “imaginary unit” for the number i fulfilling

(i^2)=-1

is attributed to René Descartes. In general a complex number is written as linear superposition of the real unit 1 and the imaginary unit i,

z = x+iy, {x,y}∈ℝ.

Here x and y are referred to as real and imaginary part of z, respectively, and one writes x=Re[z], y=Im[z]. The ensemble of all possible complex numbers of the above form is denoted ℂ. Usually, the ensemble ℂ is referred to as “complex plane”.

Elementary calculus

The rule  for adding / subtracting two complex numbers is

(Subscript[z, 1]) = (Subscript[x, 1])+(Subscript[iy, 1]), (Subscript[z, 2]) = (Subscript[x, 2])+(Subscript[iy, 2]) ⇒ (Subscript[z, 1])±(Subscript[z, 2])=((Subscript[x, 1])±(Subscript[x, 2]))+i((Subscript[y, 1])±(Subscript[y, 2])).

As for the multiplication, (Subscript[z, 1])and (Subscript[z, 2]) are treated as normal numbers, simply using the distributive law and (i^2)=-1,

(Subscript[z, 1]) = (Subscript[x, 1])+(Subscript[iy, 1]), (Subscript[z, 2]) = (Subscript[x, 2])+(Subscript[iy, 2]) ⇒ (Subscript[z, 1])(Subscript[z, 2])=((Subscript[x, 1])(Subscript[x, 2])-(Subscript[y, 1])(Subscript[y, 2]))+i((Subscript[x, 1])(Subscript[y, 2])+(Subscript[x, 2])(Subscript[y, 1])).

One verifies that the multiplication of complex numbers is commutative

(Subscript[z, 1])*(Subscript[z, 2])=(Subscript[z, 2])*(Subscript[z, 1]).

The rule for addition and subtraction is “component by component”, treating (Subscript[z, 1]) and (Subscript[z, 2])as two-dimensional vectors in a Cartesian frame. The vector interpretation suggests the following definition of a the modulus of a complex number:

z = x+iy ⇒|z|=√((x^2)+(y^2)), {x,y}∈ℝ.

We define also the complex conjugate of complex number through

z = x+iy  ⇒(SuperStar[z])=x-iy, {x,y}∈ℝ,

which shows that

((|z|)^2) =(SuperStar[z])z.

The inverse of a complex number  is defined through  (z^(-1))z=(z^(-1))z=1 and one sees that

(z^(-1))=((SuperStar[z])/((|z|)^2)).

With the above definition, the division (Subscript[z, 1])/(Subscript[z, 2]) of two complex numbers is defined through

(Subscript[z, 1])/(Subscript[z, 2]) ≡(Subscript[z, 1])((Subscript[z, 2])^(-1))=(((Subscript[z, 1])(SuperStar[(Subscript[z, 2])]))/((|(Subscript[z, 2])|)^2)) .

Polar form and complex exponential function

The polar form of a complex number is defined by writing

z = x+iy, x=r Cos[ϕ], y= r Sin[ϕ], where r=√((x^2)+(y^2) )=|z| and ϕ∈[0,2π].

The polar form of a complex number can be written in the alternative form

z = r (e^(iϕ)), where r=√((x^2)+(y^2) )=|z| and ϕ∈[0,2π].

One uses here the formula by Euler and Moivre

(e^(iϕ)) = Cos[ϕ] + i Sin[ϕ], ϕ∈[0,2π],

which can be proven by using the Taylor series of the functions which are involved:

(e^(iϕ)) = (∑_(n=0)^∞)(((iϕ))^n)/n!= (∑_(k=0)^∞)(((iϕ))^(2k))/(2k)!+(∑_(k=0)^∞)(((iϕ))^(2k))/(2k+1)!=(∑_(k=0)^∞)(((-1))^k)((ϕ^(2k))/(2k)!)+i{(∑_(k=0)^∞)(((-1))^k)((ϕ^(2k+1))/(2k+1)!)}=Cos[ϕ]+i Sin[ϕ].

For two complex numbers, (Subscript[z, 1])and (Subscript[z, 2]), we have thus

(Subscript[z, 1])(Subscript[z, 2]) = (Subscript[r, 1])(Subscript[r, 2]) (e^(i((Subscript[ϕ, 1])+(Subscript[ϕ, 2]))))

(Subscript[z, 1])/(Subscript[z, 2]) = (Subscript[r, 1])/(Subscript[r, 2]) (e^(i((Subscript[ϕ, 1])-(Subscript[ϕ, 2])))), (Subscript[r, 2])≠0.

The polar form of complex numbers can be used to prove elementary theorems for the addition and subtraction of arguments in trigonometric functions. Using that (e^(Subscript[iϕ, 1]))(e^(Subscript[iϕ, 2])) =(e^(i((Subscript[ϕ, 1])+(Subscript[ϕ, 2])))), it follows that

(Cos[(Subscript[ϕ, 1])] + i Sin[(Subscript[ϕ, 1])])(Cos[(Subscript[ϕ, 2])] + i Sin[(Subscript[ϕ, 2])])=Cos[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])] + i Sin[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])].

Performing the multiplication on the left-hand side of the above equation we obtain

Cos[(Subscript[ϕ, 1])] Cos[(Subscript[ϕ, 2])]-Sin[(Subscript[ϕ, 1])] Sin[(Subscript[ϕ, 2])]+i(Cos[(Subscript[ϕ, 2])] Sin[(Subscript[ϕ, 1])]+Cos[(Subscript[ϕ, 1])] Sin[(Subscript[ϕ, 2])])==Cos[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])] + i Sin[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])].

Therefore

Cos[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])]=Cos[(Subscript[ϕ, 1])] Cos[(Subscript[ϕ, 2])]-Sin[(Subscript[ϕ, 1])] Sin[(Subscript[ϕ, 2])],

Sin[(Subscript[ϕ, 1])+(Subscript[ϕ, 2])]=Cos[(Subscript[ϕ, 2])] Sin[(Subscript[ϕ, 1])]+Cos[(Subscript[ϕ, 1])] Sin[(Subscript[ϕ, 2])].

Roots of a complex number

We consider now an algebraic equation of the form

(z^n)=a, a∈ℂ,

and we ask for all solutions of that equation. Writing a in polar form we have

(z^n)=|a|(e^(iϕ)).

We use now that

a=|a|(e^(iϕ))=|a|(e^(i(ϕ+2πk))), k∈ℤ,

since Exp[i2πk]=1.Therefore

(Subscript[z, k])=|a|^(1/n)(e^(i((ϕ+2πk)/n))), k=0,...,n-1,

are the n different solutions of the equation. Geometrically speaking they are arranged on a circle of radius ((|a|)^(1/n)), starting with (Subscript[z, 0])=((|a|)^(1/n))(e^(i(ϕ/n)))and separated by an angle 2π /n.

GraphicsColumn[{Animate[Show[Graphics[Circle[{0,0},r^(1/n)]],Graphics[Table[{If[k==0,Red,Black],Arrow[rootArray[n,r,ϕ][[k+1]]]},{k,0,n-1}]],Axes->True,Frame->True,FrameLabel->{"Re[z]","Im[z]"},PlotRange->{{-Sqrt[10],Sqrt[10]},{-Sqrt[10],Sqrt[10]}}],{{r,0.1,"\[LeftBracketingBar]a\[RightBracketingBar]"},0,10,Appearance->"Labeled"},{{ϕ,0,"ϕ"},0,2π,0.05,Appearance->"Labeled"},{{n,2,"n"},2,12,1,Appearance->"Labeled"},AnimationRunning->False],Style["Roots of a complex number a=\[LeftBracketingBar]a\[LeftBracketingBar](e^(iϕ)). Slide the rulers for \[LeftBracketingBar]a\[RightBracketingBar] ∈ [0,10], ϕ ∈ [0,2π], n=2,3,...,12.",14]},ImageSize->900]

Matrix form of complex numbers

While complex numbers can in general not be represented by real numbers, they can be represented by 2×2 matrices of real numbers. Introducing the matrix representations

1≡(({{1, 0}, {0, 1}})), I≡(({{0, -1}, {1, 0}})), with  I·I=-1,

of the real and imaginary units, the matrix representation of a complex number  z = x + i y ({x,y} ∈ ℝ) is

z=x1+yI = (({{x, -y}, {y, x}})).

Applying the rules of matrix addition and multiplication, one finds straightforwardly that

(Subscript[z, 1])±(Subscript[z, 2])=(({{(Subscript[x, 1])±(Subscript[x, 2]), -((Subscript[y, 1])±(Subscript[y, 2]))}, {(Subscript[y, 1])±(Subscript[y, 2]), (Subscript[x, 1])±(Subscript[x, 2])}}))

represents the sum / difference of two complex numbers (Subscript[z, 1]),(Subscript[z, 2]) and ((Subscript[x, 1])(Subscript[x, 2])-(Subscript[y, 1])(Subscript[y, 2]))+i((Subscript[x, 1])(Subscript[y, 2])+(Subscript[x, 2])(Subscript[y, 1]))

(Subscript[z, 1])·(Subscript[z, 2])=(({{(Subscript[x, 1])(Subscript[x, 2])-(Subscript[y, 1])(Subscript[y, 2]), -((Subscript[x, 1])(Subscript[y, 2])+(Subscript[x, 2])(Subscript[y, 1]))}, {(Subscript[x, 1])(Subscript[y, 2])+(Subscript[x, 2])(Subscript[y, 1]), (Subscript[x, 1])(Subscript[x, 2])-(Subscript[y, 1])(Subscript[y, 2])}}))

their product. Logically, the inverse matrix (z^(-1)) should be the representation of 1/z, and this is, indeed, the case,

(z^(-1))=(1/((x^2)+(y^2)))(({{x, y}, {-y, x}})), if det[z]=(x^2)+(y^2)≠0.

Noting that

det[z]=|z|

is the squared modulus of the complex number z which is represented by the matrix z, Eq. () is thus the matrix representation of in the inverse (z^(-1))of a complex number given in Eq. ().  In this context it is important to note that the matrix multiplication (Subscript[z, 1])·(Subscript[z, 2]) defined in Eq. () is commutative,

(Subscript[z, 1])·(Subscript[z, 2])=(Subscript[z, 2])·(Subscript[z, 1]),

reflecting the fact that the multiplication of complex numbers is commutative. Therefore (Subscript[z, 1])·((Subscript[z, 2])^(-1))=((Subscript[z, 2])^(-1))·(Subscript[z, 1])represents a normal division (Subscript[z, 1])/(Subscript[z, 2]) and not a “left or right division”, which would have to be distinguished in non-commutative algebras.  We close the section by the remark that the complex number ζ=(e^(iϕ))=Cos[ϕ] + i Sin[ϕ] is represented by the orthogonal matrix

ζ=(({{Cos[ϕ], -Sin[ϕ]}, {Sin[ϕ], Cos[ϕ]}})), det[ζ]=1.

which represents rotations in (Subscript[ℝ, 2]). The norm preserving property of the multiplication of an arbitrary complex number z with ζ, which guarantees that|ζz|= |z|, is here translated by the multiplication theorem for determinants, det(ξ·z)=det(ξ)det(z)=det(z).

Fourier series

Periodic functions as Fourier series

In the following we will consider functions depending on one real parameter, t  (“time”), which have the periodicity property

f[t]=f[t+k T], T>0,k∈ℤ,

where T is the period of the function. If f[t] fulfills a certain number of conditions (the Dirichlet conditions),

f[t] is defined everywhere in (0,T) except at a finite number of points,

Between two discontinuities both f[t] and its derivative, f'[t], are continuous functions,

it can be written in form of trigonometric series of the form

f[t]=(∑_(n=-∞)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])t)), where (Subscript[ω, 0])=((2π)/T).

The series () fulfills clearly the periodicity condition (),  since (e^((Subscript[inω, 0])(t+k T)))=(e^((Subscript[inω, 0])t))(e^(ink 2π))=(e^((Subscript[inω, 0])t)). The Fourier coefficients (Subscript[f, n]) can be obtained by using the orthogonality relation

(1/T)(∫_0^T)dt (e^(-(Subscript[ikω, 0])t))(e^((Subscript[ilω, 0])t))=(Subscript[δ, kl])≡{({{1, k=l,}, {0, k≠l.}})

With the general form () of a Fourier series we can write

(1/T)(∫_0^T)dt (e^(-in (Subscript[ω, 0])t))f[t]=(1/T)(∫_0^T)dt (e^(-in (Subscript[ω, 0])t)){(∑_(k=-∞)^∞)(Subscript[f, k])(e^((Subscript[ikω, 0])t))}=(∑_(k=-∞)^∞)(Subscript[f, k]){(1/T)(∫_0^T)dt (e^(-in (Subscript[ω, 0])t))(e^((Subscript[ikω, 0])t))}=(∑_(k=-∞)^∞)(Subscript[f, k])(Subscript[δ, nk])=(Subscript[f, n]).

The Fourier coefficients are thus given by the integrals

(Subscript[f, n])=(1/T)(∫_0^T)dt (e^(-(Subscript[inω, 0])t))f[t],

and since f[t] and (e^(-(Subscript[inω, 0])t))are periodic with T, any shift of the upper and lower integration bound will not change the integral. One can thus in particular write

(Subscript[f, n])=(1/T)(∫_α^(T+α))dt (e^(-(Subscript[inω, 0])t))f[t],    α∈ℝ.

One can thus in particular chose α=-T/2 and write

(Subscript[f, n])=(1/T)(∫_(-T/2)^(T/2))dt (e^(-(Subscript[inω, 0])t))f[t].

Fourier series for real functions

The Fourier coefficients are in general complex and fulfill the symmetry condition

(SuperStar[(Subscript[f, k])])=(Subscript[f, -k])   for f[t]∈ℝ

for real functions. In this case one has

f[t]=(∑_(n=0)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])t))+(∑_(n=-1)^(-∞))(Subscript[f, n])(e^((Subscript[inω, 0])t))=(∑_(n=0)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])t))+(∑_(n=1)^∞)(Subscript[f, -n])(e^(-(Subscript[inω, 0])t))=(∑_(n=0)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])t))+(∑_(n=1)^∞)(SuperStar[(Subscript[f, n])])(e^(-(Subscript[inω, 0])t))=(∑_(n=0)^∞)(Re[(Subscript[f, n])]+iIm[(Subscript[f, n])])(Cos[(Subscript[nω, 0])t]+iSin[(Subscript[nω, 0])t])+(∑_(n=1)^∞)(Re[(Subscript[f, n])]-iIm[(Subscript[f, n])])(Cos[(Subscript[nω, 0])t]-iSin[(Subscript[nω, 0])t])=(Subscript[f, 0])+(∑_(n=1)^∞)2Re[(Subscript[f, n])]Cos[(Subscript[nω, 0])t]+(∑_(n=1)^∞)-2Im[(Subscript[f, n])]Sin[(Subscript[nω, 0])t]

Therefore f[t] can we written in the completely real form

f[t]=((Subscript[a, 0])/2)+(∑_(n=1)^∞)(Subscript[a, n])Cos[(Subscript[nω, 0])t]+(∑_(n=1)^∞)(Subscript[b, n])Sin[(Subscript[nω, 0])t]

where the real coefficients are given by

(Subscript[a, n])=2Re[(Subscript[f, n])]=(2/T)(∫_0^T)dt Cos[(Subscript[nω, 0])t]f[t].

(Subscript[b, n])=-2Im[(Subscript[f, n])]=(2/T)(∫_0^T)dt Sin[(Subscript[nω, 0])t]f[t].

The inverse relation reads

(Subscript[f, n])=((Subscript[a, n])-(Subscript[ib, n]))/2

and is interesting in so far as the coefficients (Subscript[a, n]) and  (Subscript[b, n]) reflect the symmetry of f[t],

(Subscript[a, n])=(2/T)(∫_0^T)dt Cos[(Subscript[nω, 0])t](f^((+)))[t], where (f^((+)))[t]=(f[t]+f[-t])/2,

(Subscript[b, n])=(2/T)(∫_0^T)dt Sin[(Subscript[nω, 0])t](f^((-)))[t], where (f^((-)))[t]=((f[t]-f[-t])/2).

Examples

Sawtooth function

Define the “sawtooth function” in the interval [0,2π).

f0[t_]=t/(2π);FigSawTooth=Plot[ {f0[Mod[t,2π]]},{t,-4π-1,4π+1},Frame->True,GridLines->{{-4π,-2π,0,2π,4π},{0.5}},Exclusions->None,GridLinesStyle->Directive[Gray,Dotted,Thickness[0.003]],FrameLabel->{Style["t",12],Style["f[t]",12]},PlotLabel->"Saw tooth function"]

Compute the Fourier coefficients as integral over the symmetric interval  [-π,π).

Print["(Subscript[f, n]) = ",f1[n_]=FourierCoefficient[f0[Mod[t,2π]],t,n,FourierParameters->{1, 1}]]

Show Fourier coefficients

Ref1coeff=Table[{n,Re[f1[n]]},{n,-10,10}];Imf1coeff=Table[{n,Im[f1[n]]},{n,-10,10}];ListPlot[{Ref1coeff,Imf1coeff},Frame->True,Filling->Axis,PlotRange->All,FrameLabel->{Style["n",14],Style["(Subscript[f, n])",14]},PlotLegends->{"Re{(Subscript[f, n])}","Im{(Subscript[f, n])}"}]

Fourier sum approximation

In case that only a finite number of Fourier coefficients are considered on obtains an approximation.

f[t]≈(∑_(n=-m)^m)(Subscript[f, n])(e^((Subscript[inω, 0])t)).

f0approx=Function[{t,m},ComplexExpand[Re[Sum[f1[n]Exp[I n t],{n,-m,m}]]]];Animate[Show[{Plot[ {f0[Mod[t,2π]]},{t,-4π-1,4π+1},Frame->True,GridLines->{{-4π,-2π,0,2π,4π},{0.5}},Exclusions->None,GridLinesStyle->Directive[Gray,Dotted,Thickness[0.003]],FrameLabel->{Style["t",12],Style["f[t]",12]},PlotLabel->"Saw tooth function"],Plot[f0approx[t,n],{t,-4π-1,4π+1},PlotStyle->Red,PlotLegends->{"Fourier sum"}]}],{n,0,15,1,Appearance->"Labeled"},AnimationRunning->False]

Show the Fourier sum approximation of Parseval’s theorem

(1/T)(∫_0^T)dt ((|f[t]|)^2)≈(∑_(k=-m)^m)((|(Subscript[f, k])|)^2).

This is the exact value of the integral on the left-hand side:

Print["(1/T)(∫_0^T)dt \[LeftBracketingBar]f[t](\[RightBracketingBar]^2) = ",Integrate[f0[t]^2,{t,0,2π}]/(2π)]

and this the approximation as a function of m

Parsevalf0approx=Function[{m},Sum[Abs[f1[n]]^2,{n,-m,m}]];Show[{ListPlot[Table[{m,Parsevalf0approx[m]},{m,0,15}],Frame->True,PlotRange->{{-0.1,15.1},{0,0.4}},FrameLabel->{Style["m",12],Style["Integral",12]}],Plot[1/3,{m,-0.1,15.1},PlotStyle->Dashed]}]

Trigonometric functions Cos[t] and Sin[t]

Compute the Fourier coefficients

We consider as simple cases the trigonometric functions Cos[t] and Sin[t]. Writing

Cos[t]=(((e^it)+(e^(-it)))/2) =(∑_(n=-m)^m)(Subscript[f, n])(e^int),       Sin[t]=(((e^it)-(e^(-it)))/(2i)) =(∑_(n=-m)^m)(Subscript[g, n])(e^int)

we infer

(Subscript[f, n])=({{{, ({{(1/2), n=1}, {(1/2), n=-1}, {0, otherwise}})}}) and     (Subscript[g, n])=({{{, ({{-(i/2), n=1}, {(i/2), n=-1}, {0, otherwise}})}})

Check by symbolic calculus:

f2[n_]=FourierCoefficient[Cos[t],t,n,FourierParameters->{1, 1}];Print["(Subscript[f, n]) = ",f2[n]]

g2[n_]=FourierCoefficient[Sin[t],t,n,FourierParameters->{1, 1}];Print["(Subscript[g, n]) = ",g2[n]]

These are the Fourier coefficients :

GraphicsRow[{ListPlot[{Table[{n,Re[f2[n]]},{n,-5,5}],Table[{n,Im[f2[n]]},{n,-5,5}]},Filling->Axis,Frame->True,FrameLabel->{"n","(Subscript[f, n])"},PlotLegends->{"Re{(Subscript[f, n])}","Im{(Subscript[f, n])}"},PlotRange->{-0.6,0.6},ImageSize->400],ListPlot[{Table[{n,Re[g2[n]]},{n,-5,5}],Table[{n,Im[g2[n]]},{n,-5,5}]},Filling->Axis,Frame->True,FrameLabel->{"n","(Subscript[g, n])"},PlotLegends->{"Re{(Subscript[g, n])}","Im{(Subscript[g, n])}"},PlotRange->{-0.6,0.6},ImageSize->400]},ImageSize->Full]

Compute the convolution of Cos[t] and Sin[t]

Print["(1/(2  π))(∫_0^(2  π))dτ f[t-τ]g[τ] = ",Integrate[Cos[t-τ]*Sin[τ],{τ,0,2π}]/(2π)]

Check the result by the convolution theorem for Fourier series

Print["(∑_(k = -∞)^∞) (Subscript[f, k])(Subscript[g, k]) = ",FullSimplify[Sum[f2[n]*g2[n]*Exp[I*n*t],{n,-Infinity,Infinity}],Assumptions->Element[t,Reals]]]

Compute the correlation of Cos[t] and Sin[t]

Print["(1/(2  π))(∫_0^(2  π))dτ f[t+τ](SuperStar[g])[τ] = ",Integrate[Cos[t+τ]*Sin[τ],{τ,0,2π}]/(2π)]

Check the result by the correlation theorem for Fourier series

Print["(∑_(k = -∞)^∞) (Subscript[f, k])(Subscript[g, k]^( * )) = ",FullSimplify[Sum[f2[n]*Conjugate[g2[n]]*Exp[I*n*t],{n,-Infinity,Infinity}],Assumptions->Element[t,Reals]]]

More complicated functions leading to Fourier sums

Definitions

Define f[t]=((Cos[t])^2)Sin[t] and g[t]=Cos[t]((Sin[t])^2). One sees easily that

f[t]=(((((e^it)+(e^(-it)))/2))^2)(((e^it)-(e^(-it)))/(2i))=(((e^i2t)+2+(e^(-i2t)))/4)(((e^it)-(e^(-it)))/(2i))=(1/8) i (e^(-i t))-(1/8) i (e^(i t))+(1/8) i (e^(-3 i t))-(1/8) i (e^(3 i t)),

and

g[t]=(((e^it)+(e^(-it)))/2)(((((e^it)-(e^(-it)))/(2i)))^2)=(((e^it)+(e^(-it)))/2)(((e^i2t)-2+(e^(-i2t)))/(-4))=(1/8)(e^(-i t))+(1/8)(e^(i t))-(1/8)(e^(-3 i t)) -(1/8)(e^(3 i t)) .

Fourier coefficients

Print["(Subscript[f, n]) = ",f3[n_]=FourierCoefficient[Cos[t]^2*Sin[t],t,n,FourierParameters->{1, 1}]]

Print["(Subscript[g, n]) = ",g3[n_]=FourierCoefficient[Cos[t]*Sin[t]^2,t,n,FourierParameters->{1, 1}]]

GraphicsRow[{ListPlot[{Table[{n,Re[f3[n]]},{n,-5,5}],Table[{n,Im[f3[n]]},{n,-5,5}]},Filling->Axis,Frame->True,FrameLabel->{"n","(Subscript[g, n])"},PlotLegends->{"Re{(Subscript[g, n])}","Im{(Subscript[g, n])}"},PlotRange->{-0.2,0.2},ImageSize->400],ListPlot[{Table[{n,Re[g3[n]]},{n,-5,5}],Table[{n,Im[g3[n]]},{n,-5,5}]},Filling->Axis,Frame->True,FrameLabel->{"n","(Subscript[f, n])"},PlotLegends->{"Re{(Subscript[f, n])}","Im{(Subscript[f, n])}"},PlotRange->{-0.2,0.2},ImageSize->400]},ImageSize->Full]

Basic calculation rules

The following list resumes the basic calculation rules for Fourier series, which follow directly from the definition

Linearity: Given  f[t] ⟷(Subscript[f, n])and g[t] ⟷(Subscript[g, n]) ⟹ αf[t]+βf[t] ⟷(Subscript[αf, n])+(Subscript[βg, n]).

Complex conjugate: Given  f[t] ⟷(Subscript[f, n])⟹ (SuperStar[f])[t]⟷(Subscript[f, -n]^( * )). Proof:

(SuperStar[f])[t]=(SuperStar[((∑_(n=-∞)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])t)))])=(∑_(n=-∞)^∞)(Subscript[f, n]^( * ))(e^(-(Subscript[inω, 0])t))=(∑_(n=+∞)^(-∞))(Subscript[f, -n]^( * ))(e^((Subscript[inω, 0])t))=(∑_(n=-∞)^∞)(Subscript[f, -n]^( * ))(e^((Subscript[inω, 0])t)).

Reflection in time: Given  f[t] ⟷(Subscript[f, n])⟹ f[-t]⟷(Subscript[f, -n]). Proof:

f[-t]=(∑_(n=-∞)^∞)(Subscript[f, n])(e^(-(Subscript[inω, 0])t))=(∑_(n=+∞)^(-∞))(Subscript[f, -n])(e^((Subscript[inω, 0])t))=(∑_(n=-∞)^∞)(Subscript[f, -n])(e^((Subscript[inω, 0])t)).

Shift in time: Given  f[t] ⟷(Subscript[f, n])⟹ f[t+τ]⟷(Subscript[f, n])(e^((Subscript[inω, 0])τ)). Proof:

f[t+τ]=(∑_(n=-∞)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])(t+τ)))=(∑_(n=-∞)^∞)(Subscript[f, n])(e^((Subscript[inω, 0])τ))(e^((Subscript[inω, 0])t)).

Convolution, correlation, Parseval’s theorem

Convolution theorem

The convolution of two periodic functions is defined through the integral

(Subscript[(f*g), T])[t]≡(1/T)(∫_0^T)dτ f[t-τ]g[τ]

and the convolution theorem states that

(Subscript[(f*g), T])[t]⟷ (Subscript[f, k])(Subscript[g, k]).

Proof: Defining h(t)≡(Subscript[(f*g), T])[t], we have

(Subscript[h, n])=(1/T)(∫_0^T)dt (e^(-(Subscript[inω, 0])t)){(1/T)(∫_0^T)dτ f[t-τ]g[τ]}=(1/(T^2))(∫_0^T)dτg[τ](∫_0^T)dt (e^(-(Subscript[inω, 0])t))f[t-τ](Overscript[=, u=t-τ])(1/(T^2))(∫_0^T)dτ g[τ](∫_τ^(T+τ))du (e^(-(Subscript[inω, 0])(u+τ)))f[u]=(1/(T^2))(∫_0^T)dτ g[τ](e^(-(Subscript[inω, 0])τ))(∫_(-τ)^(T-τ))du (e^(-(Subscript[inω, 0])u)))f[u](Overscript[=, use periodicity])(1/(T^2))(∫_0^T)dτ g[τ](e^(-(Subscript[inω, 0])τ))(∫_0^T)du (e^(-(Subscript[inω, 0])u))f[u]=(Underscript[{(1/T)(∫_0^T)dτ g[τ](e^(-(Subscript[inω, 0])τ))}, (Subscript[g, n])])(Underscript[{(1/T)(∫_0^T)du (e^(-(Subscript[inω, 0])u))f[u]}, (Subscript[f, n])])=(Subscript[f, n])(Subscript[g, n]).

The convolution of two periodic signals, f[t] and g[t], describes a filtering process which transforms the input signal, f[t], into a periodic output signal, h[t]=(f*g)[t], where g[t] acts as filter or transfer function,

f[t]⟶ (f*g)[t]=(∑_(n=-∞)^∞)(Subscript[f, n])(Subscript[g, n])(e^((Subscript[inω, 0])t)), where (Subscript[ω, 0])=((2π)/T).

Schematically one could write

(f*g)[t]⟵ g[t] ⟵f[t].

Since the convolution is obviously commutative, the role of input signal and transfer function can be commuted without changing the result.

Example

Consider as an example for f[t] the sawtooth function presented earlier and as filter function

g[t]=1+Cos[t].

The corresponding Fourier coefficients are

(Subscript[f, n])=({{{, ({{(1/2), n==0}, {(i/(2 n π)), otherwise}})}})    and    (Subscript[g, n])=({{{, ({{(1/2), n==-1∨n==1}, {1, n==0}, {0, otherwise}})}}),

and show that the filter g[t] picks the the frequencies n(Subscript[ω, 0]) for n=-1,0,1 and (Subscript[ω, 0])=1. Therefore

(Subscript[h, n])=(Subscript[f, n])(Subscript[g, n])=({{{, ({{(1/2), n==0}, {(i/(4 π n)), n==-1∨n==1}})}}).

The corresponding filtered sawtooth function is thus

h[t]=(1/2)-((Sin[t])/(2 π)).

Plot[ {f0[Mod[t,2π]],gPeriodic[t],hPeriodic[t]},{t,-4π-1,4π+1},Frame->True,GridLines->{{-4π,-2π,0,2π,4π},{0.5}},Exclusions->None,GridLinesStyle->Directive[Gray,Dotted,Thickness[0.003]],FrameLabel->{Style["t",12],Style[" ",12]},PlotLabel->"Filtered saw tooth function",PlotLegends ->{"f[t]","g[t]","h[t]=(f*g)[t]"}]

Correlation theorem

The correlation integral for periodic functions is similarly defined as the convolution integral,

(Subscript[(f∘g), T])[t]≡(1/T)(∫_0^T)dt f[t+τ](SuperStar[g])[t],

and the correlation theorem for Fourier series states that

(Subscript[(f∘g), T])[t]⟷ (Subscript[f, k])(SuperStar[(Subscript[g, k])]).

Proof: Defining h[t]=(1/T)(∫_0^T)dτ f[t+τ](SuperStar[g])[τ]we have

(Subscript[h, n])=(1/T)(∫_0^T)dt (e^(-(Subscript[inω, 0])t)){(1/T)(∫_0^T)dτ f[t+τ](SuperStar[g])[τ]}=(1/(T^2))(∫_0^T)dτ(SuperStar[g])[τ](∫_0^T)dt (e^(-(Subscript[inω, 0])t))f[t+τ](Overscript[=, u=t+τ])(1/(T^2))(∫_0^T)dτ(SuperStar[g])[τ](∫_(-τ)^(T-τ))du (e^(-(Subscript[inω, 0])(u-τ)))f[u]=(1/(T^2))(∫_0^T)dτ(SuperStar[g])[τ](e^((Subscript[inω, 0])τ))(∫_(-τ)^(T-τ))du (e^(-(Subscript[inω, 0])u)))f[u](Overscript[=, use perodicity])(1/(T^2))(∫_0^T)dτ(SuperStar[g])[τ](e^((Subscript[inω, 0])τ))(∫_0^T)du (e^(-(Subscript[inω, 0])u))f[u]=(Underscript[{(1/T)(∫_0^T)dτ(SuperStar[g])[τ](e^((Subscript[inω, 0])τ))}, (SuperStar[(Subscript[g, n])])])(Underscript[{(1/T)(∫_0^T)du (e^(-(Subscript[inω, 0])u))f[u]}, (Subscript[f, n])])=(Subscript[f, n])(SuperStar[(Subscript[g, n])]).

Parseval’s theorem

Parseval’s theorem follows directly from the correlation theorem, considering the autocorrelation of two periodic functions at t=0. It follows from Eq. () that

(1/T)(∫_0^T)dt f[t+τ](SuperStar[f])[t]=(∑_(k=-∞)^∞)(e^((Subscript[ikω, 0])t))((|(Subscript[f, k])|)^2).

Setting t=0 we thus obtain

(1/T)(∫_0^T)dt ((|f[t]|)^2)=(∑_(k=-∞)^∞)((|(Subscript[f, k])|)^2),

The physical interpretation of Parseval’s theorem is that the average energy of a signal equals the sum of its squared Fourier coefficients and ((|(Subscript[f, k])|)^2) can be considered as contribution of the (angular) frequency (Subscript[kω, 0]).  A concrete example is the velocity v[t] of a particle, which is supposed to be periodic with the period T. The average kinetic energy of the particle  is then

(1/T)(∫_0^T)dt (m/2)(((v[t]))^2)=(∑_(k=-∞)^∞)(e^((Subscript[ikω, 0])t)){(m/2)((|(Subscript[v, k])|)^2)}

and (m/2)((|(Subscript[v, k])|)^2)is contribution of the angular frequency (Subscript[kω, 0])(k-th harmonic) to this average energy.

Dirac comb

Definition

We ask now the question if there is a periodic function (Subscript[δ, T])[t] whose Fourier coefficients are all identically to 1, such that

(Subscript[δ, T])[t]=(∑_(n=-∞)^∞)(e^((Subscript[inω, 0])t)) ,with (Subscript[ω, 0])=2π/T.

Consequently, according to the convolution theorem for Fourier series, we have for an arbitrary periodic function f[t] with period T

f[t]=(∫_0^T)dτ (Subscript[δ, T])[t-τ]f[τ] ⟷ (OverTilde[f])[ω]= (Underscript[(Subscript[(OverTilde[δ]), T])[ω], =1]) (OverTilde[f])[ω].

Visibly such a function would be ill-defined without further specifications, since (Subscript[δ, T])[k×T]= ∞ for k ∈ ℤ. One can, however, define a function (Subscript[Δ, T])[t,η]such that its Fourier coefficients are arbitrarily close to 1 as η approaches some limit (Subscript[η, 0]),

(Subscript[δ, T])[t]=(lim_(ϵ→(Subscript[ϵ, 0])))(Subscript[Δ, T])[t,η], where  (lim_(η→(Subscript[η, 0])))(Subscript[Δ, T,n])[η]=1,

where  (Subscript[Δ, T])[t,η]is some appropriately defined “normal” function, depending parametrically on η, and (Subscript[Δ, T,n])[η] are its Fourier coefficients. Two examples will be discussed in the following.

Representation 1

We define a periodic repetition of boxcar functions,

(Subscript[Δ, T])[t,ϵ]=(∑_(n=-∞)^(+∞))(Subscript[W, T])[t-n T,ϵ]

where

(Subscript[W, T])[t,ϵ]={({{(T/(2ϵ)) for |t|< ϵ}, {0 otherwise.      }})

Here is an example for the periodic repetition for T=2π and ϵ varying between 0.05 and π-0.1:

Animate[Plot[ W[Mod[t-π,2π]-π,2*π,ϵ],{t,-4π-1,4π+1},Frame->True,GridLines->{{-4π,-2π,0,2π,4π},{0.5}},Exclusions->None,Filling->Axis,GridLinesStyle->Directive[Gray,Dotted,Thickness[0.003]],PlotRange->{0,65},PlotPoints->2000,FrameLabel->{Style["t",12],Style["(Subscript[Δ, T])[t, ϵ]",12]},PlotLabel->"Dirac comb for ϵ→0"],{ϵ,π-0.1,0.05,Appearance->"Labeled"},AnimationRunning->False]

The Fourier coefficients of  (Subscript[Δ, T])[t,ϵ] are then given by

(Subscript[Δ, T,n])[ϵ]=(1/T)(∫_(-T/2)^(T/2))dt (e^(-(Subscript[inω, 0])t))(Subscript[δ, T])[t]=(1/T)(∫_(-ϵ)^ϵ)dt ((T/(2ϵ)))(e^(-(Subscript[inω, 0])t))=(Subscript[1/(2ϵ)[(e^(-(Subscript[inω, 0])t))/(-(Subscript[inω, 0]))], 0]^T)=(1/(2ϵ))((e^(-(Subscript[inω, 0])ϵ))/(-(Subscript[inω, 0]))+(e^((Subscript[inω, 0])ϵ))/(Subscript[inω, 0])),

which can be simplified to

(Subscript[Δ, T,n])[ϵ]=((Sin[n ϵ (Subscript[ω, 0])])/(n ϵ (Subscript[ω, 0]))).

In the limit ϵ → 0  one obtains thus

(lim_(ϵ→0))(Subscript[Δ, T,n])[ϵ]=1,

and therefore

(lim_(ϵ→0))(Subscript[Δ, T])[t,ϵ]=(Subscript[δ, T])[t].

Representation 2

We consider now the analytical periodic function

(Subscript[Δ, T])[t,r]=(∑_(n=-∞)^∞)(r^(|n|)) (e^((Subscript[inω, 0])t)) ,with (Subscript[ω, 0])=2π/T, 0<r<1,

whose Fourier coefficients are by construction

(Subscript[Δ, T,n])[t]=(r^(|n|)).

The function  (Subscript[Δ, T])[t,r]has a simple analytical form,

(Subscript[Δ, T])[t,r]=(∑_(n=-∞)^∞)(r^(|n|)) (e^((Subscript[inω, 0])t))=(∑_(n=0)^∞)(r^n) (e^((Subscript[inω, 0])t))+(∑_(n=-∞)^(-1))(r^(-n)) (e^((Subscript[inω, 0])t))=(∑_(n=0)^∞)(r^n) (e^((Subscript[inω, 0])t))+(∑_(n=0)^∞)(r^n) (e^(-(Subscript[inω, 0])t))-1=(∑_(n=0)^∞)(((r (e^((Subscript[iω, 0])t))))^n)+(∑_(n=0)^∞)(((r (e^(-(Subscript[iω, 0])t))))^n)-1.

Noting that 0 < |r (e^(±(Subscript[iω, 0])t))| < 1 and the identity  (∑_(n=0)^∞)(q^n)  =1/(1-q)  (0 < |q| < 1) for a geometrical series, we can sum up the two series above to obtain

(Subscript[Δ, T])[t,r]=1/(1-r (e^((Subscript[iω, 0])t)))+1/(1-r (e^(-(Subscript[iω, 0])t)))-1,

which can be simplified to

(Subscript[Δ, T])[t,r]=((1-(r^2))/(1+(r^2)-2 r Cos[t (Subscript[ω, 0])])).

By construction

(lim_(r→1))(Subscript[Δ, T])[t,r]=(Subscript[δ, T])[t].

The figure below illustrates the form of (Subscript[Δ, T])[t,ϵ] for T= 2π, varying r between 0.01 and 0.99. Clearly (lim_(r→1))(Subscript[Δ, T])[t,ϵ]=(Subscript[δ, T])[t].

Animate[Plot[ G[Mod[t-π,2π]-π,2*π,r],{t,-4π-1,4π+1},Frame->True,GridLines->{{-4π,-2π,0,2π,4π},{0.5}},Exclusions->None,Filling->Axis,GridLinesStyle->Directive[Gray,Dotted,Thickness[0.003]],PlotRange->{0,All},PlotPoints->2000,FrameLabel->{Style["t",12],Style["(Subscript[Δ, T])[t, r]",12]},PlotLabel->"Dirac comb for r→1"],{r,0.01,0.99,Appearance->"Labeled"},AnimationRunning->False]

Noting that Cos[t (Subscript[ω, 0])]= 1 for t = n×T (n ∈ ℤ) and Cos[t (Subscript[ω, 0])]=-1 for  t = (2n+1)×T/2 we have

((1-r)/(1+r))<(Subscript[Δ, T])[t,r]< ((1+r)/(1-r)).

Solution of differential equations with an example from physics

Periodically forced damped harmonic oscillator

Equation of motion

We consider now an application of Fourier series to the solution of a damped harmonic oscillator which is driven by a periodic force. Newton’s equation of motion for such a system reads

M(Overscript[x, ¨])[t]=-K x[t]-γ (OverDot[x])[t] + f[t]

where M is the mass of the oscillator and x[t] is the time-dependent deviation from its equilibrium position. The forces on the right-hand side are the elastic force, -K x, where K is Hooke’s constant, the friction force, -γ (OverDot[x]),  where γ>0 is the friction constant, and the driving force, f(t). Supposing that the period of the latter is T and defining (Subscript[Ω, 0])=2π/T, we write f[t] as Fourier series

f[t]=(∑_(n=-∞)^∞)(Subscript[f, n])(e^(in (Subscript[Ω, 0])t )).

Particular solution

A simple way to solve the equation of motion is to use the ansatz that the solution can be written as a Fourier series as well, with the same period,

x[t]=(∑_(n=-∞)^∞)(Subscript[x, n])(e^(in (Subscript[Ω, 0])t )).

To obtain this solution, it is convenient to rewrite Eq. () in the form

(Overscript[x, ¨])[t]+ η(OverDot[x])[t]+(Subscript[ω, 0]^2) x[t]= ((f[t])/M)≡a[t],

where (Subscript[ω, 0])=(Sqrt[K/M])is the (angular) frequency of the oscillator and η=γ/M is the damping constant, which has the dimension of a frequency. Inserting now the Fourier series () into the equation of motion () we obtain

(∑_(n=-∞)^∞)(-(((n (Subscript[Ω, 0])))^2)+i η n (Subscript[Ω, 0])+(Subscript[ω, 0]^2))(Subscript[x, n])(e^(in (Subscript[Ω, 0])t ))=(∑_(n=-∞)^∞)((Subscript[f, n])/M)(e^(in (Subscript[Ω, 0])t )).

The Fourier series on the left- and righthand side are equal if all coefficients are equal, i.e.

(Subscript[x, n])=(1/(M((Subscript[ω, 0]^2)-(((n (Subscript[Ω, 0])))^2)+i η n (Subscript[Ω, 0]))))(Subscript[f, n])

Introducing the abbreviation

(Subscript[h, n])=(1/(M((Subscript[ω, 0]^2)-(((n (Subscript[Ω, 0])))^2)+i η n (Subscript[Ω, 0]))))

we may thus write

(Subscript[x, n])=(Subscript[h, n])(Subscript[f, n])

where the products (Subscript[h, n])(Subscript[f, n]) are the Fourier coefficients of a periodic convolution integral. Wee may therefore write

x[t]==(∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n])(e^(in (Subscript[Ω, 0])t ))=(Subscript[(h*f), T])[t],

where

h[t]=(∑_(n=-∞)^∞)(Subscript[h, n])(e^(in (Subscript[Ω, 0])t ))

is the transfer function describing the physical system, which is independent of the driving force. The name derives from the fact that h[t] transforms the driving force, f[t], which can be considered as “input” of a physical system (the oscillator), into a displacement, x[t], which can be considered as “output”
(see figure below).

Graphics[{Arrow[{{-4,0},{-2,0}}],{LightGray,EdgeForm[Thick],{Rectangle[{-2,-1},{2,1}]}},Arrow[{{2,0},{4,0}}],Text[Style["f[t]",14],{-3,0.3}],Text[Style["h[t]",14],{0,0}],Text[Style["x[t]",14],{3,0.3}]}]

The system represented by h[t] is linear since

x[t]=(Subscript[(h*(αf+βg)), T])[t]=α(Subscript[(h*f), T])[t]+β(Subscript[(h*g), T])[t].

Initial conditions and general solution

A point that has not yet been addressed is that the output of the system depends also on the initial conditions. The initial conditions enter through the solution of the homogenous differential equation, (Overscript[x, ¨])[t]+ η(OverDot[x])[t]+(Subscript[ω, 0]^2) x[t]=0 . To see this we define the linear differential operator

(OverHat[L])=((d^2)/(dt^2))+η(d/dt)+(Subscript[ω, 0]^2)

and write the original differential equation in the form (OverHat[L]) x[t] = f[t]/M. Defining (Subscript[x, p])[t] to fulfill (OverHat[L]) (Subscript[x, p])[t] = f[t]/M (the index “p” stands here for “particular”) and (Subscript[x, h])[t] to be the solution of the homogeneous differential equation, (OverHat[L]) (Subscript[x, h])[t] = 0, it follows immediately that (OverHat[L]) {(Subscript[x, p])[t]+(Subscript[x, h])[t]} = (OverHat[L]) (Subscript[x, p])[t]+(OverHat[L]) (Subscript[x, h])[t] = f[t]/M, which shows that (Subscript[x, p])[t]+(Subscript[x, h])[t] is a solution as well. The point is now that (Subscript[x, h])[t] has two free coefficients which can be fixed to fulfill the initial conditions for the solution of the linear differential equation (OverHat[L]) x[t] = f[t]/M. Inserting the ansatz  (Subscript[x, h])[t]=exp(λ t) in the differential equation (OverHat[L]) (Subscript[x, h])[t] = 0 leads then to

((λ^2) + η λ + (Subscript[ω, 0]^2))(e^(λ t))=0  ⇒ (Underscript[P(λ) ≡ (λ^2) + γ λ + (Subscript[ω, 0]^2) , characteristic polynomial])=0

We have thus two solutions for λ, which are the two  roots of the characteristic polynomial P(λ),

(Subscript[λ, 1,2])=-(η/2)±(Sqrt[((η^2)/4)-(Subscript[ω, 0]^2)]).

Since (OverHat[L]) is a linear differential operator the general solution of the homogeneous differential equation is

(Subscript[x, h])[t] =(Subscript[C, 1])(e^((Subscript[λ, 1])t)) + (Subscript[C, 2])(e^((Subscript[λ, 2])t)) ,

if the two roots of P(λ) are different. If they are equal, i.e. if  η/2 = (Subscript[ω, 0]), we write

(Subscript[x, h])[t] =(Subscript[C, 1])(e^((Subscript[λ, 1])t)) + (Subscript[C, 2])t (e^((Subscript[λ, 1])t)) .

This case corresponds to the critically damped limit of the motional regime of the oscillator, whereas η/2 < (Subscript[ω, 0]) defines the underdamped and η/2 > (Subscript[ω, 0]) the overdamped regime. With these preliminary considerations, the general form of the solution of the driven harmonic oscillator with friction is thus

x[t] =(Subscript[x, h])[t] + (∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n])(e^(in (Subscript[Ω, 0])t )),

and the free constants (Subscript[C, 1,2]) of the homogeneous solution (Subscript[x, h])[t],  which have the dimension of a position, can be fixed by solving

x[0] = (Subscript[C, 1])+(Subscript[C, 2])+(∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n]) = (Subscript[x, 0]),

(OverDot[x])[0] = (Subscript[λ, 1])(Subscript[C, 1])+(Subscript[λ, 2])(Subscript[C, 2])+(∑_(n=-∞)^∞)in (Subscript[Ω, 0]) (Subscript[h, n])(Subscript[f, n]) = (Subscript[v, 0]).

In the critically damped case one would have

(OverDot[x])[0] = (Subscript[λ, 1])(Subscript[C, 1])+ (Subscript[C, 2])+(∑_(n=-∞)^∞)in (Subscript[Ω, 0]) (Subscript[h, n])(Subscript[f, n]) = (Subscript[v, 0])

Stationary solution and stability

Since η > 0 it follows that

(lim_(t→∞)) x[t] =(∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n])(e^(in (Subscript[Ω, 0])t ))≡ (Subscript[x, s])[t],

which is referred to as stationary solution. The stationary solution is thus the one that is asymptotically established, independently of the initial conditions. Asymptotically means here for t ≫ τ, where τ=1/η is the characteristic time scale for the transition to the asymptotic regime. Mathematically speaking the existence of an asymptotic or stationary  regime results from the fact that the roots of the characteristic polynomial fulfill

ℜ[(Subscript[λ, k])] >0, k=1,2

Concrete example

The solution of the equation of motion discussed above with Fourier series is useful if the driving force, f[t], is described by few Fourier coefficients. In this case the Fourier series for the partial solution becomes a Fourier sum which can be explicitly calculated. We consider an example,

f[t]=M a Cos[(Subscript[Ω, 0])t].

Here we have only two non-vanishing Fourier coefficients

(Subscript[f, n])=((M a)/2)((Subscript[δ, n,1])+(Subscript[δ, n,-1])),

such that

(Subscript[x, p])[t]=(∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n])(e^(in (Subscript[Ω, 0])t ))=((a (e^(i(Subscript[Ω, 0])t )))/(2((Subscript[ω, 0]^2)-(Subscript[Ω, 0]^2)+i η  (Subscript[Ω, 0]))))+((a (e^(-i(Subscript[Ω, 0])t )))/(2((Subscript[ω, 0]^2)-(Subscript[Ω, 0]^2)-i η  (Subscript[Ω, 0]))))=Re[((a (e^(i(Subscript[Ω, 0])t )))/((Subscript[ω, 0]^2)-(Subscript[Ω, 0]^2)+i η  (Subscript[Ω, 0])))]

This reads explicitly

XpSol[t_]=ComplexExpand[Re[((a (E^(I(Subscript[Ω, 0])t )))/((Subscript[ω, 0]^2)-(Subscript[Ω, 0]^2)+I η (Subscript[Ω, 0])))]]//FullSimplify;Print["(Subscript[x, p])[t]=",XpSol[t]]

The solution of the homogeneous equation of motion reads

λ1=-η/2+Sqrt[η^2/4-(Subscript[ω, 0]^2)];λ2=-η/2-Sqrt[η^2/4-(Subscript[ω, 0]^2)];XhSol[t_]=(Subscript[C, 1])*Exp[λ1*t]+(Subscript[C, 2])*Exp[λ2*t];Print["(Subscript[x, h])[t] = ",XhSol[t]]

and therefore

XSol[t_]=XpSol[t]+XhSol[t];Print["x[t] = ",XpSol[t]+XhSol[t]]

Requiring now that x[0]=(Subscript[x, 0]) and (OverDot[x])[0]=(Subscript[v, 0]) leads to two linear equations for the coefficients (Subscript[C, 1])and (Subscript[C, 2]),

Eq1=XSol[0]==(Subscript[x, 0])

Eq2=(D[XSol[t],t]/.t->0)==(Subscript[v, 0])

the solutions of which are

Csolve=FullSimplify[Solve[{Eq1,Eq2},{(Subscript[C, 1]),(Subscript[C, 2])}]];C1sol=Csolve[[1,1,2]];C2sol=Csolve[[1,2,2]];Print["(Subscript[C, 1]) = ",C1sol,"\n(Subscript[C, 2]) = ",C2sol]

One checks that with these coefficients

XhSol[t_]=C1sol*Exp[λ1*t]+C2sol*Exp[λ2*t];XSol[t_]=XpSol[t]+XhSol[t];XOcheck=XSol[0]//Simplify;V0check=D[XSol[t],t]/.t->0//Simplify;Print["x[0] = ",XOcheck,"\n(OverDot[x])[0] = ",V0check]

The animated figure below illustrates the convergence of solution x[t] towards the the stationary solution  (Subscript[x, s])[t]=(Subscript[x, p])[t], where the oscillator frequency is set to (Subscript[ω, 0])=1 and the amplitude of the driving force to a=1. Resonance is observed as (Subscript[Ω, 0])→(Subscript[ω, 0]).

XpPlot[t_,Ω0_,η_]=XpSol[t]/.{a->1,(Subscript[Ω, 0])->Ω0,(Subscript[ω, 0])->1};XPlot[t_,Ω0_,η_]=XSol[t]/.{a->1,(Subscript[Ω, 0])->Ω0,(Subscript[ω, 0])->1,(Subscript[x, 0])->0,(Subscript[v, 0])->0};Animate[Plot[{XPlot[t,Ω0,η],XpPlot[t,Ω0,η]},{t,0,100},PlotStyle -> {ColorData[97,1],ColorData[97,2]},Frame->True,FrameLabel->{Style["t",14],Style["x[t]",14]},PlotRange->All,FrameStyle->14,ImageSize->400,PlotLegends->{"x[t]","(Subscript[x, s])[t]"}],{Ω0,0.1,1,0.1,Appearance->"Labeled"},{η,0.01,1,0.01,Appearance->"Labeled"},AnimationRunning->False]

The same example using “DSolve” :

Define the differential equation

Deq=x''[t]+η*x'[t]+((Subscript[ω, 0]))^2*x[t]==a*Cos[Ω0*t]

Solve the differential equation with x[0]=0 and x’[0]=0

dsol=DSolve[{Deq,x[0]==0,x'[0]==0},x[t],t]//Simplify

Plot the solution with fixed parameters (Subscript[ω, 0])=1 and a = 1

XPlot2[t_,Ω0_,η_]=(x[t]/.dsol[[1]])/.{(Subscript[ω, 0])->1,a->1};Animate[Plot[{XPlot2[t,Ω0,η],XpPlot[t,Ω0,η]},{t,0,100},Frame->True,FrameLabel->{Style["t",14],Style["x[t]",14]},PlotRange->All,PlotStyle -> {ColorData[97,1],ColorData[97,2]},FrameStyle->14,ImageSize->400,PlotLegends->{"x[t]","(Subscript[x, s])[t]"}],{Ω0,0.1,1,0.1,Appearance->"Labeled"},{η,0.01,1,0.01,Appearance->"Labeled"},AnimationRunning->False]

General periodically forced linear systems

Differential equation

We consider an ordinary linear differential equation of order r with constant coefficients,

(Subscript[a, r])((d^r)/(dt^r))x[t]+(Subscript[a, n-1])((d^(r-1))/(dt^(r-1)))x[t]+ … + (Subscript[a, 2])((d^2)/(dt^2))x[t]+(Subscript[a, 1])(d/dt)x[t]+(Subscript[a, 0])x[t]=f[t],

where (Subscript[a, r])≠0. Introducing the differential operator

(OverHat[L]) = (Subscript[a, r])((d^r)/(dt^r)) +(Subscript[a, n-1])((d^(r-1))/(dt^(r-1))) + … + (Subscript[a, 2])((d^2)/(dt^2)) +(Subscript[a, 1])(d/dt) +(Subscript[a, 0])

we can write

(OverHat[L])x[t]=f[t].

We will again assume that f[t] is a periodic function with period T and angular frequency (Subscript[Ω, 0]) =2π/T,

f[t]=(∑_(n=-∞)^∞)(Subscript[f, n])(e^(in (Subscript[Ω, 0])t )).

Particular solution

A particular solution of the differential equation () is obtained by inserting the ansatz

x[t]=(∑_(n=-∞)^∞)(Subscript[x, n])(e^(in (Subscript[Ω, 0])t )).

With the characteristic polynomial

P[λ]= (Subscript[a, r])(λ^r) +(Subscript[a, n-1])(λ^(r-1)) + … + (Subscript[a, 2])(λ^2) +(Subscript[a, 1])λ+(Subscript[a, 0])

we obtain

(∑_(n=-∞)^∞)(Subscript[x, n])P[in (Subscript[Ω, 0])](e^(in (Subscript[Ω, 0])t ))= (∑_(n=-∞)^∞)(Subscript[f, n])(e^(in (Subscript[Ω, 0])t ))

such that

(Subscript[x, n])= (1/(P[in (Subscript[Ω, 0])]))(Subscript[f, n]).

The transfer function has here the general form

h[t]=(∑_(n=-∞)^∞)(Subscript[h, n])(e^(in (Subscript[Ω, 0])t )) where   (Subscript[h, n])=(1/(P[in (Subscript[Ω, 0])])),

such that

x[t]=
 (∑_(n=-∞)^∞)(Subscript[h, n])(Subscript[f, n])(e^(in (Subscript[Ω, 0])t ))= (Subscript[(h * f), T])(t).

Fourier transforms

Heuristic introduction of the 1d Fourier transform in the time domain

As Fourier series represent periodic functions in time, arbitrary functions can be represented by Fourier integrals if certain conditions are met. To derive the corresponding representation, we consider a non-periodic function to be a periodic function with an arbitrarily long period, T, and write

f[t]=(∑_(n=-∞)^∞){(1/T)(∫_(-T/2)^(T/2))dt' (e^(-(Subscript[inω, 0])t'))f[t']}(e^((Subscript[inω, 0])t)), where (Subscript[ω, 0])=((2π)/T).

where the term in in curly brackets is the Fourier coefficient (Subscript[f, n])of f[t]. Using the expression for the angular frequency (Subscript[ω, 0]), we write

f[t]=(1/(2π))(∑_(n=-∞)^∞)(Subscript[ω, 0]){(∫_(-T/2)^(T/2))dt' (e^(-(Subscript[inω, 0])t'))f[t']}(e^((Subscript[inω, 0])t))

and consider T → ∞, (Subscript[ω, 0]) → 0, such that (Subscript[nω, 0]) →ω and  (Subscript[ω, 0])→ dω. The function f[t] becomes then a Fourier integral

f[t]=(1/(2π))(∫_(-∞)^∞)dω (e^(iω t))(OverTilde[f])[ω]

where the function (OverTilde[f])[ω]takes the role of the Fourier coefficients,

(OverTilde[f])[ω]=(∫_(-∞)^∞)dt' (e^(-iω t'))f[t'].

Renaming t'→t, this defines the Fourier transform pair ,

(OverTilde[f])[ω]=(∫_(-∞)^∞)dt (e^(-iω t))f[t],

f[t]=(1/(2π))(∫_(-∞)^∞)dω (e^(iω t))(OverTilde[f])[ω],

where (OverTilde[f])[ω]is the Fourier transform of f[t] and f[t] the inverse Fourier transform of (OverTilde[f])[ω].

Examples for 1d Fourier transforms

Damped cosine function

Definition

A periodic function, such as f[t] = cos[Ωt], becomes non-periodic if one considers the damped version,

f[t]=(e^(-γ|t|))cos[Ωt].

Here is a plot of the function:

Animate[Plot[{Exp[-γ*Abs[t]]Cos[Ω*t],Exp[-γ*Abs[t]],-Exp[-γ*Abs[t]]},{t,-2π,2π},PlotRange->{-1.1,1.1},Frame->True,PlotStyle-> {ColorData[97,1],{ColorData[97,1],Dashed},{ColorData[97,1],Dashed}},FrameLabel->{"t","f[t]"}],{Ω,0,3,0.01,Appearance->"Labeled"},{γ,0.01,3,0.001,Appearance->"Labeled"},AnimationRunning->False]

Fourier transform

Let us calculate the Fourier transform of the above example:

(OverTilde[f])[ω]=(∫_(-∞)^(+∞))dt (e^(-iωt))Exp[-γ|t|]Cos[Ωt]=(∫_(-∞)^(+∞))dt (Cos[ωt]-iSin[ωt])Exp[-γ|t|]Cos[Ωt]
=(∫_(-∞)^(+∞))dt Cos[ωt]Exp[-γ|t|]Cos[Ωt]=2(∫_0^∞)dt Re{(e^(-iωt))}Exp[-γt]Cos[Ωt]=2Re{(∫_0^∞)dt Exp[-γt](e^(-iωt))(((e^(iΩt))+(e^(-iΩt)))/2)}=Re{(∫_0^∞)dt (e^(-(i(ω-Ω)+γ)t))+(e^(-(i(ω+Ω)+γ)t)))}=Re{(Subscript[[(e^(-(i(ω-Ω)+γ)t))/(-(iω-Ω+γ))], 0]^∞)}+Re{(Subscript[[(e^(-(i(ω+Ω)+γ)t))/(-(iω-Ω+γ))], 0]^∞)}=Re{1/(i(ω-Ω)+γ)}+Re{1/(i(ω+Ω)+γ)}=γ/((γ^2)+(((ω-Ω))^2))+(γ/((γ^2)+(((ω+Ω))^2))).

Check by symbolic calculus:

Map[FullSimplify,FTf4[ω_,Ω_,γ_]=Apart[FourierTransform[Exp[-γ*Abs[t]]Cos[Ω t],t,ω,FourierParameters->{1, -1},Assumptions->{γ>0,Element[Ω,Reals]}],ω]]

We have thus the correspondence

Exp[-γ|t|]Cos[Ωt]⟷(γ/((γ^2)+(((ω-Ω))^2)))+(γ/((γ^2)+(((ω+Ω))^2))).

Animate[Plot[FTf4[ω,Ω,γ],{ω,-4,4},Frame->True,FrameLabel->{"ω","(OverTilde[f])[ω]"},PlotRange->All,PlotPoints->200],{Ω,0,3,0.1,Appearance->"Labeled"},{γ,0.01,3,0.001,Appearance->"Labeled"},AnimationRunning->False]

The special case Ω = 0 leads to

Exp[-γ|t|]⟷((2γ)/((γ^2)+(ω^2))).

and the limit γ → 0 will be discussed later in the context of the Dirac distribution.

Boxcar function

Definition

The boxcar function is defined through

Π[t,σ]={({{1  if |t|≤σ}, {0  if |t|≻σ}})

Animate[Plot[If[Abs[t]<σ,1,0],{t,-5,5},Frame->True,Exclusions->None,FrameLabel->{"t","Π[t]"}],{σ,0.1,3,Appearance->"Labeled"},AnimationRunning->False]

Fourier transform

The Fourier transform is easily obtained

(OverTilde[Π])[ω,σ]=(∫_(-∞)^(+∞))dt (e^(-iωt)) Π[t,σ]=(∫_(-σ)^(+σ))dt (e^(-iωt)) =(Subscript[[(( (e^(-iωt)))/(-iω))], -σ]^(+σ))=(((e^(-iωσ))-(e^(iωσ)))/(-iω))=(2/ω)(((e^(iωσ))-(e^(-iωσ)))/(2i))=((2Sin[ωσ])/ω).

Animate[Plot[FTΠ[ω,σ],{ω,-10π,10π},Frame->True,FrameLabel->{"ω","(OverTilde[Π])[ω]"},PlotRange->{-2,6}],{σ,0.1,3,0.1,Appearance->"Labeled"},AnimationRunning->False]

We have thus the correspondence

Π[t,σ] ⟷ ((2 Sin[σ ω])/ω).

Dirac distribution

Definition

Similarly as one can define a periodic function whose Fourier coefficients are “almost equal to 1”, one can define a non-periodic function whose Fourier transform is “almost equal to 1” - the Dirac distribution. Setting formally (OverTilde[δ])[ω]≡1, the function has the integral representation

δ[t]≡(1/(2π))(∫_(-∞)^(+∞))dω (e^(iωt)),

such that

f[t]=(∫_(-∞)^∞)dτ δ[t-τ]f[τ] ⟷ (OverTilde[f])[ω]= (Underscript[(OverTilde[δ])[ω], =1]) (OverTilde[f])[ω],

according to the convolution theorem of the Fourier transform. Strictly speaking δ[t] is not well defined, noting that δ[0] does not exist if one uses straightforwardly  the definition (). One may , however, consider that the Dirac distribution is the limit of some “normal” function, Δ[t,α],  where α is some parameter, whose Fourier transform fulfills (lim_(α→(Subscript[α, 0])))(OverTilde[Δ])[ω,α]=1.

Representation 1

The simplest representation of the Dirac distribution uses a normalized boxcar function

(Subscript[Δ, 1])[t,σ]=(1/(2σ))Π[t,σ]={({{(1/(2σ))  if |t|≤σ}, {0  if |t|≻σ}})

Animate[Plot[Π[t,σ]/(2σ),{t,-4,4},Frame->True,FrameLabel->{"t","Δ[t,σ]"},PlotRange->{0,20},PlotPoints->200,Filling->Axis],{σ,0.03,3,0.01,Appearance->"Labeled"},AnimationRunning->False]

The function fulfills trivially

(∫_(-∞)^(+∞))dt (Subscript[Δ, 1])[t,σ]=1

and the Fourier transform is given by

FTΔ1[ω_,σ_]=FourierTransform[Π[t,σ]/(2σ), t, ω, FourierParameters -> {1, -1},  Assumptions -> {σ > 0}]

We have thus the Fourier transform pair

(Subscript[Δ, 1])[t,σ] ⟷ ((Sin[σ ω])/(σ ω))≡(Subscript[(OverTilde[Δ]), 1])[ω,σ].

Animate[Plot[FTΠ[ω,σ]/(2σ),{ω,-10π,10π},Frame->True,FrameLabel->{"ω","(Subscript[(OverTilde[Δ]), 1])[ω,σ]"},PlotRange->{-0.3,1.1}],{σ,0.03,3,0.01,Appearance->"Labeled"},AnimationRunning->False]

Since the Fourier transform (Subscript[(OverTilde[Δ]), 1])[ω,σ] fulfills

(lim_(σ→0))(Subscript[(OverTilde[Δ]), 1])[ω,σ]=1.

we can conclude

(lim_(σ→0))(Subscript[Δ, 1])[t,σ]=δ[t].

Representation 2

Here we chose the Dirac distribution to be represented by the limit of Gaussian distribution, which is well known from mathematical statistics,

(Subscript[Δ, 2])[t,σ]=((e^(-((t^2)/(2 (σ^2)))))/((Sqrt[2 π]) σ)).

Animate[Plot[Quiet[Δ2[t,σ]],{t,-1,1},Frame->True,FrameLabel->{"t","(Subscript[Δ, 2])[t,σ]"},PlotRange->{0,45},PlotPoints->2000,Filling->Axis,WorkingPrecision->300],{σ,0.01,0.3,0.01,Appearance->"Labeled"},AnimationRunning->False]

By construction

(∫_(-∞)^(+∞))dt (Subscript[Δ, 2])[t,σ]=1

and here the Fourier transform is again a Gaussian

(Subscript[(OverTilde[Δ]), 2])[ω,σ]=(e^(-(1/2)(σ^2)(ω^2)))

Animate[Plot[Quiet[FTΔ2[ω,σ]],{ω,-10,10},Frame->True,FrameLabel->{"ω","(Subscript[(OverTilde[Δ]), 2])[ω,σ]"},PlotRange->{0,All},PlotPoints->2000,Filling->Axis,WorkingPrecision->300],{σ,0.01,0.3,0.01,Appearance->"Labeled"},AnimationRunning->False]

The Fourier transform (Subscript[(OverTilde[Δ]), 2])[ω,σ] fulfills

(lim_(σ→0))(Subscript[(OverTilde[Δ]), 1])[ω,σ]=1.

and we can again conclude

(lim_(σ→0))(Subscript[Δ, 2])[t,σ]=δ[t].

Representation 3

We will now derive a third version of the Dirac distribution, which arises naturally in the frequency domain by looking at the Fourier transform pair () for the damped cosine function. In the Limit γ → 0 we have

(lim_(γ→0))f[t]=Cos[Ωt]

and with the general definition of the Dirac distribution

(lim_(γ→0))(OverTilde[f])[ω]=(∫_(-∞)^(+∞))dt (e^(-iωt))(((e^(iΩt))+(e^(-iΩt)))/2)=π{(1/(2π))(∫_(-∞)^(+∞))dt (e^(-i(ω-Ω)t))+(1/(2π))(∫_(-∞)^(+∞))dt (e^(-i(ω+Ω)t))}=πδ(ω-Ω)+πδ(ω+Ω).

Cos[Ωt]⟷ πδ(ω-Ω)+πδ(ω+Ω)

Therefore

(lim_(γ→0))(OverTilde[f])[ω]=(lim_(γ→0)){γ/((γ^2)+(((ω-Ω))^2))+γ/((γ^2)+(((ω+Ω))^2))}=πδ[ω-Ω]+πδ[ω+Ω],

and we can conclude

(lim_(γ→0)){γ/(π ((γ^2)+(((ω∓Ω))^2)))}=πδ(ω∓Ω).

A third form for the Dirac distribution is thus

δ[ω]=(lim_(ϵ→0)){(ϵ/(π ((ϵ^2)+(ω^2))))},    ω∈ℝ,

noting that

L[ω]=(ϵ/(π ((ϵ^2)+(ω^2))))

is a normalized Lorentzian function, fulfilling (∫_(-∞)^(+∞))dω L[ω]=1. Here the Dirac distribution is naturally introduced in the frequency domain, but the definition is, of course, valid for any real variable ω → x. Use in particular ω → t

(Subscript[Δ, 3])[t,σ]=(1/(π (1+((t^2)/(σ^2))) σ)).

Animate[Plot[Quiet[Δ3[t,σ]],{t,-1,1},Frame->True,FrameLabel->{"t","(Subscript[Δ, 3])[t,σ]"},PlotRange->{0,35},PlotPoints->2000,Filling->Axis,WorkingPrecision->300],{σ,0.01,0.3,0.01,Appearance->"Labeled"},AnimationRunning->False]

Here the Fourier transform is thus

(Subscript[(OverTilde[Δ]), 3])[ω,σ]=(e^(-σ|ω|)).

Animate[Plot[Quiet[FTΔ3[ω,σ]],{ω,-10,10},Frame->True,FrameLabel->{"ω","(Subscript[(OverTilde[Δ]), 3])[ω,σ]"},PlotRange->{0,All},PlotPoints->2000,Filling->Axis,WorkingPrecision->300],{σ,0.01,0.3,0.01,Appearance->"Labeled"},AnimationRunning->False]

Basic calculation rules

The following list resumes the basic calculation rules for Fourier series, which follow directly from the definition

Linearity: Given  f[t] ⟷(OverTilde[f])[ω] and g[t] ⟷(OverTilde[g])[ω] ⟹ αf[t]+βg[t] ⟷α(OverTilde[f])[ω]+β(OverTilde[g])[ω].

Complex conjugate: Given  f[t] ⟷(OverTilde[f])[ω]⟹ (SuperStar[f])[t]⟷(SuperStar[(OverTilde[f])])[-ω] . Proof:

(SuperStar[f])[t]=(SuperStar[((1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(iωt))(OverTilde[f])[ω])])=(1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(-iωt))(SuperStar[(OverTilde[f])])[ω](Overscript[=, ω→-ω]) (1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(iωt))(SuperStar[(OverTilde[f])])[-ω].

Reflection in time/frequency: Given  f[t] ⟷(OverTilde[f])[ω] ⟹  f[-t]⟷(OverTilde[f])[-ω] . Proof:

f[-t]=(1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(-iωt))(OverTilde[f])[ω](Overscript[=, ω→-ω]) (1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(iωt))(OverTilde[f])[-ω].

Shift in time: Given   f[t] ⟷ (OverTilde[f])[ω] ⟹ f[t+τ]⟷(OverTilde[f])[ω](e^(iωτ)). Proof:

f[t+τ]=(1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(iω(t+τ)))(OverTilde[f])[ω]=(1/(2π))(Underoverscript[∫ , (-∞), (+∞)]) dω (e^(iωt))(e^(iωτ))(OverTilde[f])[ω].

Shift in frequency: Given   f[t] ⟷ (OverTilde[f])[ω] ⟹ f[t](e^(-(Subscript[iω, 0])t))⟷(OverTilde[f])[ω+(Subscript[ω, 0])]. Proof:

(OverTilde[f])[ω+(Subscript[ω, 0])]=(Underoverscript[∫ , (-∞), (+∞)]) dt (e^(-i(ω+(Subscript[ω, 0]))t))f[t] = (Underoverscript[∫ , (-∞), (+∞)]) dt (e^(-iω t))(e^(-(Subscript[iω , 0])t))f[t]

Convolution, correlation, and Parseval’s theorem

As for Fourier series, one can introduce convolution and correlation integrals for non-periodic functions, as well as corresponding Parseval theorem. These will be discussed in the following.

Convolution theorem

The convolution integral of two non-periodic functions, f[t] and g[t], is defined as

(f*g)[t]≡(∫_(-∞)^(+∞))dτ f[t-τ]g[τ]

and the convolution theorem states that

(f*g)[t]⟷ (OverTilde[f])[ω](OverTilde[g])[ω],

supposing that the Fourier transform pairs f[t]⟷ (OverTilde[f])[ω]and g[t]⟷ (OverTilde[g])[ω]exist. Relation () is thus the equivalent of the convolution theorem () for periodic functions. The proof is straightforward. Defining h[t]≡(f*g)[t] we have

(OverTilde[h])[ω]=(∫_(-∞)^(+∞))dt (e^(-iωt)){(∫_(-∞)^(+∞))dτ f[t-τ]g[τ]}=(∫_(-∞)^(+∞))dτ g[τ](∫_(-∞)^(+∞))dτ (e^(-iωt))f[t-τ](Overscript[=, u=t-τ])(∫_(-∞)^(+∞))dτ g[τ](∫_(-∞)^(+∞))dτ (e^(-iω(u+τ)))f[u]=(Underscript[(∫_(-∞)^(+∞))dτ g[τ](e^(-iωτ)), (OverTilde[g])[ω]])(Underscript[(∫_(-∞)^(+∞))dτ (e^(-iωu))f[u], f[ω]])=f[ω](OverTilde[g])[ω].

Convolution theorem in frequency space

The convolution theorem can also be formulated in (angular) frequency space. Defining

((OverTilde[f])*(OverTilde[g]))[ω]≡(∫_(-∞)^(+∞))dΩ (OverTilde[f])[ω-Ω](OverTilde[g])[Ω]

we find

((OverTilde[f])*(OverTilde[g]))[ω] ⟷(1/(2π))f[t]g[t].

The proof is analogous o the one for the convolution theorem in the time domain and is left to the reader. Apart from the factor 1/2π , the convolution the convolution theorem in time and (angular) frequency are thus symmetric: Products in time/frequency translate into convolutions in frequency/time. This is, of course, the result of the almost symmetric definitions of the Fourier transform and its inverse.

Correlation theorem

The correlation integral for non-periodic functions is defined through

(f∘g)[t]≡(∫_(-∞)^(+∞))dτ f[t+τ](SuperStar[g])[τ]

and the correlation theorem of the Fourier transform states that

(f∘g)[t]⟷(OverTilde[f])[ω](SuperStar[(OverTilde[g])])[ω].

Proof:

(OverTilde[h])[ω]=(∫_(-∞)^(+∞))dt (e^(-iωt)){(∫_(-∞)^(+∞))dτ f[t+τ](SuperStar[g])[τ]}=(∫_(-∞)^(+∞))dτ (SuperStar[g])[τ](∫_(-∞)^(+∞))dτ (e^(-iωt))f[t+τ](Overscript[=, u=t+τ])(∫_(-∞)^(+∞))dτ (SuperStar[g])[τ](∫_(-∞)^(+∞))dτ (e^(-iω(u-τ)))f[u]=(Underscript[(∫_(-∞)^(+∞))dτ (SuperStar[g])[τ](e^(iωτ)), (SuperStar[(OverTilde[g])])[ω]])(Underscript[(∫_(-∞)^(+∞))dτ (e^(-iωu))f[u], f[ω]])=f[ω](SuperStar[(OverTilde[g])])[ω].

Parseval’s theorem

It follows from the correlation theorem that the autocorrelation function of a signal can be written as

(f∘f)[t]≡(∫_(-∞)^(+∞))dτ f[t+τ](SuperStar[f])[τ] = (∫_(-∞)^(+∞))dω (e^(iωt)) ((|(OverTilde[f])[ω]|)^2)

and setting t=0 one obtains Parseval’s theorem of the Fourier transform

(∫_(-∞)^(+∞))dt ((|f[t]|)^2)=(1/(2π))(∫_(-∞)^(+∞))dω ((|(OverTilde[f])[ω]|)^2).

Three-dimensional Fourier spatial transform

Definition

So far the Fourier series and Fourier transforms concerned the time/frequency domains and have applications in signal processing, but Joseph Fourier invented Fourier series in order to solve a boundary  value problem for heat conduction. Here we briefly present the the continuous Fourier transform in n-dimensional space.  The Fourier transform of f[x] is here defined as

(OverTilde[f])[k]=(∫_(-∞)^∞)(d^n)x (e^(-ik·x))f[x]=(∫_(-∞)^∞)(Subscript[dx, 1])⋯(∫_(-∞)^∞) (Subscript[dx, n])(e^(-i((Subscript[k, 1])·(Subscript[x, 1])+⋯+(Subscript[k, n])·(Subscript[x, n]))))f[(Subscript[x, 1]),⋯,(Subscript[x, n])],

where n =1,2,... is the dimension of the problem, and the inverse Fourier transform is correspondingly

f[x]=(1/(((2π))^n))(∫_(-∞)^∞)(d^n)k  (e^(ik·x))(OverTilde[f])[k]=(1/(((2π))^n)) (∫_(-∞)^∞)(Subscript[dk, 1])⋯(∫_(-∞)^∞) (Subscript[dk, n])(e^(-i((Subscript[k, 1])·(Subscript[x, 1])+⋯+(Subscript[k, n])·(Subscript[x, n]))))(OverTilde[f])[(Subscript[k, 1]),⋯,(Subscript[k, n])].

The generalization of the Dirac distribution ()  is

δ[x]=(1/(((2π))^n))(∫_(-∞)^∞)(d^n)k  (e^(ik·x))=(1/(((2π))^n)) (∫_(-∞)^∞)(Subscript[dk, 1])⋯(∫_(-∞)^∞) (Subscript[dk, n])(e^(i((Subscript[k, 1])·(Subscript[x, 1])+⋯+(Subscript[k, n])·(Subscript[x, n])))).

Solving the diffusion equation

The time evolution of the concentration of particles which are diffusing freely in a liquid at rest is described by the diffusion equation, which has been derived by Adolf Fick in the middle of the 19th century [] and was later given a random walk interpretation by Albert Einstein [],

(∂/(∂t))c[x,t]=(Subscript[D, t]) Δc[x,t],   where  Δ = (((∂^2)/(∂(Subscript[x, 1]^2)))+ ⋯+((∂^2)/(∂(Subscript[x, n]^2))))

and (Subscript[D, t]) is the translational diffusion constant whose physical dimension is (m^2)/s.  In order to solve this partial differential equation we transform it into an ordinary differential equation by applying the Fourier transform,

(∂/(∂t))(OverTilde[c])[k,t]=-(Subscript[D, t])  ((|k|)^2) (OverTilde[c])[k,t].

This differential equation is easily solved and leads to

(OverTilde[c])[k,t]=(OverTilde[c])[k,0](e^(-(Subscript[D, t]) ((|k|)^2) t)),

where (OverTilde[c])[k,0]is the Fourier transform of the initial concentration. If we impose that all particles are initially concentrated in the origin of the coordinate system,  x = 0, i.e.

c[x,0]=δ[x] ⟹   (OverTilde[c])[k,t]=1,

it follows that

(OverTilde[c])[k,t]=(e^(-(Subscript[D, t]) ((|k|)^2) t)).

In one dimension the the inverse Fourier transform is given by

c[x_,t_]=FullSimplify[InverseFourierTransform[Exp[-(Subscript[D, t])k^2 t],k,x,FourierParameters->{1,-1}],Assumptions->{(Subscript[D, t])>0,t>0}]

such that in n dimensions

c[x,t]=((e^(-(((|x|)^2)/(4(Subscript[D, t])t))))/(((2(Sqrt[π D t])))^n)).

In both normal and Fourier space the solution of the diffusion equation has thus Gaussian form and its width is given by

σ[t]=(Sqrt[2(Subscript[D, t]) t]).

This shows that the diffusion constant, (Subscript[D, t]), describes the speed of the diffusion. The following animation illustrates the solution in two dimensions ( (Subscript[D, t])=1):

c2d[x_,y_,t_]=(c[x,t]c[y,t])/.(Subscript[D, t])->1;Animate[Plot3D[c2d[x,y,t],{x,-5,5},{y,-5,5},PlotRange->{0,5},PlotPoints->100,AxesLabel->{Style["x",14],Style["y",14],Style["c[x,y,t]",14]}],{t,0.02,0.2,Appearance->"Labeled"},AnimationRunning->False]

Elements of complex analysis

Complex functions

General definition

Complex functions map the complex plane to the complex plane

z∈ℂ ⟶ f[z]∈ℂ

where z = x + I y (x,y ∈ ℝ), and maybe decomposed into a real and imaginary part, which are each a function of the real and imaginary part of its argument, z:

f[z]=u[x,y] + I v[x,y]     (u,v ∈ ℝ).

Complex functions are typically defined by “analytical continuation” of real functions, such as the exponential function, which are defined by series expansions.  The basic example is the exponential function, the complex version of which is

(e^z)=(∑_(n=0)^∞)((z^n)/(n!)).

The result of setting here z = x + I y  is though not obvious. Here it is helpful to resort to the matrix representation of complex numbers, which leads to the evaluation of  the matrix exponential

(e^z)=(e^(x1+yI)), where  1=(({{1, 0}, {0, 1}})) and      I=(({{0, -1}, {1, 0}})).

One may then use a theorem from matrix calculus that says that  (e^(A+B))=(e^A)·(e^B)  if  [A,B]=A·B - B·A = 0. Since [1,I] = 0 it follows that

(e^z)=(e^(x1))·(e^(y1)), where (e^(x1))=(∑_(n=0)^∞)((((x1))^n)/(n!)) and  (e^(yI))=(∑_(n=0)^∞)((((yI))^n)/(n!)).

Obviously

(e^(x1))=(∑_(n=0)^∞)((((x1))^n)/(n!)) =( (∑_(n=0)^∞)((x^n)/(n!)))1 = (e^x)1 = (({{ (e^x), 0}, {0,  (e^x)}}))

and splitting the the power series for (e^(yI)) into even and odd powers we obtain

(e^(yI))=(∑_(n=0)^∞)((((yI))^n)/(n!))= (∑_(n=0)^∞)((((yI))^n)/(n!))= (∑_(k=0)^∞)((((yI))^(2k))/((2k)!))+(∑_(k=0)^∞)((((yI))^(2k+1))/((2k+1)!))= (∑_(k=0)^∞)((y^(2k))/((2k)!))(((-1))^k)+ I(  (∑_(k=0)^∞)((y^(2k+1))/((2k+1)!))(((-1))^k))=( (∑_(k=0)^∞)(((-1))^k)((y^(2k))/((2k)!)))1+(  (∑_(k=0)^∞)(((-1))^k)((y^(2k+1))/((2k+1)!)))I =Cos[y]1+Sin[y]I=(({{Cos[y], -Sin[y]}, {Sin[y], Cos[y]}})).

We obtain thus

(e^z)=(e^(x1))·(e^(y1))=(({{ (e^x), 0}, {0,  (e^x)}}))·(({{Cos[y], -Sin[y]}, {Sin[y], Cos[y]}}))=(({{(e^x)Cos[y], -(e^x)Sin[y]}, {(e^x)Sin[y], (e^x)Cos[y]}}))

and therefore

(e^z)=(e^x)Cos[y]+I (e^x)Sin[y].

This example shows what one unconsciously does when one works with complex numbers as with real ones, using simply (I^2)=-1. The fact that we can write (e^z)= (e^(x+Iy))=(e^x)(e^Iy)=(e^x)(Cos[y]+I Sin[y])is consequence of the fact that the multiplication of complex numbers is commutative, which is, in turn, the result of [1,I]=0.

Examples

Here are some examples :

fTable={Exp[x+I*y],Cos[x+I*y],Sin[x+I*y],Cosh[x+I*y],Sinh[x+I*y],Tan[x+I*y],Cot[x+I*y],Tanh[x+I*y],Coth[x+I*y],Exp[(x+I y)^2]};gTable=Table[ComplexExpand[fTable[[i]]],{i,1,Length[fTable]}];TableForm[Table[ToString[TraditionalForm[fTable[[i]]]]<>" = "<>ToString[TraditionalForm[gTable[[i]]]],{i,1,Length[fTable]}],TableAlignments->Center]

Holomorphic functions

Definition

Real functions are said to be differentiable at a given point on the real axis if the left and right derivative in that point are identical,

f'[(Subscript[x, 0])]=(lim_(x→(Subscript[x, 0]^(+))))((f[x]-f[(Subscript[x, 0])])/(x-(Subscript[x, 0])))=(lim_(x→(Subscript[x, 0]^(-))))((f[x]-f[(Subscript[x, 0])])/(x-(Subscript[x, 0]))),

where x→(Subscript[x, 0]^(+)) means “approaching from the right”, i.e. x>(Subscript[x, 0]), and x→(Subscript[x, 0]^(-)) means “approaching from the left”, i.e. x<(Subscript[x, 0]). In the complex plane a function is said to be differentiable, or “holomorphic”, if the limit

f'[(Subscript[z, 0])]=(lim_(z→(Subscript[z, 0])))((f[z]-f[(Subscript[z, 0])])/(z-(Subscript[z, 0])))

does not depend on the direction of z→(Subscript[z, 0]). The choice of a particular direction can be mathematically expressed by writing

z = (Subscript[z, 0])+ϵ (e^(iϕ))

where ϵ → 0 and ϕ is fixed. This is illustrated in the animation below, which shows the derivative of exp[z] at (Subscript[z, 0])=3/2 + I.

Exp[1.5+1.*I]

Animate[Graphics[{Arrow[{{-3,0},{3,0}}],Arrow[{{0,-3},{0,3}}],{PointSize[0.015],Red,Point[{1.5,1}]},{PointSize[0.015],Point[{1.5+ϵ*Cos[ϕ],1+ϵ*Sin[ϕ]}]},{Arrowheads[0.02],Dashed,Red,Arrow[{{1.5+ϵ*Cos[ϕ],1+ϵ*Sin[ϕ]},{1.5,1}}]},Text[Style["Re{z}",14],{3.5,0}],Text[Style["Im{z}",14],{0,3.3}],Text[Style["(Subscript[z, 0])",12,Red],{1.4,0.9}],Text[Style["z",12],{1.7+ϵ*Cos[ϕ],1+ϵ*Sin[ϕ]}],Text[ "((Δf)/(Δz)) = "<>ToString[(Exp[(1.5+ϵ*Cos[ϕ])+I*(1+ϵ*Sin[ϕ])]-Exp[1.5+I])/(ϵ*Exp[I ϕ])],{1.5,2.5}]},Background->LightGray],{ϵ,0.9,0.0000001,Appearance->"Labeled"},{ϕ,0,2π,Appearance->"Labeled"},AnimationRunning->False]

With these definitions the differentiability of a function f[z] in a point (Subscript[z, 0]) can then be expressed in the form

(d/(dϕ)){(lim_(ϵ→0))((f[z]-f[(Subscript[z, 0])])/(z-(Subscript[z, 0])))}=0.

Cauchy-Riemann equations

We will now discuss necessary and sufficient conditions for the differentiablity of a complex function. For this purpose we use again the matrix representation of complex numbers and write the differential of df = du + I dv of a function f[x,y]=u[x,y]+ I v[x,y] in matrix form:

df =du×1 + dv×I = (dx (Subscript[u, x])((Subscript[x, 0]),(Subscript[y, 0]))+dy (Subscript[u, y])((Subscript[x, 0]),(Subscript[y, 0])))(({{1, 0}, {0, 1}}))+(dx (Subscript[v, x])((Subscript[x, 0]),(Subscript[y, 0]))+dy (Subscript[v, y])((Subscript[x, 0]),(Subscript[y, 0])))(({{0, -1}, {1, 0}}))

precisely

Df=((Subscript[u, x])[(Subscript[x, 0]),(Subscript[y, 0])]dx+(Subscript[u, y])[(Subscript[x, 0]),(Subscript[y, 0])]dy)*Id+((Subscript[v, x])[(Subscript[x, 0]),(Subscript[y, 0])]dx+(Subscript[v, y])[(Subscript[x, 0]),(Subscript[y, 0])]dy)*Imag;Print["df = ",MatrixForm[Df]//TraditionalForm]

where we use the abbreviations (Subscript[u, x])[(Subscript[x, 0]),(Subscript[y, 0])]≡∂(Subscript[u, x])[x,y](Subscript[|, ((Subscript[x, 0]), (Subscript[y, 0]))]) etc. Writing now the differential dz of the argument in matrix form

Print["dz = ",MatrixForm[dz=dx*Id+dy*Imag]]

the derivative in a point ((Subscript[x, 0]),(Subscript[y, 0])) becomes

(Subscript[(df/dz)|, (Subscript[z, 0])])⟶ df[(Subscript[x, 0]),(Subscript[y, 0])]·(dz^(-1))

which is explicitly

Print["df = ",TraditionalForm[MatrixForm[fPrime=Df.Inverse[dz]]]]

The differentiability is defined by the property that the above expression does not depend on the direction of dz = dx + I dy. Writing dz in polar form,

Print["dz = ",TraditionalForm[MatrixForm[dzϕ=dz/.{dx->ϵ*Cos[ϕ],dy->ϵ*Sin[ϕ]}//FullSimplify]]]

we obtain for the matrix df[(Subscript[x, 0]),(Subscript[y, 0])]·(dz^(-1))

Print["df[(Subscript[x, 0]),(Subscript[y, 0])]·(dz^(-1)) = ",TraditionalForm[MatrixForm[fPrimeϕ=fPrime/.{dx->ϵ*Cos[ϕ],dy->ϵ*Sin[ϕ]}//FullSimplify]]]

Computing now d(df[(Subscript[x, 0]),(Subscript[y, 0])]·(dz^(-1)))/dϕ leads to

Print["(d/(dϕ))(df[(Subscript[x, 0]),(Subscript[y, 0])]·(dz^(-1))) =  ",TraditionalForm[MatrixForm[D[fPrimeϕ,ϕ]//FullSimplify]]]

and the necessary and sufficient conditions for this matrix to vanish for any ϕ are therefore

(Subscript[(Subscript[(Subscript[(Subscript[((∂u)/(∂x))|, ((Subscript[x, 0]),(Subscript[y, 0]))])= ((∂v)/(∂y))|, ((Subscript[x, 0]),(Subscript[y, 0]))])         and        ((∂u)/(∂y))|, ((Subscript[x, 0]),(Subscript[y, 0]))])= -((∂v)/(∂x))|, ((Subscript[x, 0]),(Subscript[y, 0]))]).

This set of partial differential equations are referred to as  Cauchy - Riemann equations. One verifies easily that the complex function given in the last section are differentiable.

References

A. Fick. Ueber Diffusion. Annalen der Physik, 170(1):59–86, 1855.

A. Einstein. Über die von der molekularkinetischen Theorie Der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der Physik, 322(8):549–560, 1905.