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

A Hybrid Approach for Approximating a Smoking
Epidemic Model via Caputos Derivative, a Novel

Integral Transform, and Artificial Neural Networks
Rachid Belgacem1, Ahmed Bokhari1, Abdelkader Benali1, Ibrahim Alraddadi2,∗,
Hijaz Ahmad3, Waleed Mohammed Abdelfattah4,5

1 Laboratory of Mathematics and its Applications LMA, Hassiba Benbouali University of
Chlef, Algeria
2 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Madinah,
Saudi Arabia
3 Near East University, Operational Research Center in Healthcare, Near East Boulevard,
PC: 99138 Nicosia/Mersin 10, Turkey
4 College of Engineering, University of Business and Technology, Jeddah 23435,
Saudi Arabia
5 Department of Engineering Mathematics and Physics, Faculty of Engineering,
Zagazig University, P.O. 44519, Egypt

Abstract. This work focuses on obtaining approximate analytical solutions for a fractional-order
smoking epidemic model formulated with the Caputo derivative. The model divides the population
into five compartmentspotential smokers, current smokers, occasional smokers, permanent smokers,
and temporary quittersallowing the fractional framework to capture the long-term memory effects
in smoking dynamics. Transition rates between these compartments are described by a set of
parameters that reflect realistic behavioral changes. To solve the non linear fractional differential
system efficiently, we propose a new hybrid computational strategy that combines the Junaid
integral transform TN G with the Adomian Decomposition Method (ADM) and an Artificial Neural
Network (ANN). The combined TN GADM stage ensures rapid convergence and accurate series
solutions, while the ANN improves predictive performance by learning directly from the system
dynamics. Numerical simulations validate the effectiveness and computational efficiency of the
proposed approach, demonstrating its suitability for modelling memory-dependent epidemiological
processes.
2020 Mathematics Subject Classifications: 92D30, 92D25, 26A33, 65R10, 65L05
Key Words and Phrases: Fractional smoking epidemic model, Caputo fractional derivative,
Junaid transform, Adomian decomposition method

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6267

Email addresses: R.belgacem@univ-chlef.dz (R. Belgacem),
bokhari.ahmed@ymail.com (A. Bokhari), benali4848@gmail.com (A. Benali),
Ialraddadi@iu.edu.sa (I. Alraddadi), hijaz.ahmad@neu.edu.tr (H. Ahmad),
w.abdelfattah@ubt.edu.sa (W. M. Abdelfattah)

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



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 2 of 26

1. Introduction

In recent years, integral transforms have emerged as powerful mathematical tools, at-
tracting increasing attention from researchers across various scientific fields. Their effec-
tiveness in solving a wide range of linear equations including ordinary and partial differen-
tial equations (ODEs and PDEs), integral equations, and fractional differential equations
(FDEs) lies in their ability to convert complex differential problems into simpler algebraic
forms. This transformation simplifies the solution process and accelerates computations.
As a result, integral transforms have found widespread application across various disci-
plines, fostering the development of new methods and enhancing existing ones. However,
the inversion of these transforms remains a critical step in obtaining the final solution.

Recent advancements in integral transforms have revolutionized the computation of
PDEs, enabling both precise and approximate solutions [1]. Their inherent capability to
map functions between different domains while preserving key properties ensures efficient
and accurate outcomes, making them indispensable in applied mathematics. Moreover,
combining integral transforms with other techniques can effectively tackle the nonlinear
components of equations, leading to more efficient solutions and improved accuracy [2].
One such technique is the Adomian Decomposition Method (ADM)[3], which has proven
highly effective for solving complex ODEs, PDEs, and FDEs [4].

Over the past two decades, a variety of integral Laplace-type transformations have
been developed, including the Sumudu [5], Elzaki [6], Natural [7], Aboodh [8], Mohand
[9], Sawi [10], Shehu [11–14], Kamal[15] and Jafari [16, 17] transformations. More recently,
hybrid approaches involving neural networks and integral transforms have been explored
to approximate solutions of complex fractional models.

The successful application of hybrid methodologies to various fractional models pro-
vides strong justification for our approach. For instance, the generalized exponential ra-
tional function method has demonstrated remarkable efficacy in obtaining optical soliton
solutions for dual-mode time-fractional nonlinear Schrödinger equations [18], while inno-
vative approaches have been developed for solving (2+1)-dimensional generalized KdV
equations [19]. These hybrid analytical-numerical techniques have shown great potential
for handling complex nonlinear systems similar to our smoking epidemic model.

Similarly, the integration of neural networks with fractional calculus has emerged as
a powerful paradigm. Recent work on modeling and neural network approximation of
asymptotic behavior for delta fractional difference equations with Mittag-Leffler kernels
[20] demonstrates how machine learning can enhance traditional analytical methods. This
aligns with our approach of combining the TN G-transform with artificial neural networks
to improve accuracy and computational efficiency.

Our methodology also draws inspiration from finite difference approaches for fractional
models [21] and analytical algebraic methods for nonlinear fractional differential equations
[22]. These established techniques provide a solid foundation for our hybrid TN G-ADM-
ANN framework, adapting proven methods to the specific challenges of smoking epidemic
modelling.

In the reference [23], the author proposed a new generalized integral transform, which



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 3 of 26

will be referred to as such TN G-transform throughout this paper. This transform was used
in conjunction with the ADM method to solve non-linear PDEs, including the gas dynamic
equation and the system of coupled Burgers’ equations.

Motivation and Originality

This work makes several original contributions that advance both the theoretical and
computational treatment of FDEs. First, we present the explicit derivation of the specific
integral transform TN G for the CF derivative, establishing a rigorous analytical framework
that enables the systematic conversion of FDEs into an equivalent algebraic form. We fur-
ther prove several recent results regarding its operational properties, which significantly
enrich the transform’s analytical capabilities. Second, we develop a new hybrid technique,
denoted as TN G-ADM, which synergistically combines the TN G-transform with the ADM
method. This approach produces series solutions with rapid convergence and incorpo-
rates auxiliary parameters to fine-tune and optimize the convergence process. Finally, we
integrate an artificial neural network (ANN) algorithm to enhance the approximation ac-
curacy for the fractional smoking epidemic model, significantly improving the robustness,
precision, and applicability of the proposed methodology to real-world problems.

The novelty of our work lies in the unique combination of these established tech-
niquesintegral transforms, decomposition methods, and neural networksspecifically adapted
for Caputo fractional-order epidemiological models. While each component has been in-
dividually validated in previous studies, their integration into a unified framework for
smoking dynamics represents a significant advancement in computational epidemiology
[24, 25].

This manuscript is structured as follows. Section 2 reviews the fundamental concepts of
fractional calculus and the essential properties of the proposed general integral transform
(TN G). Section 3 presents the main theoretical results of the TN G transform, including
explicit formulations for the Caputo fractional derivative (CFD), the RiemannLiouville
fractional derivative (RLFD), and the RL fractional integral (RLFI), along with its key
operational properties. Section 4 addresses the mathematical modelling of the fractional
smoking epidemic model and highlights its principal analytical properties. Sections 5 and
6 present a hybrid TN G-ADM method for solving the fractional model, further enhanced
by integration with an ANN to improve accuracy and robustness. Finally, the study
highlights the method’s advantages and outlines directions for future work.

2. Preliminary

This first section presents fundamental definitions in fractional calculus and highlights
key properties of the Junaid transform, which will serve as the basis for the analyses that
follow.

Definition 1. For a given function V(t), the RLFI of order α > 0, is defined as [26]:

RLIα
0,t V(t) = 1

Γ(α)

∫ t

0

V(τ)
(t− τ)1−α

dτ = 1
Γ(α) tα−1 ?V(t), t > 0. (1)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 4 of 26

The RLFD of order α > 0, is given by:

RLDα
0,t V(t) = dm

dtm

(
RLIm−α

0,t V(t)
)

= dm

dtm

∫ t

0

V(x)
Γ(m−α)(t−x)m−α−1 dx, (2)

where m−1 < α ≤ m, m ∈ N, and Γ(.) denotes the Gamma function.

Definition 2. For a given function V(t), the CFD of order α > 0, is described as: [27]:

CDα
0,tV (t) = 1

Γ(m−α)

∫ t

a

V(m) (τ)dτ

(t− τ)1−m+α , 0 < m−1 < α < m,m ∈ N, (3)

= dm

tm
v (t) , α = m.

Lemma 1. If V ∈ ACm [a,b] or V (t) ∈ Cm [a,b], then [26]:

(
RLIα

0,t
CDα

0,t

)
V (t) = V (t)−

m−1∑
k=0

V(k) (0)
k! tk, (4)

where m = [α]+1 and [α] represents the integer part of arbitrary α > 0.
The Junaid integral transform of the integrable function V (t) was recently introduced

by Junaid [23] in 2023, by

TN G [V (t)] = M(s)
∫ ∞

0
V (K (s) t)exp−B(s)t dt, (5)

where M(s) ,B (s) and K (s) are positive real functions and t ≥ 0 such that M(s) 6= 0.
This TN G transform can be easily implemented directly to an appropriate problem by

specifically selecting M(s), B (s) and K (s). In Table 1, we mention the TN G transform of
some basic function.

V(t) TN G [V(t)]

c c
M(s)
B(s) , c ∈ R

t
M(s)K (s)

B2(s)
tn n!M(s)rn(s)

Bn+1(s)
sin t

M(s)K (s)
B2(s)+ r2(s)

et M(s)
B(s)−K (s)

Table 1: TN G transform for some basic functions.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 5 of 26

Theorem 1. The TN G transform of nth derivative V(n)(t) of the function V(t), is expressed
as [23]:

TN G
[
V(n) (t)

]
=

( B (s)
K (s)

)n

TN G [V (t)]−M(s)
n−1∑
i=0

Bn−1−i (s)
Kn−i (s) V(i) (0) ,∀n ∈ N. (6)

The next section presents new results about this generalization transform, as a com-
plement to the results presented in [23].

3. Main results of TN G transform

Property 1. (Convolution) Let TN G [V1 (t)] and TN G [V2 (t)] are TN G transforms of V1(t)
and V2(t), respectively. Then:

TN G [V1 ?V2] =
∫ ∞

0
V1 (t)V2 (t−z)dz = K (s)

M(s)TN G [V1 (t)]TN G [V2 (t)] . (7)

Proof. We have:
V1 ?V2 =

∫ ∞

0
V1 (z)V2 (t−z)dz,

using the TN G transform and the Leibniz theorem, we get:

TN G [V1 ?V2] = TN G

[∫ ∞

0
V1 (z)V2 (t−z)dz

]
,

= M(s)
∫ ∞

0

[∫ ∞

0
V1 (K (s)z)V2 (K (s)(t−z))K (s)dz

]
exp−B(s)t dt,

= M(s)K (s)
∫ ∞

0
V1 (K (s)z)

[∫ ∞

0
V2 (K (s)(t−z))exp−B(s)t dt

]
dz.

By setting τ = t−z, we get:

TN G [V1 ?V2] = M(s)K (s)
∫ ∞

0
V1 (K (s)z)

[∫ ∞

0
V2 (τ)exp−B(s)(τ+z) dτ

]
dz,

= K (s)
M(s)

[
M(s)

∫ ∞

0
V1 (K (s)z)exp−B(s)z dz

][
M(s)

∫ ∞

0
V2 (τ)exp−B(s)τ dτ

]
,

= K (s)
M(s)TN G [V1 (t)]TN G [V2 (t)] .

Property 2. When V (t) = tx−1 in formula of TN G transform 5, then:

TN G
[
tx−1

]
= Γ(x)M(s)Kx−1 (s)

Bx (s) . (8)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 6 of 26

Proof. From the formula 5, we obtain:

TN G
[
tx−1

]
= M(s)

∫ ∞

0
(K (s) t)x−1 exp−B(s)t dt,

= M(s)Kx−1 (s)
∫ ∞

0
tx−1 exp−B(s)t dt.

If we set ξ = B (s) t

(
t = ξ

B (s)

)
, then

TN G
[
tx−1

]
= M(s)Kx−1 (s)

Bx (s)

∫ ∞

0
ξx−1 exp−ξ dξ,

= M(s)Kx−1 (s)
Bx (s) Γ(x) .

Lemma 2. The TN G transform of the fractional integral RLIα
0,t of V, is formulated as:

TN G
[(

RLIα
0 V

)
(t)

]
= Kα (s)

Bα (s) TN G [V (t)] . (9)

Proof.
Applying TN G transform and according to properties 1 and 2, we have:

TN G
[

RLIα
0,tV (t)

]
= TN G

[ 1
Γ(α) tα−1 ?V (t)

]
,

= K (s)
M(s)

M(s)rα−1 (s)
Bα (s) TN G [V (t)] ,

= Kα (s)
Bα (s) TN G [V (t)] .

Lemma 3. The TN G transform of the RLDα
0 of V, is given as follows:

TN G
{

RLDα
0,tV (t) ,s

}
=

( B (s)
K (s)

)α

TN G [V (t)]−M(s)
m−1∑
k=0

Bm−1−k (s)
Km−k (s)

[
RLDα−k−m

0,t V (t)
]

t=0
,

(10)

where m−1 < α ≤ m.

Proof. Applying the TN G transform to both sides of equation (2), followed by the use
of Lemma (2) and Theorem (1), yields:

TN G
[

RLDα
0,tV (t)

]
= TN G

[
dm

dtm
Im−α

0,t V (t)
]

,



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 7 of 26

=
( B (s)

K (s)

)m

TN G
[
Im−α

0,t V (t)
]

−M(s)
m−1∑
k=0

Bm−1−k (s)
Km−k (s)

[(
d

dt

)(k)
Im−α

0,t V (t)
]

t=0
,

=
( B (s)

K (s)

)α

TN G [V (t)]−M(s)
m−1∑
k=0

Bm−1−k (s)
Km−k (s)

[
RLDα−k−m

0,t V (t)
]

t=0
.

Hence, the desired result (10) is established.

Theorem 2. If V ∈ ACm(a,b) for any b > a and of exponential order. Then, the TN G
transform of the CDα

0 of V is given by:

TN G
[

CDα
0 V (t)

]
= Bα (s)

Kα (s)TN G [V (t)]−
m−1∑
k=0

M(s)Bα−k−1 (s)
Kα−k (s) V(k) (0) . (11)

Proof. Since
(

RLIα
0,t

CDα
0,t

)
V (t) = V (t) −

∑m−1
k=0

V(k) (0)
k! tk. According to Lemma 1, we

have:

TN G
[(

RLIα
0,t

CDα
0,t

)
V (t)

]
= TN G

[
V (t)−

m−1∑
k=0

V(k) (0)
k! tk

]
, (12)

therefore, the Eq 12 becomes,

Kα (s)
Bα (s) TN G

[
CDα

0 V (t)
]

= TN G [V (t)]−
m−1∑
k=0

V(k) (0)
k!

k!M(s)Kk(s)
Bk+1 (s) .

Finally, we obtain:

TN G
[

CDα
0,tV (t)

]
= Bα (s)

Kα (s)TN G [V (t)]−
m−1∑
k=0

M(s)Bα−k−1 (s)
Kα−k (s) V(k) (0) .

This concludes the proof.

4. Main results and Mathematical Modeling of Smoking epidemic
Model

Many mathematical models are employed to grasp biological phenomena through an
examination of their dynamic behavior. In our research, we consider a system comprising
five non linear FDEs that characterize the smoking epidemic model, incorporating the
CDα

0,t, where 0 < α < 1. This model is formulated in [28] as follows:

CDα
0,t(P (t)) = λ−βP (t)S(t)−µP (t),

CDα
0,t(O (t)) = βP (t)S(t)−α1O(t)−µO(t),

CDα
0,t(S (t)) = α1O(t)+α2S(t)Q(t)− (µ+γ)S(t),

CDα
0,t(Q(t)) = −α2S(t)Q(t)−µQ(t)+γ (1− δ)S(t),

CDα
0,t(L(t)) = δγS(t)−µL(t),

(13)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 8 of 26

subject to the following initial conditions

P (0) = b1,O (0) = b2,S (0) = b3,Q(0) = b4,L(0) = b5. (14)

The behavior of the model dynamics under various parameter settings is illustrated in
Figure 1.

Figure 1: A graphical illustration of the envisaged mathematical model concerning the
smoking epidemic

For α = 1, the system reduces to ordinary differential equations. In this model 13, we
divide the total population size at time t denoted by N(t) into five subgroups: potential
smokers represented by P (t), current smokers denoted as S(t), occasional smokers labeled
O(t), individuals who have permanently quit smoking indicated by L(t), and those who
have temporarily quit smoking represented by Q(t).The parameters used in model 13,
along with their detailed descriptions, are listed in Table 2 below.

Table 2: Parameters used in the model and their descriptions.

Parameter Description
λ Recruitment rate in population P
β Effective contact rate between S and P
µ Natural death rate
α1 Rate at which occasional smokers convert to regular smokers
α2 Contact rate between smokers and temporary quitters who revert back to smoking
γ Rate of quitting smoking
1− δ Fraction of smokers who temporarily quit smoking (at rate γ)
δ Fraction of smokers who permanently quit smoking

Lemma 4. When α = 1 in system (13), all possible solutions are bounded and confined
within the region defined by:

Ω =
{

(P,O,S,Q,L) ∈ R5
+ : P +O +S +Q+L 6

λ

µ

}
. (15)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 9 of 26

Proof. Let (P,O,S,Q,L) ∈ R5
+ be any solution with not negative subject to initial

conditions, then:

dN(t)
d(t) = λ−µN(t),

Consequently,
0 6 N(t) 6 λ

µ
+N(0)e−µt,

where N(0) represents the initial value of N(t). Thus,

0 6 N(t) 6 λ

µ
as t → ∞.

Therefore, all possible solutions of the system lie within the region Ω. Hence, Ω is positively
invariant.

Combining the equations from the (13), and considering the linearity of CDα
0,t, we

derive:
CDα

0,t(N (t)) = λ−µN(t). (16)

When the initial conditions for the model (13) are positive, it has a unique positive so-
lution (see the work of [28]). We also address in the following that system (13) has two
equilibrium, the Smoking-Free Equilibrium (SFE) and the Smoking-Present Equilibrium
(SPE). Then, the local and global stability results of the equilibrium are also obtained.
First, employing the next-generation matrix method as formulated in [? ], we analyze the
existence of the point (SFE), Let Y = (O,S,Q,L,P )T , whereby the model can be expressed
as:

CDα
0,t(Y) = F(Y)+H(Y), (17)

where

F(Y) =


βP (t)S(t)

0
0
0
0


, (18)

and

H(Y) =


−(α1 +µ)O

α1O +α2SQ− (µ+γ)S
−α2SQ−µQ+γ − δγS

δγS −µL
λ−βPS −µP

 . (19)

By replacing S = 0 in (13), we obtian the following result.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 10 of 26

Proposition 1. For the our model (13), an SEP exists, denoted by E0, with E0 =
(λ

µ
,0,0,0,0).

Proposition 2. The Reproduction Number, R0, is determined as:

R0 = ρ(−FV −1) = α1βP0
(µ+γ)(µ+α1) , (20)

where:

F = DF(E0) =


0 0 βP0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0

 , (21)

V = DH(E0) =


0 −(α1 +µ) 0 0 0
0 α1 −(µ+γ) 0 0
0 0 γ(1− δ) −µ 0
0 0 γδ 0 −µ

−µ 0 −βP0 0 0

 . (22)

Here, DH(E0) and DF(E0) denote the Jacobian matrices of H(Y) and F(Y) evaluated at
E0 respectively.

Proof. A simple calculation gives: 21 and 22 (see [? ]). So,

V = DV(E0) =

−α1 −µ 0 0
α1 −µ−γ 0
0 γ −γδ −µ

 , F = DF(E0) =

0 βP0 0
0 0 0
0 0 0

 ,

V −1 =


− 1

µ+α1
0 0

− α1
(µ+γ)(µ+α1) − 1

(µ+γ) 0

− α1γ(1− δ)
µ(µ+γ)(µ+α1) − γ(1− δ)

µ(µ+γ) − 1
µ

 ,

and

−FV −1 =


α1βP0

(µ+γ)(µ+α1
βP0
µ+γ 0

0 0 0
0 0 0

 ,

then:
R0 = ρ(−FV −1) = α1βP0

(µ+γ)(µ+α1) .



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 11 of 26

Proposition 3. For the model(13) there exists E∗ the SPE, with E∗ = (P ∗,O∗,S∗,Q∗,L∗),
where:

P ∗ = λ

µ+βS∗ .

O∗ = βλS∗

(µ+βS∗)(µ+α1) .

Q∗ = γS∗ −S∗γδ

α2S∗ +µ
.

L∗ = δγS∗

µ
.

(23)

with S∗ satisfies the equation:

AS∗2 +BS∗ +C = 0, (24)

where
A = βγα2(α1 +µ)(δ −1)+α2β(α1 +µ)(µ+γ).
B = µα2γ(δ −1)(α1 +µ)−α1βλα2 +(α1 +µ)(µ+γ)(βµ+µα2).
C = µ2(µ+γ)(α1 +µ)−α1βλµ.

Proof. To assess the existence of a positive E∗ of (13), consider S∗ > 0 and
CDα

0,t(P ∗) =C Dα
0,t(O∗) =C Dα

0,t(S∗) =C Dα
0,t(Q∗) =C Dα

0,t(L∗) = 0.

This gives

λ−βP ∗S∗ −µP ∗ = 0 (25)

βP ∗S∗ −α1O∗ −µO∗ = 0 (26)

α1O∗ +α2S∗Q∗ − (µ+γ)S∗ = 0 (27)

−α2S∗Q∗ −µQ∗ +γ (1− δ)S∗ = 0 (28)

δγS∗ −µL∗ = 0 (29)

From Equations (25), (27), (28), (29), we obtain 23. Finally, substituting O∗ and Q∗ in
Equation (26) gives:

S∗
[

α1βλ

(α1 +µ)(µ+βS∗) − α2γ(δ −1)S∗

α2S∗ +µ
−µ−γ

]
= 0.

Since S∗ 6= 0, we obtain AS∗2 +BS∗ +C = 0, with:

A = βα2(α1 +µ)(δγ +µ) > 0,

C = µ2(µ+γ)(α1 +µ)−α1βλµ = µ

(µ+γ)(α1 +µ)(1−R0),

B = −µα2γ(1− δ)(α1 +µ)+µβ(α1 +µ)(µ+γ)+ α2
µ(α1 +µ)(µ+γ)(1−R0).



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 12 of 26

Theorem 3. (i) If R0 = 1 , and µ

γ
> α2 −1, there is no current positive E∗,

(ii) if R0 < 1, and B > 0, there is no current positive equilibrium point E∗,

(iii) If R0 > 1, a unique positive endemic equilibrium E∗ exists.

Proof. From Equation 24, we deduce that:

(i) When R0 = 1, it follows that C = 0 yielding the solutions S∗
1 = 0 and S∗

2 = −B

A
.

Consequently, no positive solution exists if B > 0 ( since A is always positive) .

(ii) If R0 < 1 and B > 0, then A,C > 0. There is no positive solution to the equation.

(iii) If R0 > 1, then there A,C > 0 (There is one change in the sign of the terms, so

there is a positive solution) and ∆ = B2 −4AC > 0. which gives S1 = −B −
√

∆
2A

and

S2 = −B +
√

∆
2A

. Note that:
√

B2 −4AC > −B. So: S2 = −B +
√

∆
2A

it is a positive
solution to the equation.

The stability of the point E0 is studied in the following theorem, by using the result
proven in [29, 30].

Theorem 4. The SFE point E0 is locally asymptotically stable for R0 < 1, and unstable
for R0 > 1.

Proof. Evaluating the Jacobian matrix of (13) at E0(λ

µ
,0,0,0,0), yields:

J(E0) =


−µ 0 −βP0 0 0
0 −(α1 +µ) βP0 0 0
0 α1 −µ−γ 0 0
0 0 γ −γδ −µ 0
0 0 δγ 0 −µ

 ,

then:

J(E0) =


−µ 0 −βP0 0 0
0 −(α1 +µ) βP0 0 0
0 0 βP0 − (µ+γ)(α1+µ)

α1
0 0

0 0 γ −γδ −µ 0
0 0 δγ 0 −µ

 ,

The characteristic polynomial takes the form:

P (Λ) = (Λ+µ)3(Λ+α1 +µ)
(

Λ−βP0 + (µ+γ)(α1 +µ)
α1

)
.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 13 of 26

Accordingly, the eigenvalues of J(E0) are:

Λ1 = −µ, Λ2 = −(α1 +µ), Λ3 = (µ+γ)(α1 +µ)
α1

(R0 −1) .

If R0 < 1, then Λ3 < 0. As a result, all eigenvalues of J(E0) are negative and satisfy the
required stability conditions:

|arg(Λi)| > α
π

2 , i = 1,2,3,

ensuring that E0 is locally asymptotically stable (see [29, 30]). However, when R0 = 1, we
have Λ3 = 0, implying that E0 is locally stable. If R0 > 1, then Λ3 > 0 and

|arg(Λ3)| ≤ α
π

2 ,

which implies that E0 is unstable.

The local stability analysis of the equilibrium point E∗ is examined by by evaluating
the Jacobian matrix of system (13) at E∗, yielding:

J (E∗) =


−(µ−βS∗) 0 −βP ∗ 0 0

βS∗ −(α1 +µ) βP ∗ 0 0
0 α1 α2Q∗ − (µ+γ) α2S∗ 0
0 0 −α2Q∗ +γ − δγ −(α2S∗ +µ) 0
0 0 γδ 0 −µ


We note that Λ1 = −µ is eigenvalue then the local stability of E∗ of depends of sign of
reale of eigenvalue of the following matrix:

J1 (E∗) =


−βS∗ −µ 0 −βP ∗ 0

βS∗ −α1 −µ βP ∗ 0
0 α1 α2Q∗ − (µ+γ) α2S∗

0 0 −α2Q∗ +γ(1− δ) −α2S∗ −µ

 ,

the characteristic polynomial is represented by:

P (Λ) = (Λ+(βS∗ +µ))(Λ3 +A1Λ2 +A2Λ+A3).

Such that,

A1 = (βS∗ +µ)(α1µ)
βS∗ −α2Q∗ +(µ+γ)+α2S∗ +µ,

A2 = (βS∗ +µ)(α1 +µ)
βS∗ (α2Q∗ − (µ+γ)−α2S∗ −µ)+(α2Q∗ − (µ+γ))(−α2S∗ −µ)− µα1βP ∗

βS∗

+α2S∗(α2Q∗ −γ(1− δ)),

A3 = (µ+βS∗)(γ +µ−α2Q∗)(α1 +µ)(α2S∗ +µ)−µα1βP ∗(α2 +µ)
βS∗



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 14 of 26

+ α2S∗(α2Q∗ −γ +γδ))(βS∗ +µ)(+α1 +µ)
βS∗ .

The eigenvalues consist of
Λ1 = −(βS∗ +µ),

and the roots of the
Λ3 +A1Λ2 +A2Λ+A3 = 0.

All eigenvalues have negative real parts if and only if these conditions are verified:

A1 > 0, A2 > 0, A3 > 0, and A1A2 −A3 > 0.

Hence, we present the following result:

Theorem 5. The equilibrium E∗ of system (13) is locally asymptotically stable under the
condition that:

A1 > 0, A2 > 0, A3 > 0, and A1A2 > A3.

5. A hybrid New general transform Adomian Decomposition Method
for solving the Fractional Model

The Junaid Integral Transform-Adomian Decomposition Method (TN G-ADM) is a hy-
brid approach that combines the Junaid integral transform with the ADM technique. This
method leverages the strengths of both techniques to simplify and efficiently solve complex
equations.

5.1. Numerical Solution for the Model 13

We can express the model (13) as:(
CDα

0,tZ
)

(t) = R(t,Z (t)) = Gi (t,P,O,S,Q,L) , i = 1, ...,5. (30)

Z (0) = Z0, (31)

where Z (t) = (P,O,S,Q,L), Z (0) = (b1, b2, b3, b4, b5).
Given that Gi(0) = 0, the problem has a solution. Upon applying the Iα

a+ to both sides
of 30, we obtain:

Y (t)−Z (0) = Iα
a+(R(t,Z (t))). (32)

The application of the TN G transform to both sides of 32, yields the following relation:

TN G [Z (t)]−TN G [Z (0)] =
(K (s)

B(s)

)α

TN G [R(t,Z (t))] , (33)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 15 of 26

With initial conditions 30. Hence:

TN G [(P (t))] = M(s)
B(s) b1 +

(K (s)
B(s)

)α

TN G [λ−βP (t)S(t)−µP (t)] ,

TN G [(O (t))] = M(s)
B(s) b2 +

(K (s)
B(s)

)α

TN G [βP (t)S(t)−α1O(t)−µO(t)] ,

TN G [(S (t))] = M(s)
B(s) b3 +

(K (s)
B(s)

)α

TN G [α1O(t)+α2S(t)Q(t)− (µ+γ)S(t)] ,

TN G [(Q(t))] = M(s)
B(s) b4 +

(K (s)
B(s)

)α

TN G [−α2S(t)Q(t)−µQ(t)+γ (1− δ)S(t)] ,

TN G [(L(t))] = M(s)
B(s) b5 +

(K (s)
B(s)

)α

TN G [δγS(t)−µL(t)] .

(34)

Assuming the method yields the solution as:

P(t) =
∞∑

n=0
Pn(t), O(t) =

∞∑
n=0

On(t), S(t) =
∞∑

n=0
Sn(t), Q(t) =

∞∑
n=0

Qn(t), L(t) =
∞∑

n=0
Ln(t).

(35)
The expressions for PS and SQ are:

P(t)S(t) =
∞∑

n=0
Xn(t), S(t)Q(t) =

∞∑
n=0

Zn(t), (36)

where Xn and Zn is as:

Xm(t) = 1
m!

dm

dλm

[
m∑

n=0
λnPn(t)

m∑
n=0

λnSn(t)
]

λ=0
, (37)

Zm(t) = 1
m!

dm

dλm

[
m∑

n=0
λnSn(t)

m∑
n=0

λnQn(t)
]

λ=0
. (38)

Substituting (35) and (36) into (34) results in:

TN G [
∑∞

n=0 Pn(t)] = M(s)
B(s) b1 +

(K (s)
B(s)

)α

TN G [λ−β
∑∞

n=0 Xn(t)−µ
∑∞

n=0 Pn(t)] ,

TN G [
∑∞

n=0 On(t)] = M(s)
B(s) b2 +

(K (s)
B(s)

)α

TN G [β
∑∞

n=0 Xn(t)−α1
∑∞

n=0 On(t)−µΣ∞
n=0On(t)] ,

TN G [
∑∞

n=0 Sn(t)] = M(s)
B(s) b3 +

(K (s)
B(s)

)α

TN G [α1
∑∞

n=0 On(t)+α2
∑∞

n=0 Zn(t)− (µ+γ)
∑∞

n=0 Sn(t)] ,

TN G [
∑∞

n=0 Qn(t)] = M(s)
B(s) b4 +

(K (s)
B(s)

)α

TN G [−α2
∑∞

n=0 Zn(t)−µ
∑∞

n=0 Qn(t)+γ (1− δ)
∑∞

n=0 Sn(t)] ,

TN G [
∑∞

n=0 Ln(t)] = M(s)
B(s) b5 +

(K (s)
B(s)

)α

TN G [δγ
∑∞

n=0 Sn(t)−µ
∑∞

n=0 Ln(t)] .

Now by comparing the same terms on both sides we can get:

TN G [P0(t)] = M(s)
B(s) b1,TN G [O0(t)] = M(s)

B(s) b2,TN G [S0(t)] = M(s)
B(s) b3, (39)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 16 of 26

TN G [Q0(t)] = M(s)
B(s) b4,TN G [L0(t)] = M(s)

B(s) b5. (40)

Similarly, we have:

TN G [P1(t)] =
(K (s)

B(s)

)α

TN G [λ−βX0(t)−µP0(t)]

TN G [O1(t)] =
(K (s)

B(s)

)α

TN G [βX0(t)− (α1 +µ)O0(t)]

TN G [S1(t)] =
(K (s)

B(s)

)α

TN G [α1O0(t)+α2Z0(t)− (µ+γ)S0(t)]

TN G [Q1(t)] =
(K (s)

B(s)

)α

TN G [−α2Z0(t)−µQ0(t)+γ (1− δ)S0(x)]

TN G [L1(t)] =
(K (s)

B(s)

)α

TN G [δγS0(t)−µL0(t)]

(41)



TN G [P2(t)] =
(K (s)

B(s)

)α

TN G [λ−βX1(t)−µP1(t)]

TN G [O2(t)] =
(K (s)

B(s)

)α

TN G [βX1(t)−α1O1(t)−µO1(t)]

TN G [S2(t)] =
(K (s)

B(s)

)α

TN G [α1O1(t)+α2Z1(t)− (µ+γ)S1(t)]

TN G [Q2(t)] =
(K (s)

B(s)

)α

TN G [−α2Z1(t)−µQ1(t)+γ (1− δ)S1(t)]

TN G [L2(t)] =
(K (s)

B(s)

)α

TN G [δγS1(t)−µL1(x)]

(42)



TN G [P3(t)] =
(K (s)

B(s)

)α

TN G [λ−βX2(t)−µP2(t)]

TN G [O3(t)] =
(K (s)

B(s)

)α

TN G [βX2(t)− (α1 +µ)O2(t)]

TN G [S3(t)] =
(K (s)

B(s)

)α

TN G [α1O2(t)+α2Z2(t)− (µ+γ)S2(t)]

TN G [Q3(t)] =
(K (s)

B(s)

)α

TN G [−α2Z2(t)−µQ2(t)+γ (1− δ)S2(t)]

TN G [L3(t)] =
(K (s)

B(s)

)α

TN G [δγS2(t)−µL2(t)]

(43)



TN G [Pk+1(t)] =
(K (s)

B(s)

)α

TN G [λ−βXk(t)−µPk(t)]

TN G [Ok+1(t)] =
(K (s)

B(s)

)α

TN G [βXk(t)−α1Ok(t)−µOk(t)]

TN G [Sk+1(t)] =
(K (s)

B(s)

)α

TN G [α1Ok(t)+α2Zk(t)− (µ+γ)Sk(t)]

TN G [Qk+1(t)] =
(K (s)

B(s)

)α

TN G [−α2Zk(t)−µQk(t)+γ (1− δ)Sk(t)]

TN G [Lk+1(t)] =
(K (s)

B(s)

)α

TN G [δγSk(t)−µLk(t)]

(44)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 17 of 26

by using the inverse transformation of TN G , we have the initial conditions, and X0 =
b1b3 , Z0 = b3b4 

P1 (t) = [λ−βb1b3 −µb1] tα

Γ (α +1)
O1 (t) = [βb1b3 − (µ+α1)b2] tα

Γ (α +1)
S1 (t) = [α1b2 +α2b3b4 − (µ+γ)b3] tα

Γ (α +1)
Q1 (t) = [−α2b3b4 −µb4 +γ(1−σ)b3] tα

Γ (α +1)
L1 (t) = [σγb3 −µb5] tα

Γ (α +1) .

(45)

Consequently: X1(t) = P1(t)S0(t) + S1(t)P0(t) , Z1(t) = Q1(t)S0(t) + S1(t)Q0(t), we
pose:

wP = λ−βb1b3 −µb1
wO = βb1b3 − (µ+α1)b2
wS = α1b2 +α2b3b4 − (µ+γ)b3
wQ = −α2b3b4 −µb4 +γ(1−σ)b3
wL = σγb3 −µb5

(46)

So: 

P2 (t) = λ

Γ(α +1) tα + 1
Γ(2α +1) [−βwSb1 − (βb3 +µ)wP ] t2α,

O2 (t) = 1
Γ(2α +1) [βwSb1 +βwpb3 − (µ+α1)wO] t2α,

S2 (t) = 1
Γ(2α +1) [α1wO +α2wQb3 +α2wSb4 − (µ+γ)wS ] t2α,

Q2 (t) = 1
Γ(2α +1) [(−α2b4 +γ(1−σ))wS +(−α2b3 −µ)wQ] t2α,

L2 (t) = 1
Γ(2α +1) [σγwS −µwL] t2α,

(47)

we pose:
wP P = −βwSb1 − (βb3 +µ)wP ,
wOO = βwSb1 +βwpb3 − (µ+α1)wO,
wSS = α1wO +α2wQb3 +α2wSb4 − (µ+γ)wS ,
wQQ = (−α2b4 +γ(1−σ))wS +(−α2b3 −µ)wQ,
wLL = σγwS −µwL,

(48)

we have:
X2(t) = P2(t)S0(t)+S2(t)P0(t)+S1(t)P1(t),

Z2(t) = Q2(t)S0(t)+S2(t)Q0(t)+S1(t)Q1(t).



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 18 of 26

So:

P3 (t) = λ

Γ(α +1) tα − λ(βb3 +µ)
Γ(2α +1) t2α + 1

Γ(3α +1)

[
−β(wP P b3 +wSSb1)−µwP P − βwP wSΓ(2α +1)

(Γ(α +1))2

]
t3α,

O3 (t) = βλb3
Γ(2α +1) t2α + 1

Γ(3α +1)

[
βb3wP P +βb1wSS − (µ+α1)wOO + βwP wSΓ(2α +1

(Γ(α +1))2

]
t3α,

S3 (t) = 1
Γ(3α +1)

[
α1wOO +α2wQQb3 +α2wSSb4 − (µ+γ)wSS + α2wQwSΓ(2α +1)

(Γ(α +1))2

]
t3α,

Q3 (t) = 1
Γ(3α +1)

[
−α2wQQb3 −α2wSSb4 −µwQQ +γ(1− δ)wSS − α2wQwSΓ(2α +1)

(Γ(α +1))2

]
t3α,

L3 (t) = 1
Γ(3α +1) [(δγwSS −µwLL)] t3α.

(49)
On solving the above system using initial conditions:

(b1, b2, b3, b4, b5) = (40, 10, 20, 10, 5).

P1 (t) = λ−800β −40µ

Γ (α +1) tα

O1 (t) = 800β −10µ−10α1
Γ (α +1) tα

S1 (t) = 10α1 +200α2 −20(µ+γ)
Γ (α +1) tα

Q1 (t) = −200α2 −10µ+20γ(1−σ)
Γ (α +1) tα

L1 (t) = 20σγ −5µ

Γ (α +1) tα

(50)

We define:
wP = λ−800β −40µ,

wO = 800β −10µ−10α1,

wS = 10α1 +200α2 −20(µ+γ) ,

wQ = −200α2 −10µ+20γ(1−σ),
wL = 20σγ −5µ.

So 

P2 (t) = λ

Γ (α +1) tα + 1
Γ(2α +1) [−40βwS − (20β +µ)wP ] t2α,

O2 (t) = 1
Γ(2α +1) [40βwS +20βwp − (µ+α1)wO] t2α,

S2 (t) = 1
Γ(2α +1) [α1wO +20α2wQ +10α2wS − (µ+γ)wS ] t2α,

Q2 (t) = 1
Γ(2α +1) [−10α2 +γ(1−σ)wS +(−20α2 −µ)wQ] t2α,

L2 (t) = 1
Γ(2α +1) [σγwS −µwL] t2α,

(51)



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 19 of 26

we pose:
wP P = −40β wS − (20β +µ)wP ,

wOO = 40β wS +20β wP − (µ+α1)wO,

wSS = α1wO +20α2wQ +10α2wS − (µ+γ)wS ,

wQQ = (−10α2 +γ(1−σ))wS +(−20α2 −µ)wQ,

wLL = σγ wS −µwL.

P3 (t) = λ

Γ(α +1) tα − λ(20β +µ)
Γ(2α +1) t2α + 1

Γ(3α +1)

[
−β(20wP P +40wSS)−µwP P − βwP wSΓ(2α +1)

(Γ(α +1))2

]
t3α,

O3 (t) = 20βλ

Γ(2α +1) t2α + 1
Γ(3α +1)

[
20βwP P +40βwSS − (µ+α1)wOO + βwP wSΓ(2α +1

(Γ(α +1))2

]
t3α,

S3 (t) = 1
Γ(3α +1)

[
α1wOO +20α2wQQ +10α2wSS − (µ+γ)wSS + α2wQwSΓ(2α +1)

(Γ(α +1))2

]
t3α,

Q3 (t) = 1
Γ(3α +1)

[
−20α2wQQ −10α2wSS −µwQQ +γ(1− δ)wSS − α2wQwSΓ(2α +1)

(Γ(α +1))2

]
t3α.

L3 (t) = 1
Γ(3α +1) [(δγwSS −µwLL)] t3α.

(52)
In particular, The global ADM solution is written as the sum of the series:

P (t) = 40+P1(t)+P2(t)+P3(t)+ · · · ,

O(t) = 10+O1(t)+O2(t)+O3(t)+ · · · ,

S(t) = 20+S1(t)+S2(t)+S3(t)+ · · · ,

Q(t) = 10+Q1(t)+Q2(t)+Q3(t)+ · · · ,

L(t) = 5+L1(t)+L2(t)+L3(t)+ · · · ,

where the terms P1, . . . ,L3 above provide a global approximation of order 3α:

(P,O,S,Q,L)(t) ≈ (b1, b2, b3, b4, b5)+
3∑

n=1
(Pn(t),On(t),Sn(t),Qn(t),Ln(t)) .

5.2. Simulations and Discussion

We perform a numerical investigation of the fractional smoking epidemic model with
a nonlinear incidence rate to assess the influence of key parameters, particularly the frac-
tional order α. Simulations are conducted for a range of integer, fractional, and mixed
integerfractional values of α, enabling a comparative analysis of their effects on the model
dynamics. The results indicate that fractional orders introduce only subtle modifications to
the overall epidemic evolution, yet they enhance the accuracy of the numerical approxima-
tion across different fractional frameworks. Notably, for potential smokers (non-smokers),
the simulations reveal an increase in their population size when α takes fractional values.
These findings are illustrated in Figure 2, while Figure 3 provides a detailed depiction of
the variation of the R0.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 20 of 26

From the perspective, one of the infected compartments, S(t), remains nonzero and
exhibits convergence for both integer and fractional orders of α, whereas the other infected
compartment approaches zero over time. Furthermore, the infection rate is observed to
decrease progressively as the fractional order α decreases. The numerical values of the
parameters λ,β,ζ,σ,γ,α1, and α2 used in these simulations are summarized in Table 3.

Parameters λ β µ α1 α2 γ δ

Values 1 0.14 0.05 0.002 0.0025 0.8 0.1

Table 3: specific values of parameters used in system 13.

Figure 2: Numerical solutions for P (t) ,O (t) ,Q(t) ,P (t) and L(t)

6. Numerical Approximation via the Hybrid TN G-ADM–ANN Method

We now introduce a new hybrid method that combines the TN G-ADM method with an
Artificial Neural Network (ANN) algorithm to approximate the solutions S, P, L, O, Q of



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 21 of 26

Figure 3: (a): Variation of R0 with α1 (Fixed β), (b): Variation of R0 with β (Fixed α1),
(c):3D Variation of R0 as a function of α1 and β

the proposed model over the time domain t ∈ [0,T ]. The proposed method is implemented
through a two-step strategy:

(i) First, a reference solution is obtained by applying the TN G-ADM technique to the
model.

(ii) Next, an ANN is trained to learn the dynamics of the variables by minimizing the
gap between its predictions and the TN G-ADM reference solution.

The ANN-based algorithm is designed to approximate the solution of the fractional
smoking model by taking time t as input and producing a vector output [S,P,L,O,Q]. It
consists of a feedforward neural network with an input layer containing a single neuron,
two hidden layers with 50 neurons each using non linear activation functions, and an
output layer with five neurons employing a linear activation function.

Prior to training, both the input data (t) and the reference solutions (obtained via
TN G-ADM) are normalized to enhance numerical stability. The training process mini-
mizes the Mean Squared Error (MSE) between the ANN predictions and the reference
values, employing Bayesian regularization (trainbr) to prevent overfitting. The dataset



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 22 of 26

is divided into three distinct subsets: 70% allocated for training the model, 15% reserved
for validation to fine-tune hyperparameters, and the remaining 15% set aside for evalu-
ating the models performance on unseen data. The network parameters are iteratively
adjusted until convergence or until the maximum number of epochs is reached.

Once trained, the ANN provides rapid and accurate predictions of the system’s evo-
lution over the entire time domain t ∈ [0,T ]. By integrating the fractional differential
equations into its loss function, this hybrid approach leverages the reliability of the TN G-
ADM method and the generalization capability of the ANN, offering a powerful tool for
solving complex differential systems.

Figure 4: Numerical solutions for P (t), O(t), Q(t), P (t), and L(t) obtained using the ANN
algorithm for α = 1.

In this study, we evaluate the performance of the ANN by comparing its predictions
with the reference solution obtained via the TN G-ADM method.

Figure 4 illustrates the plots for each variable. In these plots, the TN G-ADM solutions
are represented by solid lines, while the ANN predictions are depicted as dashed lines. The
close overlap of the curves indicates that the ANN has accurately captured the systems
dynamics.

Table 4 summarizes the global results for each variable. For each variable, the table lists



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 23 of 26

the mean value computed using TN G-ADM (denoted as ADM (Mean)), the corresponding
ANN-predicted mean (denoted as ANN (Mean)), and the Mean Absolute Error (MAE),
defined by:

MAE = |XADM −XANN| ,

where X represents each variable.

Table 4: Comparison of global results for variables S, P , O, Q, and L obtained using the
TN G-ADM and ANN methods.

Variable TN G-ADM (Mean) ANN (Mean) MAE

S 0.85 0.84 0.01
P 0.65 0.66 0.01
O 0.90 0.89 0.01
Q 1.10 1.08 0.02
L 1.20 1.18 0.02

The low Mean Absolute Error values, ranging from 0.01 to 0.02, indicate a strong
agreement between the ANN predictions and the reference TN G-ADM solutions, with an
error of approximately 1%. This level of accuracy is generally sufficient for many scientific
and engineering applications. If higher accuracy is needed, enhancements to the network
architecture or adjustments to the training parameters may be considered.

The plotted trajectories for S, P , L, O, and Q indicate that the ANN effectively
replicates the system dynamics as captured by the TN G-ADM method. The near-perfect
alignment between the solid curves (representing TN G-ADM solutions) and the dashed
curves (ANN predictions) throughout the time domain demonstrates the ANNs high ap-
proximation accuracy. This strong agreement underscores the ANNs potential to serve as
a fast and reliable surrogate model, provided it is properly trained.

7. Conclusion

This study proposed an efficient and accurate hybrid framework, the TN G-ADM–ANN
scheme, for analyzing a fractional-order smoking epidemic model formulated with the
Caputo fractional derivative. The model captures the intrinsic memory effects and hered-
itary properties of smoking dynamics, which are often neglected in classical integer-order
models. By integrating the newly developed TN G integral transform with the Adomian De-
composition Method (ADM) and Artificial Neural Networks (ANN), the approach provides
rapidly convergent analytical–numerical solutions with low computational complexity.

The obtained numerical and graphical results confirm that both the fractional order
and the epidemiological parameters have a substantial impact on the systems qualitative
behavior and stability. These findings highlight the suitability of fractional calculus in
modeling real-world processes that depend on historical states, particularly those involving



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 24 of 26

behavioral or biological memory such as smoking spread, addiction persistence, or cessation
dynamics.

The incorporation of ANN into the ADM–transform framework significantly improves
convergence speed and approximation accuracy by leveraging data-driven learning capa-
bilities. This hybridization demonstrates a powerful synergy between analytical decom-
position and machine learning, offering a generalizable strategy for nonlinear fractional
differential systems.

Overall, the TN G-ADM–ANN framework represents a reliable and flexible tool for solv-
ing nonlinear fractional models across various scientific and engineering domains. Future
research may focus on extending this methodology to stochastic or time-delayed systems,
exploring alternative fractional operators such as Atangana–Baleanu or Prabhakar types,
and developing high-performance or deep-learning-based implementations to address large-
scale and multidimensional fractional models. These extensions would further consolidate
the relevance and applicability of the proposed approach in complex epidemiological and
dynamical systems.

References

[1] R. Selvaraj, A. Kumar, and P. Prakash. Recent advancements in integral transforms
for solving partial differential equations. Applied Mathematics and Computation,
459:128312, 2025.

[2] Y. Feng, Z. Li, and H. Zhang. Hybrid approaches combining integral transforms and
decomposition methods for nonlinear differential equations. Journal of Computational
and Applied Mathematics, 439:115623, 2025.

[3] G. Adomian. Solving Frontier Problems of Physics: The Decomposition Method.
Kluwer Academic Publishers, Boston, 1994.

[4] M. Alhazmi, A.F. Aljohani, N.E. Taha, S. Abdel-Khalek, M. Bayram, and S. Saber.
Application of a fractal fractional operator to nonlinear glucose-insulin systems: Ado-
mian decomposition solutions. Comput. Biol. Med., 196(Pt A):110453, Sep 2025.

[5] G.K. Watugala. Sumudu transform: A new integral transform to solve differen-
tial equations and control engineering problems. Int. J. Math. Educ. Sci. Technol.,
24:35–43, 1993.

[6] T.M. Elzaki. The new integral transform ’elzaki transform’. Glob. J. Pure Appl.
Math., 7:57–64, 2011.

[7] M. Khan, T. Salahuddin, M.Y. Malik, M.S. Alqarni, and A.M. Alkahtani. Numerical
modeling and analysis of bioconvection on mhd flow due to an upper paraboloid
surface of revolution. Physica A, 553:124231, 2020.

[8] K.S. Aboodh. The new integral transform. Global J. Pure Appl. Math., 9(1):35–43,
2013.

[9] M.M. Abdelrahim Mahgoub. The new integral transform mohand transform. Adv.
Theoret. Appl. Math., 12(2):113–120, 2017.

[10] M.A. Mahgoub and M. Mohand. The new integral transform “sawi transform”. Adv.
Theoret. Appl. Math., 14(1):81–87, 2019.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 25 of 26

[11] S. Maitama and W. Zhao. New integral transform: Shehu transforma generalization
of sumudu and laplace transform for solving differential equations. Int. J. Anal. Appl.,
17(2):167–190, 2019.

[12] R. Belgacem, D. Baleanu, and A. Bokhari. Shehu transform and applications to
caputo-fractional differential equations. Int. J. Anal. Appl., 17:917–927, 2019.

[13] R. Belgacem, A. Bokhari, and B. Sadaoui. Shehu transform of hilferprabhakar frac-
tional derivatives and applications on some cauchy type problems. Adv. Theory Non-
linear Anal. Appl., 5(2):203–214, 2021.

[14] A. Bokhari, D. Baleanu, and R. Belgacem. Application of shehu transform to atan-
ganabaleanu derivatives. J. Math. Comput. Sci., 20(2):101–107, 2020.

[15] H. Kamal and A. Sedeeg. The new integral transform kamal transform. Adv. Theoret.
Appl. Math., 11(4):451–458, 2016.

[16] H. Jafari. A new general integral transform for solving integral equations. J. Adv.
Res., 2020.

[17] R. Belgacem, A. Bokhari, D. Baleanu, and S. Djilali. New generalized integral trans-
form via dzherbashiannersesian fractional operator. IJOCTA, 14(2):90–98, 2024.

[18] Muhammad Amin S. Murad, Salim S. Mahmood, Homan Emadifar, Wael W. Mo-
hammed, and Karim K. Ahmed. Optical soliton solution for dual-mode time-fractional
nonlinear schrödinger equation by generalized exponential rational function method.
Results in Engineering, 2025.

[19] I. Alraddadi, F. Alsharif, S. Malik, H. Ahmad, T. Radwan, and K.K. Ahmed. Inno-
vative soliton solutions for a (2+1)-dimensional generalized kdv equation using two
effective approaches. AIMS Mathematics, 9(12):34966–34980, 2024.

[20] P. O. Mohammed, M. R. Alharthi, M. A. Yousif, A. A. Lupas, and S. M. Azzo. Mod-
eling and neural network approximation of asymptotic behavior for delta fractional
difference equations with mittag-leffler kernels. Fractal and Fractional, 9(7):452, 2025.

[21] Majeed Ahmad Yousif, Dumitru Baleanu, Mohamed Abdelwahed, Shrooq Mo-
hammed Azzo, and Pshtiwan Othman Mohammed. Finite difference β-fractional
approach for solving the time-fractional fitzhughnagumo equation. Alexandria Engi-
neering Journal, 125:127–132, 2025.

[22] Karim K. Ahmed, Muhammad Bilal, Javed Iqbal, Majeed Ahmad Yousif, Dumitru
Baleanu, and Pshtiwan Othman Mohammed. An analytical algebraic method for solv-
ing nonlinear fractional differential equations with conformable fractional derivatives.
Contemporary Mathematics, 6(5):5925–5954, 2025.

[23] I.M. Junaid. The generalization of integral transforms combined with he’s polynomial.
Eur. J. Pure Appl. Math., 16(2):1024–1046, 2023.

[24] Shengqiang Zhang, Yanling Meng, Amit K. Chakraborty, and Hao Wang. Controlling
smoking: A smoking epidemic model with different smoking degrees in deterministic
and stochastic environments. Mathematical Biosciences, page 109132, 2023.

[25] R. Ullah, M. Khan, and G. Zaman. Dynamical features of a mathematical model on
smoking. J. Appl. Environ. Biol. Sci., 6(1):92–96, 2016.

[26] A.A. Kilbas, H.M. Srivastava, and J.J. Trujillo. Theory and Applications of Fractional
Differential Equations. Elsevier, Amsterdam, 2006.



R. Belgacem et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6267 26 of 26

[27] M. Caputo and M. Fabrizio. A new definition of fractional derivative without singular
kernel. Progress in Fractional Differentiation and Applications, 73:1–13, 2015.

[28] M. Abdullah, A. Ahmad, N. Raza, M. Farman, and M.O. Ahmad. Approximate
solution and analysis of smoking epidemic model with caputo fractional derivatives.
Int. J. Appl. Comput. Math., 4:1–16, 2018.

[29] D. Matignon. Stability results for fractional differential equations with applications
to control processing. Comput. Eng. Syst. Appl., 2:963, 1996.

[30] C. Li, A. Muhammadhaji, L. Zhang, and Z. Teng. Stability analysis of a fractional-
order predatorprey model incorporating a constant prey refuge and feedback control.
Adv. Differ. Equ., 2018:325, 2018.


	Introduction
	Preliminary
	Main results of  TNG  transform
	Main results and Mathematical Modeling of Smoking epidemic Model
	A hybrid New general transform Adomian Decomposition Method for solving the Fractional Model
	Numerical Solution for the Model 13 
	Simulations and Discussion

	Numerical Approximation via the Hybrid  TNG -ADM–ANN Method
	Conclusion