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:

1. Find the steady state solution to un-damped / damped systems with pure / near resonance.
2. 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.

 

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 > 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 <- Var("t") n <- 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 Anti-derivatives

Let’s aim at computing the Fourier coefficients for the 2π-periodic triangle wave by using the coefficients of the 2π-periodic square wave, using the anti-derivative property.

 

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 > 0, ifelse(i %% 2 == 0, t - i*pi, -t + (i+1)*pi),                 ifelse(i %% 2 == 0, -t + i*pi, t - (i-1)*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,                  TeX('$\\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"

The next animation shows how the superposition of first few terms in the Fourier series approximates the triangle wave function.

 

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 <- 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,                               TeX('$\\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², -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>0,-1,1)) / 2)     ifelse(t > 0, (t-2*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.

The next section shows how to find the largest gain corresponding to the following system with damping.

 

Let’s plot the amplitudes of the terms to see which term has the largest gain.

1 2 3 4 5 6

n = seq(1,15,2) f = (1/sqrt((49-n^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

As expected, since the right end of the beam is free, it will bend down when more loads are applied on the beam to the transverse direction of its length. The next animation shows how the beam bending varies with the value of q.

 

Numerically solving a linear system to obtain the solution of the beam-bending system represented by the 4th-order differential equation in R

 First create a near-tri-diagonal 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)*(-1e-7) # 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.