EP4210 - Computational Physics

Kirit Makwana
Department of Physics, IIT Hyderabad

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 1 / 15
 Boundary value problems

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 2 / 15
 Motivation
Consider the Poisson equation
ρ
∇2 ϕ = − (1)
ϵ0
The uniqueness theorem for the Poisson equations states that there is
a unique solution which satisfies the boundary condition if ϕ is
specified on the boundary
Consider the 1D version of this equation

d 2ϕ ρ(x)
2
=− (2)
dx ϵ0
The uniqueness theorem implies that if ϕ is specified at two points,
x = a and x = b, then ϕ is unique in the interval [a, b]
This defines a boundary value problem - specifically a “Dirichlet”
boundary value problem

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 3 / 15
 Shooting method
Consider a more general ODE
f ′′ (x) + p(x)f ′ (x) + q(x)f (x) = 0 (3)
The prime indicates derivative, f ′ = df /dx, f ′′ = d 2 f /dx 2
This is a second order ODE. The functions p(x) and q(x) are given
We need to solve for f (x) in the domain a < x < b given the
boundary conditions f (x = a) = fa and f (x = b) = fb
This second order ODE can be converted into two coupled first order
ODEs by defining
f′ =y (4)
′
y = −py − qf (5)
This can be written in matrix form
 ′   
f 0 1 f
= (6)
y −q −p y

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 4 / 15
 Shooting method

This equation acn be solved by the RK method by identifying the

function G as
  
0 1 f
G= (7)
−q −p y

We require the initial conditions for applying RK methods

 
fa
(8)
ya

We know f (x = a) and f (x = b), but we do not know

y (x = a) = f ′ (x = a)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 5 / 15
 Shooting method
We can assume some value for ya , denote it as ỹa
Now by knowing fa and ỹa , the RK method can be applied to solve
the ODE
Integrate up to x = b and obtain f˜(x = b)
Now this f˜b need not match with the boundary condition specified fb ,
because we have probably guessed a wrong starting value for ya
But now, we can try different starting values ỹa and obtain the values
f˜b at the other boundary, thereby obtaining a function

f˜b (ỹa ) (9)

Thus, given an initial condition ỹa , we can find where the solution of
the ODe lies at the other boundary f˜b
Now we need to find the root of the equation

f˜b (ỹa ) = fb (10)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 6 / 15
 Shooting method

This can be done by bisection method

Guessing two different initial values ỹa which give f˜b on two sides of
the the desired solution fb , we can then use the bisection method to
find the root of the equation f˜b (ỹa ) − fb = 0
Once the correct parameter ỹa is found, the solution is also found
This method is called ”shooting” method because we try different
initial derivative values to “shoot” the solution like a projectile and we
keep on adjusting the starting value to hit the “target” which is the
boundary condition at the other end

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 7 / 15
 Shooting method

It can be visualized as below

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 8 / 15
 Finite Difference methods
Consider a general second order ODE
f ′′ (x) + p(x)f ′ (x) + q(x)f (x) = s(x); (11)
The boundary values f (x = 0) = f0 and f (x = a) = fa are given and
so are the functions p, q and s. We have to solve for f (x) in the
domain 0 ≤ x ≤ a
Divide the space into N + 1 intervals of width ∆x

The grid points are labelled n = 0, 1, 2, ..., N + 1

n = 0 and n = N + 1 are the boundary points where the value of the
solution is known
There are N internal points where we need to find the solution fn
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210 - Computational Physics 9 / 15
 Finite difference technique
We can develop finite difference approximations for the derivatives in
the equation f ′′ (x) + p(x)f ′ (x) + q(x)f (x) = s(x)
More specifically, by Taylor series expansion we can show that
fn+1 − fn−1
= f ′ (xn ) + O(∆x 2 ) (12)
2∆x
fn+1 − 2fn + fn−1
= f ′′ (xn ) + O(∆x 2 ) (13)
∆x 2
H.W. - check this
Substituting these expressions for f ′ and f ′′ in the equation, we get
fn+1 + fn−1 − 2fn pn (fn+1 − fn−1 )
+ + qn fn = sn (14)
∆x 2 2∆x
Reorganizing this we get
     
1 pn 1 pn 2
fn+1 + + fn−1 − + f n qn − = sn
∆x 2 2∆x ∆x 2 2∆x ∆x 2
(15)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210 - Computational Physics 10 / 15
 Finite difference equations
This equation can be written in matrix form
      "fn−1 #
1 pn 2 1 pn
∆x 2 − 2∆x qn − ∆x 2 ∆x 2 + 2∆x
fn = sn (16)
fn+1
Lets write down these equations explicitly for n = 1, 2, ..., N
      "f0 #
1 p1 p1
∆x 2 − 2∆x q1 − ∆x2 2 1
∆x 2 + 2∆x
f1 = s1 (17)
f2
      "f1 #
1 p2 2 1 p2
∆x 2 − 2∆x q2 − ∆x 2 ∆x 2 + 2∆x f2 = s2 (18)
f3
. (19)
. (20)
      "fN−1 #
1 pN 2 1 pN
∆x 2 − 2∆x qN − ∆x 2 ∆x 2 + 2∆x
fN = sN (21)
fN+1
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210 - Computational Physics 11 / 15
 Matrix formulation
These N coupled linear-algebra equations can be compactly written as
 2
 
1 p
 
q1 − 1
∆x 2 2 + 2∆x 0 . . . 0
  ∆x     f1 
1 p2 2 1 p2
 − 2∆x q2 − + 0 . . 0
 ∆x 2  ∆x 2
 ∆x
2 2∆x
    f2

 0 1 − p3
2
∆x 2∆x
q3 − 2
∆x 2
1
∆x 2
3
+ 2∆x
p
0 .

0  ff34 
 .  .
 .  .
 .  .
     
 1 pN−1 1 + pN−1
 .
0... − 2∆x qN−1 − 2 2
∆x 2 ∆x 2 2∆x fN
 ∆x 

  
0... 0 1 − pN qN − 2 2
2
∆x 2∆x ∆x
N×N
(22)
 s1   
 1 p1 
− 2∆x f0
s2 ∆x 2
 ss34   0
0

=.−
   .


 ..   .
0 


. 1 Np
+ 2∆x fN+1
sN ∆x 2

(23)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 12 / 15
 Matrix formulation
This equation can be written as a matrix equation
DF = S − B (24)
 2
 
1 p1
 
q1 − + 2∆x 0 . . . 0
 ∆x 2  ∆x
2
  
 1 − p2 q2 − 2 2 1 +
p2
0 . . 0
 ∆x 2 2∆x
 ∆x  ∆x
2 2∆x
   
 0 1 − p3 q3 − 2 1 3p
+ 2∆x 0 .

0
 2
∆x 2∆x ∆x 2 ∆x 2
D= . 
 . 
 .      

 1 pN−1 pN−1 
0... − 2∆x qN−1 − 2 2 1 + 2∆x
∆x 2 ∆x 2
 ∆x

   
0... 0 1 − pN qN − 2 2
2
∆x 2∆x ∆x

(25)
 f1   s1   
 1 p1 
− 2∆x f0
f2 s2 ∆x 2
 ff34   ss34   0
0

F = . S =. B=
     .


 ..   ..   .
0 


. . 1 pN
+ 2∆x fN+1
fN sN ∆x 2

(26)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210 - Computational Physics 13 / 15
 Matrix inversion

The matrices D, S, and B are known

We need to solve for F . This can be done by inverting the matrix D

F = D −1 (S − B) (27)

D is a non-singular N × N matrix, so it can be inverted

The finite difference method converts differential equations into
algebraic equations
The matrix inversion can be done by standard linrary routines which
implement one of the several inversion methods known
For ex. the Gauss-Seidel method

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 14 / 15
 Exercise

Solve the differential equation

f ′′ (x) + sin(x)f ′ (x) + f (x) = 0 (28)

The boundary conditions are f (x = 0) = 0 and f (x = 4π) = 1

Solve with the dinite difference method, take N = 100
The functions p(x) = sin(x), q(x) = 1, s(x) = 0
Plot the solution for f

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210 - Computational Physics 15 / 15