Santanu Manna MA 204N: Finite Difference Method

Elliptic Equations:
♦ When a system represented by a parabolic equation reaches in steady state then the ultimate state is
defined by an elliptic equation. Initial condition does not play any role in an elliptic equation as the
solution is fully dependent on the boundary conditions. The most familiar elliptic equations are:

∂2u ∂2u
Laplace equation: ∇2 u ≡
+ 2 =0 (1)
∂x2 ∂y
2 2
∂ u ∂ u
Poisson equation: ∇2 u ≡ 2
+ = f (x, y) (2)
∂x ∂y 2
defined over D in the x-y plane and appropriate boundary conditions are prescribed on boundary curve
C; f (x, y) is a known function in D.

For solving the elliptic equation we subdivide the domain into rectangular mesh and discretize the p.d.e.
at the mesh points. This leads to a set of linear system of equations which can be solved by Gaussian
elimination or Gauss-Seidel method.

1 Solution of Laplace equation:

Suppose, we have ∂u ∂u
∂x and ∂y
Discretizing the continuous variable x by xi and y by yi where representing u(x, y) as ui,j , we may write
 
∂u ui+1,j − ui−1,j
= + O(h2 ) (3)
∂x xi 2h
where xi+1 − xi = h
Similarly, we may write  
∂u ui,j+1 − ui,j−1
= + O(k 2 ) (4)
∂y yj 2k
Now adding Eqs. (3) and (4) again, we get

∂2u u(x + h, y) − 2u(x, y) + u(x − h, y)

2
= + O(h2 )
∂x h2
or

∂2u ui+1,j − 2ui,j + ui−1,j

2
= + O(h2 ) (5)
∂x h2

1
Santanu Manna MA 204N: Finite Difference Method

and similarly,

∂2u ui,j+1 − 2ui,j + ui,j−1

2
= + O(k 2 ) (6)
∂y k2
2
∂ u
We may derive the expression for ∂x∂y in finite difference as discussed below:
 2 
∂ u 1
= [(ui+1,j+1 − ui+1,j−1 ) − (ui−1,j+1 − ui−1,j−1 )] + O(hk) (7)
∂x∂y xi ,yj 4hk

1.1 Representation of Derivatives by Diagram

Let xi = xi−1 + h and yj = yj−1 + k, then

   
∂u (C − A) ∂u (D − E)
= , =
∂x x,y 2h ∂y x,y 2k
 2   2

∂ u C − 2B + A ∂ u D − 2B + E
2
= , =
∂x x,y h2 ∂y 2 x,y k2

1.2 Solution of Laplace Equation in Two Dimensions

The Laplace equation in two dimensions is:

∂2u ∂2u
+ 2 =0 (8)
∂x2 ∂y
We solve this in a bounded region R, in a square or a rectangle under:
• (a) Prescribed values of u on the boundary (Dirichlet’s problem)
• (b) Derivative boundary conditions (Neumann problem)
• (c) Mixed boundary conditions (Cauchy problem)
The basic principle of the finite difference solution is to discretize as well as the region at the points uij , as
given in the figure below:
In Dirichlet’s problem, the values of uij on the boundary of the square are given, and our aim is to find
the values of uij at internal points, called points satisfying the Laplace Eq (8) and the boundary values.
Taking interval size h on both x and y side the Laplace Eq (8) may be discretized as

2
Santanu Manna MA 204N: Finite Difference Method

ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1

2
+ =0 (9)
h h2
which gives the recurrence relation for different i and j from Eq. (8) as

1
ui,j = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) (10)
4
Equation (10) is called the standard five-point formula which is shown in above figure. This may
also be called the mid-point formula.
Since Laplace equation is invariant when the coordinate axes are rotated through 45°, the five-point
formula may also be written using diagonal positions as given below:

1
ui,j = (ui−1,j−1 + ui+1,j+1 + ui−1,j+1 + ui+1,j−1 ) (11)
4

3
Santanu Manna MA 204N: Finite Difference Method

The recurrence relations of Eqs. (10) or (11) may be solved by various methods such as:
• (a) Gauss–Jacobi method
• (b) Gauss–Seidel method
• (c) Successive over-relaxation method

• (d) Tridiagonal method

• (e) Liebmann’s method, and so on

Example : 1
Solve the Laplace equation
∂2u ∂2u
+ 2 =0
∂x2 ∂y
on the square with boundary values as prescribed.

Solution
Using the mid-point Eq. (11), we get:
1
u11 = (0 + u12 + 40 + u21 )
4
1
u21 = (0 + u22 + u11 + u31 )
4
1
u22 = (u12 + u32 + u21 + u23 )
4
1
u32 = (150 + u22 + u31 + u33 )
4
1
u13 = (80 + u23 + u12 + 100)
4
1
u23 = (100 + u22 + u13 + u33 )
4
1
u33 = (100 + 200 + u23 + u32 )
4
1
u31 = (0 + u32 + 100 + u21 )
4
1
u12 = (60 + u22 + u11 + u13 ) (12)
4
We have nine unknowns uij and nine linear equations which may be solved by:
• (a) Gauss elimination method

• (b) Gauss–Jacobi or Gauss iteration method using computer

4
Santanu Manna MA 204N: Finite Difference Method

2 Solution of Poisson equation:

Consider eq. (2) defined over a rectangular domain

D = [0 ≤ x ≤ a] × [0 ≤ y ≤ b]

with Dirichlet condition u = uc at x = 0, y = 0, y = b and mixed condition ∂u

∂x = αu + β at x = a.
Let the domain be subdivided into square mesh with side h, i.e., ∆x = ∆y = h so that M h = a and
N h = b. If xi = ih, i = 0(1)M and yj = jh, j = 0(1)N , then the mesh point (i, j) corresponds to the point
(xi , yj ) in the domain.
Discretizing eq. (2) at mesh point (i, j), we get:
ui−1,j − 2ui,j + ui+1,j ui,j−1 − 2ui,j + ui,j+1
2
+ = f (xi , yj )
(∆x) (∆y)2
or

ui−1,j + ui,j−1 − 4ui,j + ui+1,j + ui,j+1 = h2 fi,j , i = 1(1)M, j = 1(1)N − 1 (2.1)

The set of simultaneous equations (2.1) can be solved by Gauss-Seidel method in the following scheme:

(n+1) 1 h (n+1) (n+1) (n) (n)

i
ui,j = ui−1,j + ui,j−1 + ui+1,j + ui,j+1 − h2 fi,j (2.2)
4
Start the computations taking the initial values of u, i.e., ui,j = 0 varying i as i = 1(1)M for fixed j,
j = 1(1)N − 1. In the formula (2.2), the most recent values of u’s are used as soon as they are computed.
The boundary condition prescribed at x = a can be approximated as:
uM +1,j − uM −1,j
= αuM,j + β
2h
or
uM +1,j = uM −1,j + 2h(αuM,j + β) (2.3)
where (M + 1, j) is a fictitious point.

Example: 2 A p.d.e.
∂2u ∂2u
+ 2 = 0.5
∂x2 ∂y
is defined over a rectangular domain

0 ≤ x ≤ 0.6, 0 ≤ y ≤ 0.6

with boundary conditions: u = 1 at three sides x = 0, y = 0, y = 0.6; and ∂u

∂x = u on x = 0.6. Solve the
equation by taking h = 0.2 and approximating the derivative boundary condition by CD.

5
Santanu Manna MA 204N: Finite Difference Method

Solution:
Given ∆x = ∆y = h = 0.2, and f (x, y) = 0.5.
We have to compute ui,j , i = 1(1)3, j = 1(1)2.
Approximating the derivative boundary condition by CD at mesh point (3,1):
u4,1 − u2,1
= u3,1
0.4
or
u4,1 = u2,1 + 0.4u3,1 (1)
Similarly at mesh point (3,2):
u4,2 = u2,2 + 0.4u3,2 (2)

Discretizing the p.d.e. at various mesh points (i, j), i = 1(1)3, j = 1(1)2:

At (1,1):
4u1,1 = u0,1 + u1,0 + u2,1 + u1,2 − h2 f1,1
we have h2 fi,j = 0.02 ∀i, j, therefore

1 n
un+1 u1,2 + un2,1 + 1.98

1,1 = (3)
4
At (2,1):
4u2,1 = u1,1 + u2,0 + u3,1 + u2,2 − h2 f2,1
or
1  n+1
un+1 u1,1 + un3,1 + un2,2 + 0.98

2,1 = (4)
4
At (3,1):
4u3,1 = u2,1 + u3,0 + u4,1 + u3,2 − h2 f3,1
or
5  n
un+1 2u2,1 + un3,2 + 0.98

3,1 = (5)
18
At (1,2):
4u1,2 = u0,2 + u1,1 + u2,2 + u1,3 − h2 f1,2
or
1  n+1
un+1 u1,1 + un2,2 + 1.98

1,2 = (6)
4
At (2,2):
4u2,2 = u1,2 + u2,1 + u3,2 + u2,3 − h2 f2,2
or
1  n+1
un+1 u1,2 + un+1 n

2,2 = 2,1 + u3,2 + 0.98 (7)
4
At (3,2):
4u3,2 = u2,2 + u3,1 + u4,2 + u3,3 − h2 f3,2
or
5  n
un+1 2u2,2 + un+1

3,2 = 3,1 + 0.98 (8)
18
Solving the above equations from (3) to (8) by substituting the most recent values of u. The convergence
is achieved after 14 iterations. The iterated values are shown in the table after skipping the alternate iterative
values:

6
Santanu Manna MA 204N: Finite Difference Method

j\i 0 1 2 3
0 1 1 1 1
1 0.495 0.3688 0.4771 1
0.8642 0.8656 1.0141
0.9756 1.0126 1.1486
1.0046 1.053 1.1830
1.0120 1.0599 1.1878
1.0137 1.0622 1.1939
1.0142 1.0628 1.1944
2 0.6188 0.4919 0.6788 1
0.9070 0.5230 0.0667
0.9870 1.0274 1.1620
1.0089 1.0544 1.1866
1.0126 1.0608 1.1922
1.0139 1.0624 1.1941
1.0142 1.0628 1.1944

3 Successive Over Relaxation (SOR) Method:

The iterative methods for solving linear system of equations, namely Gauss-Jacobi and Gauss-Seidel can be
modified to increase their rate of convergence. These improved methods are known as SOR methods. They
are also called Liebmann’s method.

The scheme for solution of Poisson equation by Gauss-Seidel method is discussed in previous lecture and
is given as:
(n+1) 1  (n+1) (n+1) (n) (n)

ui,j = ui−1,j + ui,j−1 + ui+1,j + ui,j+1 − h2 fi,j (3.1)
4
where i = 1(1)M for fixed j, j = 1(1)N − 1.

(n+1)
The Successive Over Relaxation scheme for Gauss-Seidel method for computing the value of ui,j from
(3.1) is expressed as:
(n+1) (n)
ui,j = ui,j + ωRi,j
where,
1  (n+1) (n+1) (n) (n)

(n)
Ri,j = ui−1,j + ui,j−1 + ui+1,j + ui,j+1 − h2 fi,j − ui,j
4
(value from Gauss-Seidel) - (value at previous iteration)
or,
(n+1) (n) ω  (n+1) (n+1) (n) (n)

ui,j = (1 − ω)ui,j + ui−1,j + ui,j−1 + ui+1,j + ui,j+1 − h2 fi,j (3.2)
4

The optimum value of ω for fastest convergence of the SOR scheme (3.2) is given by:
p
ωopt = 2 − 2 1 − 4t2 (3.3)
where π π
t = cos + cos
M N
For large M and N , (3.3) may be approximated as:
√
 
1 1
ωopt = 2 − 2π + (3.4)
M2 N2

7
Santanu Manna MA 204N: Finite Difference Method

If M = N then  π
ωopt = 2 1 − (3.5)
M
For ω < 1, the scheme is known as under-relaxed and for 1 < ω < 2, scheme is over-relaxed.

Example: A p.d.e.
∂2u ∂2u
+ 2 = 0.5
∂x2 ∂y
is defined over a rectangular domain [0 ≤ x ≤ 0.6] × [0 ≤ y ≤ 0.6] with boundary conditions, u = 1 at three
sides x = 0, y = 0, y = 0.6; and ∂u
∂x = u on x = 0.6. Solve the equation by taking h = 0.2 and using SOR
method corresponding to Gauss-Seidel method by taking ω = 1.2.

Solution:
∆x = ∆y = h = 0.2,
and f (x, y) = 0.5.
Approximating the boundary condition by CD, we have

u4,j = u2,j + 0.4u3,j (1)

The SOR scheme at different points is defined as:

ω n+1
un+1 n
ui−1,j + un+1 n n 2


i,j = (1 − ω)ui,j + i,j−1 + ui+1,j + ui,j+1 − h fi,j (2)
4
At (1,1):
ω n+1
un+1 n
u0,1 + un+1 n n 2


1,1 = (1 − ω)u1,1 + 1,0 + u2,1 + u1,2 − h f1,1
4
or
un+1 n n n
1,1 = −0.2u1,1 + 0.3[u2,1 + u1,2 + 1.98] (3)
At (2,1):
ω n+1
un+1 n
+ un+1 n n 2


2,1 = (1 − ω)u2,1 + u 2,0 + u3,1 + u2,2 − h f2,1
4 1,1
or
un+1 n n+1 n n
2,1 = −0.2u2,1 + 0.3[u1,1 + u3,1 + u2,2 + 0.98] (4)
(3,1):
ω n+1
un+1 n
u2,1 + un3,0 + un4,1 + un3,2 − h2 f3,1


3,1 = (1 − ω)u3,1 +
4
or
1
un+1 n
2un+1 n


3,1 = −0.2u3,1 + 2,1 + u3,2 + 0.98 (5)
3

(1,2):
ω n
un+1 n
u0,2 + un+1 n n 2


1,2 = (1 − ω)u1,2 + 1,1 + u2,2 + u1,3 − h f1,2
4
or
un+1 n n+1 n


1,2 = −0.2u1,2 + 0.3 u1,1 + u2,2 + 1.98 (6)

(2,2):
ω n+1
un+1 n
u1,2 + un+1 n n 2


2,2 = (1 − ω)u2,2 + 2,1 + u3,2 + u2,3 − h f2,2
4
or
un+1 n n+1 n+1 n


2,2 = −0.2u2,2 + 0.3 u1,2 + u2,1 + u3,2 + 0.98 (7)

8
Santanu Manna MA 204N: Finite Difference Method

(3,2):
ω n+1
un+1 n
+ un+1 n n 2


3,2 = (1 − ω)u3,2 + u 3,1 + u4,2 + u3,3 − h f3,2
4 2,2
or
1
un+1 n
2un+1 n+1


3,2 = −0.2u3,2 + 2,2 + u3,1 + 0.98 (8)
3

Solving the above equations from (3) to (8) by substituting the most recent values of u. The convergence is
achieved after 7 iterations which are half those in the Gauss-Seidel method. The iterative values are shown
in the table:

j\i 0 1 2
1 0.594 0.4722 0.6415
0.8485 0.8467 1.0913
0.9466 1.0296 1.1767
1.0114 1.0591 1.1937
1.0135 1.0632 1.1951
1.0145 1.0632 1.1947
1.0145 1.0630 1.1946
2 0.7722 0.6673 0.9854
0.8943 0.9785 1.1457
0.9927 1.0487 1.1889
1.0135 1.0627 1.1947
1.0142 1.0632 1.1948
1.0145 1.0632 1.1947
1.0144 1.0630 1.1946