Subscribe to DSC Newsletter

Fourier Series and Differential Equations with some applications in R (Part 1)

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.

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

f1

The above formulae can be simplified as below, for even and odd functions

 
 

f2

A few popular 2π2π periodic functions and their Fourier series are shown below:

f3

Computing the Fourier coefficients for the Square wave

 
 

f4

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))

sq

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.

Sq

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.

 

f5

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))

trh

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.

 

Tr

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))

 

st

The following animation shows how the superposition of first few terms in the Fourier series of the Sawtooth wave approximate the function.

St

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) = t2, -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))

Pr

Here is another function which is neither even nor odd, with period 2π

f21

The following animation shows how the Fourier series of the function approximates the above function more closely as more and more terms are added.

Ev

 

Solution of ODE with ERF and Resonance

f6

 

Let’s solve the following differential equation (for a system without damping) and find out when the pure resonance takes place.

f7

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

 

f8

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)))

gain

Boundary Value Problems

f9
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.
beam

 

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

f10 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.

f11

 

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())
beam 

As can be seen from the above figure, the output of the numerical solution agrees with the one obtained above theoretically.


Views: 771

Tags: Differential, Equation, Fourier, Math, Transform

Comment

You need to be a member of Data Science Central to add comments!

Join Data Science Central

Comment by Sandipan Dey on July 8, 2020 at 7:51am

Thank you very much Vincent, i shall try it next time.

Comment by Vincent Granville on July 8, 2020 at 7:38am

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!

Comment by Sandipan Dey on July 8, 2020 at 1:25am

Due to formatting issue some parts of code are replaced by junk html elements. For example, 

  • "a >0" is to be read as "a > 0"
  • "t &amp;amp;amp;amp;amp;lt;- Var('t')" is to be read as "t <- Var('t')"

Resources

Videos

  • Add Videos
  • View All

© 2020   Data Science Central ®   Powered by

Badges  |  Report an Issue  |  Privacy Policy  |  Terms of Service