EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 16, No. 3, 2023, 1956-1969
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Multigrid Methods for The Solution of Nonlinear
Variational Inequalities

Nesba Nour El Houda1, Beggas Mohammed1, Belouafi Mohammed Essaid1,
Imtiaz Ahmad2, Hijaz Ahmad3,4,5,∗, Sameh Askar6

1 Faculty of Exact Sciences, Echahid Hamma Lakhdar University of El Oued, Operator
Theory, PDE and Applications Laboratory, Algeria
2 Institute of Informatics and Computing in Energy (IICE), Universiti Tenaga Nasional,
Kajang, Selangor, Malaysia
3 Section of Mathematics, International Telematic University Uninettuno, Corso Vittorio
Emanuele II, 39,00186 Roma, Italy
4 Operational Research Centre in Healthcare, Near East University, TRNC Mersin 10,
Nicosia, 99138, Turkey
5 Department of Computer Science and Mathematics, Lebanese American University, Beirut,
Lebanon
6 Department of Statistics and Operations Research, College of Science,
King Saud University, P.O. Box 2455, Riyadh 11451, Saudi Arabia

Abstract. In this research, we investigate the numerical solution of second member problems
that depend on the solution obtained through a multigrid method. Specifically, we focus on
the application of multigrid techniques for solving nonlinear variational inequalities. The main
objective is to establish the uniform convergence of the multigrid algorithm. To achieve this,
we employ elementary subdifferential calculus and draw insights from the convergence theory of
nonlinear multigrid methods.

2020 Mathematics Subject Classifications: 65M55, 65N30, 49N05, 65K15, 35J86

Key Words and Phrases: Multigrid method, nonlinear variational inequality, finite element,
iterative method, HJB equation

1. Introduction

The recent literature offers a wide range of diverse computational methods [6, 9, 20,
21, 24, 29, 31, 32, 34, 35] that are employed to solve complex real-world problems across

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v16i3.4835

Email addresses: ahmad.hijaz@uninettuno.it hijaz.ahmad@neu.edu.tr (H. Ahmad),
nour17nesba@gmail.com (N. N. E. Houda), beggasmr@yahoo.fr (B. Mohammed),
gogoo.said@gmail.com (B. M. Essaid), imtiaz.ahmad@uniten.edu.my (I. Ahmad),
saskar@ksu.edu.sa (S. Askar)

https://www.ejpam.com 1956 © 2023 EJPAM All rights reserved.



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1957

various fields of science and engineering [1, 2, 5, 8, 23, 33]. These methods have been de-
veloped and applied to tackle challenging problems and provide effective solutions in their
respective domains. Researchers have explored and utilized these computational tech-
niques to address a multitude of problems, advancing our understanding and facilitating
progress in numerous scientific and engineering disciplines [3, 4, 7, 27, 28, 30].

The numerical methods commonly employed to solve boundary problems typically
lead, after discretization, to the resolution of a system of large algebraic equations. These
numerical algorithms, including iteration methods such as Jacobi, Gauss-Seidel iteration,
and relaxation methods, are often used due to their customary nature. However, they
may exhibit slow convergence for small mesh sizes and complexity when adapted to gen-
eral elliptical problems. In contrast, multigrid methods offer a distinct advantage. They
are algorithms with a linear cost based on the number of discretization points, irrespec-
tive of the problem’s dimension. These methods are particularly efficient in solving both
linear and nonlinear partial differential equations (PDEs) and linear variational inequal-
ities (VIs) [16, 19]. Their linear complexity makes them a powerful tool for large-scale
problems, significantly reducing computational efforts while providing accurate solutions.

Multigrid techniques are widely acknowledged as rapid methods for solving various
types of variational equations and inequalities [22], especially for elliptic problems with
discretization leading to an M-matrix [15].
In the second section, we present an overview of nonlinear variational inequality (VI)
problems and their discretization using a conforming finite element method P1 [13]. Addi-
tionally, we introduce an algorithm that formulates the VI as stationary Hamilton-Jacobi-
Bellman equations, drawing inspiration from the Hoppe multigrid method [18, 25]. We
refer to this algorithm as the multigrid Hierarchy Jacobi (MGHJB) and provide the iter-
ation matrices associated with the algorithm.
Firstly, we present the original results concerning the approximation and smoothing prop-
erties in the L∞ norm. Subsequently, we demonstrate the uniform convergence of the
MGHJB algorithm. Finally, we apply the numerical methods to a specific application
where the operator is linear, and the second member is nonlinear, dependent on the solu-
tion. In this context, we implement the Gauss-Seidel method and the multigrid methods V
and W-cycle. The numerical experiments aim to assess the effectiveness and performance
of these methods in solving the given problem.

2. Multigrid method

2.1. Symbols and assumptions

Let Ω be an open set in RN with sufficiently regular bounds ∂Ω. For u, v ∈ H1(Ω), we
define second-order operators

A =
∑

1≤j,k≤N

ajk(x)
∂2

∂xj∂xk
+

N∑
k=1

bk(x)
∂

∂xk
+ a0(x).



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1958

The coefficients ajk(x), bk(x), and b0(x) are required to possess sufficient regularity, en-
suring:

akj(x) = ajk(x); b0(x) ≧ β > 0, x ∈ Ω,∑
1≤j,k≤N

akj(x)ξkξj ≧ α|ξ|2;
(
x ∈ Ω̄, ξ ∈ RN , α > 0

)
.

We also define the associated coercive continuous bilinear form:

a(u, v) =

∫
Ω

 ∑
1≤j,k≤N

ajk
∂u

∂xj

∂v

∂xk
+

N∑
k=1

bk
∂u

∂xk
v + a0(x)uv

 dx.

In addition, we consider the second member f is l-Lipschitzian such that:

f(u) ∈ L∞(Ω); f ≥ 0 and
l

β
< 1.

An obstacle ψ ∈W 2,∞, such that ψ ≥ 0.

2.2. Continuous problem

The objective is to find the solution u for the given variational inequalities:
Find u as the solution to:{

a(u, v − u) ≥ (f(u), v − u), ∀v ∈ H1(Ω)
u ≤ ψ; v ≤ ψ

(2.3)

It is established that this problem possesses a unique solution, as proven by the fixed point
theorem under the preceding assumption (see [10]).

2.3. Discretize

To build a multigrid loop, we consider a sequence of discretization steps, denoted as
0 < hk+1 < hk < 1, such that the grids are nested hk+1 =

hk
2 .

Next, we define Ωk = Ωhk
, Vk = Vhk

, and Ak = Ahk, and establish a series of uniform
regular triangulations denoted as {Tk, k ∈ N0} . For all Tk, we have the following estimates:

Ωk ⊂ Ωk+1 ⊂ Ω,
dist (∂Ωk, ∂Ω) ≤ c0h

2
k,

hkhk+1 ≤ c1.

We introduce Vhk
=

{
vhk

∈ C(Ω) ∩H1 So vhk
/T ∈ P1} , we simply write:

Vk =
{
vk ∈ C(Ω) ∩H1 this vk/r ∈ P1

}
.



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1959

The function of the usual basis φi
k, i ∈ (1, . . . .,m (hk)) is defined as: φi

k

(
xjk

)
= δij , x

j
k is

a node of the triangulation Tk.
We define the normal restriction operator as follows:

rkv(x) =

m(hk)∑
i=1

v
(
M i

k

)
φi
k(x). (2.4)

If we write Uk = Rmk , rk is a bijection from Uk to Vk .
In Uk, we define the scalar product:

< u, v >= h2k

m(hk)∑
i=1

uivi, ∥u∥k =< u, v >
1/2
k .

The maximum norm ∥ ·∥∞ (in Uk) and ∥ ·∥L∞ (in Vk) are equivalent standards, we denote
them ∥ · ∥∞. We have the following lemma (see [11]).
Lemma 1: There exists C1, C2 independent of k so:

∥rk(u)∥L∞ = ∥u∥L∞ ,∀u ∈ Uk.

C1∥v∥L∞ ≤ ∥r∗k(v)∥L∞ ≤ C2∥ v∥L∞ ,∀v ∈ Vk.
(2.5)

2.4. Discrete problem

In a natural progression, we introduce the discretization matrices Ak and the generic
coefficient matrices a

(
φk1, φs

k

)
, where φs, s = 1, 2, . . .m (hk) represent the customary

basic functions. With these definitions in place, we can now define the discrete problem
as follows:
Find uk ∈ Vk, which serves as the solution of:{

< Akuk, vk − uk >≥< fk(uk), vk − uk >, ∀vk ∈ Vk.
uk ≤ rkψ; vk ≤ rkψ.

(2.6)

We assume that the matrices Ak are M-matrices (see [12]).

2.5. Formulation in HJB:

The equivalence between the finite-dimensional variational inequality (2.3) and a for-
mulation in Hamilton-Jacobi-Bellman (HJB) can be readily observed (see [19]). We de-
scribe the numerical algorithm chosen to solve the stationary HJB equations.
In the classic framework, we recall some convergence results that will be instrumental in
establishing the convergence of the MGHJB algorithm described in the subsequent section.
Iterative diagram:
Step 1: Choose an initial vector u0k ∈ Rnk .

Step 2: Let u
(ν)
k ∈ Rnk , ν ≥ 0, and calculate u

(ν+1)
k ∈ Rnk as the solution of the equation:

Aν
ku

ν+1
k − Zν

k = 0, (2.7)



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1960

such that:
Zν
k = F ν

k (u
ν
k),

where:

Aν
k,i =

{
Ak,i (uk) if Ak,iu

ν
k,i − Zk,i > uνk,i − ψk,i.

uk,i if 1 ≤ i ≤ N
. (2.8a)

Zν
k,i =

{
Zk,i if Ak,iu

ν
k,i − Zk,i > uνk,i − ψk,i

uk,i if 1 ≤ i ≤ N
(2.8b)

Let u∗k be the unique solution of the discrete HJB equation

max
1≤i≤N

(
Ak,iu

∗
k − Zk,i, u

∗
k,i − ψk,i

)
= 0. (2.9)

We shall state the following theorem and present our problem formulated based on the
Hamilton-Jacobi-Bellman (HJB) equation, which is an adaptation of Hoppe’s work [19].
Theorem 2: Let uνk be the iteration defined, and it satisfies the HJB equation. Addi-
tionally, we assume that Ak is continuously differentiable. Under these conditions, the
sequence (uνk)ν≥0 converges and approaches u∗k.
Before proceeding to present the results, it is pertinent to recall the following theorem.
Theorem 3: ( [10] , [14] ) Under the previous assumptions and notation, we have:

∥u− u∗k∥∞ ≤ Ch2k |Loghk|
2 ∥f(u)∥∞. (2.10)

2.6. Description of the multigrid method for VIs

Choosing an iteration uνk, ν > 0 for the multigrid method, we obtain ūνk, by applying
an iterative method to solve the system (2.7) by α, expressing for

ūνk = Sα
k (uνk) (2.11)

where Sk represents the iteration or smoothing operator, and α denotes the number of
iterations performed. We denote the solution to (2.7) by u∗k. The error setting e

ν
k = ūνk−u∗k,

and the residual d
(ν)
k = Zν

k −Aν
kū

ν
k, we can write the equation (2.7) as

Aν
k (ū

ν
k + eνk) = Zν

k .

This leads to the residual equation

Aν
ke

ν
k = Zν

k −Aν
kū

ν
k = d

(ν)
k .

On the fine grid, after relaxation on Aν
kū

ν
k = Zν

k , the error will exhibit smooth behavior.
However, on the coarse grid, the error appears to be more oscillatory, leading to efficient
relaxation. To completely determine eνk, we need to compute eνk−1 at the (k − 1) level, as
it serves as the solution for the coarse grid system.

Aν
k−1e

ν
k−1 = d

(ν)
k−1. (2.12)



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1961

We can interpret eνk−1 (and its corresponding operators Aν
k−1, d

(ν)
k−1) as the approximation

operator at level k − 1, while eνk (and its corresponding operators Aν
k, d

(ν)
k ) represent the

approximation operator at level k. Additionally, we have the restriction operator Rk and
its inverse Pk.
Therefore, we identify an improved iteration at the k level

uν+1
k = ūνk + Pk

(
eνk−1

)
. (2.13)

Due to the nestedness, we use the well-defined identity operator

Π : Vk−1 −→ Vk,

Πv = v,

to define prolongation and restriction operators, that is:

Pk = r−1
k rk−1, Rk = Pt

k. (2.14)

Note 4: The preceding algorithm describes a loop of two mesh iterations to solve (2.7)
for two mesh levels Ωk−1. It is clear that the coarse grid system (2.12) has the same form
as the system (2.7). Therefore, we can approximate the solution of the system (2.12) by
recursively doing two-grid iterations at all grid levels {Ωk, k = 0, . . . ,mk}.

2.7. Matrix associated with the MGHJB algorithm:

The matrix of iterations for the two-grid method with α1 presmoothing and α2 post-
smoothing at level k is expressed as follows:

TGk (α1, α2) = Sα2
k

(
(Aν

k)
−1 − Pk

(
Aν

k−1

)−1Rk

)
(Aν

k)S
α1
k . (2.15)

Theorem 5: ([26]) The multigrid approach constitutes a linear iterative method, with
the iteration matrix denoted as MGk.

MG0 = 0,

MGk = Sα2
k

(
Ik − Pk (Ik −MGk−1)

(
Aν

k−1

)−1Rk

)
(Aν

k)S
α1
k ,

= TGk + Sα2
k PkMGk−1

(
Aν

k−1

)−1Rk (Aν
k)S

α1
k , k = 1, 2, ..

(2.16)

3. Convergence analysis of multigrid method in L∞ norm

In this section, we provide a comprehensive convergence analysis of the multigrid algo-
rithm. The algorithm was previously described using the maximum norm in the preceding
section.



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1962

3.1. Approximation property

Theorem 6: ([17]) The matrix Υk =
[
(Aν

k)
−1 − Pk

(
Aν

k−1

)−1Rk

]
has the approxi-

mation property
∥Υk∥∞ ≤ Ch2k |lnhk|

2 . (3.1)

Proof. The proof of the approximation property was introduced by Arnold in [25],
relying on Theorem (2).

3.2. Smoothing property

To demonstrate the smoothness property, we decompose Aν
k = Ek − Nk and use the

following assumptions:

Ek is regular and
∥∥E−1

k Nk

∥∥
∞ ≤ 1, for all k. (3.2)

∥Ek∥∞ ≤ Ch−2
k , for all k, with C independent of k. (3.3)

For the smoothing process, we employ a relaxation method with an iterative matrix.

Sk = Ik − ωE−1
k Nk, ω ∈ (0, 1).

In the following theorem, Arnold Reusken applied a concept that is relevant to our work.
The theorem’s application is appropriate as it is closely linked to the operator used in our
study.
Theorem 7: ( [25] ) Under the given assumptions and notations, it can be inferred that
there exists a constant C, which is independent of both k and α. With this constant, the
following smoothness properties can be established:

∥(Aν
k)S

α
k ∥∞ ≤ C

1√
α
h−2
k . (3.4)

By passing to the norm in (2.13), and taking into account (3.1) and (3.4), we have to
prove the the following stability limit:

∃Cs : ∥Sα
k ∥∞ ≤ Cs, for all k and α. (3.5)

The convergence analysis is performed based on the split of the following two lattices: the
iterate matrix with α2 = 0.

∥TGk(α1, 0)∥∞ =
∥∥∥((Aν

k)
−1 − Pk

(
Aν

k−1

)−1Rk

)
(Aν

k)S
α1
k

∥∥∥
∞

≤
∥∥∥((Aν

k)
−1 − Pk

(
Aν

k−1

)−1Rk

)∥∥∥
∞

∥∥(Aν
k)S

α1
k

∥∥
∞ .

Typically, a hierarchy of more than two grids is selected. In this scenario, iterative matrices
(2.16) can be defined using recursion of iterative matrices (2.15) for all levels k. If we



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1963

assume that (3.5) holds, then L∞-convergence results can be readily deduced from the
previous findings.
Theorem 8: ([26]) Consider a multigrid method for a given iterative matrix (2.16). Then
under the previous assumption, for the parameter value α2 = 0, α1 = α > 0, τ ≥ 2. For
each ζ ∈ (0, 1) there is α, α∗such that for all α ≥ α∗

∥MGk∥∞ ≤ ζ, k = 0, 1, ... (3.7)

hold.
Proof. If the approximation and smoothness properties are combined with (3.5), then

we can apply the same parameters as in [26].

The following theorem represents the main result of our work.
Theorem 9: Under the previous assumptions and notations the iterated uνk, ν ≥ 0 for
two meshes k and k − 1 satisfy:

∥∥uν+1
k − u∗k

∥∥
∞ ≤

(
C√
α
|Loghk|

2

)
∥uνk − u∗k∥∞ . (3.6)

Proof. We have∥∥uν+1
k − u∗k

∥∥
∞ =

∥∥∥((Ik − Pk (Ik −MGk−1)
(
Aν

k−1

)−1Rk

)
(Aν

k)S
α1
k

)
(uνk − u∗k)

∥∥∥
∞

≤
∥∥∥(Ik − Pk (Ik −MGk−1)

(
Aν

k−1

)−1Rk

)∥∥∥
∞

∥∥Aν
kS

α1
k

∥∥
∞ ∥uνk − u∗k∥∞

≤
(
C2√
α
h−2
k

)(
C1h

2
k |log hk|

2
)
∥(uνk − u∗k)∥∞

≤
(
C1C2√
α

)
|log hk|2 ∥uνk − u∗k∥∞

4. Numerical simulation

In this section, we present a numerical example of a nonlinear variational inequality.
To apply the proposed method to the example, we assume that the data of the problem
is sufficiently smooth. We then employ Bellman’s principle of dynamic programming to
tackle the problem and proceed with solving (2.3), as discussed previously, using the fol-
lowing specified data:

Au ≤ f(u), in Ω =
{
(x, y) | x2 + y2 ≤ 1

}
,

⟨Au− f(u), u− ψ⟩ = 0,

u ≤ ψ,

u = 0, in ∂Ω,

(4.1)



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1964

where
Au = −∆u,

f(u) = cosu,

ψ = 0.

We are confine ourselves to the FEM discretization with a uniform triangulation and P1
nested finite element function spaces. For the discretization of the domain, we have used
the PDE toolbox in MATLAB (R2017b) to generate the mesh and then the multigrid
FEM can be used to efficiently solve 4.1 as mentioned above.

Figure 1: Domain of our problem with 2048 triangle and 1089 nodes.

This numerical example is conducted to demonstrate the high efficiency of the multi-
grid method. In this study, we utilize the Gauss-Seidel method for pre/post smoothing
within the multigrid code. With respect to recursion in the multigrid method, the recursive
algorithm is terminated when the degrees of freedom, represented by the number of inte-
rior grid points, become less than 5. Figure 2 illustrates the convergence behavior of the

Figure 2: Comparison between the convergence behavior of Multigrid method and Gauss-Seidel methods.

multigrid solver (green and red curves for the maximum norm of the residual of multigrid
(V and W cycle)) with respect to the number of iterations performed. For comparison, the
convergence behavior of Gauss-Seidel ( blue curves) are included. The multigrid V -cycle
is carried out on the finest grid with 1089 nodes and 4 nodes on the coarsest one, and then
we have applied thethe Matlab backslash-operator and Gauss-Seidel on this finest grid to



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1965

get the solutions in figures 3-6.
Norm of residual obtained by Gauss Seidel method after 100 iteratons := 0.1303.

Norm of residual obtained by multigrid V-cycle after 100 iterations := 3.7937e− 11.

Norm of residual obtained by multigrid W-cycle after 100 iterations := 1.6209e− 14.

Notting that, if we perform more than 20 iterations, the multigrid solution is better than
the Matlab backslash-operator (MBO) solution.

Norm of residual obtained by Gauss Seidel method after 100 iteratons := 0.1303.

Norm of residual obtained by multigrid V-cycle after 20 iterations := 0.0049.

Norm of residual obtained by multigrid W-cycle after 20 iterations := 4.5777e− 5.

Figure 3: Solution of the problem 4.1 on fine grid with 1089 DOFs using Matlab backslash operator solver after
100 iterations.

Figure 4: Solution of the problem 4.1 on fine grid with 1089 DOFs using Gauss Seidel Method after 100
iterations.

5. Conclusion

In this study, we utilize the algebraic multigrid method, specifically the efficient itera-
tive solutions for discretizing elliptic variational inequalities, to handle the discretization



H. Ahmad et al. / Eur. J. Pure Appl. Math, 16 (3) (2023), 1956-1969 1966

Figure 5: Solution of the problem 4.1 on fine grid with 1089 DOFs using Multigrid Method V-cycle after 100
iterations.

Figure 6: Solution of the problem 4.1 on fine grid with 1089 DOFs using Multigrid Method W-cycle after 100
iterations.

of a loop domain using adaptive finite element approximation. After discretization, we
implement the multigrid method to efficiently solve the discrete problems. We investi-
gate the uniform convergence of our approach and demonstrate that the multigrid method
exhibits a significant reduction in the number of iterations compared to the maximum
norm. In our numerical experiments, we present an example of a variational inequality.
The results indicate that the Gauss-Seidel method does not perform well even after a
large number of iterations. In contrast, the multigrid method, with its debug function
that reduces high-frequency errors through relaxation and maps low-frequency errors to
the coarse grid for reduction, achieves convergence in a small number of iterations. We
acknowledge the potential for many extensions of these techniques. An interesting future
direction could involve applying a parallel full multigrid method to solve unconstrained
elliptical inequalities. This approach may further enhance the efficiency and scalability of
our numerical solution for a broader range of problems.

Acknowledgements

This project is funded by King Saud University, Riyadh, Saudi Arabia. Research
Supporting Project number (RSP2023R167), King Saud University, Riyadh, Saudi Arabia.



REFERENCES 1967

References

[1] M Adel, NH Sweilam, MMKhader, SM Ahmed, H Ahmad, and T Botmart. Numerical
simulation using the non-standard weighted average FDM for 2Dim variable-order
Cable equation. Results in Physics, 39:105682, 2022.

[2] H Ahmad AE Abouelregal and S-W Yao. Functionally graded piezoelectric medium
exposed to a movable heat flow based on a heat equation with a memory-dependent
derivative. Materials, 13(18):3953, 2020.

[3] H Ahmad, , and TA Khan. Variational iteration algorithm I with an auxiliary pa-
rameter for the solution of differential equations of motion for simple and damped
mass–spring systems. Noise & Vibration Worldwide, 51(1-2):12–20, 2020.

[4] H Ahmad, TA Khan, and S-W Yao. An efficient approach for the numerical solution
of fifth-order KdV equations. Open Mathematics, 18(1):738–748, 2020.

[5] I Ahmad, AR Seadawy, H Ahmad, P Thounthong, and F Wang. Numerical study of
multi-dimensional hyperbolic telegraph equations arising in nuclear material science
via an efficient local meshless method. International Journal of Nonlinear Sciences
and Numerical Simulation, 23(1):115–122, 2022.

[6] M Ahsan, AA Khan, S Dinibutun, I Ahmad, H Ahmad, N Jarasthitikulchai, and
W Sudsutad. The Haar wavelets based numerical solution of Reccati equation with
integral boundary condition. Thermal Science, 27(Spec. issue 1):93–100, 2023.

[7] A Akgül and H Ahmad. Reproducing kernel method for Fangzhu’s oscillator for water
collection from air. Mathematical Methods in the Applied Sciences, 2020.

[8] R Alharbi, R Jan, S Alyobi, A Yousif, and Z Khan. Mathematical modeling and
stability analysis of the dynamics of monkeypox via fractional-calculus. Fractals,
30(10):2240266, 2022.

[9] B Almutairi, I Ahmad, B Almohsen, H Ahmad, and DU Ozsahin. Numerical sim-
ulations of time-fractional PDEs arising in mathematics and physics using the local
Meshless differential quadrature method. Thermal Science, 27(Spec. issue 1):263–272,
2023.

[10] M Boulbrachene and M Haiour. The finite element approximation of Hamilton-Jacobi-
Bellman equations. Computers & Mathematics with Applications, 41(7-8):993–1007,
2001.

[11] F Brezzi and LA Caffarelli. Convergence of the discrete free boundaries for finite
element approximations. RAIRO. Analyse numérique, 17(4):385–395, 1983.

[12] PG Ciarlet and P-A Raviart. Maximum principle and uniform convergence for the
finite element method. Computer methods in applied mechanics and engineering,
2(1):17–31, 1973.



REFERENCES 1968

[13] P. Cortey-Dumont. On the finite element approximation in the l∞ norm of variational
inequalities with nonlinear operators. Numerische Mathematik, 47:45–57, 1985.

[14] P Cortey-Dumont and S I’Analyse. Numérique des Equations de Hamilton- Jacobi-
Bellman. Mathematical Methods in the Applied Sciences, 9:198–209, 1987.

[15] W Hackbusch. Multi-grid methods and applications., volume 4. Springer Science &
Business Media, 2013.

[16] W Hackbusch and HD Mittelmann. On multi-grid methods for variational inequali-
ties. Numerische Mathematik, 42:65–76, 1983.

[17] M HAIOUR. Etude de la convergence uniforme de la méthode multigrilles appliquées
aux problèmes frontières libres. PhD thesis, Université de Annaba-Badji Mokhtar,
2004.

[18] R HW Hoppe. Multi-grid methods for Hamilton-Jacobi-Bellman equations. Nu-
merische Mathematik, 49:239–254, 1986.

[19] R HW Hoppe. Multigrid algorithms for variational inequalities. SIAM journal on
numerical analysis, 24(5):1046–1065, 1987.

[20] H Irshad, M Shakeel, I Ahmad, H Ahmad, C Tearnbucha, and W Sudsutad. Sim-
ulation of generalized time fractional Gardner equation utilizing in plasma physics
for non-linear propagation of ion-acoustic waves. Thermal Science, 27(Spec. issue
1):121–128, 2023.

[21] R Jan, A Khan, S Boulaaras, and SA Zubair. Dynamical behaviour and chaotic
phenomena of HIV infection through fractional calculus. Discrete Dynamics in Nature
and Society, 2022, 2022.

[22] D Kinderlehrer and G Stampacchia. An introduction to variational inequalities and
their applications. SIAM, 2000.

[23] J-F Li, I Ahmad, H Ahmad, DShah, Y-M Chu, P Thounthong, and M Ayaz. Numer-
ical solution of two-term time-fractional PDE models arising in mathematical physics
using local meshless method. Open Physics, 18(1):1063–1072, 2020.

[24] H Ahmad M Adel, ME Ramadan and T Botmart. Sobolev-type nonlinear Hilfer
fractional stochastic differential equations with noninstantaneous impulsive. AIMS
Mathematics, 7(11):20105–20125, 2022.

[25] A. Reusken. On maximum norm convergcnce of multigrids methods for elliptic bound-
ary value problems. SIAM Journal on Numerical Analysis, 29(6):1569–1578, 1992.

[26] A Reusken. Introduction to multigrid methods for elliptic boundary value problems.
Inst. für Geometrie und Praktische Mathematik, 2008.



REFERENCES 1969

[27] NA Shah, I Ahmad, O Bazighifan, AE Abouelregal, and H Ahmad. Multistage
optimal homotopy asymptotic method for the nonlinear Riccati ordinary differential
equation in nonlinear physics. Applied Mathematics, 14(6):1009–1016, 2020.

[28] M Shakeel, I Hussain, H Ahmad, I Ahmad, P Thounthong, Y-F Zhang, and Ying-
Fang. Meshless technique for the solution of time-fractional partial differential equa-
tions having real-world applications. Journal of Function Spaces, 2020:1–17, 2020.

[29] M Shakeel, MN Khan, I Ahmad, H Ahmad, N Jarasthitikulchai, and W Sudsutad.
Local meshless collocation scheme for numerical simulation of space fractional PDE.
Thermal Science, 27(Spec. issue 1):101–109, 2023.

[30] TA Sulaiman, A Yusuf, S Abdel-Khalek, M Bayram, and H Ahmad. Nonautonomous
complex wave solutions to the (2+ 1)-dimensional variable-coefficients nonlinear Chi-
ral Schrödinger equation. Results in Physics, 19:103604, 2020.

[31] R Jan T-Q TangZ Shah and E Alzahrani. Modeling the dynamics of tumor–
immune cells interactions via fractional calculus. The European Physical Journal
Plus, 137(3):367, 2022.

[32] T-Q Tang, Z Shah, R Jan, W Deebani, and M Shutaywi Meshal. A robust study to
conceptualize the interactions of CD4+ T-cells and human immunodeficiency virus
via fractional-calculus. Physica Scripta, 96(12):125231, 2021.

[33] F Wang, I Ahmad, H Ahmad, MD Alsulami, KS Alimgeer, C Cesarano, and TA Nofal.
Meshless method based on RBFs for solving three-dimensional multi-term time frac-
tional PDEs arising in engineering phenomenons. Journal of King Saud University-
Science, 33(8):101604, 2021.

[34] F Wang, NA Ali, I Ahmad, H Ahmad, KM Alam, and P Thounthong. Solution
of Burgers’ equation appears in fluid mechanics by multistage optimal homotopy
asymptotic method. Thermal Science, 26(1 Part B):815–821, 2022.

[35] F Wang, E Hou, I Ahmad, H Ahmad, and Y Gu. An efficient meshless method for
hyperbolic telegraph equations in (1+ 1) dimensions. CMES-Computer Modeling in
Engineering and Sciences, 128(2):687–98, 2021.