Numerical Solution of Two Dimensional Poison Equations by Using Fourth-

Order Compact Finite Difference Method

Kedir Aliyi Koroche

Department of Mathematics, College of Natural and Computational Sciences, Ambo University,
Ambo, Ethiopia
*Corresponding author

Email address: [email protected] ( Kedir A.)

Abstract

In this paper, the fourth-order compact finite difference scheme has been presented for solving
the two-dimensional posing equation. First, the given solution domain is discretized with
uniform and no uniform mesh size and then the partial derivative is replaced into functional
values at each grid point by using Taylor series expansion, and then it gives a system of algebraic
equations. Then they obtained a system of algebraic equations is solved by the Thomas method.
The stability and convergent analysis of the scheme are investigated. To validate the applicability
of the proposed method, one model example is considered and solved for different values of the
mesh sizes in both directions. Numerical results are presented in tables in terms of maximum
absolute error, L2 and L∞ norm. The numerical results presented in tables and graphs confirm
that the approximate solution is in good agreement with the exact solution.
Keywords: Elliptic equation, poison equation, compact difference scheme, convergence of
methods.
1. Introduction
A partial differential equation (PDE) is an equation containing a partial
derivative of the dependent variable [9]. These equations arise in almost all
areas of applied mathematics, physics, and some branches of engineering [5], for instance,
in fluid mechanics, elasticity, heat transfer, energy systems, environmental
flows, hydraulics, neutron diffusion in nuclear reactors, and structural
analysis [18]. This partial differential equation is classified into parabolic, hyperbolic, and
elliptic types of equations [5, 9, 11, 16]. Therefore the partial differential equations model sorts
of phenomena, display different behavior, and require different numerical techniques for their
solution [11] and simplest examples of the elliptic type of PDEs are Poisson’s equation and
Laplace equation [5]. These elliptic-type equations are generally associated with equilibrium or
steady-state problems [5].For instance Steady-state condition in a communications circuit and
electrical circuit.
The Steady-state condition in a communications circuit can be defined as a condition in which
some specified characteristic of a condition, such as a value, rate, periodicity, or amplitude,
exhibits only negligible change over an arbitrarily long period. Again the Steady-state
conditions, in an electrical circuit define as the condition that exists after all initial transients or
fluctuating conditions have damped out, and all currents, voltages, or fields remain essentially
constant or oscillate uniformly [12]. For example, the velocity potential for the steady flow of
incompressible non-viscous fluid satisfies Laplace’s equation and the electric potential
associated with a two-dimensional electron distribution of charge density is satisfies Poisson’s
equation [5]. The Poisson equation is a generalization of Laplace's equation . This equation is
named after the French mathematician geometer, and physicist Simon Denis Poisson [1].
Boundary conation of Elliptic types of PDEs equation arises in the study of steady-state or time-
independent solutions of heat equations. Because these solutions do not depend on time, initial
conditions are irrelevant and only boundary conditions are specified. Applications of Poisson
equation also include the static displacement u(x , y) of a stretched membrane fastened in space
along the boundary of a region; the electrostatic and gravitational potentials in certain force
fields; and, in fluid mechanics for an ideal fluid [17].
Poisson’s equation is also a very powerful tool for modeling the behavior of electrostatic
systems, but unfortunately may not only be solved analytically for very simplified models [1].
Because of this, these methods are based on advanced mathematical techniques [5].
Among Elliptic types of PDEs, Poisson type equations are the most practical and frequently
investigated [1]. In solving these types of partial differential equations, we are looking for a
function of more than a variable that satisfies the same relation between different partial
derivatives [11]. But poison equation may not only be solved analytically for very simplified
models. Consequently, numerical simulation must be utilized for that model problem due to them
has complex geometries behavior within their practical value [1]. Therefore, several numerical
methods are available to solve the Poison equation. The numerical methods are, in general,
simple but generate erroneous results [5].
In many application areas, such as aero acoustics and electromagnetic, the propagation of
acoustic and electromagnetic waves needs to be accurately simulated over very long periods and
far distances. Some numerical methods are not accurate for solving such types of the equation.
For instance, the finite difference method is used as the direct conversion of the partial
differential equation from continuous function and operator into their discretely sampled
counterpart. This converts the entire problem into a system of linear equations that may be
readily solved employing matrix inversion, Jacobi, Gauss-elimination, the successive over-
relaxation method [7]. The accuracy of such a method is therefore directly tied to the ability of a
finite grid to approximate a continuous system and errors may be arbitrarily reduced by simply
increasing the number of samples [6]. But the method has required high cost regarding storage
capacity in the computational domain. Nyaya Fekadie [16] solved the elliptic equation by using
the finite volume method. He used the finite volume method and the solution domain is
subdivided into a finite number of small control volumes by a grid that grid defines the
boundaries of the control volumes while the computational node lies at the center of the control
volume to solve the elliptic equation. Nodal points are used within these control volumes for
interpolating the field variable and usually, the single node at the center of the control volume is
used for each control volume. But the method gives better accuracy only for the small number of
a grid point and is difficult to compute the solution in a complex computational domain when
step length is very small. Genet M and Lemi A. [1] presented the solution of two-dimensional
poison equations using the finite difference method. This method is mathematically simple and
guarantees the necessary accuracy for a relatively small number of grid points. This confirms
that the method is not accurate for relatively a few grid points (i.e. for mesh size very large) and
is difficult to apply for high dimension geometric spaces. Hence this method does not always
converge to the exact solutions for coarser step lengths. Mohammad and Azim [20] also presented
the Numerical Solution of Poisson’s Equation Using a Combination of Logarithmic and
Multiquadric Radial Basis Function Networks. In multiquadric radial basis functions MQ-RBFs,
some parameters influence the accuracy of the solution. The solution diverges until the optimal
shape parameters are obtained. As compared to the exact solution, the approximate solution
needs further improvement.

Even though the accuracy of the aforementioned methods is promising, they require large
memory and long computational time. Besides, the methods are not suitable for higher-
dimensional and problems involving complex geometries. So, the treatment of the mesh size and
shape parameter in the applied method presents severe difficulties that have to be addressed to
ensure the accuracy of the solution for the poison equation, and efficiency of the method applied.
Therefore due to this end, the accumulation of errors is generated throughout solving the
poison equation. Thus still, the accuracy of the method needs attention; because the treatment of
the method used to solve the poison equation is not trivial distribution.

To reduce the accumulation of errors, the numerical algorithm must be highly accurate.
To accomplish this goal, high-order compact finite difference schemes have been developed to
solve PDEs types of the equation in a different application (see [8, 14, 19]). High-order finite
difference schemes can be classified into two main categories: explicit schemes and Pade-type or
compact schemes. Explicit schemes compute the numerical derivatives directly at each grid by
using large stencils, while compact schemes obtain all the numerical derivatives along a grid line
using smaller stencils and solving a linear system of equations. Experience has shown that
compact schemes are much more accurate than the corresponding explicit scheme of the same
order [14]. Therefore to this end, this paper aims to apply the fourth-order compact finite
difference method that is capable of solving the Two-Dimensional poison Equation and obtain an
innovative solution of poison Equation in the specified solution domain.

Statement of the problem

Consider that the following poison equation which is considered in [1] given by:
U xx +U yy =−f (x , y), ( x , y ) ∈ ( a , b ) ×(c , d ) (1)
with is subject to Dirichlet boundary condition.
U ( a , y )=U ( b , y ) =U ( x , c )=U ( x , d )=0 , (2)
f ( x , y )is assumed to be sufficiently smooth functions in D=[a , b]×[c , d ] for the existence and
the uniqueness of the solution.
Discretized the solution domain
Now we define a mesh size h and k and the constant grid point by drawing a horizontal and
vertical line of distance ‘ h ’ and ‘ k ’ respectively in the ‘ x ’ and‘ y ’ directions. These lines are
called gridlines and the point at which they interacting is known as the mesh point. The mesh
point that lies at end of the domain is called the boundary point. The solution to the problem that
lies at boundary points is called boundary condition. These boundary conditions are used to find
the solution of the given model problem at interior points. The points that lie inside the region
(inside the solution domain) are called interiors points. The goal is to approximate the solution ‘
u j n ’ at the interior mesh points. Hence we discretized the solution domain and generate a grid by
using both uniform and non-uniform discretize of grid point given as follows.
 A uniform Cartesian grid point can be generated as:-
b−a
a=x 0 < x 1< x 2 <…< x M =b , x j +1=x j + jh , h= ;
M
d−c
c= y 0< y 1< y 2< …< y N =d , y n +1= y n +nk , k = . (3)
N
 A non-uniform Cartesian grid point can be generated as:-
x j=h ×rand ( 1 ) , h=x j+1 −x j
y n=k ×rand (1), k = y n+1 − y n (4)
where j=0 ( 1 ) M ,n=0 ( 1 ) N . M And N are the maximum numbers of grid points respectively in
the x and y-direction. Then the present paper is organized as follows. Section two is a description
of numerical methods, Section three Stability, and convergence analysis, Section four is the
Results of numerical experiments. Section five is the Discussion and Section six is the
conclusion.
2. Formulation of the numerical scheme
Assuming that U (x , y) has continuous higher order partial derivative on the region
p
∂ U p
D=[a , b] x [c , d ]. For the sake of simplicity, we useU ( x j , y n )=U jn, p
=∂ x U j n and
∂x
p
∂ U p
p
=∂ y U j n for p ≥1 is pth order derivatives. To construct the scheme, assume that for the
∂y
approximate value of the following from the model problem in Eq.(1) as follow:

( U xx +U yy ) ( x 0 , y 0 ) ≈ a0 U 0 +a1 ( U 1 +U 3 ) + a2 ( U 2+ U 4 ) + a3 (U 5 +U 6+U 7 +U 8 ) (5)

f ( x 0 , y 0 ) ≈ b0 f 0+ b1 (f 1+ f 2 + f 3+ f 4 ) (6)
By using Taylor series expansion, we have:-
2 3 4 5
h 2 h 3 h 4 h 5 6
U 1=U 0+ h ∂ x U 0 + ∂ x U 0+ ∂x U 0 + ∂ x U 0 + ∂x U 0 +O(h ) (7)
2! 3! 4! 5!
2 3 4 5
h 2 h 3 h 4 h 5 6
U 3=U 0−h ∂ x U 0+ ∂x U 0− ∂x U 0 + ∂ x U 0 − ∂ x U 0 +O(h ) (8)
2! 3! 4! 5!
 2 3 4 5
h 2 h 3 h 4 h 5 6
U 2=U 0 +h ∂ y U 0 + ∂ y U 0 + ∂ y U 0 + ∂ y U 0 + ∂ y U 0+O (h ) (9)
2! 3! 4! 5!
2 3 4 5
h 2 h 3 h 4 h 5 6
U 4=U 0 −h ∂ y U 0 + ∂ y U 0− ∂ y U 0+ ∂ y U 0− ∂ y U 0 +O( h ) (10)
2! 3! 4! 5!
Adding Eq. (7) to Eq. (8) and Eq. (9) to Eq. (10), we obtain

( )
4
h 4 6
a ( U 1+U 3 )=a 1 2U 0+ h ∂x U 0 + ∂ x U 0 +O(h )
4!

( )
4
h 4 6
a ( U 2+U 4 ) =a1 2 U 0 +h ∂ y U 0 + ∂ U +O(h ) (11)
4! y 0

6 6
h 6 h 6
where T 1= ∂ x U 0 and T 2= ∂ x U 0 are their local truncation errors. Again using the Taylor
360 360
series expansion we have:

2
h 2
U 6 =U 0 +h ( ∂ x U 0+ ∂ y U 0 ) +
2!
(∂ x +2 ∂2xy U 0 +U 0 ∂2y U 0 )+ ¿

3
h 3
3!
( ∂x U 0 +3 ∂3xxy U 0 +3 ∂3yyx U 0 +∂3y U 0 ) +¿

4
h 4
4!
( ∂x U 0 +4 ∂xxxy U 0 +6 ∂ xxyy U 0 +4 ∂ xyyy U 0 +∂ y U 0 ) + ¿
4 4 4 4

5
h 5
5!
( ∂x U 0 +5 ∂ x y U 0 +10 ∂ x y U 0+10 ∂ x y U 0+5 ∂x y U 0 +∂ y U 0 )+ O(h )
5
4
5 5
3 2
5
2
5
3
6
4 (12)

2
h 2
U 5=U 0−h ( h ∂x U 0 +∂ y U 0 ) +
2!
( ∂x −2 ∂2xy U 0 +U 0 ∂2y U 0 )−¿

3
h 3
3!
( ∂x U 0 +3 ∂3xxy U 0 +3 ∂3yyx U 0 +∂3y U 0 ) +¿

4
h 4
4!
( ∂x U 0−4 ∂ xxxy U 0 +6 ∂ xxyy U 0−4 ∂xyyy U 0+ ∂ y U 0 ) −¿
4 4 4 4

5
h 5
5!
( ∂x U 0 +5 ∂ x y U 0 +10 ∂ x y U 0+10 ∂ x y U 0+5 ∂x y U 0 +∂ y U 0 )+ O(h )
5
4
5 5
3 2
5
2
5
3
6
4 (13)
 2
h 2
U 7 =U 0 +h ( h ∂ x U 0+ ∂ y U 0 ) +
2!
(∂ x + 2 ∂2xy U 0 +U 0 ∂2y U 0 )+ ¿

3
h 3
3!
( ∂x U 0 +3 ∂3xxy U 0 +3 ∂3yyx U 0 +∂3y U 0 ) +¿

4
h 4
4!
( ∂x U 0 +4 ∂xxxy U 0 +6 ∂ xxyy U 0 +4 ∂ xyyy U 0 +∂ y U 0 ) + ¿
4 4 4 4

5
h 5
5!
( ∂x U 0 +5 ∂ x y U 0 +10 ∂ x y U 0+10 ∂ x y U 0+5 ∂x y U 0 +∂ y U 0 )+ O(h )
5
4
5
3 2
5 5
2
5
3
6
4 (14)

2
h 2
U 8 =U 0−h ( h ∂x U 0 +∂ y U 0 ) +
2!
( ∂x −2∂ 2xy U 0 +U 0 ∂2y U 0 )−¿

3
h 3
3!
( ∂x U 0 +3 ∂3xxy U 0 +3 ∂3yyx U 0 +∂3y U 0 ) +¿

4
h 4
4!
( ∂x U 0−4 ∂ xxxy U 0 +6 ∂ xxyy U 0−4 ∂xyyy U 0+ ∂ y U 0 ) −¿
4 4 4 4

5
h 5
5!
( ∂x U 0 +5 ∂ x y U 0 +10 ∂ x y U 0+10 ∂ x y U 0+5 ∂x y U 0 +∂ y U 0 )+ O(h )
5
4
5
3 2
5 5
2
5
3
6
4 (15)

Now adding Eqs. (11-14) all together we obtain:

[ ]
4 h4 4
a 3 ( U 5+U 6+ U 7 +U 8 ) =a3 4 U 0 +2 h ( ∂ U 0 + ∂ U 0 ) + 4 ! ( ∂x U 0 +6 ∂ xxyy U 0+ ∂ y U 0 )
3 2 2 4 4
x y
(16)
6
+O (h )

6
h
T 3=
360
( ∂x +∂ y ) U 0 Is it a local truncation error? From the model problem we have:
6 6

f 0=−[ ∂2x U 0 +∂2y U 0 ]

∂2x f 0=−[ ∂4x U 0 + ∂4x y U 0 ]

4 2 2

∂2y f 0=−[ ∂4y x U 0+ ∂4y U 0 ]

2 2 4 (17)

Using Eq. (17), from Eq. (6) we obtain:

 b 0 f 0=−b0 [ ∂2x U 0 +∂2y U 0 ] (18)

[ ]
2
2h 2
b 1 ( f 1+ f 2 + f 3+ f 4 )=b 1 4 f 0 +
2
( ∂ x f 0 +∂ y f 0) +O(h )
2 4

[ ]
2
2h 4
¿−b 1 4 ( ∂x U 0 +∂ y U 0 ) +
2 2
2
( ∂x U 0 +∂ x y U 0 + ∂ y x U 0+ ∂ y U 0 ) +O(h )
4
4
4
2 2
4
2 2
4
4 (19)

Now substituting Eqs. (5) and (6) into the model problem in Eq. (1), we obtain:

a 0 U 0+ a1 ( U 1+ U 3 ) +a 2 ( U 2 +U 4 ) +a 3 ( U 5 +U 6+ U 7 +U 8 ) =b0 f 0+ b1 (f 1+ f 2 + f 3+ f 4)

(20)

Again substituting Eqs. (11),(16) (18) and (19) into Eq.(20) we obtain :

( ) ( )
4 4
h 4 h 4
a 0 U 0+ a1 2 U 0 +h ∂ x U 0+ ∂x U 0 + a1 2 U 0 +h ∂ y U 0 + ∂ y U 0 +¿
4! 4!

[ ]
4
4h 4
a 3 4 U 0 +2 h ( ∂ x U 0+ ∂ y U 0 ) +
3 2 2
4!
( ∂ x U 0+ 6 ∂4xxyy U 0 +∂4y U 0 ) =−b 0 [ ∂2x U 0 +∂2y U 0 ]−¿

[ ]
2
2h 4
−b 1 4 ( ∂x U 0 +∂ y U 0 ) +
2 2
2
( ∂x U 0 +∂ 4x y U 0 + ∂4y x U 0+ ∂4y U 0 )
4 2 2 2 2 4

This gives the system of linear equation in the form of:

a 0+ 2 a1+ 2 a2+ 4 a4 =0

2
h ( a 1+2 a 3 )=−b 0−4 b1

2
h ( a 2+2 a 3 )=−b 0−4 b1

2
h a3 =−2 b1

2
h ( a 1+2 a 3 )=−12 b1

2
h ( a 2+2 a 3 )=−12 b1 (21)
Now by solving the system of linear equation in Eq.(21) we obtain the value of arbitrary
−8 b 1 −2 b1 b1
constant given by:a 1=a2= 2
, a3 = 2
, a 0=40 2
, b 0=8 b1 .
h h h

Now after certain simplification with b 1=1 , in Eq.(20), we obtain the proposed scheme given by

2
−h
4 [ U 1 +U 2 +U 3 +U 4 ]+ [ U 5+ U 6 +U 7 +U 8 ] −20 U 0= ( 8 f 0 + f 1+ f 2 + f 3 +f 4 )
2

Implies that:-

4 [ U i+1 j+U i j+1 +U i −1 j+U i j−1 ]+ [ U i+1 j−1 +U i+1 j+ 1+U i−1 j+1 +U i−1 j−1 ] −20 U i j

2
h
¿− ( 8 f ij + f i j+ 1+ f i j +1+ f i−1 j + f i j−1 ) (22)
2

With its local truncation is:

( )
4 4 2 2 2
( ∆ y) 6 ( ∆ x ) 4 ( ∆ y ) 3 ( ∆ x )
4
T i , j =12h ∂ − ∂ + ∂− ∂ U (23)
360 y 360 x 12 y 12 y i , j

Hence from Eq.(22), we obtain tri-diagonal coefficient matrix of system of linear equation. For solve
this system of equation, we use the Thomas method. Because of to solve these types of system of
the equation the most recommended numerical method is the Thomas method. This is due to the
coefficient matrix contains several zero entries.

3. Stability analysis and Convergent of the proposed method

The Fourier analysis (Von-Neumann) stability analysis technique is applied to investigate the
stability analysis of the proposed method. Such an approach has been used by many researchers
like [3, 4, 13, 21, 22]. Now assume that the trial solution of the given problem at the points
( x i , y j )is
j pi h K a
u j n=λ e (24)

Where p= √−1 , K a =aπ / N , k ∈ R ( R set of a real number), λ ∈ C (C set of a complex

number) anda=1 ( 1 ) N . Substituting Eq. (24) into Eq. (22), we obtain:
4 [ λ j e ph K (i+1) + λ j +1 e pih K + λ j e ph K (i −1) + λ j −1 e pihK ]+
a a a a

+ [ λ j−1 e ph K (i +1)+ λ j+1 e ph K (i +1)+ λ j+1 e ph K (i−1) + λ j−1 e ph K (i−1) ]−20 λ j e pih K
a a a a a

2
h
¿− ( 8 λ e )
j pih K j+1 pih K j −1 ph K (i +1) j ph K (i−1) j+ 1 ph K (i−1 )
+λ e +λ e
a
+λ e
a
+λ e a a a

Divided both sides of this equation by λ j e pih K us obtain: a

4 [ e ph K + λ a +e− ph K + λ−1
a
a ] + [ λa e
a −1 ph K
+ λ a e ph K + λ a e− ph K + λ−1
a
a e
− ph K
]−20
a a a

2
h
¿−
2
(8+ λa +e ph K + e− ph K + λ−1
a
a )
a

−1 −1
16 cos ( h k a ) +8 λ a +8 λ a + 4 λ a cos ( h k a ) + 4 λ a cos ( h k a )−40=¿

2 2 2 2 −1
8 h +h λ a+ 2 h cos ⁡(h k a)+ h λ a

This implies that:

( 1
4 ) ( 1 2
4
1 2
a 2+cos ( h k a ) − h + 4− h cos ( h k a ) −2 ( 5+ h ) =0
λ a 2+cos ( h k a ) − h2 + λ−1
2
2
) [( ) ]
Multiplying both sides of the above equation by λ a we obtain:

2
λ −λ a
[
2 ( 8−h2 ) cos ( h k a )−4 ( 5+ h2 ) ] − (8−4 cos ( h k ) +h ) a
=0
2

a
( 8−4 cos ( h k a ) +h 2) (8−4 cos ( h k a ) +h2 )

2
λ −λ a
[
2 ( 8−h ) cos ( h k a )−4 ( 5+ h )
2 2
] −1=0
a
( 8+ 4 cos ( h k a ) +h ) 2

Let X =[ ( 8−h ) cos ( h k a ) −4 ( 5+ h ) ] and Y =( 8+ 4 cos ( h k a ) +h ) .

2 2 2

2 2 λa X
λ a− =1
y

( ) ( )
2 2
X X
By using perfect square we have: λ a− =1+
Y Y
 λ a=
X
Y √
± 1+
X 2 1
Y Y ( )
= ( X + √Y 2 + X 2 ) (25)

Since for any value of mesh-sizeh , | X|=|[ ( 8−h2 ) cos ( h k a )−4 ( 5+ h2 ) ]|≤ 1 and

|Y |=|( 8+4 cos ( h k a ) +h2 )|>1. Hence from Eq.(25) we have:

| √ ( ) |=|Y1 ( X +√ Y + X )|
2
X X 2 2
|λ a|= Y
± 1+
Y

≤ |Y1 ||X ± √ Y + X | 2 2

¿ |Y1 ||X|+√|Y |+|X | Triangular inequality

2 2

≤ |Y1 ||X|+√|Y| +|X| < 1

2 2

Hence we obtain the required criteria for stability investigation of the proposed method.
Therefore the proposed method is strictly stable for solving two-dimensional poison equations.

Theorem 2: The difference equation given in the form of Eq.(12) is stable if for which the
eigenvalues of the coefficient matrix of the system of the differential equation are satisfied
Real ( λ j ) < 0. Proof: See reference [3].

Since from the principal part of the local truncation error, the derived local
truncation error for the proposed scheme is

6
h
T i,j=
360
(2 ∂6x + 9 ∂6xxxy +14 ∂6xxyy U 0 +9 ∂6xyyy + 2∂ 6y ) U i , j

lim ‖T i , j‖=lim h6
h→ 0 h →0
‖3601 (2 ∂ +9 ∂
6
x
6
xxxy +14 ∂6xxyy U 0 +9 ∂6xyyy +2 ∂6y )U i , j → 0 ‖ Thus this implies

that,T i j → 0 ash → 0. So that, the scheme is consistent with the order of

O ( ∆ x 6 +∆ y 6 ) =0(h6) . Hence the scheme is convergent.
Criteria for Investigating the Accuracy of the Method

This section presented the criteria that the accuracy of the present method is investigated. The
accuracy of the solution will depend on how small we make the step size, ∆ x=∆ y=h . To test
the performance of the proposed method to give an accurate solution for the given model
problem, maximum absolute error, L2 and L∞ norms are calculated by using the following
formula:

√
N
L∞= max |u ( xi , y j ) −ui , j| , L = 1 ∑|u ( x i , y j )−ui j|
2
i=1 ( 1 ) M
2
1≤ n ≤N
N j=0

where N is the maximum number of step,u ( x i , y j ) is the exact solution and ui j approximation
solution of the poison equation in Eq.(1) at the grid point ( x i , y j ).

4. Numerical experiments and their Results

To test the validity of the proposed method, we have considered the following three model
problem considered in [1]. Numerical results and errors are computed and the outcomes are
represented tabularly and graphically.
Example 1: Consider the classical two-dimensional equation considered in [1]
2
u xx +u yy =−2 π sin ( πx ) sin ⁡( πy ) , ( x , y ) ∈ ( 0 ,1 ) × ( 0 ,1 )
The subjected Dirichlet boundary condition is given by:
U ( 0 , y )=u ( 1 , y )=u ( x , 0 ) =u ( x , y )=0
1
∆ x=∆ y=h= and the exact solution is given by:
40
u ( x , y )=sin ( πx ) sin ⁡(πy )

Table 1 Comparison of Pointwise maximum absolute error (L¿¿ ∞)¿ and root mean square error
(L2 ) with uniform mesh size equal ∆ x=∆ y=h=1/40

Pointwise maximum
Specified Grid absolute error obtained Pointwise maximum absolute error and
points by Genet Mekonnen and root mean square error by present
Lemi Guta in [1] methods
 x y L∞ L∞ L2
1/4 1/4 2.65E-02 2.5710e-04 4.0651e-05
1/2 1/4 3.75E-02 3.6359e-04 5.7489e-05
3/ 4 1/4 2.65E-02 2.5710e-04 4.0651e-05
1/4 1/2 3.75E-02 3.6359e-04 5.7489e-05
1/2 1/2 5.30E-02 5.1420e-04 8.1302e-05
3/ 4 1/2 3.75E-02 3.6359e-04 5.7489e-05
1/4 3/ 4 2.65E-02 2.5710e-04 4.0651e-05
1/2 3/ 4 3.75E-02 3.6359e-04 5.7489e-05
3/ 4 3/ 4 2.65E-02 3.6359e-04 2.5710e-04

1
1

0.8
0.8

0.6
0.6
U(x,y)
U(x,y)

0.4
0.4

0.2
0.2

0 0
1 1

1 1
0.8 0.8
0.5 0.5 0.6
0.6
0.4 0.4
0.2 0.2
0 0 0 0
y x y x

Figure 1. Physical Behavior of Approximate solution Figure 2. Physical Behavior of Exact solution

forgiven example on uniform mesh size for given example on uniform mesh size

∆ x=∆ y=h=1/40 ∆ x=∆ y=h=1/40

 0.08
U-ap at x=0.75
0.07 U-exact at x=0.75

0.06

0.05
U(x,y)

0.04

0.03

0.02

0.01

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Figure 3. Variation of exact versus numerical solution for given example with uniform mesh size
∆ x=∆ y=h=1/40

Table 2. Comparison of Pointwise maximum absolute error(L¿¿ ∞ )¿ and root mean square error
(L2 ) with non-uniform mesh size equal
Pointwise maximum absolute
error obtained by Genet Pointwise maximum absolute error and root mean
Specified Grid
Mekonnen and Lemi Guta in [1] square error by present methods
points
x y L∞ L∞ L2
0.127 0.0975 3.51 E−02 4.3477e-04 1.0248e-04

0.6324 0.0975 6.92 E−02 2.0690e-03 4.8766e-04

0.8147 0.0975 1.015 E−01 1.5356e-03 3.6195e-04

0.127 0.2285 9.16 E−02 2.2017e-03 5.1895e-04

0.6324 0.2785 3.532 E−01 1.6866e-03 3.9754e-04

0.8147 0.2785 1.269 E−01 6.6610e-04 1.5700e-04

0.127 0.5469 9.59 E−02 2.2450e-03 5.2914e-04

0.6324 0.5469 2.105 E−01 1.1009e-03 2.5948e-04

0.8147 0.5469 7.33 E−02 1.4153e-03 3.3358e-04

 1.4 1
1.2
0.8
1

0.8 0.6
U(x,y)

U(x,y)
0.6
0.4
0.4
0.2
0.2

0 0
1 1
1 1
0.8 0.8
0.5 0.6 0.5 0.6
0.4 0.4
0.2 0.2
y 0 0 y 0 0
x x

Figure 4. Physical Behavior of Approximate solution for Figure 5. Physical Behavior of Exact solution for

given example on non-uniform mesh size given example on non-uniform mesh size

U-ap at x=0.75
U-exact at x=0.75
-0.5
10

-0.6
10
U(x,y)

-0.7
10

-0.8
10

-0.9 -0.7 -0.5 -0.3 -0.1

10 10 10 10 10
x

Figure 6.Variation of exact versus numerical solution for given example on non-uniform mesh size
 -4
x 10
1.4

1.2
Maximum point wise Absolute error)
Max-Abs.Error with uniform grid point
Max-Abs.Error with non-uniform grid point
1

0.8

0.6

0.4

0.2

0
-2 -1 0
10 10 10
x

Figure 7. Variation of Point-wise absolute errors between exact and numerical solution forgiven

example on uniform versus non-uniform mesh size.

5. Discussions

In this paper, we presented a fourth-order compact finite difference method

to obtaining an innovative solution for two-dimensional poison equations. The
innovative solution, obtained within the fourth-order compact finite difference method, discussed
only the case of the Dirichlet boundary condition. Regarding this partial differential equation, we
note that there are two main ways of compact finite difference discretization (over a uniform grid
and non-uniform grid points). When we apply the compact finite difference method to the
continuous two-dimensional Poisson equation the equation is replaced by a "discrete"
approximation. The number of those discrete points can be selected uniformly or non-uniformly
depending on the mesh size (h). The mesh is the set of locations where the discrete solution is
computed. Two key parameters of the mesh are and the local distance between adjacent points in
space. Fourth-order Compact finite difference discretization is simple to implement by using
both equal mesh size and non-uniform mesh size as sow above in the table and graph. The full
discretization of the poison equation by the present method leads to the system of linear
equations which is solved by using the Thomas method. The convergence has been shown in the
sense of maximum pointwise absolute error norm ( ) and root mean error(L¿¿ 2)¿ their values
are given in table and graph to comparing them with pre-existing results. The stability and
convergence of the present method are also investigated by using the Von-Neumann technique.
The results presented in Tables 1 and 2 demonstrates fourth-order finite difference method gives
a more accurate numerical solution than the pre-existing method in the literature. As we see from
Figure 7 the present method is more accurate when we investigate the solution of the model
problem on non-uniform grid point discretization of the solution domain. Moreover, Figures 3
and 6 specifies that the present method gives an accurate solution for the 2D poison equation on
both uniform and non-uniform grid point discretization of solution domain and the approximate
exact solution very well.

6. Conclusion

The key purpose of this work is to formulate and investigate the fourth-order compact
finite difference method for solving two-dimensional poison equations. To
further collaborate the applicability of the proposed method; tables of point-
wise absolute error and root mean square error and graphs have been
plotted for Examples 1, for the exact solution versus the numerical solutions
at different values of x on both uniform and non-uniform grid points. Table 1, shows
the absolute errors obtained by my fourth-order compact finite difference
method have been compared with absolute errors obtained by [1] on uniform
grid points and it shows that the present method is the more convergent
method. Table 2, also shows the absolute errors obtained by the present
method have been compared with absolute errors obtained by [1] on non-
uniform grid points and then also it is showing that the present method is
accurate than the previous method. Generally, the present method is
computational: stable, effective, simple to use, convergent, and gives an
accurate solution than some previously existing methods.

Acknowledgements
The author wish to express his thanks to the authors of literatures for the provision of initial idea
for this work.
Last but not least, the author wish to express his thanks to his beloved friends Hailu Muleta
Chemeda; your kindness and helps will be a great memory for me.
References

[1] Genet Mekonnen Assefa and Lemi Guta .Solution of Two Dimensional Poisson Equation

Using Finite Difference Method with Uniform and Non-uniform Mesh Size, Advances in

Physics Theories and Applications, 79(2019), 1-10.

[2] Sukumar N. and Bolander J. E. Numerical Computation of Discrete Differential Operators on

Non-Uniform Grids Department of Civil and Environmental Engineering, One Shields

Avenue, University of California, Davis, U.S.A, CA. (2003), 95616,

[3] K. A. Koroche. Numerical Solution for One Dimensional Linear Types of Parabolic Partial

Differential Equation and Application to Heat Equation. Mathematics and Computer Science. 5(4)

(2020) 76-85.

[4] M. Dugassa, G. File, T. Aga. Numerical Solution of Second-Order One Dimensional Linear

Hyperbolic Telegraph Equation. Ethiopian Journal of Education and Sciences, 14(1) (2018)

39-52.

[5] Madhumangal Pal. Numerical Analysis for Scientists and Engineers: Theory and C
Programs.

[6] Gayaz Khakimzyanov G.K .On supra convergence phenomenon for second order centered

Finite differences on non-uniform grids, Institute of Computational Technologies,

Novosibirsk, Russia. Last modified: February 6(2017).

[7] Yousif W.S. and Martins M.M. (2008), Explicit De-couple Group AOR method for solving

elliptic partial differential equations, Neural, Parallel and Scientific Computations, no. 4,

531-542.

[8] I.M. Navon and H.A. Riphagen. SHALL4–An implicit compact fourth-order FORTRAN
 program for solving the shallow-water equations in conservation-law form. Computers &

Geosciences, 12(1986), 129–150.

[9] B Neta. Numerical Solution Of Partial Differential Equations MA 3245 Lecture Notes,(2003)

[10] Chew, C.S. A Generalized Finite-Difference (GFD) ALE scheme for in

compressible flows around moving solid bodies on hybrid mesh free Cartesian grids.

Journal of Computational Physics, 218(2006), 510-545.

[11] Randall J.Leveque. Finite Difference Methods for Differential Equation, Spring Quertes,

(1998).

[12] Profile, S., Atm, A.T.M. And Mil-Std, T.A.H.O., Notice Of Change Inch-Pound Mil-Std-

188-176 Notice 1 20 Oct 97.

[13] S. Shokofeh, J.Rashidinia. Numerical solution of the hyperbolic telegraph the Equation by cubic B-
spline collocation method. Applied Mathematics and Computation 281(2016) 28–38.

[14] M.H. Carpenter, D. Gottlieb and S. Abarbanel. Time-stable boundary conditions for finite-

difference schemes solving hyperbolic systems: methodology and application to high-

order compact system. J. Comp. Phys., 111(1994),220–236.

[15] M.K. Jain, S.R.K. Iyengar, R. K. Jain. computational methods for partial
differential
equations, New age international publishers, 2007.

[16] Eyaya Fekadie Anley. Numerical Solutions of Elliptic Partial Differential Equations by

Using Finite Volume Method, Pure and App. Math. J.; 5(4) (2016); 120-129.

[17] Abdelwahab Kharab and Ronald B. Guenther. An Introduction to Numerical Methods

A

MATLAB Approach, Fourth-edition, Taylor & Francis Group, (2019).

[18] Yogash Jaluria. Computer Methods for Engineering with MATLAB® Applications, second

edition, Taylor & Francis Group LLC.(2011).

[19] I.M. Navon and H.A. Riphagen. An implicit compact fourth order algorithm for solving the

Shallow-water equations in conservation-law form. Monthly Weather Review, 107:1107–

1127, 1979

[20] Mazarei, M.M. and Aminataei, A. Numerical solution of Poisson's equation using a

combination of logarithmic and multiquadric radial basis function networks. Journal of

Applied Mathematics, 2012.

[21] Aliyi, K. and Muleta, H., Numerical Method of the Line for Solving One Dimensional

Initial-Boundary Singularly Perturbed Burger Equation.

[22] S.Gowrisankar, and S. Natesan .Uniformly convergent numerical method for singularly

perturbed parabolic initial-boundary-value problems with equidistributed grids, Int. J.

Comput. Math. 2014b, 91 (3) , 553–577.

AUTHORS PROFILE

Kedir Aliyi Koroche is Graduate B.sc in Mathematics from Ambo University in 2016 G.C. and
M. Sc. in Mathematics (specialized in Numerical Analysis) from Jimma University in 2019 G.C.
Now he works at Ambo University as Lecturer and Researcher. He has two published research
paper.