EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 3, Article Number 6524
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Finite Elements Method for Solving Elliptic and
Hyperbolic Partial Differential Equations

Abdulkafi Mohammed Saeed1,∗, Shahad Saleh Ayad Almutairi1

1 Department of Mathematics, College of Science, Qassim University, Buraydah,
Kingdom of Saudi Arabia

Abstract. In various disciplines of engineering such as hydrodynamics, heat conduction, geome-
chanics, civil engineering, nuclear, and even biomedical, the finite element method (FEM) is known
to be an effective and sophisticated technique to tackle complex computational problems. The main
idea of FEM is to break down a complicated space or domain into a number of small, countable, and
finite elements, or “finite elements” by approximation methods whose behavioral outputs could be
predicted using simpler equations. In this work, we will analyze the comparison between the Finite
Difference Method (FDM) and FEM in order to see the effectiveness of each method. Within the
scope of this research, several types of partial differential equations (PDEs) will be solved with
FEM to estimate a solution that is as accurate as possible to the exact answer. Computations will
be completed with the MATLAB software, and a discussion will follow.

2020 Mathematics Subject Classifications: 76S05, 76D10, 80A32, 35Q79

Key Words and Phrases: Boundary value problem, Finite difference method, Galerkin Finite
element method, Poisson’s equation, Wave equation

1. Introduction

Physicists and engineers apply Poisson’s equation to solve problems in areas including
electrostatics, magnetostatics, quantum mechanics, and fluid dynamics.Poisson’s equation
explains the relationship between the electric potential and the charge distribution in elec-
trostatics. It links the current distribution and the magnetic potential in magnetostatics.
This equation is vital for scientists and engineers since it aids in comprehending how var-
ious physical systems behave [1–4].
The wave equation describes the propagation of sound, light, water waves, and other waves
through various media. That makes the equation useful in applied mathematics, physics,
and engineering. This second-order PDE, also known as the wave equation, is crucial for
analyzing the dynamic nature of physical phenomena. The wave equation is vital when
studying dynamic systems[5, 6].

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6524

Email addresses: abdulkafi.ahmed@qu.edu.sa ([Abdulkafi M. Saeed)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 2 of 21

Real-world problems come with differential equations that are complicated and often
have no explicit solution formulas. When they do exist, they involve integrals that have
to be calculated through numerical quadrature methods. Numerical analysis has become
extremely important in solving differential equations—under specific conditions of initial
or boundary values—making this approach effective for overcoming numerous challenges.
Partial differential equations (PDEs) can be solved numerically using the Finite Volume
Method (FVM), FDM, and FEM. FDM creates approximations of solutions to differential
equations via local Taylor series expansion. When dealing with increasingly complex ge-
ometries in higher dimensions, the topologically square network of lines may present certain
challenges in the discretization of PDEs in FDM. These restrictions made integral-based
formulations more desirable, which allowed for the creation of the two less sophisticated
FVM and FDM techniques [1, 5, 7].
The primary application of FEM is in the solution of differential equations with boundary
conditions. The solution is calculated by partitioning the domain into finite number of
regions, or elements, and applying basis functions to them. Hydrodynamic and structural
analysis, as well as multidisciplinary problems are some of the many fields where FEM is
implemented today. The reason for this tremendous use is due to the method’s precision in
modeling and numerically solving intricate problems. By enabling the simulation of nat-
ural phenomena and the analysis of engineering designs, both scientific and engineering
evaluation are possible. This permits the prediction of system performance under design
conditions and natural phenomena, which aids in the development of safer, cost-effective
systems. The risk of disasters can also be reduced with the aid of FEM. The quality of
life is highly influenced by the technology we have today; this is gradually advanced by a
growing understanding of nature through FEM’s versatile nature [8] .
Galerkin finite element method (GFEM) is formulated as a modification of the basic pro-
cess of FEM, construction by Galerkin’s method. It is simpler than other methods because
GFEM uses a systematic approach with piecewise trial functions for the approximation of
partial differential equations over finite elements [6, 9, 10].
This study is an extension of the work done by S. S. Almutairi, A. M. Saeed [11] to fur-
ther investigate and scrutinize alternative approximate methods for solving various types
of boundary value problems.

This study explores the applications of GFEM and compares it with FDM to deter-
mine which approach is more effective. The paper is structured into five sections. Section
2 provides a brief overview of the FDM. Section 3 presents the GFEM formulations for
solving PDEs of both elliptic and hyperbolic types. Section 4 discusses numerical exper-
iments and compares the two methods. Finally, the paper concludes with a summary of
findings and final remarks.

2. Brief Description of Finite difference Method

FDM is a numerical method that approximates the solution of differential equations
by means of discretization. Breaking continuous functions into discrete values simplifies



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 3 of 21

calculations and enhances various types of quantitative analysis. This transformation
makes it easier to work with data, optimize processes, and apply computational methods
effectively. It is quite popular in various fields of science and technologie that involve
differential equations, like physics, engineering, or even finance. There are three different
kinds of finite difference approximations: forward distinctions, backward distinctions, and
center distinctions. By evaluating the function at two future or present values, the forward
difference approach calculates the derivative. Both the values at the current location and
the previous one are of interest to backward differences. Center is more accurate than
the others, as it takes the average of forward and backward difference, making it a more
accurate form of approximation. The Forward Difference Method is one of the numerical
techniques developed over many decades [12, 13]

In two dimensions, the Poisson’s equation can be expressed as [14, 15].

∂2u

∂x2
+
∂2u

∂y2
= g(x, y), (1)

with the boundary constraint u = f(x, y) along the boundary C.
In this case, we additionally assumed that the mesh points are consistent across the x
and y dimensions. This assumption makes it possible to reduce the (1) central difference
approximation to

ui,j =
1

4
(ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − h2gi,j),

where gi,j = g(xi, yi).

Let i, j = 0, 1, 2, 3, 4 and u = 0 along the boundary C. For j = 0, 1, 2, 3, 4, then
ui,0 = 0, ui,4 = 0, and u0,j = 0, u4,j = 0. When i = 0, 1, 2, 3, 4. In Figure 1, the boundary
values (filled circles) are displayed.

Figure 1: FDM mesh points.

An example of this would be i, j = 0, 1, 2, 3. As a result, equation (1) transforms into
a system of nine equations with nine unknowns. These equations can be represented in
matrix form for a more structured and efficient notation.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 4 of 21



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





u1,1
u1,2
u1,3
u2,1
u2,2
u2,3
u3,1
u3,2
u3,3


=



−h2g1,1
−h2g1,2
−h2g1,3
−h2g2,1
−h2g2,2
−h2g2,3
−h2g3,1
−h2g3,2
−h2g3,3


Consequently, a system of N equations was created by (1). Here, n represents the number
of subintervals along the x and y axes. The coefficient matrix is both positive definite and
symmetric, with many elements being zero, making it sparse. Due to this sparsity, the
system of equations is best solved using an iterative approach rather than a direct method.

We consider a hyperbolic PDE, the numerical solution to the wave equation. The
differential equation provides the wave equation [16] :

∂2u

∂t2
(x, t)− α2∂

2u

∂x2
(x, t) = 0, 0 < x < l, t > 0, (2)

based on the terms

u(0, t) = u(l, t) = 0, for t > 0,

u(x, 0) = f(x), and
∂u

∂t
(x, 0) = g(x), for 0 ≤ x ≤ l,

where α is a constant that depends on the problem’s physical circumstances.
To define the x-axis grid points using h = l/m, choose an integer m > 0. Choose a
time-step size k > 0 as well. The definition of the mesh points (xi, tj) is

xi = ih, and tj = jk,

for each i = 0, 1, . . . ,m and j = 0, 1, . . ..
The wave equation changes to (xi, tj) at any inner mesh point.

∂2u

∂t2
(xi, tj)− α2∂

2u

∂x2
(xi, tj) = 0. (3)

The difference technique is generated using the centered-difference quotient for the
second partial derivatives given by

∂2u

∂t2
(xi, tj) =

u (xi, tj+1)− 2u (xi, tj) + u (xi, tj−1)

k2
− k2

12

∂4u

∂t4
(xi, µj) ,

where µj ∈ (tj−1, tj+1), and



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 5 of 21

∂2u

∂x2
(xi, tj) =

u (xi+1, tj)− 2u (xi, tj) + u (xi−1, tj)

h2
− h2

12

∂4u

∂x4
(ξi, tj) ,

where ξi ∈ (xi−1, xi+1). Substituting these into (3) gives

u (xi, tj+1)− 2u (xi, tj) + u (xi, tj−1)

k2
− α2u (xi+1, tj)− 2u (xi, tj) + u (xi−1, tj)

h2

=
1

12

[
k2
∂4u

∂t4
(xi, µj)− α2h2

∂4u

∂x4
(ξi, tj)

]
.

Ignoring the incorrect phrase

τi,j =
1

12

[
k2
∂4u

∂t4
(xi, µj)− α2h2

∂4u

∂x4
(ξi, tj)

]
, (4)

causes the distinction, formula

wi,j+1 − 2wi,j + wi,j−1

k2
− α2wi+1,j − 2wi,j + wi−1,j

h2
= 0.

Define λ = αk/h. The difference equation can therefore be written as

wi,j+1 − 2wi,j + wij−1 − λ2wi+1,j + 2λ2wi,j − λ2wi−1,j = 0.

Then calculate the most advanced time-step estimate, wij+1, to obtain

wi,j+1 = 2
(
1− λ2

)
wi,j + λ2 (wi+1,j + wi−1,j)− wi,j−1. (5)

For any i = 1, 2, . . . ,m − 1 and j = 1, 2, . . ., this equation is true. The boundary
conditions give

w0,j = wm,j = 0, for each j = 1, 2, 3, . . . , (6)

and the initial condition implies that

wi,0 = f (xi) , for each i = 1, 2, . . . ,m− 1.

This set of equations can be written in matrix form to give


w1,j+1

w2,j+1
...

wm−1,j+1

 =



2(1− λ2) λ2 0 · · · 0
λ2 2(1− λ2) λ2 · · · 0

0
... · · · · · ·

...
...

. . .
...

0
... λ2 2(1− λ2)


=


w1,j

w2,j
...

wm−1,j

−


w1,j−1

w2,j−1
...

wm−1,j−1


(7)



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 6 of 21

According to equations (4) and (7), values from the j th and (j − 1) st time steps are
needed for the (j + 1) st time step. Figure 2 is referred to. Because (6) gives values for
j = 0, but the starting-velocity condition has to give values for j = 1, which are needed
in (4) to compute wi,2, this leads to a minor initial problem.

∂u

∂t
(x, 0) = g(x), 0 ≤ x ≤ l.

Figure 2: Discretization of mesh points

One method is to use a forward-difference approximation to replace ∂u/∂t.

∂u

∂t
(xi, 0) =

u (xi, t1)− u (xi, 0)
k

− k

2

∂2u

∂t2
(xi, µ̃i) ,

for some µ̃i in (0, t1).The equation’s u (xi, t1) can be solved to yield

u (xi, t1) = u (xi, 0) + k
∂u

∂t
(xi, 0) +

k2

2

∂2u

∂t2
(xi, µ̃i)

= u (xi, 0) + kg (xi) +
k2

2

∂2u

∂t2
(xi, µ̃i) .

When the truncation phrase is removed, the approximate

wi,1 = wi,0 + kg (xi) , for each i = 1, . . . ,m− 1.

The truncation error in our approximation is only O(k), whereas the truncation error
in (5) is O

(
k2
)
. A better approximation to u (xi, 0) is obtained by expanding in t, u (xi, t1)

is a second Maclaurin polynomial. Then,

u (xi, t1) = u (xi, 0) + k
∂u

∂t
(xi, 0) +

k2

2

∂2u

∂t2
(xi, 0) +

k3

6

∂3u

∂t3
(xi, µ̂i) ,

for some µ̂i in (0, t1). If f
′′ exists, then

∂2u

∂t2
(xi, 0) = α2∂

2u

∂x2
(xi, 0) = α2d

2f

dx2
(xi) = α2f ′′ (xi) ,



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 7 of 21

and

u (xi, t1) = u (xi, 0) + kg (xi) +
α2k2

2
f ′′ (xi) +

k3

6

∂3u

∂t3
(xi, µ̂i) .

An estimate with error O
(
k3
)
is therefore produced:

wi1 = wi0 + kg (xi) +
α2k2

2
f ′′ (xi) .

The difference equation can be used if f ∈ C4[0, 1] but f ′′ (xi) is not readily available.

f ′′ (x0) =
1

h2
[f (x0 − h)− 2f (x0) + f (x0 + h)]− h2

12
f (4)(ξ), (8)

to write

f ′′ (xi) =
f (xi+1)− 2f (xi) + f (xi−1)

h2
− h2

12
f (4)

(
ξ̃i

)
,

for some ξ̃i in (xi−1, xi+1). This implies that

u (xi, t1) = u (xi, 0) + kg (xi) +
k2α2

2h2
[f (xi+1)− 2f (xi) + f (xi−1)] +O

(
k3 + h2k2

)
.

Because λ = kα/h, we can write this as

u (xi, t1) = u (xi, 0) + kg (xi) +
λ2

2
[f (xi+1)− 2f (xi) + f (xi−1)] +O

(
k3 + h2k2

)
=

(
1− λ2

)
f (xi) +

λ2

2
f (xi+1) +

λ2

2
f (xi−1) + kg (xi) +O

(
k3 + h2k2

)
.

Consequently, the difference equation,

wi,1 =
(
1− λ2

)
f (xi) +

λ2

2
f (xi+1) +

λ2

2
f (xi−1) + kg (xi) .

For all i = 1, 2, . . . ,m− 1, wi,1 may be found using this method. We apply the system in
(7) to determine subsequent approximates.

3. Formulation of the Finite element Method

GFEM solves differential equations numerically by dividing the domain into smaller
finite elements. Within each of these elements, the function is approximated using piece-
wise trial functions, allowing for an efficient and flexible solution approach [17, 18]. The
following steps make up the GFEM for solving a differential equation [17, 18]:



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 8 of 21

(i) The differential equation is multiplied by a weight function ω(x) they make up the
integral for the entire domain;

(ii) Integrate by parts if required to lower the highest order term’s order;

(iii) Choose the interpolation order—whether linear, quadratic, or higher—and define the
corresponding shape functions, Ni, for i = 1...m. Using the trial function p = p̃(x) =
m∑
i=1

Ni(x)pi, the function is approximated based on these selected shape functions;

(iv) Calculate all integrals within each element, either through exact analytical methods
or numerical approximation, to establish a system of equations for the unknown
values pi′s;

(v) Solve the system of equations to determine the unknown values pi′s, ensuring accu-
racy in the numerical solution process.

Let’s examine the model problem: Poisson’s equation, subject to homogeneous Dirich-
let boundary conditions [19].

−∆u = f in Ω, (9)

u = 0 on ∂Ω.

In this case, it is clear that a(·, ·) is both symmetric and bilinear, satisfying the property
a(u, u) = |u|21,Ω := ||u||2. Furthermore, when a(u, u) = 0, it implies that ∇u = 0, meaning
that u must be a constant function throughout the domain. The value of this constant
should be zero as u|Γ = 0. Consequently, the Riesz representation theorem indicates
that (9) has a unique solution since a(·, ·) forms an inner product on V . The Galerkin
methods are a class of techniques that are used to estimate the solution to (9). Let Vh ⊂ V
be a finite dimensional subspace. Limit the variational formulation to the subspace Vh,
meaning that the goal is to find uh ∈ V such that it satisfies the given conditions within
this restricted space.

a(u, v) = (f, v), for all v ∈ V, (10)

where

a(u, v) =

∫
Ω
∇u∇v dx, (f, v) =

∫
Ω
fv dx for all f ∈ L2(Ω). (11)

It is evident that in this scenario, a(·, ·) is symmetric and bilinear, and a(u, u) = |u|21,Ω :=

||u||2 Moreover, a(u, u) = 0 suggests that ∇u = 0, and as a result, u is constant. Since
u|Γ = 0, this constant’s value should be 0. Consequently, the Riesz representation theorem
indicates that (9) has a unique solution since a(·, ·) forms an inner product on V . The
Galerkin methods are a class of techniques that are used to estimate the solution to (9).
Let Vh ⊂ V be a subspace of finite dimensions. Find uh ∈ V s.t. by restricting the
variational form in the subspace Vh.

a(uh, vh) = (f, vh), for all vh ∈ V. (12)



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 9 of 21

Vh = span{φ1, φ2, ..., φN} allows us to Every function v ∈ Vh has a unique representation:
To calculate v =

∑N
i=1 viφi As a result, one can define an isomorphism. Vh ∼= RN by

v =
N∑
i=1

viφi ←→ v = (v1, ..., vN )T .

We then refer to v as the coordinate vector of v with regard to the basis {φ}Ni=1. We
present the stiffness matrix after introducing the elasticity language.

A = (aij)N X N with aij = a(φi, φj),

and the load vector f = {⟨f, φ⟩}Nk=1 ∈ RN . The following linear algebraic system can
therefore be used to formulate the variational (12) on Vh.

Au = f.

The a(·, ·) -inner product of two functions, vu ∈ Vh is realized via the matrix product by
definition.

a(uh, vh) = a(
∑
i

uiφi,
∑

vjφj) =
∑
i,j

a(φi, φj)uivj = vtAu.

Therefore, for any vector u ∈ RN , uTAu = a(u, u) ≥ 0 and equals 0 if and only if u is
zero. Namely, A is an SPD matrix, and thus the solution u = A−1f exists and is unique.
Once we get the coefficient vector u, we can use a linear combination of basis functions
to generate uh. A well-known and frequently applied example of Galerkin methods is
the FEM, which builds a finite-dimensional subspace Vh from domain triangulations τh.
The domain’s decomposition into a limited number of pieces is the source of the name.
Piecewise polynomials are commonly used to define a finite dimensional space.

The scalar wave equation that models the propagation of acoustic waves in a bounded
domain Ω3 with boundary Γ is examined : [20–22].

1

c(x)2
∂2u

∂t2
−△u = 0, in Ω× (0, T ),

u(·, 0) = 0,
∂u

∂t
(·, 0) = f, in Ω,

∂nu|Γ1
= 0, on Γ× (0, T ).

(13)

Where t is the time variable, T is a final time, and f is a load vector. u(x, t) is the pressure,
and c(x) is the wave speed based on x = (x1, x2, x3) ∈ Ω.
We now use continuous piecewise linear functions in space and time to build a FEM for
(13). The standard method for discretizing Ω× (0, T ) is to use Kh = {K} to partition the
domain Ω into tetrahedra K(h = h(x), where Jk = {J} is a partition of the time interval
(0, T ) into time intervals J = (tk−1, tk) .k = 1, . . . , N , of uniform length τ = tk − tk−1.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 10 of 21

The finite element approach for (13) is formulated by introducing the finite element space
W u

h , which is defined as follows :

W u :=
{
u ∈ H1(Ω× J) : u(·, 0) = 0, ∂nu|Γ = 0

}
,

W u
h :=

{
u ∈W u : u|K×J ∈ P1(K)× P1(J),∀K ∈ Kh,∀J ∈ Jk

}
.

(14)

The set of piecewise linear functions on K and J , respectively, is denoted by P1(K)
and P1(J).

We use CG(1) in space and time to construct a discrete scheme, and we look for a
discrete solution in the space W u

h for u spanned by the functions.

u(x, t) =
N∑
l=0

M∑
i=1

uliφi(x)ψl(t). (15)

In space and time, the typical continuous piecewise linear functions φi(x) and ψl(t)
are represented, respectively.
The system of linear equations that results from substituting this into (13) is as follows:

M
(
uk+1 − 2uk + uk−1

)
= τ2F k − τ2A

(
1

6
uk−1 +

2

3
uk +

1

6
uk+1

)
, k = 1, . . . , N − 1,

(16)
with initial conditions :

u(0) =
∂u

∂t

∣∣∣∣
t=0

= 0.

Here, F k is the load vector at time level tk, wherek = 1, 2, 3 . . ., mathbfu is the
unknown discrete field values of u, M is the mass matrix in space, A is the stiffness
matrix, and τ is the time step.

For each element e, the explicit formulas for the entries in system (16) can be provided
as follows :

M e
i,j =

(
1

c2
φi, φj

)
e

,

Ke
i,j = (∇φi,∇φj)e ,

F e
j = (f, φj)e .

Where the L2(e) scalar product is indicated by (., .)e. The contribution of element e
to the global assembled matrix in space M is represented by the matrix Me, while the
contribution of element e to the global assembled matrix K is represented by the matrix
Ke.

W u
h is a vector space with finite dimensions, dimW u

h =M ×N , with a basis consisting

of the standard continuous piecewise linear functions {φi}Mi=1 in space and hat functions

{ψl}Nl=1 in time, where:



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 11 of 21

Figure 3: Discretization of ψ(t)

ψl(t) =


0, if t /∈ [tl−1, tl+1] ,
t−tl−1

ti−tl−1
, if t ∈ [tl−1, tl] ,

tl+1−t
tl+1−tl

, if t ∈ [tl, tl+1] .

(17)

The discrete system of equations Using basis of functions {φi}Mi=1 in space and {ψl}Nl=1

in time we have: uh =
∑N

l=0

∑M
i=1 u

l
hiφi(x)ψl(t). Substituting into (13) we get:

−
N∑
l=0

M∑
i=1

ul
hl
i

∫
Ω

1

c2
φi(x)

∫ tl+1

tl−1

∂ψl(t)

∂t

∂v(x, t)

∂t
dxdt

+
N∑
l=0

M∑
i=1

ulhi

∫
Ω

∫ tl+1

tl−1

∇φi(x)∇v(x, t)ψl(t)dxdt

=

∫
Ω

∫ T

0
f(x)v(x, t)dxdt ∀v ∈W u

h .

We take v(x, t) = φj(x)ψm(t) and get:

−
N∑
l=0

M∑
i=1

ulhi

∫
Ω

1

c2
φi(x)φj(x)

∫ tl+1

tl−1

∂ψl(t)

∂t

∂ψm(t)

∂t
dxdt

+

N∑
l=0

M∑
i=1

ulhi

∫
Ω
∇φi(x)∇φj(x)

∫ tl+1

tl−1

ψl(t)ψm(t)dxdt

=

∫
Ω

∫ T

0
f(x)v(x, t)dxdt ∀v ∈W u

h .

(18)

The vector of unknown coefficients is represented by U = ulhi, M = (mij), and A =
(aij) , are stiffness and mass matrices. Accordingly, M ×M in space, with coefficients

mij =

∫
Ω
φi(x)φj(x)dx,

aij =

∫
Ω
∇φi(x)∇φj(x)dx.

K = (klm) , P = (plm) are matrices of mass and stiffness. N ×N with coefficients in time



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 12 of 21

klm =

∫ tl+1

tl−1

∂ψl(t)

∂t

∂ψm(t)

∂t
dt,

klm =

∫ tl+1

tl−1

ψl(t)ψm(t)dt,

as well as the load vector b = (bjm) with coefficients

bjm =

∫
Ω

∫ tl+1

tl−1

f(x)φj(x)ψm(t)dxdt.

K = (klm) and P = (plm) are first calculated. As you can see, klm = 0, plm = 0 unless
when l = m− 1, l = m, l = m+ 1. Using the test function specification (17), we calculate
the first diagonal elements kll :

kll =

∫ tl+1

tl−1

ψ′
l(t)ψ

′
l(t)dt

=

∫ tl

tl−1

(
1

τ

)2

dt+

∫ tl+1

tl

(
−1
τ

)2

dt =
2

τ
,

pll =

∫ tl+1

tl−1

ψl(t)ψl(t)dt

=

∫ tl

tl−1

(
t− tl−1

τ

)2

dt+

∫ tl+1

tl

(
tl+1 − t

τ

)2

dt =
2

3
τ.

Similarly, the other elements can be written as

kl,l+1 =

∫ tl+1

tl−1

ψ′
l(t)ψ

′
l+1(t)dt

=

∫ tl+1

tl

−1
τ

1

τ
dt

= −1

τ
,

kl−1,l =

∫ tl+1

tl−1

ψ′
l−1(t)ψ

′
l(t)dt

=

∫ tl

tl−1

−1
τ

1

τ
dt

= −1

τ
.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 13 of 21

pl,l+1 =

∫ tl+1

tl−1

ψl(t)ψl+1(t)dt

=

∫ tl

tl−1

(
t− tl−1

τ

)
·
(
t− tl
τ

)
dt

=
1

6
τ,

pl−1,l =

∫ tl+1

tl−1

ψl−1(t)ψl(t)dt

=

∫ tl+1

tl

(
tl+1 − t

τ

)
·
(
tl − t
τ

)
dt

=
1

6
τ.

(19)

Verification: in (19), write (17) for ψl+1 and ψl−1.
Using the same method, we calculate the coefficients of b to obtain:

bjm =

∫
Ω

∫ tm+1

tm−1

f(x)φj(x)ψm(t)dxdt

=

∫
Ω
f(x)φj(x)

∫ tm

tm−1

t− tm−1

τ
dxdt

+

∫
Ω
f(x)φj(x)

∫ tm+1

tm

tm+1 − t
τ

dxdt

≈ τ
∫
Ω
f (xj)φj(x)dx.

The system of linear equations (16) is obtained by substituting estimated coefficients
to (18):

M
(
uk+1 −2uk + uk−1

)
= τ2F k

− τ2A
(
1

6
uk−1 +

2

3
uk +

1

6
uk+1

)
, k = 1, . . . , N − 1.

(20)

To get an explicit approach, M is approximated using the lumped mass matrix ML

in space, the diagonal approximation obtained by taking the row sum of M , and mass
lumping in time by replacing the terms 1

6u
k−1+ 2

3u
k+ 1

6u
k+1 by uk and (20) is multiplied

by
(
ML

)−1
to produce an effective explicit formulation:

uk+1 = τ2F k + 2uk − τ2
(
ML

)−1
Kuk − uk−1, k = 1, . . . , N − 1. (21)



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 14 of 21

4. Numerical Experiments

Example 1. This experiment have been done by S. S. Almutairi, A. M. Saeed [11]
but we extended it to additional points to ensure the accuracy of the method.

∂2u

∂x2
+
∂2u

∂y2
= xy. (22)

Depending on the boundary conditions u(x, 0) = 0, u(x, 2) = 0, u(0, y) = 0, u(2, y) = 0.
With the exact solution of the PDE is

u =
∑
m≥1

∑
n≥1

−16
π2

(−1)n+m

nmλnm
sin(

nπx

2
) sin(

mπy

2
). (23)

Table 1: A comparison of the approximation solutions obtained using FDM and GFEM, alongside the analytical
solution of the given problem.

NODES FDM GFEM Analytical Solution

1 0 0 0
2 −0.0684455218855219 −0.066074304273125 −0.0679026913
3 −0.119349193976091 −0.122668821548822 −0.1226832533
4 −0.11477288471148 −0.117143703703704 −0.1176222268
5 −0.119349193976091 −0.122668821548822 −0.1226832533
6 −0.146756094276094 −0.14343667851387 -0.1473700856
7 −0.216687033077673 −0.221429225589226 −0.222990496
8 -0.216574276094276 -0.213255170322128 -0.2194532299
9 −0.213255170322128 −0.216574276094276 −0.2194532299
10 -0.240081885521885 -0.237711277602432 -0.2512047163
11 -0.267567407407407 -0.262825492159407 -0.2707333253
12 -0.262825492159407 -0.267567407407407 -0.2707333253
13 -0.273941548821549 -0.270622555204864 -0.2809825828
14 -0.329065589225589 -0.324323898288515 -0.3354518272

Table 1 presents the approximate solutions obtained using GFEM and FDM, alongside
the analytical solution of the proposed problem. The results indicate that GFEM performs
better than FDM, as its approximations are closer to the analytical solution. This is
because GFEM relies on a solid mathematical foundation, using weighted residuals and
carefully chosen basis functions to minimize errors. It’s especially effective for handling
complex geometries and localized phenomena, making it ideal for simulations that demand
precision. On the other hand, FDM is straightforward and computationally efficient, but
its accuracy depends heavily on maintaining a structured mesh. This can lead to larger
discretization errors, especially in cases where the problem involves irregular domains or
varying material properties. Because of this limitation, FDM may not be the best choice
when high precision is required. Simply put, GFEM stands out because of its adaptability



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 15 of 21

and accuracy, making it a better choice for solving differential equations in situations
where precision matters most.

Figure 4: Graph representing the numerical solution of the Poisson equation computed using FDM.

Figure 5: Graph representing the numerical solution of the Poisson equation computed using GFEM.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 16 of 21

Figure 6: Graph representing the analytical solution.

Figures 4, 5 and 6 clearly show the numerical solution of Poisson’s equation using
FDM, FEM, and analytical solution.

Example 2. Consider wave equation

∂2u

∂t2
= c2(

∂2u

∂x2
+
∂2u

∂y2
),

with boundary conditions

u(0, y, t) = u(a, y, t) = 0, 0 ≤ y ≤ b, t > 0,

u(x, 0, t) = u(x, b, t) = 0, 0 ≤ x ≤ a, t > 0.

And initial conditions

u(x, y, 0) = f(x, y), (x, y) ∈ R,

u(x, y, 0) = g(x, y), (x, y) ∈ R.

The general solution is

u(x, y, t) =
∞∑
n=1

∞∑
m=1

sin(
mπ

a
x) sin(

nπ

b
y)(Bmn cos(c

2

√
(
mπ

a
)2 + (

nπ

b
)2t)+B∗

mn sin(c
2

√
(
mπ

a
)2 + (

nπ

b
)2t)).

(24)



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 17 of 21

Table 2: Comparison results of the approximation solution FDM and analytical solution.

X Y Analytical Solution FDM

0.15789 0.15789 −0.091701 −0.095024
0.78947 1.1053 -3.117 -3.0759
0.94737 0.63158 -2.3209 -2.2537
1.1053 1.4211 -4.3357 -4.3143
1.5789 0.47368 -2.1181 -2.0426
1.7368 2.5263 -2.0561 -1.9891
2.0526 2.3684 -2.3209 -2.2537
2.3684 1.8947 -2.5588 -2.4973
2.8421 1.2632 -0.64947 -0.66045

The error of FDM are L1 error: 0.031623, L2 error: 0.041968 and L∞ error: 0.083343.

Figure 7: Comparison of the analytical solution and the FDM approximation solution for wave equation

Table 3: Comparison results of the approximation solution GFEM and analytical solution.

X Y Analytical Solution GFEM

0.15 1.95 -0.56165 -0.56281
0.225 0.075 -0.061833 -0.062814
0.9 2.85 -0.50533 -0.50546
1.2 2.7 -1.2571 -1.2247

1.575 0.075 -0.30975 -0.31918
1.95 2.85 -0.56165 -0.56281
2.025 0.825 -2.9906 -2.9721
2.325 1.425 -2.9808 -2.9531
2.775 1.125 -0.89784 -0.88465

The errors of FEM are L1 error: 2.722e-05, L2 error: 0.00099633 and L∞ error:



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 18 of 21

0.056977.

Figure 8: Comparison of the analytical solution and the GFEM approximation solution for wave equation

The GFEM exhibits lower error rates than the FDM, which suggests that the GFEM is
more reliable for solving these kinds of equations, according to tables 2 and 3. In response,
this study suggests using GFEM as a substitute for elliptic and hyperbolic type PDEs.
Also figures 7 and 8 clearly demonstrate that the GFEM solution agrees more closely with
the analytical solution than FDM.

Figure 9: Graph of the numerical solution of wave equation Using FDM.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 19 of 21

Figure 10: Graph of the numerical solution of wave equation Using GFEM.

Figure 11: Graph of the analytical solution of wave equation.

Figures 9, 10, and 11 display the solutions of the wave equation obtained using the
FDM, GFEM, and the analytical solution. The results show that the GFEM achieves
higher accuracy compared to the FDM when applied to the wave equation. Based on these
findings, we conclude that GFEM provides more precise and stable solutions, making it a
more effective choice for solving elliptic and hyperbolic PDEs.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 20 of 21

5. Conclusion

In this paper, Poisson’s equation and the Wave equation have been solved using the
GFEM. To assess its effectiveness relative to the FDM, several numerical experiments .
The results consistently indicate that GFEM achieves superior accuracy, particularly in
capturing fine wave structures and handling boundary conditions with greater stability.
This improved performance is most evident in scenarios involving complex geometries and
high-frequency wave propagation, where FDM tends to suffer from numerical dispersion
and phase errors. Based on these findings, we conclude that GFEM presents a reliable and
efficient alternative for solving both elliptic and hyperbolic partial differential equations.
Furthermore, given its flexibility in mesh adaptation and compatibility with higher-order
basis functions, GFEM holds strong potential for broader applications in engineering, fluid
dynamics, material science, and other computational fields that require robust and precise
numerical modeling.

Acknowledgements

The authors gratefully acknowledge Qassim University, represented by the Deanship of
Graduate Studies and Scientific Research, on the financial support for this research under
the number (QU-J-PG-2-2025- 53051) during the academic year 1446 AH / 2024 AD.

References

[1] A Quarteroni. Numerical approximation of partial differential equations. Berlin ; New
York : Springer-Verlag, 1994.

[2] G Smith. Numerical Solution of Partial Differential Equations: Finite Difference
Methods. Clarendon Press, 1985.

[3] Kh. Lotfy. Two temperature generalized magneto-thermoelastic interactions in an
elastic medium under three theories. Applied Mathematics and Computation, 227:871–
888, 2014.

[4] Kh. Lotfy, A.A. El-Bary, W. Hassan, A.R. Alharbi, and M.B. Almatrafi. Electromag-
netic and thomson effects during photothermal transport process of a rotator semi-
conductor medium under hydrostatic initial stress. Results in Physics, 16:102983,
2020.

[5] B Szabo and I babuska. Finite Element Analysis: Method, Verification and Valida-
tion. Wiley, 2021.

[6] T Hughes. The finite element method : linear static and dynamic finite element
analysis. Mineola, NY : Dover Publications, 2000.

[7] E Cheney and D Kincaid. Numerical Mathematics and Computing. Pacific Grove,
CA : Brooks/Cole Pub. Co, 1999.

[8] O Zienkiewicz. The finite element method. Singapore : McGraw-Hill, 1989.
[9] O Zienkiewicz, R Taylor, and J Zhu. The Finite Element Method: Its Basis and

Fundamentals. Elsevier Science, 2005.



Abdulkafi M. Saeed, Shahad S. Almutairi / Eur. J. Pure Appl. Math, 18 (3) (2025), 6524 21 of 21

[10] S Brenner and R Scott. The Mathematical Theory of Finite Element Methods.
Springer New York, 2007.

[11] A Shahad and S Abdulkafi. A comparative study of finite difference and galerkin
finite element methods for solving boundary value problems. European Journal of
Pure and Applied Mathematics, 18(2):5895, 2025.

[12] P Causon and P Mingham. Introductory Finite Difference Methods for PDEs. Book-
boon, 2010.

[13] U Yashkun and M Chandio. Finite difference method with dirichlet problems of 2d
laplace’s equation in elliptic domain. Pakistan Journal of Engineering Technology
and Science (PJETS), 6(2):136–144, 2016.

[14] S Abdulkafi. Improved rotated finite difference method for solving fractional elliptic
partial differential equations. American Scientific Research Journal for Engineering,
Technology, and Sciences (ASRJETS), 26(1):261–270, 2016.

[15] P Zhou. Finite Difference Method, pages 63–94. Springer Berlin Heidelberg, Berlin,
Heidelberg, 1993.

[16] R. Burden and J. Faires. Numerical Analysis. PWS-Kent Publishing Company, 1989.
[17] H Grosshans, S Gopireddy, R Humza, and E Gutheil. Modeling and simulation of

single particle and spray drying of pvp- and mannitol-water in hot air. Springer,
Cham, pages 309–339, 2016.

[18] J West and Linster. The Evolution of Fuzzy Rules in Two-Player Games. Southern
Economic Journal, 69(3):705–717, 2003.

[19] P Ciarlet. The Finite Element Method for Elliptic Problems. Studies in Mathematics
and its Applications. North Holland, 1978.

[20] G Michael B Wolfgang and R Rolf. Adaptive galerkin finite element methods for the
wave equation. Computational Methods in Applied Mathematics, 10(1):3–48, 2010.

[21] B Klaus-Jurgen H Seounghyun. A finite element method enriched for wave propaga-
tion problems. Computers and structures, 94(1):1–12, 2012.

[22] J Kim and K Bai. A finite element method for two-dimensional water-wave problems.
International journal for numerical methods in fluids, 30(1):105–122, 1999.