EP4210/EP24040 - Computational Physics

Kirit Makwana
Department of Physics, IIT Hyderabad

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 1 / 28
 2D PDEs

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 2 / 28
 Finite difference in 2 dimensions
Now lets consider the Poisson equation in 2 dimensions
ρ
∇2 V = − (1)
ϵ0
∂2V ∂2V
=⇒ + = s(x, y ) (2)
∂x 2 ∂y 2
Let the space be discretized on a grid of step size ∆x and ∆y

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 3 / 28
 Discretization

The second order partial derivative can be written as

∂ 2 V (xn , ym ) V (xn + ∆x, ym ) + V (xn − ∆x, ym ) − 2V (xn , ym )

2
=
∂x ∆x 2
(3)
+O(∆x 2 )
(4)

We will get a similar expression for ∂ 2 V /∂y 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 4 / 28
 Jacobi Iteration method
The Poisson equation in discrete form becomes
Vn+1,m + Vn−1,m − 2Vn,m Vn,m+1 + Vn,m−1 − 2Vn,m
+ = sn,m (5)
∆x 2 ∆y 2
This can be re-arranged to give

∆y 2 (Vn+1,m + Vn−1,m ) ∆x 2 (Vn,m+1 + Vn,m−1 )

Vn,m = + (6)
2(∆x 2 + ∆y 2 ) 2(∆x 2 + ∆y 2 )
sn,m ∆x 2 ∆y 2
− (7)
2(∆x 2 + ∆y 2 )

Assume ∆x = ∆y = ∆ we have

Vn+1,m + Vn−1,m + Vn,m+1 + Vn,m−1 − sn,m ∆2

Vn,m = (8)
4

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 5 / 28
 Jacobi iteration method
This form of the finite difference approximation leads to an iterative
procedure for finding the solution
0 at all the grid points satisfying the
Assume a certain function for Vn,m
boundary condition
Then obtain a new solution
0
Vn+1,m 0
+ Vn−1,m 0
+ Vn,m+1 0
+ Vn,m−1 − sn,m ∆2
1
Vn,m = (9)
4
and iterate this for the i th step keeping the boundary values fixed
i
Vn+1,m i
+ Vn−1,m i
+ Vn,m+1 i
+ Vn,m−1 − sn,m ∆2
i+1
Vn,m = (10)
4
It can be shown mathematically that this iterative process will
converge to the correct solution

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 6 / 28
 Convergence

Till how long should we continue this cycle?

For this we define a global error by making use of the l 2 norm
sX
2 i+1 i |2
l = |Vn,m − Vn,m (11)
n,m

If this norm reduces below a certain threshold, then we consider the

iteration converged
However, the convergence with respect to N, M still needs to be done
In the Jacobi method, the new solution at every point is created by
considering the neighboring values of the older solution
A modification is the Gauss-Seidel method, where the newer neighbor
values are utilized if available - however this requires running a loop
over each point, which may be slower

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 7 / 28
 Exercise

Solve the Laplace equation to find the electric potential

∇2 V = 0 (12)

in a square domain x, y ∈ [0, 1] with the boundary conditions

V (x = 0) = sin(2πy ), and the potential on the other boundaries is
zero
Plot the solution with a color plot
Take N = M = 100 (100 × 100 internal points)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 8 / 28
 Grid
Lets consider a 2D grid as shown below

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 9 / 28
 Grid ordering
We have to solve the Laplace equation inside this grid. The boundary
condition is given to us, i.e., the potential V is known on the
boundary grid points
We have labelled the internal grid points as i = 1, 2, 3, ... with the
ordering (1, 1), (2, 1), (3, 1), ..., (n, m), ...(N, M), where N = 3 in this
case
n is the index number along x direction and m is the index along y
direction
We are taking a rectangular domain for simplicity but this method
can work for any domain shape
The relation between i and (n, m) is

i = (m − 1) ∗ N + n (13)

The inverse relation to obtain (n, m) from i can also be found

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 10 / 28
 Creation of matrices

We want to find the solution matrix, which is a 1-dimensional column

vector of the function values Vn,m with the ordering of index i
 
V1,1
 
V1  V2,1 
 
 V2   . 
   
 .   VN,1 
 .  =  V1,2 
F =    (14)
   
 .   . 
 
VN×M  . 
VN,M

Now lets create a matrix Dx that can act on this vector F and it can
give us the second order partial derivative w.r.t. x, as defined in Eq. 3

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 11 / 28
 Matrix creation
Lets write down a few derivatives explicitly

∂ 2 V1,1 V2,1 + V0,1 − 2V1,1

2
= (15)
∂x ∆x 2

∂ 2 V2,1 V3,1 + V1,1 − 2V2,1

2
= (16)
∂x ∆x 2

∂ 2 V3,1 V4,1 + V2,1 − 2V3,1

= (17)
∂x 2 ∆x 2
Or more generally

∂ 2 Vn,m Vn+1,m + Vn−1,m − 2Vn,m

2
= (18)
∂x ∆x 2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 12 / 28
 Dx matrix

When forming the Dx matrix, we will not consider the boundary grid
points because they will be taken care of separately
−2 1 0 ... 0 0 0 ... 0 0 0 0   V1,1 
1 −2 1 ... 0 0 0 ... 0 0 0 0 V2,1
.  . 
 ..0 0 0 ... 1 −2 0 ... 0 0 0 0   VN,1 
.   V1,2 
1  .  . 
Dx F = 2
.  . 
∆x  0.

0 0 ... 0 ... 0 −2 1 0 ...
 VN,M−1 
0   V1,M 
0 ... 0 0 0 ... 0 1 −2 1 ... 0   V2,M 
.  . 
. .
0 .. 0 0 ... 0 0 0 ... 0 1 −2 VN,M
(19)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 13 / 28
 Algorithm

define N and M
define dx
define NXM array D with zeros
loop n over 1 to N
loop m over 1 to M
define the 1D label i = (m-1)*N+n
Make the diagonal element D[i-1,i-1]=-2
if we are not on the left boundary (n>1)
Make the left neighbor D[i-1,i-2]=1
if we are not on the right boundary (n<N)
Make the right neighbor D[i-1,i]=1

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 14 / 28
 Dy

Similarly, the partial derivatives in y direction are

−2 0 0 ... 1 0 0 ... 0 0 0 0
 V1,1

0 −2 0 ... 0 1 0 ... 0 0 0 0 V2,1
 .  .
 .  VN,1

 ..0 1 0 ... −2 0 ... 1 ... 0 0 0  V1,2

1  .  .
Dy F =  .  .
∆y 2 
 .
0 0 0 ... 1 ... 0 −2 0 ... 0
 VN,M−1 
0   V1,M 
 .   V2,M 
 .  . 
. .
0 .. 0 0 ... 0 0 1 ... 0 0 −2 VN,M
(20)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 15 / 28
 Algorithm to create Dy

define N and M
define dx
define NXM array D with zeros
loop n over 1 to N
loop m over 1 to M
define the 1D label i = (m-1)*N+n
Make the diagonal element D[i-1,i-1]=-2
if we are not on the bottom boundary (m>1)
index of the bottom neighbor j=(m-2)*N+n
D[i-1,j-1]=1
if we are not on the top boundary (m<M)
index of the top neighbor j=(m)*N+n
D[i-1,j-1]=1

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 16 / 28
 Boundary conditions

The laplacian matrices will act on F (Dx + Dy )F to give the finite

difference form of the Laplacian operator
However, the boundary conditions are still missing, we need to provide
those
The boundary points are (n, m) = (1, :), (N, :), (:, 1), and (:, M)
For (1, :) we have to add the left boundary point while for (N, :) we
have to add the right boundary points
For (:, 1) we have to add lower boundary points and for (:, M) we
have to add upper boundary points
Define 4 column vectors for taking care of these 4 boundaries
The values of the potential are known on the boundaries V (0, :),
V (N + 1, :), V (:, 0), and V (:, M + 1)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 17 / 28
 Boundary matrices

Define four boundary matrices

 V (0, 1)   0 
0 .
 .   . 
. V (N + 1, 1)
1  V (0, 2)  Bright = 1 
  .

Bleft =  . 0
 (21)
∆x 2 
 .


2
∆x  V (N + 1, 2)


V (0, M) .
. .
0 V (N + 1, M)

 V (1, 0)   0 
V (2, 0) 0
. .
1  V (N, 0)  Btop = 1 
  .

Bbot =  0 V (1, M + 1)
 (22)
∆y 2  .
 2
∆y  V (2, M + 1)

. .
0 V (N, M + 1)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 18 / 28
 Boundary condition handling

4 m a t r i c e s B l e f t , B r i g h t , Bbottom , Btop o f s i z e NXM

l o o p i from 1 t o N
l o o p j from 1 t o M
row num=( j = 1) *N+i =1
i f ( i ==1) B l e f t [ row num ]= V [ 0 , j ] / d e l t a x ˆ2
i f ( i==N) B r i g h t [ row num ] = V [ N+1, j ] / d e l t a x ˆ2
i f ( j ==1) Bbottom [ row num ] = V [ i , 0 ] / d e l t a y ˆ2
i f ( j==M) Btop [ row num ] = V [ i ,M+1]/ d e l t a y ˆ2

B = Bleft + Bright + Bbot + Btop (23)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 19 / 28
 Solution setup

Similarly we will setup the source matrix S

 S(1, 1) 
S(2, 1)
 S(N,. 1) 
 S(1, 2) 
S = 
 S(N,. 2)  (24)
 . 
.
S(N, M)

Then the finite difference equations become

Dx F + Dy F + B = S (25)

F = (Dx + Dy )−1 (S − B) (26)

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 20 / 28
 Arbitrary geometry
The finite difference technique can be aplied in an arbitrary shaped
domain also

Consider the internal points (solid circles) surrounded by boundary

points (empty circles)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 21 / 28
 Mapping matrix

There is no clear relationship between the 2D indices (n, m) and the

1D label i = 1, 2, 3, ...
Instead, lets create a 2D matrix which is a replica of the rectangular
shaped region that contains the arbitrarily shaped domain
 
0 0 0 0 0
0 0 6 0 0
 
M= 0 3 4 5 0  (27)
0 0 1 2 0
0 0 0 0 0

This is an N × M matrix which is the smallest rectangular matrix that

encompasses the region of interest
This matrix has to created manually because it depends on the
arbirtrary shape of the region of interest

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 22 / 28
 Creation of D matrices
The F vector will contain K number of points (K = 6 in this case)
which are the internal points where we want to find the solution
So we will create a K element array to store the solution
The matrices Dx and Dy will be K × K matrices
To create the Dx matrix
loop n over 0 to N-1
loop m over 0 to M-1
if M[n,m]!=0 (so we are at an internal point)
rownum=M[n,m]-1
Dx[rownum,rownum]=-2
if M[n+1,m]!=0 (we are not at the right boundry)
colnum=M[n+1,m]-1
Dx[rownum,colnum]=1
if M[n-1,m]!=0 (not on left boundary)
colnum=M[n-1,m]-1
Dx[rownum,colnum]=1
Dx=Dx/deltax**2
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/EP24040 - Computational Physics 23 / 28
 Creation of D matrices

Similarly for Dy
loop n over 0 to N-1
loop m over 0 to M-1
if M[n,m]!=0 (so we are at an internal point)
rownum=M[n,m]-1
Dy[rownum,rownum]=-2
if M[n,m-1]!=0 (we are not at the top boundry)
colnum=M[n,m-1]-1
Dy[rownum,colnum]=1
if M[n,m+1]!=0 (not on bottom boundary)
colnum=M[n,m+1]-1
Dx[rownum,colnum]=1
Dy=Dy/deltay**2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 24 / 28
 Handling boundaries

The boundary conditions will be a K × 1 column vector

For this we can again create a map of the boundary conditions in a
2D matrix
For the example given above, this will look like
 
0 0 V3,5 0 0
 0 V2,4 6 V4,4 0 
 
B2 = V1,3
 3 4 5 V5,3 
 (28)
 0 V2,2 1 2 V5,2 
0 0 V3,2 V4,1 0

Again, this is something that will have to be created manually, or via

some input data

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 25 / 28
 Boundary handling

The boundary matrix B can be initialized with zeros and then can be
filled
loop n over 0 to N-1
loop m over 0 to M-1
if M[n,m]!=0 (so we are at an internal point)
rownum=M[n,m]-1
if (M[n-1,m]==0)(left boundary)
B[rownum]=B[rownum]+B2[n-1,m]/deltax**2
if (M[n+1,m]==0)(right boundary)
B[rownum]=B[rownum]+B2[n+1,m]/deltax**2
if (M[n,m+1]==0)(bottom boundary)
B[rownum]=B[rownum]+B2[n,m+1]/deltay**2
if (M[n,m-1]==0)(top boundary)
B[rownum]=B[rownum]+B2[n,m-1]/deltay**2

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 26 / 28
 Full solution

The finite difference version of the Poisson equation in 2D will be

(Dx + Dy )F + B = S (29)

So the solution becomes

F = (Dx + Dy )−1 (S − B) (30)

The solution F can be mapped back to 2D space by using the

mapping matrix M
This method can be generalized to 3 dimensions also and to arbitrary
shapes

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 27 / 28
 Applications of finite difference methods

Solution of heat equation in steady state

∇2 T = 0 (31)

Solving this equation can be used in modeling engines, heaters,

cooling systems in electronics
The fluid pressure in incompressible hydrodynamics obeys this
equation and is used in aerodynamic design
Gravitational potential also obeys this equation, so it has applications
in astrophysics
More generally, many other PDEs can be solved by finite difference
formulation, including time dependent PDEs

Kirit MakwanaDepartment of Physics, IIT Hyderabad

EP4210/EP24040 - Computational Physics 28 / 28