EP4210/EP24040 - Computational Physics

Kirit Makwana
Department of Physics, IIT Hyderabad

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 1 / 34
 Time dependent PDEs

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 2 / 34
 Motivation

Consider the continuity equation

∂ρ
+ ∇ · (ρv ) = 0 (1)
∂t
In a one dimensional space

∂ρ(x, t) ∂(ρu(x, t))

+ =0 (2)
∂t ∂x
Lets assume the flow field is given u(x, t), we have to solve for ρ(x, t)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 3 / 34
 Discretization
Discretize space and time in points xj and t n respectively and denote

ρnj = ρ(xj , t n ) (3)

where xj = j∆x and t n = n∆t

j = 0, 1, 2, ..., N + 1, n = 0, 1, 2, ...
Consider u to be a constant so that
∂ρ(x, t) ∂ρ(x, t)
+u =0 (4)
∂t ∂x
The finite difference form is
ρn+1
j − ρnj ρnj+1 − ρnj−1
= −u (5)
∆t 2∆x
This is taking the forward difference in time and central difference in
space
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 4 / 34
 Forward Euler method

We can find the error using Taylor expansion

∂ρ(xj , t n )
ρn+1
j = ρnj + ∆t + O(∆t 2 ) (6)
∂t
and
ρnj+1 − ρnj−1 ∂ρ(xj , t n )
= + O(∆x 2 ) (7)
2∆x ∂x
Combining these two equations we get

ρn+1 − ρnj ρnj+1 − ρnj−1 ∂ρ(xj , t n ) ∂ρ(xj , t n )

 
j
+u = +u (8)
∆t 2∆x ∂t ∂x
+ O(∆t) + O(∆x 2 ) (9)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 5 / 34
 Forward Euler method

So the forward Euler method is

∆t n
ρn+1 = ρnj − u (ρ − ρnj−1 ) (10)
j
2∆x j+1
Arrange the ρnj in a column vector

ρn1


 ρn2 
n
 
 . 
ρ = (11)

 . 
ρnN

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 6 / 34
 Matrix form

In matrix form this method is

 
0 1/2 0 0 ...
−1/2 0 1/2 0 ...
 
n+1 n ∆t  0 −1/2 0 1/2 ... n
ρ =ρ −u  ρ (12)
∆x  .
 

 . 
.... 0 0 −1/2 0

This can be written as

 
n+1 ∆t
ρ = I −u Q ρn (13)
∆x
= Pρn (14)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 7 / 34
 Stability of method
Because of the simple time stepping scheme, we can say that the
solution is

ρn = P n ρ0 (15)

Lets consider the norm of the vector ρ which can be defined as

sX
|ρ| = ρ2j (16)
j

It can be shown that when n becomes large, the norm tends to

|ρn | → |S(P)|n (17)

where S(P) denotes the spectral radius of the matrix P

This is simply the largest eigenvalue (λi ) of P

S(P) = maxi (λi ) (18)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 8 / 34
 Stability
The matrix P is of the form
 
1 a 0 0 ...
−a 1 a 0 ...
 
 .
P= (19)


 . 
... 0 0 −a 1

Just take a 2X2 matrix of this form gives

    
1 a x1 x
=λ 1 (20)
−a 1 x2 x2

Solving these two equations gives

λ2 − 2λ + 1 + a2 = 0 (21)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 9 / 34
 Stability

The eigenvalue λ is
p
2± 4 − 4(1 + a2 )
λ= (22)
2
The norm of this eigenvalue will always be more than 1
This implies that the norm of ρ will continue to grow to infinity for
large step size
This means the method is ”unstable”

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 10 / 34
 Implicit Euler method
The above method is called as an ”explicit” method - because it
explicitly writes the next step in terms of previous step
The method can be made stable by taking the next time step on the
right hand side of the equation
∆t n+1
ρn+1
j = ρnj − u (ρ − ρn+1
j−1 ) (23)
2∆x j+1
This method can be written as

ρn = T ρn+1 + Bρn+1 (24)

where T takes the derivative and B takes care of boundary terms

This is called as an implicit method because it ρn+1 implicitly
depends on itself
To solve this, we have to calculate the inverse (T + B)−1
It can be shown that this method is stable
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 11 / 34
 Implicit Euler method
Define
u∆t
a≡ (25)
2∆x
Take the usual grid of xn = n∆x with n = 0, 1, 2, 3, ..., N + 1
x0 and xN+1 are the boundary points and there are N internal points
The implicit equation is

ρnj = ρn+1
j + aρn+1 n+1
j+1 − aρj−1 (26)

This is written in matrix form as

 n    n+1 
ρ1 1 a 0 0 .... 0 ρ1
 ρn2  −a 1 a 0 ...  ρn+1 
 2 
  .  + Bρn+1
  
 . = . (27)
    
 .   .  . 
ρnN ... 0 0 −a 1 ρn+1
N

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 12 / 34
 Periodic boundary conditions
Lets apply periodic boundary conditions, i.e.,

ρn+1
0 = ρn+1
N (28)
ρn+1 n+1
N+1 = ρ1 (29)

This comes from imagining multiple copies of the same domain next
to each other
  n+1 
−aρn+1
  
N 0 0 ... 0 −a ρ1
 0  0 0 ... 0 0  ρn+1 
 2 
Bρn+1 = 
  
 .  = .
 
 . 
  (30)
 .  .  . 
aρn+1
1 a 0 ... 0 0 ρn+1
N

Now that T and B are defined, the solution becomes

ρn+1 = (T + B)−1 ρn (31)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 13 / 34
 Exercise

Take the domain x ∈ [0, 1]

Take speed u = 0.1, choose ∆t = 0.001 and N = 1000
The initial density is
  2 
0 x − 0.5
ρ (x) = exp − (32)
0.15

Solve the density continuity equation up to T = 10

Plot the density at intervals of ∆T = 1, i.e., at t = 0, 1, 2, .., 10

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 14 / 34
 Leapfrog method

There are a variety of different combinations of finite differences that

can be tried
These can also lead to stability of the method

ρn+1
j − ρjn−1 ρnj+1 − ρnj−1
= −u (33)
2∆t 2∆x
This is taking the central difference in both space and time
This is the leapfrog method
∆t n
ρn+1 = ρn−1 −u (ρ − ρnj−1 ) (34)
j j
∆x j+1
This depends on both the current and the previous time step and it is
stable

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 15 / 34
 Lax-Friedrichs scheme
Another popular method is the Lax-Friedrichs scheme
1 ∆t n
ρn+1 = (ρnj−1 + ρnj+1 ) − u (ρ − ρnj−1 ) (35)
j
2 2∆x j+1
This is like the explicit Euler method
However, the current step and position is taken as an average of the
neighboring positions
It’s matrix representation is
1 u∆t
2 − 2∆x
 
0 0 0 .... 0
 1 + u∆t 1 u∆t
 2 2∆x 0 2 − 2∆x 0 .... 0
n+1
 n
ρ =  .  ρ (36)

 . 
1 u∆t
0 ... 0 0 2 + 2∆x 0
Calculate eigenvalues of the reduced 2 × 2 matrix
1 u∆t
 
−λ 2 − 2∆x
det 1 u∆t =0 (37)
2 + 2∆x −λ
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 16 / 34
 CFL condition

The eigenvalues are

r
1 u 2 ∆t 2
λ=± 1− (38)
2 ∆x 2
u∆t
We will get |λ| < 1 if −3 < ∆x <5
A more detailed analysis shows that generally the method is stable if
∆x
|u| < (39)
∆t
In other words
∆x
∆t < (40)
|u|

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 17 / 34
 CFL condition

This condition is known as the Courant-Friedrichs-Lewy (CFL)

condition
It physical meaning is that the time step should be such that the
velocity doesn’t carry material (or information) across more than one
spatial grid interval
Practically, there is a safety factor put in while choosing ∆t
∆x
∆t < CN (41)
|u|

where CN < 1 is the Courant number

This condition is quite general across various explicit time step
methods

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 18 / 34
 Nonlinear scheme

We have considered the equation

∂ρ(x, t) ∂(ρ(x, t)u(x, t))

+ =0 (42)
∂t ∂x
We have simplified this into a linear problem by taking u(x, t) to be a
constant
If we consider the full non-linear problem, then the whole ρu term has
to be considered
Modifying the Lax-Friedrichs scheme
1 ∆t n n
ρn+1 = (ρnj−1 + ρnj+1 ) − (u ρ n
− uj−1 ρnj−1 ) (43)
j
2 2∆x j+1 j+1

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 19 / 34
 Burgers equation
Navier Stokes equation
∂v ∇P
+ v · ∇v = − + ν∇2 v (44)
∂t ρ
Ignoring the viscous term gives the Euler equation
Ignoring the pressure term gives the general Burgers equation, and
ignoring both viscosity and pressure gives the Burgers equation
∂v
+ v · ∇v = 0 (45)
∂t
In one spatial dimension

∂v ∂v ∂v 1 ∂(vv )
+v = + =0 (46)
∂t ∂x ∂t 2 ∂x
This is kind of like the nonlinear version of the continuity equation
So the same methods can also be applied to Burgers equation
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 20 / 34
 Boundary conditions

Various boundary conditions are frequently used

The simplest is to set the function at the boundaries ρn0 and ρnN+1
Another boundary condition is to set the derivatives to be zero at the
boundary
ρn0 = ρn1 and ρnN+1 = ρnN
This is useful when the boundaries are far away from region of
interest and can be assumed to be unchanging
Another very frequently used boundary in study of wave phenomena
and self-similar systems is periodic boundary condition
This means ρn0 = ρnN and ρnN+1 = ρn1
This works as well in 2D and 3D

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 21 / 34
 Lax-Wendroff
We can write
∂ρ ∆t 2 ∂ 2 ρ
ρ(x, t + ∆t) = ρ(x, t) + ∆t + + O(∆t 3 ) (47)
∂t 2 ∂t 2
Since the original equation (with constant velocity) is

∂ρ ∂ρ
= −u (48)
∂t ∂x
we can write
∂2ρ ∂2ρ ∂2ρ
 
∂ ∂ρ
2
= −u = −u −u = u2 (49)
∂t ∂t∂x ∂x ∂x ∂x 2

Thus we can write

∂ρ ∆t 2 2 ∂ 2 ρ
ρ(x, t + ∆t) = ρ(x, t) + ∆t + u + O(∆t 3 ) (50)
∂t 2 ∂x 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 22 / 34
 Lax Wendroff

This equation becomes

∂ρ ∆t 2 2 ∂ 2 ρ
ρ(x, t + ∆t) = ρ(x, t) − u∆t + u + O(∆t 3 ) (51)
∂x 2 ∂x 2
This can be discretized as
u∆t n u 2 ∆t 2 n
ρn+1 = ρnj − (ρj+1 − ρnj−1 ) + (ρ − 2ρnj + ρnj−1 ) (52)
j
2∆x 2∆x 2 j+1
This is the Lax-Wendroff scheme which has higher accuracy
To take it to the nonlinear form where u(x, t) is a function of x and t
can be done in the following 2 step way

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 23 / 34
 Richtmeyer 2-step Lax-Wendroff

First define a half-way point solution

n+ 1 1 ∆t n n
ρj+ 12 = (ρnj + ρnj+1 ) − (ρ u − ρnj ujn ) (53)
2 2 2∆x j+1 j+1
Then the full step is taken using this half step
The second step is then
∆t n+ 21 n+ 12 n+ 1 n+ 1
ρn+1
j = ρnj − (uj+ 1 ρj+ 1 − uj− 12 ρj− 12 ) (54)
∆x 2 2 2 2

It can be shown that this is the nonlinear analogue of the

Lax-Wendroff scheme
The 2-step process makes it easier to implement

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 24 / 34
 Exercise

Implement the Richtmeyer 2-step Lax-Wendroff scheme for solving

the Burger’s equation
Take domain of x ∈ [0, 1] and an initial velocity of 0.01 sin(2πx), up
to time T = 20
Take periodic boundary conditions, add an extra ghost point to take
care of the periodic boundary condition
Map the continuity equation to the Burgers equation and identify ρ
and u

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 25 / 34
 Implicit-Explicit methods
The drawback of stable methods is that they are dissipative
This means that the solution decays slowly, for ex. the norm of the
solution is not conserved
To overcome this, a combination of implicit and explicit methods
Explicit method
∆t n
ρn+1 = ρnj − u (ρ − ρnj−1 ) (55)
j
2∆x j+1
Implicit method
∆t n+1
ρn+1
j = ρnj − u (ρ − ρn+1
j−1 ) (56)
2∆x j+1
Combination of implicit-explicit
∆t n+1
ρn+1
j = ρnj − u (ρ + ρnj+1 − ρn+1 n
j−1 − ρj−1 ) (57)
4∆x j+1
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 26 / 34
 Schroedinger equation
At small enough scales, particles are described by a wavefunction
ψ(x, t)
The probability of finding a particle at a particulat location (x, t) is
given by the absolute value square of this wavefunction |ψ(x, t)|2
Therefore, the wavefunction should satisfy
Z
dx|ψ(x, t)|2 = 1 (58)
x

The dynamics of the wavefunction is given by the time- dependent

Schroedinger equation
∂ψ
iℏ = Ĥψ(x, t) (59)
∂t

Ĥ is the Hamiltonian operator, it gives the energy of the particle

This is a time dependent PDE
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 27 / 34
 Normalization
Schroedinger equation

∂ψ(x, t)
iℏ = Ĥψ(x, t) (60)
∂t
∂ψ(x, t) ℏ2 ∂ 2 ψ(x, t)
iℏ =− + V (x)ψ(x, t) (61)
∂t 2m ∂x 2
This equation can be normalized to some suitable units. Lets take
x = aq, where q is dimensionless and a is some typical length scale of
the problem
Define a characteristic time as t0 = 2ma2 /ℏ. Then the time can be
written as t = st0 , where s is dimensionless
The potential can also be made dimensionless by writing

ℏ2
V = Ṽ (62)
2ma2

where Ṽ is dimensionless
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 28 / 34
 Normalized SE
The normalized Schrodinger equation is

∂ 2 ψ(q, s)
 
∂ψ(q, s)
= −i − + Ṽ (q)ψ(q, s) (63)
∂s ∂q 2

The wavefunction satisfies the probability

Z Z
|ψ|2 dx = 1 =⇒ |ψ|2 adq = 1 (64)
Z
=⇒ |ψ̃|2 dq = 1 (65)

√
where ψ̃ = aψ

∂ 2 ψ̃(q, s)
 
∂ ψ̃(q, s)
= −i − + Ṽ (q)ψ̃(q, s) (66)
∂s ∂q 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 29 / 34
 Discretization

Lets say we have to solve this in some domain of q and at the end
points of this domain Ṽ = ∞
Like in finite difference technique, the domain is divided into N + 1
intervals of ∆q
qj , j = 0, ..., N + 1
From the boundary conditions, the wavefunction at the boundary has
to be zero at all times, ψ(q0 , s) = 0 and ψ(qN+1 , s) = 0
Let us also discretize time in steps of ∆s, so s n = n∆s,
n = 0, 1, 2, 3, ...

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 30 / 34
 Discretization
The wave function can thus be discretized as

ψjn = ψ(qj , s n ) (67)

The wave function can be written in the form of a column vector

 n
ψ1
 ψ2n 
 n
n
ψ 
Ψ = 3 (68)
 . 
 
 . 
ψNn

∂ 2 ψjn n
ψj+1 − 2ψjn + ψj−1
n
= + O(∆q 2 ) (69)
∂q 2 ∆q 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 31 / 34
 Matrix form
In matrix form we write the differential form
  n
−2 1 0 0 .... ψ1
1 −2 1 0 .... ψ n
∂ 2 Ψn
  
1   2
=  .   .  = DΨn (70)
∂q 2 ∆q 2 
 .
 
 . 
.... 0 0 1 −2 ψN n

The boundary terms are not needed because ψ is zero there

Similarly, the potential term can be written as
  n
V (q1 ) 0 0 .... ψ1
 0 V (q2 ) 0 ....   ψ2n 
 
n   .  = PΨn
 
Ṽ (qj )ψj =  .

  (71)
 .  . 
.... 0 0 V (qN ) ψN n

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 32 / 34
 Crank Nicolson

The time derivative can similarly be written as

∂ψjn ψjn+1 − ψjn

= (72)
∂s ∆s
Or in matrix form as
∂Ψn Ψn+1 − Ψn
= (73)
∂s ∆s
Now the Schroedinger equations can be written as

Ψn+1 − Ψn
= −i(−D + P)Ψn (74)
∆s

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 33 / 34
 Crank Nicolson method
However, here there is the Crank-Nicolson method which takes an
average on the spatial derivative also

Ψn+1 − Ψn (Ψn + Ψn+1 )

= −i(−D + P) (75)
∆s 2
This method provides some advantages
Then the time evolution becomes
   
I i n+1 I i
+ (−D + P) Ψ = − (−D + P) Ψn (76)
∆s 2 ∆s 2
=⇒ F Ψn+1 = G Ψn (77)

where I is just the identity matrix

Given an initial Ψ0 , we can solve for later times by multiplying F −1
on both sides

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 34 / 34