American Scientific Research Journal for Engineering, Technology,  and Sciences  (ASRJETS) 
ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 

© Global Society of Scientific Research and Researchers  

http://asrjetsjournal.org/  
 

Solution of Laplace`s Equation Using Numerical Methods 

Eltigani Ismail Hassan* 

Department of Mathematics and Statistics, College of Science, Al Imam Mohammad Ibn Saud Islamic 

University, Riyadh 11642, KSA. 

Department of Mathematics, College of applied science and industrial, University of Bahri 

Email: ttfsr555@gmail.com 

 

 

Abstract 

This paper is simple review of the solution of Laplace`s equation in rectangular coordinates system, cylindrical 

polar coordinates system and spherical polar coordinate system. It also covers numerical method in the solution 

of Laplace`s equation. 

Keywords: Laplace equation; Numerical method; Laplacian difference; Uniform grid. 

1. Introduction 

During the last two centuries several methods have been advanced for solving partial differential equations. 

Among these we consider only two techniques known as the method of separation of variables and Laplace 

transformation , the method of separation of variables is perhaps the oldest systematic method  for solving 

partial differential equations. It has been considerably refined and generalized in meantime and remains a 

method of great importance to day. In this study we will review how the method of separation of variables. In 

solving this problem it is eventually necessary to consider the equations of whether an essentially arbitrary 

function can be expressed as an infinite series of sine and cosine functions. Series of this kind are called Fourier 

series, then we will come across to in some typical problems as the Laplace`s equation solving it by separation 

of variables and Fourier series .By the other hand we will explain in following: 

 

 

 

 

------------------------------------------------------------------------ 

*Corresponding author. 

178 

http://asrjetsjournal.org/


American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

In first we define rectangular coordinate system 

∇2𝑈 =
𝜕2𝑢
𝜕𝑥2

+
𝜕2𝑢
𝜕𝑦2

= 0 

In second we transform the rectangular coordinates system to cylindrical coordinates system is which comes as: 

𝑈(𝑥,𝑦. 𝑧) = (𝜌𝜌,𝜃𝜃, 𝑧) 

and its solution given by: 

𝜕2𝑣
𝜕𝜌𝜌2

+
1
ρ
𝜕𝑣
∂ρ

+
1
ρ2
𝜕2𝑣
∂θ2

+
𝜕2𝑣
∂z2

= 0 

In third we also transform the rectangular coordinates system to spherical polar coordinates system which comes 

as: 

𝑈 = (𝑥,𝑦, 𝑧) = (𝑟,𝜃𝜃,𝜙𝜙) 

∇2𝑈 =
𝜕2u
∂r2

+
2
r
𝜕u
∂r

+
1

r2 sin θ
∂
∂θ

(sin𝜃𝜃
𝜕𝑢
𝜕𝜃𝜃

) +
1
𝑟2

1
sin2θ

∂2u
∂θ2

= 0 

In     fourth we explain computerized computation method to solve general solution of Laplace`s equation by 

finite different equation. 

2.  Two dimensional heat flow 

Consider the flow of heat in a metal plate of uniform thickness 𝛼𝛼(𝑐𝑚), density (𝜌𝜌 𝑔𝑟⁄ 𝑐𝑚3⁄ ) ,specific heat 

𝑆(𝑐𝑎𝑙. 𝑔𝑟⁄ . deg)and thermal conductivity 𝐾(𝑐𝑎𝑙 𝑐𝑚⁄ . deg) . Let 𝑋 ∘ 𝑌plane be taken in face the point [1]. If 

temperature at any point is independent of the 𝑍coordinate and depends only on 𝑋,𝑌 and time 𝑡. Then the flow 

is said to be two dimensional .In this case, the heat flow is in the 𝑋𝑌-plane only and is zero along the normal to 

the 𝑋𝑌-plane. Consider a rectangular element𝐴𝐵𝐶𝐷 of the plane with sides as shown in figure. By 𝐴on the  

amount of heat entering the element, from the side 

AB = −kα∂y �∂u
∂x
�
x
 . 

The quantity of heat flowing out through the side 𝐶𝐷  per second = −𝑘𝛼𝛼𝜕𝑥 �𝜕𝑢
𝜕𝑦
�
𝑦+𝜕𝑦

         (1) 

and the quantity of heat flowing out through the side 𝐵𝐶  per second = −𝑘𝛼𝛼𝜕𝑦 �𝜕𝑢
𝜕𝑥
�
𝑥+𝜕𝑥

   (2) 

Hence the total gain of heat by rectangular element per second 

179 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

= −𝑘𝛼𝛼𝜕𝑥 �𝜕𝑢
𝜕𝑦
�
𝑦
− 𝑘𝛼𝛼𝜕𝑦 �𝜕𝑢

𝜕𝑥
�
𝑥

+ 𝑘𝛼𝛼𝜕𝑥 �𝜕𝑢
𝜕𝑦
�
𝑦+𝜕𝑦

+ 𝑘𝛼𝛼𝜕𝑦 �𝜕𝑢
𝜕𝑥
�
𝑥+𝜕𝑥

= 𝑘𝛼𝛼𝜕𝑥 �
�𝜕𝑢𝜕𝑥�𝑥+𝜕𝑥

−�𝜕𝑢𝜕𝑥�𝑥
𝜕𝑥

+
�𝜕𝑢𝜕𝑦�𝑦+𝜕𝑦

−�𝜕𝑢𝜕𝑦�𝑦
𝜕𝑦

� 

Also the rate of gain of heat by the element = 𝜌𝜌𝜕𝑥𝜕𝑦𝛼𝛼𝑠 𝜕𝑢
𝜕𝑡

 

Thus equating (1) and (2)                                                                        = 𝑘𝛼𝛼𝜕𝑥𝜕𝑦 �
�𝜕𝑢𝜕𝑥�𝑥+𝜕𝑥

−�𝜕𝑢𝜕𝑦�𝑥
𝜕𝑥

+

�𝜕𝑢𝜕𝑥�𝑦+𝜕𝑦
−�𝜕𝑢𝜕𝑦�𝑦

𝜕𝑦
� =  𝜌𝜌𝜕𝑥𝜕𝑦𝛼𝛼𝑠 𝜕𝑢

𝜕𝑡
(3) 

Dividing both sides by 𝛼𝛼𝜕𝑥𝜕𝑦 and taking limit as → 0 , we get: 

𝑘 �
𝜕2𝑢
𝜕𝑥2

+
𝜕2𝑢
𝜕𝑦2

� = 𝜌𝜌𝑠
𝜕𝑦
𝜕𝑡

 

𝜕𝑢
𝜕𝑡

= 𝐶2 �𝜕
2𝑢
𝜕𝑥2

+ 𝜕2𝑢
𝜕𝑦2

�                                                                                                    (4) 

where 𝐶2 = 𝑘/𝜌𝜌𝑠 is the diffusivity. Hence the equation (4) gives the temperature distribution of plane in the 

transit state [1]. 

 

3. The diffusion equation in two-dimensions 

When cylindrical co-ordinates 𝜌𝜌,𝜃𝜃, 𝑧 are used see [2],we have 

𝑥 = 𝜌𝜌 cos 𝜃𝜃 , y = ρ sin𝜃𝜃  , z = zand the equation of the conduction of heat comes                                 

   𝜕
2𝑣

𝜕𝜌2
+ 1

ρ
𝜕𝑣
∂ρ

+ 1
ρ2

𝜕2𝑣
∂θ2

+ 𝜕2𝑣
∂z2

= 0                                                                                            (5) 

Let us first consider solutions which are independent of 𝑧 

𝐶(𝑥 + 𝜕𝑥,𝑦) +  𝜕𝑦) 

𝐵(𝑥 + 𝜕𝑥,𝑦 ) 

𝐷(𝑥, 𝑦)
  

𝐴(𝑥,𝑦) 

Figure1  

180 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

𝜕2𝑣
𝜕𝜌𝜌2

+
1
ρ
𝜕𝑣
∂ρ

+
1
ρ2
𝜕2𝑣
∂θ2

 = 0                                                                                                 (6) 

the equation for 𝑅  is 𝜕
2𝑅
𝜕𝜌2

+ 1
ρ
𝜕𝑅
∂ρ

+ 1
ρ2

𝜕2𝑅
∂θ2

 = 0                                                              (7) 

By separation of variables 

𝑅∖∖𝜓𝜓 +
1
𝜌𝜌
𝑅∖𝜓𝜓 +

1
𝜌𝜌2
𝑅𝜓𝜓∖∖ = 0                                                                                          (8) 

𝑅𝜓𝜓
𝑅∖∖

𝑅
+

1
𝜌𝜌
𝑅∖

𝑅
+
𝜓𝜓∖∖

𝜓𝜓
= −𝜆                                                                                                (9) 

𝜌𝜌2
𝑅∖∖

𝑅
+ ρ

𝑅∖

𝑅
+ λ2ρ2 −

𝜓𝜓∖∖

𝜓𝜓
= 𝜇𝜇2                                                                                    (10) 

from (10),we 𝜇𝜇2𝜓𝜓 = 0                                                                                                           (11) 

𝜌𝜌2𝑅∖∖ + 𝜌𝜌𝑅∖ + (λ2ρ2 − 𝜇𝜇2)𝑅 = 0                                                                                     (12) 

The solution of (11) and (12) given respectively by : 

𝜓𝜓(𝜃𝜃) = 𝐴1 cos𝜇𝜇𝜃𝜃 + 𝐴2 sin 𝜇𝜇𝜃𝜃                                                                                              (13) 

𝑅(𝜌𝜌) =  𝐵1𝐽𝜇(𝜆𝜌𝜌) + 𝐵2𝑌𝜇(𝜆𝜌𝜌)                                                                                              (14) 

The general solution  is 

𝑉(𝜌𝜌,𝜃𝜃) = [𝐴1 cos 𝜇𝜇𝜃𝜃 + 𝐴2 sin 𝜇𝜇𝜃𝜃]�𝐵1𝐽𝜇(𝜆𝜌𝜌) + 𝐵2𝑌𝜇(𝜆𝜌𝜌)�                                            (15) 

Since 𝑉 is bounded 𝜌𝜌 = 0   then 𝐵2 = 0       gives 

𝑉(𝜌𝜌,𝜃𝜃) = [𝐴1 cos 𝜇𝜇𝜃𝜃 + 𝐴2 sin 𝜇𝜇𝜃𝜃]�𝐵1𝐽𝜇(𝜆𝜌𝜌)�                      (16) 

An initial distribution of concentration, let 𝑉 = 𝑓(𝜌𝜌,𝜃𝜃) 

when  𝑡 = 0 ,we may try to satisfy the condition ,when 𝜇𝜇 = 𝑛 

and V must have period 2𝜋𝜋  in the variable 𝜃𝜃. 

(17) 𝑉(𝜌𝜌,𝜃𝜃) = [𝐴1 cos 𝜇𝜇𝜃𝜃 + 𝐴2 sin 𝜇𝜇𝜃𝜃]�𝐵1𝐽𝜇(𝜆𝜌𝜌)� 

Now by substituting boundary condition  𝑉(1, 𝜃𝜃) =0 

181 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

 𝑉(𝜌𝜌,𝜃𝜃) = [𝐴1 cos𝑛𝜃𝜃 + 𝐴2 sin𝑛𝜃𝜃][𝐵𝐽𝑛(𝜆𝜌𝜌)]                                                      (18) 

Now equation (17)  can be expressed us 

𝑉(𝜌𝜌,𝜃𝜃) = [𝐴 cos𝑛𝜃𝜃 + 𝐵 sin𝑛𝜃𝜃][𝐽𝑛(𝜆𝑚𝑛𝜌𝜌)]                                                                     (19) 

By using super position 

𝑉(𝜌𝜌,𝜃𝜃) = ∑ ∑ [Amncos𝑛𝜃𝜃 + Bmn sin𝑛𝜃𝜃]𝐽𝑛(𝜆𝑚𝑛𝜌𝜌)∞
m=0

∞
𝑛=0                                   (20) 

We let 𝑉 = 𝑓(𝜌𝜌,𝜃𝜃)  gives 

𝑓(𝜌𝜌,𝜃𝜃) = ∑ 𝐷𝑛∞
n=0 cos𝑛𝜃𝜃 + 𝐸𝑛 sin𝑛𝜃𝜃(21)Formula (21) represent Fourier Series of𝑓(𝜌𝜌,𝜃𝜃) with period 2𝜋𝜋 

𝐷𝑛 = 1
𝜋 ∫ 𝑓(𝜌𝜌,𝜃𝜃) cos𝑛𝜃𝜃 𝑑𝜃𝜃2𝜋

0 and𝐸𝑛 = 1
𝜋 ∫ 𝑓(𝜌𝜌,𝜃𝜃) sin𝑛𝜃𝜃 𝑑𝜃𝜃2𝜋

0    for 𝑛 = 1,2, … 

We have 𝐴𝑚𝑛 = 2
𝐽𝑛+1
2 (𝜆𝑚𝑛)∫ 𝜌𝜌𝐽𝑛(𝜆𝑚𝑛𝜌𝜌)𝐷𝑛𝑑𝜌𝜌

1
0  

And    𝐵𝑚𝑛 = 2
𝐽𝑛+1
2 (𝜆𝑚𝑛)∫ 𝜌𝜌𝐽𝑛(𝜆𝑚𝑛𝜌𝜌)𝐸𝑛𝑑𝜌𝜌

1
0  

4. The elementary solution 

If we solve 

∇2𝑈 = 0 in spherical coordinates (𝑟,𝜃𝜃,𝜙𝜙) 

∇2𝑈 =
1

ℎ1ℎ2ℎ3
[
𝜕
𝜕𝑞1

[
ℎ2ℎ3
ℎ1

𝜕𝑢
𝜕𝑞1

] +
𝜕
𝜕𝑞2

[
ℎ1ℎ3
ℎ2

𝜕𝑢
𝜕𝑞2

] +
𝜕
𝜕𝑞3

[
ℎ2ℎ1
ℎ3

𝜕𝑢
𝜕𝑞3

] 

whereℎ1 = 1, ℎ2 = 𝑟   and ℎ3 = 𝑟 sin𝜃𝜃 and   𝑞1 = 𝑟, 𝑞2 = 𝜃𝜃 ,  and 𝑞3 = 𝜙𝜙.If 𝑈 is independent of𝜙𝜙 with 

boundary condition 𝑈(1,𝜃𝜃) = 𝑓(𝜃𝜃) 

the solution of ∇2𝑈 = 0 in spherical coordinates which is independent of 𝜙𝜙 expressed as: 

𝑟2
𝜕2𝑢
𝜕𝑟2

+ 2𝑟
𝜕𝑢
𝜕𝑟

+
1

sin 𝜃𝜃
𝜕
𝜕𝜃𝜃

�sin𝜃𝜃
𝜕𝑢
𝜕𝜃𝜃
� = 0                                                      ( 22) 

Take 𝑢(𝑟,𝜃𝜃) = 𝑅(𝑟)𝜓𝜓(𝜃𝜃)(23) 

We get 

𝑟2𝑅\\𝜓𝜓 + 2𝑟𝑅\𝜓𝜓 +
𝑅

sin𝜃𝜃
𝜕
𝜕𝜃𝜃

�sin 𝜃𝜃 𝜓𝜓\� = 0 

182 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

Dividing the above equation by 𝑅𝜓𝜓 ,then we get 

𝑟2
𝑅\\

𝑅
+ 2𝑟

𝑅\

𝑅
= −

1
𝜓𝜓 sin𝜃𝜃

𝜕
𝜕𝜃𝜃

�sin 𝜃𝜃 𝜓𝜓\� = 𝜆2 

𝑟2𝑅\\ + 2𝑟𝑅\ + 𝜆2𝑅 = 0                                                                                          (24) 

And 

𝜕
𝜕𝜃𝜃

�sin𝜃𝜃 𝜓𝜓\� − 𝜆2ψ sin θ = 0 

sin𝜃𝜃
𝜕2𝜓𝜓
𝜕𝜃𝜃2

+ cos 𝜃𝜃
𝜕𝜓𝜓
𝜕𝜃𝜃

− 𝜆2𝜓𝜓 sin𝜃𝜃 = 0                                                                  (25) 

the solution of (24) is given by 

𝑅(𝑟) = 𝐴𝑟𝑛 + 𝐵
𝑟𝑛+1

                                                                                                     (26) 

where𝜆2 = −𝑛(𝑛 + 1)                                                                                               (27) 

from (27) and equation(25) 

sin𝜃𝜃
𝜕2𝜓𝜓
𝜕𝜃𝜃2

+ cos 𝜃𝜃
𝜕𝜓𝜓
𝜕𝜃𝜃

+ 𝑛(𝑛 + 1) sin 𝜃𝜃 𝜓𝜓 = 0                                                   (28) 

Equation (28) is  a Legendre ODE, with general solution 

𝜓𝜓(𝜃𝜃) = 𝐶𝑃𝑛(cos 𝜃𝜃) + 𝐷𝑄𝑛(cos 𝜃𝜃)                                                                       (29) 

𝑢(𝑟,𝜃𝜃) = �𝐴𝑟𝑛 +
𝐵

𝑟𝑛+1
� [𝐶𝑃𝑛(cos 𝜃𝜃) + 𝐷𝑄𝑛(cos 𝜃𝜃)]                                        (30) 

If 𝑢(𝑟,𝜃𝜃) represent temperature on sphere radius with center at origin, there the heat must be bounded. When 

𝜃𝜃 = 0  or𝜋𝜋 (along 𝑍 - axis in spherical coordinates) when 

𝑄𝑛(1) → ∞ this provided that 𝐷 = 0                                                                      ( 31) 

To avoid infinitely temperature at center of sphere 𝑟 = 0  we take 

𝐵 = 0                                                                                                                                (32) 

Now from (31) and (32) the solution is given by (30) 

𝑢(𝑟,𝜃𝜃) = 𝐸𝑟𝑛𝑃𝑛(cos 𝜃𝜃)                                                                                  (33) 

183 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

By the supper position (33) can be written as 

𝑢(𝑟,𝜃𝜃) = �𝐸𝑛𝑟𝑛𝑃𝑛(cos𝜃𝜃)                                                                                    (34)
∞

𝑛=0

 

By applying the boundary condition  𝑢(1,𝜃𝜃) = 𝑓(𝜃𝜃) into (34). 

𝑓(𝜃𝜃) = � EnPn(cos θ)                                                                                              (35)
∞

𝑛=0

 

Let 𝜇𝜇 = cos 𝜃𝜃,then𝜃𝜃 = 𝑐𝑜𝑠−1𝜇𝜇 

𝑓(𝑐𝑜𝑠−1𝜇𝜇) = � EnPn(µ)                                                                                          (36)
∞

𝑛=0

 

Now 

𝐸𝑛 =
(2𝑛 + 1)

2
�𝑓(𝑐𝑜𝑠−1𝜇𝜇)𝑃𝑛(𝜇𝜇)𝑑𝜇𝜇                                                                  (37)
1

−1

 

From (34), we have 

𝑢(𝑟,𝜃𝜃) = ��
(2𝑛 + 1)

2
�𝑓(𝜃𝜃)𝑃𝑛(cos 𝜃𝜃) sin𝜃𝜃𝑑𝜃𝜃
𝜋

0

� [𝑟𝑛𝑃𝑛(cos 𝜃𝜃)]
∞

𝑛=0

. 

5. The Laplacian Difference Equation 

Central differences based on the grid and scheme used for the finite - difference solution in two independent 

variables such as the Laplace equation [3]. 

𝜕2𝑇
𝜕𝑥2

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

∇𝑥2
                                                                                       (38) 

And 

𝜕2𝑇
𝜕𝑦2

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

∇𝑦2
                                                                                        (39) 

respectively   which have errors of O[∇(𝑥)2]  and O[∇(𝑦)2] 

Substituting these expressions into equation into equation 

184 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

𝜕2𝑇
𝜕𝑥2

+
𝜕2𝑇
  𝜕𝑦2

 = 0 

Gives 

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

∇𝑥2
+
𝑇𝑖,𝑗+1 − 2𝑇𝑖,𝑗 + 𝑇𝑖,𝑗−1

∇𝑦2
= 0                                                      (40) 

For the square grid, ∇𝑥 = ∇𝑦 and by collection of terms ,the equation becomes 

𝑇𝑖+1,𝑗 + Ti−1,j + Ti,j+1 + Ti,j−1 − 4Ti,j = 0                                                                       (41) 

This relationship, which holds for all interior points on the plate, is referred to as the Laplacian different 

equation. 

 

 

 

 

 

A heated plate where boundary temperature are hold at const balance 

ant levels. This called adirichlet  boundary condition. Where the edges are hold at constant temperature for the 

case illustrated in figure 2 balance figurec1 is ,according to equation (41) 

𝑇21 + T01 + T12 + T10 − 4T11 = 0                                                                               (42) 

 and  T10 = 0However, T01 = 75 

therefore equation  𝜕
2𝑇

𝜕𝑥2
+ 𝜕2𝑇

  𝜕𝑦2
 = 0, can be expressed as −4𝑇11 + 𝑇12 + 𝑇21 = −75 

Similar equation can be developed for the other interior points. the result is the following set of nine 

simultaneous equations with nine unknowns. 

4𝑇11    − 𝑇21    − 𝑇12                                                                                = 75 

−𝑇114𝑇21−𝑇11−𝑇22                                                                   = 0 

−𝑇214𝑇31−𝑇32                                                              = 50 

(1,3)              (2,3) 

(1,2)              (2,2)             (3,2) 

(1,1)              (2,1)            (3,1) 

100 

  

75C

  

50C ْ◌ 

0C ْ◌ 
Figure 2  

 

185 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

−𝑇11  4𝑇12−𝑇22−𝑇13                                        = 75 

−𝑇21−𝑇124𝑇22−𝑇32−𝑇23                            = 0 

−𝑇31−𝑇224𝑇32𝑇33                        = 50 

−𝑇124𝑇13−𝑇23 = 175 

−T22−T134𝑇23−T33      = 100 

−T32−𝑇234𝑇33  = 150 

5.1 Example 

Temperature of heated plate with fixed boundary conditions problem statement [4]. Use Lobmann's   method 

(Gauss-Seidel) to solve for temperature of the heated in figure (4.2).Employ over relaxation with a value of 1.5 

for the   weighting factor and iterate  𝜀𝜀𝑎 = 1%. 

solution 

From equation 𝑇𝑖,𝑗 =
𝑇𝑖+1,𝑗+𝑇𝑖−1,𝑗+𝑇𝑖,𝑗+1+𝑇𝑖,𝑗−1

4
 

at𝑖 = 1, 𝑗 = 1 

is𝑇11 = 0+75+0+0   
4

  =18.75  and applying relaxation yield 

𝑇11 = 1.5(18.75) + (1 − 1.5)0 = 28.125 

𝑇21 =
0 + 28.125 + 0 + 0

4
= 7.03125 

𝑇21 = 1.5(7.03125) + (1 − 1.5)0 = 10.54688 

𝑇31 =
50 + 10.54688 + 0 + 0

4
= 15.13672 

𝑇31 = 1.5(15.13672) + (1 − 1.5)0 = 22.70508 

The computation is repeated for other rows to give 

𝑇12 = 38.67188          𝑇22 = 18.45703      𝑇32 = 34.18579 

𝑇13 = 80.12696          𝑇23 = 74.46900             𝑇33 = 96.99554 

186 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

Because all the 𝑇𝑖,𝑗`𝑠 are initially zero, all 𝜀𝜀𝑎  for the first iteration will be 100%. 

For the second iteration 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      𝑇33 = 67.68736 

The error for𝑇11  can be estimated as equation −𝜕𝑞
𝜕𝑥
− 𝜕𝑞

𝜕𝑞
= 0 

�(𝜀𝜀𝑎)𝑖,𝑗� = �
32 − 1953 − 28 − 12500

32 − 51953
� 100% = 13.5% 

Because this value is  the stopping criterion of 1% the computation is continued. 

The ninth iteration gives the result: 

𝑇11 = 43.0061           𝑇21 = 33.29755        𝑇31 = 33.88506 

𝑇12 = 63.21152      𝑇22 = 56.11238          𝑇32 = 52.33999 

𝑇13 = 78.58718        𝑇23 = 76.06402        𝑇33 = 69.71050 

where the maximum error is 0.71%. 

6. Explicit Solution for uniform grid increments 

Consider  a rectangular domain where the increments in both 𝑥 and 𝑦 are uniform [5]. The appropriate equation 

to use for an explicit solution to the Laplace equation is  

𝑇𝑖𝑗 =
𝑇𝑖−1,𝑗+𝑇𝑖+1,𝑗+𝛽2�𝑇𝑖,𝑗+1+𝑇𝑖,𝑗+1�

2(1+𝛽2)
                                                                 (43)   

Where 𝛽𝛽 = Δ𝑥
Δ𝑦

 

The solution  will start by loading  the boundary conditions, and then calculating  the values of 𝑇𝑖𝑗  in the interior 

points of domain. While we are initially temped to calculate𝑇𝑖𝑗   only once with equation (43), it should be 

mentioned that these values are only a fist approximation to the solution .We should, therefore , add additional 

index ,𝑘, representing the current iteration, to each solution value .The solution values will now  be referred to 

as 𝑇𝑖𝑗𝑘 , and equation (43) will  be modified to read:   

187 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

𝑇𝑖,𝑗𝑘+1 =
𝑇𝑘𝑖−1,𝑗+𝑇𝑘𝑖+1,𝑗+𝛽2�𝑇𝑘𝑖,𝑗−1+𝑇𝑘𝑖,𝑗+1�

2(1+𝛽2)
                                                                             (44) 

The iterative process should be repeated until convergence  is achieved in every interior point  of the domain , or 

until a maximum number of iterations , say 100 , have been performed . Convergence can be achieved ,for 

example ,if ,given a tolerance value 𝜀𝜀 maximum difference two consecutive iterations is less than the tolerance , 

,i.e., if 

max
𝑖,𝑗

�𝑇𝑖,𝑗𝑘+1 − 𝑇𝑘𝑖,𝑗� ≤ 𝜀𝜀. 

Consider , as an example ,a rectangular domain of length 𝐿 = 5𝑐𝑚 , and height 𝐻 = 3.5𝑐𝑚 , with increments 

∆𝑥 = 1𝑐𝑚, 𝑎𝑛𝑑 ∆𝑦 = 0.5𝑐𝑚 , as illustrated in the figure blow. 

There will be 𝑛 = 𝐿
∆𝑥�   sub-intervals in 𝑥 , and𝑚 = 𝐻

∆𝑦�  sub-intervals in  , with 

𝑥𝑖 = (𝑖 − 1)∆𝑥, 𝑓𝑜𝑟 𝑖 = 1,2, … ,𝑛 + 1, 

and 

𝑦𝑗 = (𝑗 − 1)∆𝑦 , 𝑓𝑜𝑟 𝑗 = 1,2, … ,𝑚 + 1 

The boundary conditions are given as follows: 𝑇𝑖𝑗 = 5 along the left and right sides of the domain , while the 

temperature are given by the function 𝑇𝑏(𝑥) = 5. 𝑥. (1 − 𝑥) for the top and bottom sides of the domain 

,respectively [5]. 

Solution is achieved by using function LaplaceExplicit.min Matlab : 

188 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

function [x,y,T]= LaplaceExplicit(n,m,Dx,Dy) 

echo off; 

numgrid(n,m); 

R = 5.0; 

T = R*ones(n+1,m+1); % All T(i,j) = 1 includes all boundary conditions 

x = [0:Dx:n*Dx];y=[0:Dy:m*Dy]; % x and y vectors 

for i = 1:n % Boundary conditions at j = m+1 and j = 1 

6 

T(i,m+1) = T(i,m+1)+ R*x(i)*(1-x(i)); 

T(i,1) = T(i,1) + R*x(i)*(x(i)-1); 

end; 

TN = T; % TN = new iteration for solution 

err = TN-T; 

% Parameters in the solution 

beta = Dx/Dy; 

denom = 2*(1+beta^2); 

% Iterative procedure 

epsilon = 1e-5; % tolerance for convergence 

imax = 1000; % maximum number of iterations allowed 

k = 1; % initial index value for iteration 

% Calculation loop 

while k<= imax 

189 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

for i = 2:n 

for j = 2:m 

TN(i,j)=(T(i-1,j)+T(i+1,j)+beta^2*(T(i,j-1)+T(i,j+1)))/denom; 

err(i,j) = abs(TN(i,j)-T(i,j)); 

end; 

end; 

T = TN; k = k + 1; 

errmax = max(max(err)); 

if errmax< epsilon 

[X,Y] = meshgrid(x,y); 

figure(2);contour(X,Y,T',20);xlabel('x');ylabel('y'); 

title('Laplace equation solution - Dirichlet boundary conditions 

- Explicit'); 

figure(3);surfc(X,Y,T');xlabel('x');ylabel('y');zlabel('T(x,y)'); 

title('Laplace equation solution - Dirichlet boundary conditions 

- Explicit'); 

fprintf('Convergence achieved after %i iterations.\n',k); 

fprintf('See the following figures:\n'); 

fprintf('==========================\n'); 

fprintf('Figure 1 - sketch of computational grid \n'); 

fprintf('Figure 2 - contour plot of temperature \n'); 

fprintf('Figure 3 - surface plot of temperature \n'); 

190 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

return 

end; 

end; 

fprintf('\n No convergence after %i iterations.',k); 

 

 

To activate the function for the case illustrated in the figure above we use: 

≫ [𝑋,𝑌,𝑇] = 𝐿𝑎𝑝𝑙𝑎𝑐𝑒𝐸𝑥𝑝𝑙𝑖𝑐𝑖𝑡(5,7,1,0.5) 

The solution is returned in the vectors 𝑥 and 𝑦 ,and in matrix 𝑇 .The function produces three plots: a sketch of 

the grid (similar to the figure above),the solution as a contours, and the solution as a surface .The last two 

figures are shown next: 

Laplace equation solution –Dirichletboundary  conditions-Explicit 

 

Laplace equation solution - Dirichletboundary conditions -  Explicit 

 

 

 

 

191 



American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2016) Volume 16, No  1, pp 178-192 

7. Conclusion and Recommendation  

This paper showed that Numerical methods are usually easier to use in the solution of Laplace's equation. We 

notice that our results agree with other works cited in article. It is aim to continue higher order for two 

dimensional.    

References 

[1] Partial Differential equations Of Mathematical Physics.1969. pp  120-126. 

[2] Paul W. Berg& James L. MC Gregor Elementary differential equation 1960, pp55. 

[3] Steven C. Chapra &Roymond P- Canale Numerical methods for Engineers 2002.pp 40. 

[4] Grewal M.A. ph. D. Higher Engineering Methematics 1990.pp 86-99. 

[5] Numerical Solution of Laplace Equation By Gilberto. Urroz,October 2004. pp 4-7.    

 

 

192