MEC503 Computer Techniques February 16, 2022

Lecture 4: FINITE DIFFERENCE OF ELLIPTIC EQUATION

University of Mosul Younis Najim Dep. of Mechanical Engineering

Typical examples of elliptic PDE are;

★ Laplace equation (𝑢 𝑥𝑥 + 𝑢 𝑦𝑦 = 0) which can be used to model a variety of problems involving the potential
of an unknown variable, steady two-dimensional heat conduction.
★ The Poisson’s equation (𝑢 𝑥𝑥 + 𝑢 𝑦𝑦 = 𝑓 (𝑥, 𝑦)) which can be a model to potential field caused by a given
electric charge or mass density distribution or steady, 2D heat conduction with heat source or sink, stress
equation for torsion of an elastic rod, and the displacement of an membrane under constant pressure.
★ The Helmholtz equation (𝑢 𝑥𝑥 + 𝑐𝑢 = 𝑓 (𝑥, 𝑦)) which represents a time-independent, 2D form of the wave
equation.

4.1 BOUNDARY CONDITIONS FOR ELLIPTIC PDE EQUATION

The most common boundary value conditions from for two-dimensions domain bounded by the boundary
curve C are:
1. Dirichlet-condition 𝑢(𝑥, 𝑦) = 𝐺 1 (𝑥, 𝑦) on the boundary cure C.
2. Neumann-condition 𝜕𝑢 𝜕𝑛 = 𝐺 2 (𝑥, 𝑦) on the boundary cure C.
3. Robin (mixed) condition𝑎 𝑢(𝑥, 𝑦) + 𝑏 𝜕𝑢 𝜕𝑛 = 𝐺 3 (𝑥, 𝑦) on the boundary cure C.

4.2 THE LAPLACE EQUATION

Use finite difference formulation to solve the steady two-dimensional heat conduction problem described as:

𝜕 2𝑇 𝜕 2𝑇
+ =0 (1)
𝜕𝑥 2 𝜕𝑦 2

4-1
4-2

Central differences based on the grid scheme:

𝜕 2𝑇 𝑇𝑖+1, 𝑗 − 2𝑇𝑖, 𝑗 + 𝑇𝑖−1, 𝑗 𝜕 2𝑇 𝑇𝑖, 𝑗+1 − 2𝑇𝑖, 𝑗 + 𝑇𝑖, 𝑗−1

2
= 2
+ O(Δ𝑥) 2 , 2
= 2
+ O(Δ𝑦) 2
𝜕𝑥 Δ𝑥 𝜕𝑦 Δ𝑦
substitute into (1):
𝑇𝑖+1, 𝑗 − 2𝑇𝑖, 𝑗 + 𝑇𝑖−1, 𝑗 𝑇𝑖, 𝑗+1 − 2𝑇𝑖, 𝑗 + 𝑇𝑖, 𝑗−1
+ =0
Δ𝑥 2 Δ𝑦 2
For the square grid Δ𝑥 = Δ𝑦 ⇒ 𝑇𝑖+1, 𝑗 + 𝑇𝑖−1, 𝑗 + 𝑇𝑖, 𝑗+1 + 𝑇𝑖, 𝑗−1 − 4𝑇𝑖, 𝑗 = 0 (2)
★ This relationship, which holds for all interior points on the plate, is referred to as the Laplacian difference
equation.
★ Boundary conditions must be specified to obtain a unique solution. The simplest case is where the
temperature at the boundary is set at a fixed value.

For the case illustrated in figure above: To solve this system in Matlab (review Matlab lectures)
[ 𝐴] [𝑇] = [𝐵] ⇒ [ 𝐴] −1 [ 𝐴] [𝑇] = [ 𝐴] −1 [𝐵]. Use inv() to get matrix inverse in matlab.
Syntax
 4-3
node (i,j) difference equation
1,1 4𝑇11 − 𝑇21 − 𝑇12 =75
2,1 −𝑇11 + 4𝑇21 − 𝑇31 − 𝑇22 =0
3,1 −𝑇21 + 4𝑇31 − 𝑇32 =50
1,2 −𝑇11 + 4𝑇12 − 𝑇22 − 𝑇13 =75
2,2 −𝑇21 − 𝑇12 + 4𝑇22 − 𝑇32 − 𝑇23 =0
3,2 −𝑇31 − 𝑇22 + 4𝑇32 − 𝑇33 =50
3,1 −𝑇12 + 4𝑇13 − 𝑇23 =175
3,2 −𝑇22 − 𝑇13 + 4𝑇23 − 𝑇33 =100
3,3 −𝑇32 − 𝑇23 + 4𝑇33 =150

 4 −1 0 −1 0 0 0 0 0  𝑇11   75 
    
−1 4 −1 0 −1 0 0 0 0  𝑇21   0 
    
 0 −1 4 0 0 −1 0 0 0  𝑇31   50 
    
−1 0 0 4 −1 0 −1 0 0  𝑇12   75 
    
𝐴 =  0 −1 0 −1 4 −1 0 −1 0  𝑇22  =  0  (3)
 0 0 −1 0 −1 4 0 0 −1 𝑇32   50 
    
 0 0 0 −1 0 0 4 −1 0  𝑇13  175
    
 0 0 0 0 −1 0 −1 4 −1 𝑇23  100
    
 0 0 0 0 0 −1 0 −1 4  𝑇33  150
    
A=[4,-1,0,-1,0,0,0,0,0;-1,4,-1,0,-1,0,0,0,0;0,-1,4,0,0,-1,0,0,0;
-1,0,0,4,-1,0,-1,0,0;0,-1,0,-1,4,-1,0,-1,0;
0,0,-1,0,-1,4,0,0,-1;0,0,0,-1,0,0,4,-1,0;
0,0,0,0,-1,0,-1,4,-1;0,0,0,0,0,-1,0,-1,4];
B=[75;0;50;75;0;50;175;100;150];
T=inv(A)*B;

For larger-sized grids and since significant number of the terms will be zero. When applied to such
large systems, full-matrix elimination methods (such as Gauss Elimination Method, Matrix inversion,
Cramer’s rule) waste great amounts of computer memory storing these zeros. For this reason, approximate
methods provide a viable approach for obtaining solutions for elliptical equations. The most commonly
employed approach is the successive iteration method of Gauss-Seidel, which when applied to PDEs is also
referred to as Liebmann’s method.

4.2.1 Neumann boundary conditions

Consider a heat conduction problem where we impose Neumann boundary conditions ( 𝜕𝑇
𝜕𝑥 | 𝑥=0 = 0,
𝜕𝑇
𝜕𝑦 | 𝑦=0 =
0), as shown below.
4-4
homogeneous Neumann boundary conditions at:
𝜕𝑇 𝜕𝑇
=0 for 𝑥 = 0 and 0 < 𝑦 < 1, | 𝑦=0 = 0 for 𝑦 = 0 and 0 < 𝑥 < 1
𝜕𝑥 𝜕𝑦
Dirichlet boundary conditions for:

𝑇 =1 for 𝑦 = 1 and 0 ≤ 𝑥 ≤ 1, 𝑇 =0 for 𝑥 = 1 and 0 ≤ 𝑦 ≤ 1

𝜕𝑇 𝑇1, 𝑗 − 𝑇−1, 𝑗
Assume Δ𝑥 = Δ𝑦, 𝑎𝑡 = 0, at 𝑖 = 0 yields: =0 ⇒ 𝑇−1, 𝑗 = 𝑇1, 𝑗
𝜕𝑥 2Δ𝑥
where we have introduced ghostcells with negative indices outside the physical domain to express the
derivative at the boundary.
Or you can use forward differences approximation for ( 𝜕𝑇
𝜕𝑥 = 0) of similar order such that:

𝜕𝑇 −3𝑇𝑖, 𝑗 + 4𝑇𝑖+1, 𝑗 − 𝑇𝑖+2, 𝑗 𝜕𝑇 4𝑇1, 𝑗 − 𝑇2, 𝑗

|𝑖, 𝑗 = at𝑥 = 0, 𝑖 = 0, =0 ⇒ 𝑇0, 𝑗 =
𝜕𝑥 2Δ𝑥 𝜕𝑥 3

Once all 𝑇0, 𝑗 and 𝑇𝑖,0 are determined, the solution can be found using similar steps in previous example. Do
it at home!.

4.2.2 JACOBI & GAUSS-SEIDEL (LIEBMANN) METHODS

The Gauss-Seidel method is the most commonly used iterative method.
" 𝑖−1
#
1 ∑︁
Jacobi : 𝑥𝑖(𝑘+1) = − 𝑎𝑖, 𝑗 𝑥 (𝑘)
𝑗 + 𝑏𝑖 , 𝑖 = 1, 2..𝑛, 𝑘 ≥ 0
𝑎𝑖𝑖 𝑗=1, 𝑗≠𝑖

With the Gauss-Seidel method, we use the new values 𝑥𝑖𝑘+1 as soon as they are known. For example, once
we have 𝑥 1𝑘+1 computed from the first equation, its value is then used in the second equation to obtain the
new 𝑥 2𝑘+1 and so on. Hence, the Gauss-Seidel is faster than Jacobi method.
" 𝑖−1 𝑛
#
1 ∑︁ ∑︁
Gauss-Seidel : 𝑥𝑖(𝑘+1) = − 𝑎𝑖, 𝑗 𝑥 (𝑘+1)
𝑗 − 𝑎𝑖 𝑗 𝑥 (𝑘)
𝑗 + 𝑏𝑖 , 𝑖 = 1, 2..𝑛, 𝑘 ≥ 0
𝑎𝑖𝑖 𝑗=1 𝑗=𝑖+1
 4-5
Assume that we are given a set of 3 equations:

Liebmann’s method:
Eq(2) can be written as:
𝑇𝑖+1, 𝑗 + 𝑇𝑖−1, 𝑗 + 𝑇𝑖, 𝑗+1 + 𝑇𝑖, 𝑗−1
𝑇𝑖, 𝑗 =
4
Solved iteratively for 𝑗 = 1 𝑡𝑜 𝑛 and 𝑖 = 1 𝑡𝑜 𝑚.
★ the system will eventually converge on a stable solution when the coefficient matrix (A) is diagonal
dominant (4 on diagonal and 0, -1 everywhere).
Example 4.2.1 Consider steady two-dimensional heat transfer in a square cross section (3 cm × 3 cm) with
the prescribed temperatures at the top, right, bottom, and left surfaces to be 100°C, 200°C, 300°C, and
500°C, respectively. Using a uniform mesh size Δ𝑥 = Δ𝑦. (a) Derive the finite difference equations and (b)
Write Matlab code to find the nodal temperatures with the Gauss-Seidel iterative method.

%MEC503: This program to solve the 2D, steady heat conduction problem
clear all; close all; clc
4-6
n=4;
w=3;h=3;
x=linspace(0,w,n);y=linspace(0,h,n);
T=zeros(n,n);
%Assign boundary conditions
T(1:n,1)=500; %left wall
T(1:n,n)=200; %right wall
T(1,1:n)=100; %top wall
T(n,1:n)=300; %bottom wall
%Start
errorn=1;tol=1e-5;k=1;
while errorn>tol
k=k+1;
T_old=T;
for i=2:n-1;
for j=2:n-1;
T(i,j)=0.25*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1));
end
end
errorn=max(max(abs(T-T_old)));
end
figure1 = figure(’Color’,[1 1 1]);
contourf(x,y,T,’LevelStep’,30);xlabel(’x-axis’);ylabel(’y-axis’);colorbar;
title(’MEC503: Steady 2D Heat Conduction’)

4.2.3 Successive Over Relation (SOR) METHODS

★ Overrelaxation is sometimes employed to accelerate the rate of convergence by applying the following
formula after each iteration:
𝑇𝑖,new new old
𝑗 = 𝜔𝑇𝑖, 𝑗 + (1 − 𝜔)𝑇𝑖, 𝑗

" 𝑖−1 𝑛
#
𝜔 ∑︁ ∑︁
SOR : 𝑥𝑖(𝑘+1) = (1 − 𝜔)𝑥𝑖(𝑘) + − 𝑎𝑖, 𝑗 𝑥 (𝑘+1)
𝑗 − 𝑎𝑖 𝑗 𝑥 (𝑘)
𝑗 + 𝑏𝑖 , 𝑖 = 1, 2..𝑛, 𝑘 ≥ 0
𝑎𝑖𝑖 𝑗=1 𝑗=𝑖+1

𝑇𝑖,new new
𝑗 values from present iteration, 𝑇𝑖, 𝑗 values from previous iteration, and 1 ≤ 𝜔 ≤ 2.
★ The iterations are repeated until the absolute values of all the percent relative errors (𝜀 𝑎 )𝑖 , j fall below a
prespecified stopping criterion (𝜀 𝑠 ) These percent relative errors are estimated by:
 
 𝑇𝑖,new −𝑇𝑖,old𝑗 
(𝜀 𝑎 )𝑖, 𝑗 =  𝑇 new  100%
𝑗
𝑖, 𝑗

Revoke example in Section 4.2 and solve with initial values =0, and 𝜔 = 1.5, 𝑇𝑖, 𝑗 = (𝑇𝑖+1, 𝑗 + 𝑇𝑖−1, 𝑗 + 𝑇𝑖, 𝑗+1 +
𝑇𝑖, 𝑗−1 )/4:
𝑖 = 1, 𝑗 = 1 → 𝑇11 𝑛 = (0+75+0+0)/4 = 18.75, 𝑇11 = 𝜔(𝑇11𝑛 )+(1−𝜔)𝑇 𝑜 = 1.5(18.75)−0.5(0) = 28.125
11
𝑖 = 2, 𝑗 = 1 → 𝑇21 = (0 + 28.125 + 0 + 0)/4 = 7.03125,
𝑛 𝑛 ) + (1 − 𝜔)𝑇 𝑜 = 1.5(7.03125) −
𝑇21 = 𝜔(𝑇21 21
0.5(0) = 10.54688
𝑖 = 3, 𝑗 = 1 → 𝑇31 𝑛 = (50+10.54688+0+0)/4 = 15.13672, 𝑛 )+(1−𝜔)𝑇 𝑜 = 1.5(15.13672)−
𝑇31 = 𝜔(𝑇31 11
0.5(0) = 22.70508
The computation is repeated for the other rows to give 𝑇12 = 38.67188, 𝑇22 = 18.45703, 𝑇32 =
34.18579, 𝑇13 = 80.12696, 𝑇23 = 74.46900, 𝑇33 = 96.99554
𝑗 = 0, all (𝜀 𝑎 )𝑖, 𝑗 of first iteration are zeros. Repeat the process for the second iteration,
Since all initial 𝑇𝑖,𝑜𝑙𝑑
 4-7
the results are:
𝑇11 = 32.51953, 𝑇21 = 22.35718, 𝑇31 = 28.60108𝑇12 = 57.95288, 𝑇22 = 61.63333, 𝑇32 =
71.86833𝑇13 = 75.21973, 𝑇23 = 87.95872, 𝑇32 = 67.68736
The error for T1,1 (𝜀 𝑎 )1,1 can be estimated:
 𝑇 new −𝑇 old   
(𝜀 𝑎 )1,1 =  1,1𝑇 new1,1  100% ⇒ (𝜀 𝑎 )1,1 =  32.51953−28.12500
32.51953
 100% = 13.5%
1,1
Then continue to 𝜀 𝑎2,1 . 𝜀 𝑎3,1 , ....𝜀 𝑎3,3 and check that max(𝜀 𝑎𝑖, 𝑗 ) is less than 𝜀 𝑠

Example 4.2.2 Write Matlab program to solve the system:

10𝑥 1 − 𝑥 2 + 2𝑥 3 = 6
−𝑥1 + 11𝑥 2 − 𝑥3 + 3𝑥4 = 25
2𝑥1 − 𝑥 2 + 10𝑥3 − 4𝑥 4 = −11
3𝑥2 − 𝑥 3 + 8𝑥 4 = 15

Using Gauss-Seidel Successive iteration method.

% Solution of x in Ax=b using Gauss Seidel successive iteration Method
close all;clear;clc
A=[10,-1,2,0;-1,11,-1,3;2,-1,10,-4;0,3,-1,8];
b=[6;25;-11;15];%a column matrix b
[𝑚, 𝑛]=size(A);% dimension matrix A (m number of row and n number of col)
x=zeros(n,1);% initial guess, form matrix x with zeros, n rows and 1 col
err=Inf; %convergence criterion do iteration as long as normval<tol
tol=0.00001; itr=0;
%start iterations;
while err>tol
x_old=x;
for i=1:n
sigma=0;
for j=1:i-1 % before current x
sigma=sigma+A(i,j)*x(j);
end
for j=i+1:n %after current x
sigma=sigma+A(i,j)*x_old(j);
end
x(i)=(1/A(i,i))*(b(i)-sigma);
end
itr=itr+1;
err=norm(x_old-x);
end
%fprintf(’Solution of the system is : \ n%f\ n%f \ n%f\ n%f\ n
’Number of iteration is %d iterations’,x,itr);% print results here
4-8
4.2.4 OPTIMAL OVER RELAXATION
We have shown that:
𝑢𝑖,𝑘 𝑗 = 𝜔𝑢𝑖,𝑛 𝑗 + (1 − 𝜔)𝑢𝑖,𝑘−1
𝑗

The convergence behavior of the SOR iteration will be heavily dependent on the choice of 𝜔.
★ The optimal relaxation parameter (𝜔 𝑜 𝑝𝑡 ) value that results in the lowest possible iteration number for 𝜔
hold constant throughout the computation.
★ For Laplace and Poisson equations in a rectangular domain it is possible to calculate an optimal relaxation
parameter (𝜔 𝑜 𝑝𝑡 ) , which we will call as theoretically optimal.
suppose we have a rectangular domain of length 𝐿 𝑦 and width 𝐿 𝑥 . Taking the same space discretization in
𝐿𝑥
both directions i.e (Δ𝑥 = Δ𝑦), then the number of grid point along 𝑥-axis is 𝑛𝑥 = Δ𝑥 and along 𝑦-axis is
𝐿𝑦
𝑛 𝑦 = Δ𝑦 . The theoretical optimal 𝜔 𝑜 𝑝𝑡 is then given by:

1h 𝜋
𝜋
i
𝜌= cos + cos
2 𝑛𝑥 𝑛𝑦

2
𝜔 𝑜 𝑝𝑡 = √︃
1 + 1 − 𝜌 2𝑗

Here 𝜌( 𝐴), is the largest (in magnitude) of all the eigenvalues of 𝐴𝑛×𝑛 . Formally, 𝜌( 𝐴𝑛×𝑛 ) = 𝑚𝑎𝑥|𝜆|,
where the maximum is taken over all the eigenvalues of 𝐴𝑛×𝑛 .

If Δ𝑥 ≠ Δ𝑦 the following holds:

Δ𝑥 2
𝜋


𝜋


cos 𝑛𝑥 + Δ𝑦 cos 𝑛 𝑦
𝜌=
Δ𝑥 2


1+ Δ𝑦
Moreover, if the domain is not rectangular, one can use Garabedian’s estimate:

2
𝜔 𝑜 𝑝𝑡 =
1 + 3.014 𝐴ℎ𝑟

Where 𝐴𝑟 is the rectangular domain area and here Δ𝑥 = Δ𝑦 = ℎ

The procedure to find the optimal over relaxation factor is explained in the follwoing example:

4 3 0 
 
Example 4.2.3 Find the optimal choice of 𝜔 for the SOR method for the matrix. 𝐴 = 3 4 −1
0 −1 4 
 
The procedure:
1- The matrix
 must be tridiagonal, and positive definite.
4 3 0 
 
Det(A) = 3 4 −1 = 4(4(4) − (−1)(−1)) − 3(3(4) − 0(−1) + 0(3(−1) − 0(4) = 24
0 −1 4 
 
4 3 
Det(sub A) =   = 4(4) − (3)(3) = 7
3 4  
Det(sub sub A) = 4 = 4
∴ the matrix A is positive definite.
2- Find 𝑇 = 𝐷 −1 (𝐿 + 𝑈) where L is the lower triangular matrix and U is the upper triangular matrix and D
is the diagonal matrix of A.
 1 0 0   0 −3 0  0 −0.75 0 
−1
4 1   
𝑇 = 𝐷 (𝐿 + 𝑈) = 𝐴 =  0 4 0  −3 0 1 = −0.75 0 0.25
 0 0 1   0 1 0  0 0.25 0 
 4   
 4-9
3-Find Eigenvalue
 of the T which denotes
 by 𝜆 and can be found as:
 −𝜆 −0.75 0 
  √ √
= |𝑇 − 𝜆𝐼 | = −0.75 −𝜆 0.25 = −𝜆(𝜆2 − 0.625), 𝜆1 = 0, 𝜆2 = 0.625, 𝜌(𝑇𝐽 ) = 0.625
 0 0.25 −𝜆 
4-Find the optimal 𝜔 𝑜 𝑝𝑡 : 𝜔 𝑜 𝑝𝑡 = √︃2 2 = √ 2 = 1.24
1+ 1−𝜌 𝑗 1+ 1−0.625

Example 4.2.4 Write Matlab program to find the optimal over relaxation factor for the previous example.

%MEC503: find the optimal 𝜔.

close all;clear all; clc
A=[4,-1,0,-1,0,0,0,0,0;-1,4,-1,0,-1,0,0,0,0;0,-1,4,0,0,-1,0,0,0;
-1,0,0,4,-1,0,-1,0,0;0,-1,0,-1,4,-1,0,-1,0; 0,0,-1,0,-1,4,0,0,-1;
0,0,0,-1,0,0,4,-1,0; 0,0,0,0,-1,0,-1,4,-1;0,0,0,0,0,-1,0,-1,4];
DA=diag(diag(A));
T=inv(DA)*(A-DA);
EIGA=eig(T); EIGA_MAX=max(abs(EIGA)); omega_opt=2/(1+sqrt(1-EIGA_MAX^2));
disp(’The optimal w is’)
disp(omega_opt)
run
The optimal w is
1.1716