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:
In the next part more applications on differential equation / Fourier series (e.g., heat / diffusion / wave PDEs) will be discussed.
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:
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.
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)) 
1
2
3

a_0 = 1/ pi * integrate (t, 0, pi )$value print ( paste (a_0, pi /2)) #[1] "1.5707963267949 1.5707963267949" 
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.
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.
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))) 
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.
Comment
Thank you very much Vincent, i shall try it next time.
Hi Sandipan, I also experience the same formatting problem. If you put 2 spaces around the " > " symbol, it fixes the issue. By the way, great article!
Due to formatting issue some parts of code are replaced by junk html elements. For example,
"a &amp;amp;amp;amp;gt;0" is to be read as "
a > 0"
"t &amp;amp;amp;amp;lt;
Var
('
t')
" is to be read as "t < Var('t')"
© 2020 Data Science Central ® Powered by
Badges  Report an Issue  Privacy Policy  Terms of Service
Upcoming DSC Webinar
Most Popular Content on DSC
To not miss this type of content in the future, subscribe to our newsletter.
Other popular resources
Archives: 20082014  20152016  20172019  Book 1  Book 2  More
Upcoming DSC Webinar
Most popular articles
You need to be a member of Data Science Central to add comments!
Join Data Science Central