Sturm--Liouville Problems

Fourier series is a particular case of more general orthogonal expansion with respect to eigenfunctions of differential operators. We start with classical Fourier expansion of a periodic function and then present some important expansions that are widely used in applications.

Fourier Series¶

Fourier series is a way to decompose a periodic function with finite period \( 2\ell \) into an infinite sum of its projections onto an orthonormal basis that consists of trigonometric polynomials. Therefore, Fourier series provides a periodic extension of a function initially defined on a finite interval of length \( 2\ell . \) They are named in honour of Jean-Baptiste Joseph Fourier (1768--1830), who used trigonometric series in representing solutions of partial differential equations, after preliminary investigations by Leonhard Euler (1707--1783), Jean le Rond d'Alembert (1717--1783), and Daniel Bernoulli (1700--1782).

There are two equivalent forms of Fourier series---one is complex form

\[ f(x) \,\sim\, \sum_{k=-\infty}^{\infty} \alpha_k e^{k{\bf j} \pi x/\ell} \]

where

\[ \alpha_k = \frac{1}{2\ell} \int_{-\ell}^{\ell} f(x)\, e^{-k{\bf j} \pi x/\ell} \,{\text d} x , \qquad k=0, \pm 1, \pm 2, \ldots ; \]

and another is trigonometric form

\[ f(x) \,\sim\, \frac{a_0}{2} + \sum_{k=1}^{\infty} \left[ a_k \cos \frac{k \pi x}{\ell} + b_k \sin \frac{k \pi x}{\ell} \right] , \]

where coefficients are defined according to Euler--Fourier formulas (they were discovered by Leonhard Euler many years before Fourier used them in his research):

\begin{align*} a_0 &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\,{\text d} x , \\ a_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \cos \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots ; \\ b_k &= \frac{1}{\ell} \int_{-\ell}^{\ell} f(x)\, \sin \frac{k \pi x}{\ell} \,{\text d} x , \qquad k= 1, 2, 3, \ldots . \end{align*}

The trigonometric form is based on Euler's formulas:

\[ \sin \theta = \frac{1}{2{\bf j}} \,e^{{\bf j}\theta} - \frac{1}{2{\bf j}} \,e^{-{\bf j}\theta} = \Im \,e^{{\bf j}\theta} = \mbox{Im} \,e^{{\bf j}\theta}, \qquad \cos \theta = \frac{1}{2} \,e^{{\bf j}\theta} - \frac{1}{2} \,e^{-{\bf j}\theta} = \Re \,e^{{\bf j}\theta} = \mbox{Re} \,e^{{\bf j}\theta}. \]

Here, j is the unit vector in positive vertical direction on the complex plane, so \( {\bf j}^2 =-1. \) You can use the following commands to calculate the nth partial sum of the Fourier series of the expression f on the interval

syms x k L n

The next command tells MATLAB that k is an integer. That will allow simple and simplify to evaluate \( \sin \frac{k\pi x}{\ell} \quad\mbox{and}\quad \cos \frac{k\pi x}{\ell} \) for a symbolic integer k.

evalin(symengine,'assume(k,Type::Integer)');

The kth Fourier cosine coefficient a_k of f is given by the command

a = @(f,x,k,L) int(f*cos(k*pi*x/L)/L,x,-L,L);

The kth Fourier sine coefficient b_k of f is given by the command

b = @(f,x,k,L) int(f*sin(k*pi*x/L)/L,x,-L,L);

The nth partial sum is given by

fs = @(f,x,n,L) a(f,x,0,L)/2 + ...
    symsum(a(f,x,k,L)*cos(k*pi*x/L) + b(f,x,k,L)*sin(k*pi*x/L),k,1,n);

Sometimes, it is convenient to use polar form:

\[ f(x) \,\sim\, \frac{a_0}{2} + \sum_{k=1}^{\infty} A_k \sin \left( \frac{k \pi x}{\ell} + \phi_k \right) = \frac{a_0}{2} + \sum_{k=1}^{\infty} A_k \cos \left( \frac{k \pi x}{\ell} - \varphi_k \right) , \]

where \( A_k = \sqrt{a_k^2 + b_k^2} \quad\mbox{and} \quad \varphi_k = \arctan \left( b_k / a_k \right) , \quad \phi_k = \arctan \left( a_k / b_k \right) . \)

Computing an integral in Matlab is fairly painless, and it's tempting to simply use a partial Fourier sum depending to the number of terms n. Let's consider the following half-wave rectifier:

syms x
f = piecewise(0< x< 1, 1-x, 1< x, 0)

The 10th partial sum is

pretty(fs(f,x,10,1))

We can also have matlab calculuate the general Fourier coefficients. To do this and get matlab to simplify the results, we can use simple. The following command gives the kth Fourier cosine coefficient, suppressing the results of all of the steps of simple except for the simplest.

[A,how]=simple(a(f,x,k,1))

If you don't want to see how simple found the answer, you can suppress the output, then just display the simplified answer. The following command does that for the kth Fourier sine coefficient.

[B,how]=simple(b(f,x,k,1)); B

Here are the plots of the partial sums for n=10,50. The plot also shows the function f.

ezplot(fs(f,x,2,1),-1,1)
hold on
ezplot(f,-1,1)
hold off
title('Partial sum with n=10')

So we get exactly the same graph. Now we are ready to expand this function into Fourier series by finding its Fourier coefficients:

\begin{align*} a_0 &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\,{\text d} x = \frac{2}{\pi} , \\ a_k &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\, \cos \left( k\, x \right) {\text d} x = -\frac{2}{\pi}\, \frac{1}{k^2 -1} \, \cos \frac{k\pi}{2} , \qquad k= 2, 3, \ldots ; \\ b_k &= \frac{1}{\pi} \int_{-\pi /2}^{\pi /2} \cos x\, \sin \left( k\, x \right) {\text d} x =0, \qquad k= 1, 2, 3, \ldots . \end{align*}

Then the corresponding Fourier series becomes

\[ f(t) \,=\,\frac{1}{\pi} + \frac{1}{2}\, \cos t - \frac{2}{\pi} \, \sum_{k\ge 1} \frac{1}{k^2 -1} \,\cos \left( k\, x \right) . \]

A partial sum with N = 20 terms gives a very good approximation:

We show how these commands work by examples.

Example: We consider a piecewise continuous function that cannot be extended into Taylor series because it is idenically zero o the interval (1,2). However, we can find its Fourier series.

f[x_] = Piecewise[{{1 - x, 0 < x < 1}, {0, 1 < x < 2}}]

Out[2]= { 1-x, if 0<x<1 and 0 for 1<x<2.

The standard Mathematica command FourierTrigSeries provides you the Fourier series of the function that is extended by zero outside the given interval (0<x<2 in our case):

FourierTrigSeries[f[x], x, 3]

Out[2]= 1/(4 \[Pi]) + (2 Cos[x] Sin[1/2]^2)/\[Pi] + (Cos[2 x] Sin[1]^2)/(2 \[Pi]) + (2 Cos[3 x] Sin[3/2]^2)/(
9 \[Pi]) + ((2 - 2 Sin[1]) Sin[x])/(2 \[Pi]) - ((-2 + Sin[2]) Sin[2 x])/(4 \[Pi]) - ((-3 + Sin[3]) Sin[3 x])/(9 \[Pi])

We can plot the partial sum with say 50 terms:

curve = FourierTrigSeries[f[x], x, 50];
Plot[curve, {x, -1.5, 3.5}, PlotRange -> 1]

However, it is not what we want: we need the Fourier series of the function f[x] extended periodically from the interval of length 2:

FourierTrigSeries[f[x], x, 3, FourierParameters -> {1, Pi}]

Out[3]= 1/4 + (2 Cos[\[Pi] x])/\[Pi]^2 + (2 Cos[3 \[Pi] x])/(9 \[Pi]^2) + Sin[\[Pi] x]/\[Pi] + Sin[2 \[Pi] x]/(2 \[Pi]) + Sin[3 \[Pi] x]/( 3 \[Pi])

To check, we find the Fourier coefficients manually:

Integrate[f[x]*Cos[n*Pi*x], {x, 0, 2}, Assumptions -> Element[n, Integers]]

Out[4]= (1 - Cos[n \[Pi]])/(n^2 \[Pi]^2)

\[ a_0 = \frac{1}{2} , \qquad a_n = \frac{1- (-1)^n}{n^2 \pi^2} , \quad n=1,2,\ldots . \]

Integrate[f[x]*Sin[n*Pi*x], {x, 0, 2}, Assumptions -> Element[n, Integers]]

Out[5]= (n \[Pi] - Sin[n \[Pi]])/(n^2 \[Pi]^2)

which is 1/(n*π). Hence, its Fourier series becomes

\[ f(x) = \frac{1}{4} + \sum_{n\ge 1} \left[ \frac{1- (-1)^n}{n^2 \pi^2} \, \cos \left( n \pi x \right) + \frac{1}{n\pi}\, \sin \left( n \pi x \right) \right] \]

We can also represent f(x) in other forms:

\[ f(x) = \frac{1}{4} + \sum_{n\ge 1} \sqrt{ \left( \frac{1- (-1)^n}{n^2 \pi^2} \right)^2 + \frac{1}{n^2 \pi^2}} \, \sin \left( n \pi x + \arctan \frac{1- (-1)^n}{n \pi} \right) = \frac{1}{4} + \sum_{n\ge 1} \sqrt{ \left( \frac{1- (-1)^n}{n^2 \pi^2} \right)^2 + \frac{1}{n^2 \pi^2}} \, \cos \left( n \pi x - \mbox{arccot} \frac{1- (-1)^n}{n \pi} \right) . \]

Now we plot some partial sums, starting with the given function itself:

Plot[f[x], {x, -1, 3}, PlotRange -> {0, 1}, AspectRatio -> Automatic, Axes -> {True, False}, PlotStyle -> Thick]

Then we plot some finite sums:

curve5Pi = FourierTrigSeries[f[x], x, 5];
Plot[%, {x, -1, 3}]

Then we repeat calculations with 50 terms:

Next we plot the correct Fourier series partial sums:

curve5=FourierTrigSeries[f[x], x, 50, FourierParameters -> {1, Pi}];
Plot[%, {x, -1, 3}]

The option FourierParameters has two parameters and when applied, it looks as FourierParameters->{a,b}
In trigonometric form, with setting \( {\bf FourierParameters} \,-> \,\{ a, b \} \) the following series is returned:

\[ \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \left[ \frac{a_0}{2} + \sum_{k=1}^n \left( a_k \,\cos (bkt) + b_k \,\sin (bkt) \right) \right] , \]

with

\[ a_k = \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \,\int_{-\pi/|b|}^{\pi/|b|} f(t)\,\cos (bkt) \,{\text d}t , \qquad k=0,1,2,\ldots , \]

\[ b_k = \left\vert \frac{b}{\pi} \right\vert^{(1+a)/2} \,\int_{-\pi/|b|}^{\pi/|b|} f(t)\,\sin (bkt) \,{\text d}t , \qquad k=1,2,\ldots . \]

Syntax for the FourierSeries command is: FourierSeries [ function , variable , number of terms ]

Changing the FourierParameters setting allows control over the limits of integration on the coefficients, and therefore, the base frequency of the series. Some work will go into calculating what value for “b” will give the proper limits of integration. There is ordinarily no reason to change “a” from its default setting, which is 1.

\[ \cos \left( a\,\sin x \right) = J_0 (a) + \sum_{k\ge 1} J_{2k} (a) \,\cos \left( 2kx \right) . \]