EP4210/EP24040 - Computational Physics

Kirit Makwana
Department of Physics, IIT Hyderabad

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 1 / 21
 Ordinary Differential Equations (ODEs)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 2 / 21
 Ordinary differential equations ODEs
An ODE is a differential equation in which there is only one
independent variable and all derivatives are with respect to it
For ex., let x be the independent variable and let y be a function of
x, and c1 , c2 , ... are constants
Then any equation of the form

dy d 2 y d 3 y
 
f x, y , , 2 , 3 , ...., c1 , c2 , ... = 0 (1)
dx dx dx

is an ODE
The order of the PDE is the highest order derivative that appears in it
Ex. Poisson equation in 1 dimension

d 2ϕ ρ(x)
2
=− (2)
dx ϵ0

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 3 / 21
 Verlet method for equation of motion
Consider the equation of a motion of a classical object (particle)
under the action of a force - this is provided by Newton’s second law

F
a=m =⇒
d 2x
=
F = a(x , t) (3)
dt 2 m
If the force field F (x , t), which can be a function of space and time,
is known and the mass of the particle is given, then we know the
acceleration a (x , t) also as a function of space and time
Given the starting position and velocity of the particle, it should be
possible to solve this equation to obtain the position of the particle at
any later time
However, even for simple force fields this may not be possible
analytically
This can be done numerically by breaking time and space into
discrete chunks
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 4 / 21
 Verlet method
Lets Taylor expand the position of the particle around time t

x (t + ∆t) = x (t) + ∆t ddtx +

∆t 2 d 2 x
2! dt 2
+ O(∆t 3 ) (4)
t t

x (t − ∆t) = x (t) − ∆t ddtx +

∆t 2 d 2 x
2! dt 2
− O(∆t 3 ) (5)
t t

Adding up these two equations gives

x (t + ∆t) + x (t − ∆t) = 2x (t) + ∆t 2 ddtx2

2
+ O(∆t 4 ) (6)
t

This gives us a way to approximate the position at time t + ∆t using

the information at time t and t − ∆t

x (t + ∆t) = −x (t − ∆t) + 2x (t) + ∆t 2 a(t) + O(∆t 4 ) (7)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 5 / 21
 Verlet method

However, there is a problem at the first time step itself, starting from
t=0
We only know the position and velocity of the particle at time t = 0,
but we also require its position at time −∆t to take the first step
This can again be approximated by the Taylor series
2
x (−∆t) = x (0) − ∆t v (0) + ∆t2 a(0) − O(∆t 3 ) (8)

Knowing x (t = 0), v (t = 0), and a (t = 0), we can get x (−∆t)

The Verlet method given by Eq. 7 is third order accurate, but because
of the third order error in the first step, it becomes second order
accurate
This is the “local” error (i.e., error at every time step)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 6 / 21
 Verlet algorithm

define starting position x0 and velocity v0

define an acceleration function which takes position and
time as input
define the stepsize Dt and number of steps N
Calculate x-=x0-Dt*v0+(Dt**2/2)*a0
start from t=0
Loop over N
x+=-(x-)+2*x0+Dt**2*acc(x0,t)
t=t+Dt
Store the location x+
x-=x0
x0=x+

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 7 / 21
 Exercise

Exercise
Consider a particle of unit mass moving under the action of a force

F = sin(y )x̂ + cos(x)ŷ (9)

The particle starts from rest at the origin at t = 0

Calculate and plot the trajectory of the particle up to time T = 50
using the Verlet method
Check for the convergence of the method - i.e. the solution should
converge for a small enough ∆t

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 8 / 21
 Runge-Kutta (RK) methods
Consider a first order ODE of the form
df
= G (x, f (x)) (10)
dx
It can be shown that any nth order ODE can be converted into an
equation of this type
Suppose the initial condition is given, i.e., f (x = a) = fa
Then the independent variable x can be broken in steps of size ∆x
such that xn = a + n∆x, n = 0, 1, 2, ...
The function is defined at these discrete points fn = f (xn )
By the definition of limits
fn+1 − fn
lim = G (xn , fn ) (11)
∆x→0 ∆x

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 9 / 21
 Euler method

The Euler method simply derives from the first order Taylor expansion

df
fn+1 = fn + ∆x + O(∆x 2 ) (12)
dx n

The Euler method then is given by

f¯n+1 = f¯n + ∆xG (xn , f¯n ) (13)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 10 / 21
 Runge-Kutta methods

Make Taylor series expansion to second order

dfn ∆x 2 d 2 fn
fn+1 = fn + ∆x + + O(∆x 3 ) (14)
dx 2 dx 2
∆x 2 dG (xn , fn )
= fn + ∆xG (xn , fn ) + + O(∆x 3 ) (15)
2 dx
By using the chain rule we get

∆x 2 ∂Gn ∂Gn
 
fn+1 = fn + ∆xG (xn , fn ) + + Gn + O(∆x 3 ) (16)
2 ∂x ∂f

Sometimes, the function G will be known but its derivatives might

not be known or difficult to evaluate
So the derivatives of G are evaluated in the following way

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 11 / 21
 Midpoint RK method

Now consider the Taylor series expansion of G around the following

point
 
∆x ∆x ∆x ∂G
G xn + , fn + Gn =G (xn , fn ) + (17)
2 2 2 ∂x xn
∆x ∂G
Gn + O(∆x 2 )
+ (18)
2 ∂f xn
 
∆x ∆x
=⇒ G xn + , fn + Gn − Gn = (19)
2 2
 
∆x ∂Gn ∂Gn
+ Gn + O(∆x 2 ) (20)
2 ∂x ∂f

The term in the square brackets is exactly the term required in the
Eq. 16

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 12 / 21
 Midpoint RK method

Substituting this term in Eq. 16 we get

 
∆x ∆x
fn+1 = fn + ∆xG xn + , fn + Gn + O(∆x 3 ) (21)
2 2

Ignoring the third order term, we get the midpoint RK method which
is second-order accurate
It can be broken down as

k1 = ∆xG (xn , fn ) (22)

 
∆x k1
k2 = ∆xG xn + , fn + (23)
2 2
fn+1 = fn + k2 (24)

This is also called as the RK2 method

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 13 / 21
 Algorithm

define starting position x0 and f0

define step size dx
define function G(x,f)
define number of steps N
Loop over N
k1=dx*G(x0,f0)
k2=dx*G(x0+dx/2,f0+k1/2)
f1=f0+k2
store f1
x0=x0+dx
f0=f1

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 14 / 21
 RK4 method
RK methods can be generalized to higher order by expanding f (x) to
higher orders in the Taylor series
Then we need to Taylor expand G around judiciously chosen
intermediate points, such that we get the terms in the Taylor
expansion of f (x)
The fourth order RK method is given by

k1 = ∆xG (xn , fn ) (25)

∆x k1
k2 = ∆xG (xn + , fn + ) (26)
2 2
∆x k2
k3 = ∆xG (xn + , fn + ) (27)
2 2
k4 = ∆xG (xn + ∆x, fn + k3 ) (28)
 
1 k1 k4
fn+1 = fn + + k2 + k3 + (29)
3 2 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 15 / 21
 Equation of motion (EOM)
Let’s consider the equation of motion again

d 2r
= a (r , t) (30)
dt 2
This can be converted into coupled 1st order ODEs
dr
= v (r , t) (31)
dt
dv
= a (r , t) (32)
dt
Now we have two dependent variables, r and v
So the dependent variables can be thought of as a vector, giving

d r v (r , t)
   

dt v a(r , t)
= (33)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 16 / 21
 RK2 method for EOM
This equation looks like the standard ODE for RK methods
df
= G (x, f (x)) (34)
dx
except that now f is a vector, the independent variable is t, and G
function is also a vector which depends on t and f
For example in 3D space, the f vector and the G vector will be
   
x vx
y  vy 
   
z   vz 
f = 
vx  G = 
 ax  (35)
   
vy  ay 
vz az

To start the RK method, we require the initial position r0 and initial

velocity v0
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 17 / 21
 Motion of particle in electromagnetic field
Consider the motion of a charged particle in electromagnetic field.
Suppose the electric and magnetic fields are given as a function of
space and time E (r , t) and B (r , t)
The equation of motion is

d 2r
= [E + v × B ]
q
(36)
dt 2 m
Then the G vector becomes
 
vx

 vy 

 vz 
G =
  (37)
 (q/m)(E x + vy Bz − v z By ) 

(q/m)(Ey − vx Bz + vz Bx )
(q/m)(Ez + vx By − vy Bx )

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 18 / 21
 Predictor-Corrector methods

Predictor-Corrector methods are another general class of methods for

solving ODE’s
The general idea here is to first use an Euler step to get a rough
solution and then refine the solution by using the rough solution
The predictor step will be the Euler step

f˜n+1 = fn + ∆xG (xn , fn ) (38)

Then the solution f˜n+1 is used in the corrector step

G (xn , fn ) + G (xn+1 , f˜n+1 )

 
fn+1 = fn + ∆x (39)
2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 19 / 21
 Predictor-Corrector methods

The predictor corrector steps can be repeated in order to increase the

accuracy
For example

f˜n+1 = fn + ∆xG (xn , fn ) (40)

G (xn , fn ) + G (xn+1 , f˜n+1 )
 
fˆn+1 = fn + ∆x (41)
2
G (xn , fn ) + G (xn+1 , fˆn+1 )
 
fn+1 = fn + ∆x (42)
2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 20 / 21
 Exercise

Repeat the particle trajectory exercise of Verlet method with the RK2
method
Consider a particle of unit mass moving under the action of a force

F = sin(y )x̂ + cos(x)ŷ (43)

The particle starts from rest at the origin at t = 0

Calculate and plot the trajectory of the particle up to time T = 50
using the Verlet method
Check for the convergence of the method - i.e. the solution should
converge for a small enough ∆t

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 21 / 21