In this article, a few applications of Fourier Series in solving differential equations will be described. All the problems are taken from the edx Course: MITx – 18.03Fx: Differential Equations Fourier Series and Partial Differential Equations. The article will be posted in two parts (two separate blongs)
First a basic introduction to the Fourier series will be given and then we shall see how to solve the following ODEs / PDEs using Fourier series:
 Find the steady state solution to undamped / damped systems with pure / near resonance.
 Solve the 4th order differential equation for beam bending system with boundary values, using theoretical and numeric techniques..
In the next part more applications on differential equation / Fourier series (e.g., heat / diffusion / wave PDEs) will be discussed.
Some basics
Let f(t) be a periodic function with period L. Then the Fourier series of f(t)f(t) is given by the following superposition
The above formulae can be simplified as below, for even and odd functions
A few popular 2π2π periodic functions and their Fourier series are shown below:
Computing the Fourier coefficients for the Square wave
Computing the coefficients in R
Let’s first plot the square wave function using the following R code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

library (ggplot2) library (latex2exp) Sq = function (t) { ifelse (t &amp;amp;amp;amp;gt; 0, ifelse ( as.integer (t / pi ) %% 2 == 0, 1, 1), ifelse ( as.integer (t / pi ) %% 2 == 1, 1, 1)) } t = seq (3* pi , 3* pi , 0.01) ggplot () + geom_line ( aes (t, Sq (t)), size=1.5, col= 'red' ) + ylab ( 'f(t)' ) + scale_x_continuous (breaks= seq (3* pi , 3* pi , pi ), labels= c ( TeX ( '$3\pi$' ), TeX ( '$2\pi$' ), TeX ( '$\pi$' ),0, TeX ( '$\pi$' ), TeX ( '$2\pi$' ), TeX ( '$3\pi$' ))) + scale_y_continuous (breaks= c (1,0,1), labels= c ( '1' , '0' , '1' )) + ggtitle ( TeX ( paste ( 'Square wave of period $2\pi$' ))) + theme (plot.title = element_text (hjust = 0.5)) 
1
2
3
4
5
6
7
8
9
10

for (n in 1:6) { b_n = round (2/ pi * integrate ( function (t) sin (n*t), 0, pi )$value, 15) print ( paste (b_n, ifelse (n %% 2 == 0, 0, 4/n/ pi ))) } #[1] "1.27323954473516 1.27323954473516" #[1] "0 0" #[1] "0.424413181578388 0.424413181578388" #[1] "0 0" #[1] "0.254647908947032 0.254647908947033" #[1] "0 0" 
R symbolic computation package rSymPy can be used to compute the Fourier coefficients as shown below:
1
2
3
4
5

library (rSymPy) t &amp;amp;amp;amp;lt; Var ( "t" ) n &amp;amp;amp;amp;lt; Var ( "n" ) sympy ( "integrate(2*1*sin(n*t)/pi,(t,0,pi))" ) # sq # [1] "2*cos(pi*n)/(pi*n) + 2/(pi*n)" 
The next animation shown how the first few terms in the Fourier series approximates the periodic square wave function.
Notice from the above animation that the convergence of the Fourier series with the original periodic function is very slow near discontinuities (the square wave has discontinuities at points t=kπ,where k∈Z), which is known as Gibbs phenomenon.
Computing the Fourier Coefficients of the Triangle Wave using Antiderivatives
First let’s write a few lines of code to plot the triangle wave.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

Tr = function (t) { i = as.integer (t / pi ) ifelse (t &amp;amp;amp;amp;gt; 0, ifelse (i %% 2 == 0, t  i* pi , t + (i+1)* pi ), ifelse (i %% 2 == 0, t + i* pi , t  (i1)* pi )) } t = seq (3* pi , 3* pi , 0.01) ggplot () + geom_line ( aes (t, Tr (t)), size=1.5, col= 'red' ) + ylab ( 'f(t)' ) + scale_x_continuous (breaks= seq (3* pi , 3* pi , pi ), labels= c ( TeX ( '$3\pi$' ), TeX ( '$2\pi$' ), TeX ( '$\pi$' ),0,
( '$\pi$' ), TeX ( '$2\pi$' ), TeX ( '$3\pi$' ))) + scale_y_continuous (breaks= c ( pi ,0, pi ), labels= c ( TeX ( '$\pi$' ), '0' , TeX ( '$\pi$' ))) + ggtitle ( TeX ( paste ( 'Triangle wave of period $2\pi$' ))) + theme (plot.title = element_text (hjust = 0.5)) 
h
computing the constant term
1
2
3

a_0 = 1/ pi * integrate (t, 0, pi )$value print ( paste (a_0, pi /2)) #[1] "1.5707963267949 1.5707963267949" 
Computing the Fourier Coefficients of the Sawtooth function
Similarly, let’s visualize the sawtooth wave function using the following R code snippet:
1
2
3
4
5
6
7
8
9
10
11
12
13
14

St = function (t) { i &amp;amp;amp;amp;lt; floor ((t + pi ) / (2* pi )) t  i*2* pi } t = seq (3* pi , 3* pi , 0.01) ggplot () + geom_line ( aes (t, St (t)), size=1.5, col= 'red' ) + scale_x_continuous (breaks= seq (3* pi , 3* pi , pi ), labels= c ( TeX ( '$3\pi$' ), TeX ( '$2\pi$' ), TeX ( '$\pi$' ),0,
( '$\pi$' ), TeX ( '$2\pi$' ), TeX ( '$3\pi$' ))) + scale_y_continuous (breaks= c ( pi ,0, pi ), labels= c ( TeX ( '$\pi$' ), '0' , TeX ( '$\pi$' ))) + ylab ( 'f(t)' ) + ggtitle ( TeX ( paste ( 'Sawtooth wave of period $2\pi$' ))) + theme (plot.title = element_text (hjust = 0.5)) 
The following animation shows how the superposition of first few terms in the Fourier series of the Sawtooth wave approximate the function.
Computing the Fourier series for the even periodic extension of a function
Let’s first plot the even periodic extension of the function f(t) = t^{2}, 1 ≤ t ≤ 1, with the following R code.
1
2
3
4
5
6
7
8
9
10
11
12
13

Pr = function (t) { i = as.integer ( abs (t ifelse (t&amp;amp;amp;amp;gt;0,1,1)) / 2) ifelse (t &amp;amp;amp;amp;gt; 0, (t2*i)^2, (t+2*i)^2) } t = seq (5, 5, 0.01) ggplot () + geom_line ( aes (t, Pr (t)), size=1.5, col= 'red' ) + scale_x_continuous (breaks=5:5, labels=5:5) + ylab ( 'f(t)' ) + ggtitle ( TeX ( paste ('Even periodic extension of the function $ f (t) = t^2$, $1\leq t \leq 1$'))) + theme (plot.title = element_text (hjust = 0.5)) 
Here is another function which is neither even nor odd, with period 2π
The following animation shows how the Fourier series of the function approximates the above function more closely as more and more terms are added.
Solution of ODE with ERF and Resonance
Let’s solve the following differential equation (for a system without damping) and find out when the pure resonance takes place.
1
2
3
4
5
6

n = seq (1,15,2) f = (1/ sqrt ((49n^2)^2 + (0.1*n)^2)/n) ggplot () + geom_point ( aes (n, f)) + geom_line ( aes (n, f)) + ylab ( 'c_n' ) + scale_x_continuous (breaks= seq (1,15,2), labels= as.character ( seq (1,15,2))) 
Boundary Value Problems
Numerically solving a linear system to obtain the solution of the beambending system represented by the 4thorder differential equation in R
First create a neartridiagonal matrix A that looks like the following one, it takes care of the differential coefficients of the beam equation along with all the boundary value conditions.
Now, let’s solve the linear system using the following R code.
1
2
3
4
5
6
7
8
9

#Create a vector b that is zero for boundary conditions and 0.0000001 in every other entry. b = rep (1,10)*(1e7) # Create a vector v that solves Av = b. v = solve (A, b) #Create a column vector x of 10 evenly spaced points between 0 and 1 (for plotting) x = seq (0,1,1/9) #Plot v on the vertical axis, and x on the horizontal axis. print ( ggplot () + geom_line ( aes (x,v), col= 'blue' , size=2) + geom_point ( aes (x, v), size=2) + ggtitle ( paste ( 'Solution of the linear (beam bending) system' )) + theme_bw ()) 
As can be seen from the above figure, the output of the numerical solution agrees with the one obtained above theoretically.
Credit: Data Science Central By: Sandipan Dey