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

An Iterative Fractional Operational Matrix Method for
Solving Fractional Undamped Duffing Equation with

High Nonlinearity

Muhammed I. Syam1,∗, Mwaffag Sharadga2, Ishak Hashim2,
Mohd Almie Alias2, Mohammed Abuomar1, Shaher Momani3

1 Department of Mathematical Sciences, United Arab Emirates University, Al-Ain,
United Arab Emirates
2 Department of Mathematical Sciences, Faculty of Science & Technology,
Universiti Kebangsaan Malaysia, Malaysia
3 Department of Mathematics, Faculty of Science, The University of Jordan, Amman, 11942,
Jordan

Abstract. This study presents a novel iterative fractional operational matrix method for solv-
ing the highly nonlinear fractional undamped Duffing equation. The proposed method efficiently
approximates the solution by transforming the original fractional differential equation into a sys-
tem of algebraic equations using modified operational matrices. The accuracy and effectiveness of
this approach are validated through comparisons with established numerical methods and alterna-
tive analytical techniques from the literature. The results demonstrate that the proposed method
provides a highly accurate approximation with rapid convergence. Furthermore, a rigorous con-
vergence analysis is conducted to establish the existence and uniqueness of the solution. Notably,
the study explores the impact of varying the fractional order revealing its significant influence on
the system’s dynamic behavior. As the fractional order approaches unity, the fractional model
converges to its classical counterpart, highlighting the role of fractional derivatives in capturing
memory effects and hereditary properties in physical systems.

2020 Mathematics Subject Classifications: 65L03, 34A08

Key Words and Phrases: Modified operational matrices, undamped duffing equation, high
nonlinearity

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5884

Email addresses: m.syam@uaeu.ac.ae (M. Syam), P101953@siswa.ukm.edu.my (M. Sharadga),
ishak h@ukm.my (I. Hashim), almie.alias@monash.edu (M. Alias),
mohammed m@uaeu.ac.ae (M. Abuomar), shahermm@yahoo.com (S. Momani)

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



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 2 of 18

1. Introduction

The domain of fractional calculus broadens the conventional concepts of differentiation
and integration to encompass non-integer orders, thereby offering a more comprehensive
mathematical framework for these operations. The origins of fractional differentiation and
integration date back to 1695, when Leibniz and Euler first introduced these concepts [1].
A distinguishing characteristic of fractional operators is their non-local behavior, which
provides valuable insights into the historical and dynamic significance of fractional models.

In recent years, fractional differential equations have found widespread applications
in modeling diverse real-world phenomena. For instance, it is examined their utility in
analyzing and designing control systems, see [2]. Magin [3] identified three bioengineer-
ing domains where fractional calculus principles have been instrumental in developing
innovative mathematical models. Similarly, Matlob and Jamali [4] employed fractional
differential equations to study the dynamics of viscoelastic systems. In image process-
ing, fractional orders have been utilized to improve denoising techniques [5], while other
notable applications include finance [6], solid mechanics [7], and mathematical biology
[8–12].

A significant recent development in fractional calculus is the Atangana-Baleanu frac-
tional derivative [8]. Proposed by Abdon Atangana and Dumitru Baleanu in 2016, this
derivative incorporates non-singular kernels and the Mittag-Leffler function. The applica-
tion of fractional differentiation in [8] led to the formulation of a mathematical model for
heat conduction in materials. This methodology has proven effective in addressing a wide
range of real-world problems, as evidenced in [13–23].

In 1918, George Duffing introduced the Duffing equation in his publication [24]. Within
his work, Duffing proposed a simplified mathematical model

x′′(t) + a2x(t)− βx2(t)− γx3(t) = k sin(ωt), (1)

where he calculated the first term, H sin(ωt), of the periodic solution. Duffing also pro-
posed simplified forms of this equation to describe the motion of pendulums. For a sym-
metric pendulum, the equation takes the form

x′′(t) + αx(t)− γx3(t) = 0, (2)

while for an asymmetric pendulum, it becomes:

x′′(t) + αx(t)− γx2(t) = 0. (3)

Since its introduction, differential equations with polynomial nonlinearities have been
referred to as Duffing’s equations [25]. Ganji et al. [26] formulated a nonlinear differential
equation for the cubic free undamped Duffing oscillator

x′′(t) + αx(t) + βx3(t) = 0, (4)

subject to the initial conditions:

x(0) = A, x′(0) = 0. (5)



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 3 of 18

Various researchers have addressed this problem numerically. For instance, He [27] applied
the homotopy perturbation method to solve the Duffing equation, while Beléndez et al.
[28] employed a modified version of the homotopy perturbation method. Additionally,
Ramos, Syam, and Wazwaz utilized the variational iteration method [29–34], and Ghosh
et al. [35] used the Adomian decomposition method. Other techniques include the ar-
tificial parameter decomposition method by Ramos [36] and He’s parameter expanding
method by Syam [37]. Several analytical methods have been developed for solving Duffing
equations. For weak nonlinearity, approaches such as the method of multiple scales [38],
the Krylov-Bogolubov method [14], straightforward expansions [39], and the Lindstedt-
Poincaré method [39] have been employed. Strong cubic nonlinearity has also been studied,
as discussed in [40, 41].

In this paper, we extend the problem defined in equations (4)-(5) to the fractional
domain. The generalized form is given by

D2αx(t) + βx(t) + γx3(t) = 0,
1

2
< α ≤ 1, 0 < t < T, (6)

subject to the initial conditions

x(0) = A, Dαx(0) = 0. (7)

Here, the derivative Dα is interpreted in the Caputo sense.
The Operational Matrix Method (OMM) is a robust numerical approach widely em-

ployed for solving differential equations, particularly those arising in engineering and ap-
plied mathematics [42–49]. By converting differential operators into operational matrices,
OMM effectively transforms differential equations into algebraic matrix equations [43, 44].
This transformation eliminates the necessity for explicit analytical integration or sym-
bolic computation, streamlining the process of finding solutions. OMM is particularly
advantageous in scenarios where the solutions can be approximated using polynomials or
linear combinations of polynomial functions. It has found extensive applications in fields
such as control system analysis, structural dynamics, heat transfer, and fluid mechanics.
The method’s versatility and computational efficiency have made it an essential tool in
addressing a diverse array of problems across scientific and engineering disciplines. The
effectiveness of the Operational Matrix Method lies in its ability to accurately approxi-
mate solutions by converting differential or integral equations into a system of algebraic
equations. This is achieved through the use of operational matrices, which efficiently
handle the underlying mathematical operations. The method is especially well-suited for
problems where solutions can either be directly expressed in terms of polynomial functions
or approximated as linear combinations of such functions. The OMM approach simplifies
the computational process, reduces computational overhead, and enhances the accuracy
of the results. These advantages have made it a preferred technique in numerous sci-
entific and engineering domains. Applications of OMM include, but are not limited to,
control system design, modeling of structural dynamics, heat conduction problems, and
fluid flow analysis. The extensive use of OMM in these areas highlights its adaptability
and robustness in solving a broad spectrum of mathematical and physical problems. The



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 4 of 18

Operational Matrix Method stands out as a highly efficient numerical technique for solving
both linear and nonlinear differential equations, particularly those involving polynomial
forms, [45, 50–56]. Its ability to transform complex mathematical problems into simpler
algebraic systems has solidified its position as a valuable tool in applied mathematics and
engineering. The references [50–53] provide further insights into the diverse applications
of OMM, showcasing its reliability and effectiveness in tackling real-world challenges.

The organization of this paper is as follows. Section 2 outlines essential definitions and
lemmas that serve as the foundation for the subsequent analysis. In Section 3, we present
the modified iterative formulation of the Operational Matrix Method (OMM), detailing
its theoretical framework and implementation. Section 4 is devoted to establishing key
theoretical findings, including proofs for the existence and uniqueness of solutions, error
estimation, and a comprehensive convergence analysis. In Section 5, we illustrate the
practical applicability of the proposed method through three numerical examples, demon-
strating its convergence to the unique solution of Problem (6)-(7). The paper concludes
with Section 6, where we provide a summary of our findings and offer remarks on potential
future research directions.

2. Preliminarily

This section presents several essential definitions and results that are integral to the
analysis within this paper.

Definition 1. [1] For α ∈ (0, 1) and t > 0, the Caputo derivative (CD) of the function
x(t) is defined as

Dαx(t) =
1

Γ(1− α)

∫ t

0
(t− s)−αx′(s) ds, (8)

while the fractional integral (FI) operator is expressed as:

Iαx(t) =
1

Γ(α)

∫ t

0
(t− s)α−1x(s) ds. (9)

The subsequent lemma outlines key properties of the Caputo derivative and the fractional
integral.

Lemma 1. [1] For α ∈ (0, 1) and x(t) ∈ C[0, T ], the following relations hold:

IαDαx(t) = x(t)− x(0), (10)

DαIαx(t) = x(t). (11)

In this paper, we utilize the concept of the sequential Caputo derivative, defined as
follows.

Definition 2. [43] Let α ∈ (0, 1), and for a given function x(t), the sequential Caputo
derivative of order α is defined recursively as

cDLαx(t) = cDα
(
cD(L−1)αx(s)

)
, L ∈ N, (12)



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 5 of 18

where scDLαx(t) denotes the sequential Caputo derivative of order α. For ease of notation,
we will refer to it simply as DLαx(t) when discussing the sequential derivative.

Next, we introduce the Operational Matrix Method (OMM) with the following defini-
tion.

Definition 3. [42, 49] Let ts = s∆ for s ∈ {0, 1, 2, . . . ,M−1} and ∆ = T
M , where M ∈ N.

The s-block pulse function (BPF) is defined as:

µs(t) =

{
1, ts ≤ t < ts+1,

0, otherwise,
0 ≤ s < M.

The product and orthogonality properties of BPFs follow directly from the above
definition, as given in the next theorem.

Theorem 1. [43–45] Let {t0 = 0, t1, . . . , tM = T} be a uniform partition of [0, T ], and let
{µ0(t), µ1(t), . . . , µM−1(t)} be the corresponding BPFs. Then, for any 0 ≤ i, j ≤ M − 1,
the following relations hold:

µi(t)µj(t) =

{
µi(t), if i = j,

0, otherwise.
(13)

and ∫ T

0
µi(t)µj(t) dt =

{
∆, if i = j,

0, otherwise.
(14)

We conclude this section by stating the completeness property in the following lemma.

Lemma 2. [45–47] If x ∈ L2[0, T ), then we have the expansion:

x(t) =
∞∑
s=0

xsµs(t), (15)

where

xs =
1

∆

∫ (s+1)∆

s∆
x(r) dr. (16)

For numerical purposes, we approximate x(t) by:

x(t) ≈
M−1∑
s=0

xsµs(t), M ≫ 1. (17)



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 6 of 18

3. Iterative Operational Matrix Method

Syam et al. [42–44, 46–49], introduced an inefficient iterative method which reduce
the computational time and cost, and increase the effecinet of the OMM. In this section,
we will follow this approach to solve our problem (6)-(7). First, by applying the fractional
integral operator (9) to both sides of (6), we obtain the following equation:

Dαx(t)−Dαx(0) +
β

Γ(α)

∫ t

0
x(s)(t− s)α−1ds+

γ

Γ(α)

∫ t

0
x3(s)(t− s)α−1ds = 0.

Since Dαx(0) = 0, this simplifies to:

Dαx(t) = − β

Γ(α)

∫ t

0
x(s)(t− s)α−1ds− γ

Γ(α)

∫ t

0
x3(s)(t− s)α−1ds. (18)

Next, applying the fractional integral operator (9) again to both sides of equation (18),
we obtain

x(t)− x(0) = − β

Γ2(α)

∫ t

0

∫ s

0
x(w)(t− s)α−1(s− w)α−1dwds

− γ

Γ2(α)

∫ t

0

∫ s

0
x3(w)(t− s)α−1(s− w)α−1dw ds.

Since x(0) = A, we can rewrite this as:

x(t) = A− β

Γ2(α)

∫ t

0

∫ s

0
x(w)(t− s)α−1(s− w)α−1dw ds

− γ

Γ2(α)

∫ t

0

∫ s

0
x3(w)(t− s)α−1(s− w)α−1dwds. (19)

Let {t0 = 0, t1, . . . , tM = T} represent a uniform partition of [0, T ], and let {µ0(t), µ1(t),
. . . , µM−1(t)} be the corresponding block pulse functions (BPFs). We approximate x(t)
using these BPFs as follows

x(t) =

M−1∑
i=0

xiµi(t). (20)

By applying the collocation method to equation (20) with collocation points at tj , where
1 ≤ j < M , we get

x(tj) =

M−1∑
i=0

xiµi(tj) = xj , 0 ≤ j ≤ M − 1, (21)

which leads to the following equation:

xj = A− β

Γ2(α)

∫ tj

0

∫ s

0

(
M−1∑
i=0

xiµi(w)

)
(tj − s)α−1(s− w)α−1dw ds



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 7 of 18

− γ

Γ2(α)

∫ tj

0

∫ s

0

(
M−1∑
i=0

xiµi(w)

)3

(tj − s)α−1(s− w)α−1dw ds. (22)

This equation can be solved due to the fact that

µi(tj) =

{
1, if i = j

0, if i ̸= j
.

By utilizing the properties of the Riemann integral, equation (22) simplifies to

xj = A− β

Γ2(α)

j−1∑
k=0

∫ tk+1

tk

∫ s

0

(
M−1∑
i=0

xiµi(w)

)
(tj − s)α−1(s− w)α−1dw ds

− γ

Γ2(α)

j−1∑
k=0

∫ tk+1

tk

∫ s

0

(
M−1∑
i=0

xiµi(w)

)3

(tj − s)α−1(s− w)α−1dw ds. (23)

For any 1 ≤ k, j < M , we know that

µj(t) =

{
1, if j = k

0, if j ̸= k
, t ∈ [tk, tk+1).

This leads to the following simplifications∫ tk+1

tk

∫ s

0

(
M−1∑
i=0

xiµi(w)

)
(tj − s)α−1(s− w)α−1dw ds

=

∫ tk+1

tk

∫ s

0
xk(tj − s)α−1(s− w)α−1dw ds,

∫ tk+1

tk

∫ s

0

(
M−1∑
i=0

xiµi(w)

)3

(tj − s)α−1(s− w)α−1dw ds

=

∫ tk+1

tk

∫ s

0
x3k(tj − s)α−1(s− w)α−1dw ds.

Thus, we can express the equation for xj as

xj = A− β

Γ2(α)

j−1∑
k=0

ηjkxk −
γ

Γ2(α)

j−1∑
k=0

ηjkx
3
k, (24)

where the coefficients ηjk are defined as

ηjk =

∫ tk+1

tk

∫ s

0
(tj − s)α−1(s− w)α−1dw ds. (25)



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 8 of 18

4. Error Analysis

Let x ∈ L2([0, T )) be a function, and its norm is defined by the following formula:

||x|| =

√∫ T

0
|x(t)|2 dt. (26)

According to Lemma (2), the function x(t) can be approximated as

xM (t) =

M−1∑
i=0

xiµi(t). (27)

The aim of Theorem (2) is to demonstrate that the mean square error is minimized when
the coefficients xi are chosen according to Equation (16).

Theorem 2. Let x ∈ L2([0, T )) and define xM (t) as given in Equation (27). Then, the
mean square error term

A(x0, x1, . . . , xM−1) =

∫ T

0
(x(t)− xM (t))2 dt

is minimized when xi is computed according to Equation (16) for all i = 0, 1, . . . ,M − 1.

Proof. For 0 ≤ i ≤ M − 1, by applying Theorem (1), we compute

∂A
∂xi

= 2

∫ T

0
(x(t)− xM (t))µi(t) dt = 2

(∫ T

0
x(t)µi(t) dt−∆xi

)
= 0.

Therefore, we obtain

xi =
1

∆

∫ T

0
x(t)µi(t) dt.

Additionally, Theorem (1) provides the second-order derivatives as:

∂2A
∂xi∂xj

=

{
2∆, if i = j,

0, if i ̸= j.
0 ≤ i, j ≤ M − 1.

For any 0 ≤ i ≤ M − 1, we compute the determinant of the Hessian matrix:∣∣∣∣∣∣∣∣∣∣∣

∂2A
∂x2

0

∂2A
∂x0∂x1

. . . ∂2A
∂x0∂xi

∂2A
∂x1∂x0

∂2A
∂x2

1
. . . ∂2A

∂x1∂xi

...
...

. . .
...

∂2A
∂xi∂x0

∂2A
∂xi∂x1

. . . ∂2A
∂x2

i

∣∣∣∣∣∣∣∣∣∣∣
=

∣∣∣∣∣∣∣∣∣
2∆ 0 . . . 0
0 2∆ . . . 0
...

...
. . .

...
0 0 . . . 2∆

∣∣∣∣∣∣∣∣∣ = (2∆)i+1 > 0.

Thus, the error A is minimized when xi is computed according to Equation (16) for all
i = 0, 1, . . . ,M − 1.

Next, we aim to establish the order of the mean square error in approximating x(t)
over the interval [0, T ).



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 9 of 18

Theorem 3. Let x(t) be a differentiable function on [0, T ) such that

|x′(t)| ≤ F (28)

for all t ∈ [0, T ), where F is a positive constant. Then,

||A||2 ≤ C∆2 (29)

where A(x) = x(t)− xM (t), t ∈ [0, T ), and C is a positive constant, with ∆ = T
M .

Proof. Let ti = i∆ and Ai = [ti, ti+1), where ∆ = T
M and i = 0, 1, . . . ,M − 1. By

applying the mean value theorem for integrals and Equation (16), we obtain

xM (t) = xi, t ∈ [ti, ti+1),

=
1

∆

∫ ti+1

ti

x(t) dt,

= x(νi), νi ∈ [ti, ti+1), i = 0, 1, . . . ,M − 1.

Using the mean value theorem for integrals, we have:

||A||2 =
∫ T

0
(x(t)− xM (t))2 dx

=

M−1∑
i=0

∫ ti+1

ti

(x(t)− xM (t))2 dx

= ∆
M−1∑
i=0

(x(ωi)− x(νi))
2

where ωi, νi ∈ [ti, ti+1). By the mean value theorem and Equation (28), we get:

||A||2 ≤ ∆F2
M−1∑
i=0

(ωi − νi)
2

≤ ∆F2
M−1∑
i=0

∆2

= MF2∆3 = (TF2)∆2 = C∆2.

Thus, C = F2T .

The above result confirms that the mean square error in the approximation is of order
∆2.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 10 of 18

t HPM MHPM SHPM Present Mathematica

0.5 0.762476 0.768902 0.768766 0.768802 0.768802
1.0 0.176929 0.233741 0.233680 0.233692 0.233692
2.0 −1.055110 −0.891260 −0.859323 −0.859349 −0.859349
3.5 −0.461650 −0.079433 −0.093034 −0.093013 −0.093013
5.0 2.049041 0.996472 0.947107 0.947130 0.947130

Table 1: Comparison of approximate solutions for A = 1, α = 1, β = 1, and γ = 1.

t HPM MHPM SHPM Present Mathematica

1 0.056288 0.080176 0.080519 0.080527 0.080527
2 −0.808192 −0.739174 −0.729000 −0.729018 −0.729018
3 −0.339208 −0.239413 −0.238620 −0.238626 −0.238626
4 0.891267 0.706827 0.667953 0.668022 0.668022
5 0.893003 0.395315 0.387550 0.387551 0.387551

Table 2: Comparison of approximate solutions for A = 0.75, α = 1, β = 1.5, and γ = 1.5.

5. Numerical Results

We first analyze Problem (6)-(7) for the case where α = 1. Since the exact solution for
this problem is unavailable, numerical solutions were computed using the built-in Wolfram
Mathematica routine, employing the fully explicit Runge-Kutta method. These solutions
serve as the standard reference for comparison. In Tables 1 and 2, we present a comparison
of our results with those obtained using HPM [54], MHPM [26], and SHPM [54], alongside
the numerical solution. Table 1 corresponds to the parameters A = 1, α = 1, β = 1, and
γ = 1, while Table 2 pertains to A = 0.75, α = 1, β = 1.5, and γ = 1.5.

Figures 1–5 illustrate the comparisons between the current method and numerical
solutions under varying parameters. Specifically, Figure 1 evaluates A = 1, α = 1, β = 0.5,
and γ = 2; Figure 2 examines A = 1.5, α = 1, β = 1, and γ = 0.5; Figure 3 considers
A = 1.5, α = 1, β = 0.5, and γ = 1.5; Figure 4 studies A = 2, α = 1, β = 1.5, and γ = 1;
and Figure 5 assesses A = 1.5, α = 1, β = 1.5, and γ = 1.5.

Finally, the behavior of Problem (6)-(7) is investigated for A = 1, β = 0.5, and
γ = 2 with varying values of α (α = 0.75, 0.85, 0.9, 0.95, 0.995, 1). Results are de-
picted in Figures 6 and 7. Similarly, for A = 1.5, β = 1.5, and γ = 1.5, outcomes for
α = 0.75, 0.85, 0.9, 0.95, 0.995, 1 are illustrated in Figures 8 and 9. As α approaches 1, the
results converge toward the solution of the classical second-order differential equation.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 11 of 18

0 1 2 3 4 5

-1.0

-0.5

0.0

0.5

1.0

t

F
Ht
L

Figure 1: A = 1, α = 1, β = 0.5, and γ = 2. Numerical (solid) vs. present method (dashed).

0 1 2 3 4 5

-1.0

-0.5

0.0

0.5

1.0

t

F
Ht
L

Figure 2: A = 1.5, α = 1, β = 1, and γ = 0.5. Numerical (solid) vs. present method (dashed).

0 2 4 6 8 10

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

t

F
Ht
L

Figure 3: A = 1.5, α = 1, β = 0.5, and γ = 1.5. Numerical (solid) vs. present method (dashed).



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 12 of 18

2 4 6 8 10

-2

-1

1

2

Figure 4: A = 2, α = 1, β = 1.5, and γ = 1. Numerical (solid) vs. present method (dashed).

2 4 6 8 10

-1.5

-1.0

-0.5

0.5

1.0

1.5

Figure 5: A = 1.5, α = 1, β = 1.5, and γ = 1.5. Numerical (solid) vs. present method (dashed).

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=0.75

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=0.85

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=0.9

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=0.95

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=995

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Α=1

Figure 6: A = 1, β = 0.5, and γ = 2 for α = 0.75, 0.85, 0.9, 0.95, 0.995, 1.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 13 of 18

1 2 3 4 5

-1.0

-0.5

0.5

1.0

Figure 7: Results for A = 1, β = 0.5, γ = 2, α = 0.995, 1.

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=0.75

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=0.85

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=0.9

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=0.95

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=995

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Α=1

Figure 8: Results for A = 1.5, β = 1.5, γ = 1.5, α = 0.75, 0.85, 0.9, 0.95, 0.995, 1.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 14 of 18

1 2 3 4 5

-1.5

-1.0

-0.5

0.5

1.0

1.5

Figure 9: Results for A = 1.5, β = 1.5, γ = 1.5, α = 0.995, 1.

6. Conclusion

In this work, we developed and analyzed an iterative fractional operational matrix
method, based on the method described in [13], for solving the nonlinear fractional un-
damped Duffing equation. The proposed method efficiently reduces computational com-
plexity while maintaining high accuracy, as demonstrated through comparisons with ex-
isting analytical and numerical techniques. The results confirm that the method provides
an excellent approximation, outperforming classical perturbation-based approaches such
as the Homotopy Perturbation Method (HPM), Modified HPM (MHPM), and Simplified
HPM (SHPM).

A key contribution of this study is the investigation of the effect of varying the frac-
tional order α beyond the provided range. Our findings reveal that decreasing α leads to a
significant deviation from the classical Duffing oscillator, introducing memory-dependent
effects and altering the system’s response dynamics. As α approaches unity, the system
gradually transitions to the standard integer-order Duffing equation, confirming the ro-
bustness of fractional modeling in capturing real-world physical behaviors. This insight
is crucial for applications in nonlinear dynamics, where fractional-order models provide a
more generalized framework for analyzing oscillatory systems with memory effects.

Future research directions include extending the method to more complex nonlin-
ear oscillatory systems, exploring higher-dimensional fractional models, and investigating
real-world applications such as mechanical vibrations and bioengineering systems. Addi-
tionally, further refinement of the numerical scheme may improve computational efficiency
for large-scale problems.

Author Contributions:

All authors have the same contribution.

Data Availability Statement:

Not applicable.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 15 of 18

Conflicts of Interest:

The authors declare no conflict of interest.

References

[1] Keith Oldham and Jerome Spanier. The fractional calculus theory and applications
of differentiation and integration to arbitrary order. Elsevier, 1974.

[2] Radek Matuš. Application of fractional order calculus to control theory. International
Journal of Mathematical Models and Methods in Applied Sciences, 5(7):1162–1169,
2011.

[3] R. L. Magin. Fractional calculus models of complex dynamics in biological tissues.
Comput. Math. Appl., 59(5):1586–1593, 2010.

[4] M. A. Matlob and Y. Jamali. The concepts and applications of fractional order
differential calculus in modeling of viscoelastic systems: A primer. Crit. Rev. Biomed.
Eng., 47(4), 2019.

[5] J. Yu, L. Tan, S. Zhou, L. Wang, and M. A. Siddique. Image denoising algorithm
based on entropy and adaptive fractional order calculus operator. IEEE Access,
5:12275–12285, 2017.

[6] K. B. Kachhia. Chaos in fractional order financial model with fractal–fractional
derivatives. Partial Differ. Equ. Appl. Math., 7:100502, 2023.

[7] Y. A. Rossikhin and M. V. Shitikova. Applications of fractional calculus to dynamic
problems of linear and nonlinear hereditary mechanics of solids, 1997.

[8] A. Atangana and D. Baleanu. New fractional derivatives with nonlocal and non-
singular kernel: Theory and application to heat transfer model. arXiv preprint, 2016.

[9] B. Ghanbari, H. Günerhan, and H. M. Srivastava. An application of the Atangana-
Baleanu fractional derivative in mathematical biology: A three-species predator-prey
model. Chaos Solitons Fractals, 138:109910, 2020.

[10] J. F. Gómez-Aguilar, R. F. Escobar-Jiménez, M. G. López-López, and V. M.
Alvarado-Mart́ınez. Atangana-Baleanu fractional derivative applied to electromag-
netic waves in dielectric media. J. Electromagn. Waves Appl., 30(15):1937–1952,
2016.

[11] E. F. D. Goufo, M. Mbehou, and M. M. K. Pene. A peculiar application of Atan-
gana–Baleanu fractional derivative in neuroscience: Chaotic burst dynamics. Chaos
Solitons Fractals, 115:170–176, 2018.

[12] B. Ghanbari and A. Atangana. A new application of fractional Atangana–Baleanu
derivatives: Designing ABC-fractional masks in image processing. Physica A: Stat.
Mech. Appl., 542:123516, 2020.

[13] S. M. Syam, Z. Siri, S. H. Altoum, and R. Md. Kasmani. Analytical and numerical
methods for solving second-order two-dimensional symmetric sequential fractional
integro-differential equations. Symmetry, 15(6):1263, 2023.

[14] N. N. Moiseev. Asimptoticheskie metodi nelinejnoj mehaniki. Nauka, Moscow, 1981.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 16 of 18

[15] M. Al-Refai and D. Baleanu. On an extension of the operator with mittag-leffler
kernel. Fractals, 30(05):2240129, 2022.

[16] Mahmmoud M. Syam, Farah Morsi, Ayaha Abu Eida, and Muhammed I. Syam. In-
vestigating convective darcy–forchheimer flow in maxwell nanofluids through a com-
putational study. Partial Differential Equations in Applied Mathematics, 11:100863,
2024.

[17] M. M. Syam, S. Cabrera-Calderon, K. A. Vijayan, V. Balaji, P. E. Phelan, and J. R.
Villalobos. Mini containers to improve the cold chain energy efficiency and carbon
footprint. Climate, 10:76, 2022.

[18] A. H. Nayfeh and D. T. Mook. Nonlinear Oscillations. Wiley, New York, 1979.
[19] M. I. Syam, M. N. Y. Anwar, A. Yildirim, and et al. The modified fractional power

series method for solving fractional non-isothermal reaction–diffusion model equations
in a spherical catalyst. Int. J. Appl. Comput. Math., 5:38, 2019.

[20] S. Al Omari, A. M. Ghazal, M. Syam, R. Al Najjar, and M. Y. Selim. An investigation
on the thermal degradation performance of crude glycerol and date seeds blends using
thermogravimetric analysis (tga). In 5th Int. Conf. Renewable Energy: Generation
and Appl. (ICREGA 2018), pages 102–106, 2018.

[21] M. I. Syam, M. A. Raja, M. M. Syam, and et al. An accurate method for solving
the undamped duffing equation with cubic nonlinearity. Int. J. Appl. Comput. Math.,
4:69, 2018.

[22] A.-H. I. Mourad, A. M. Ghazal, M. M. Syam, O. D. Al Qadi, and H. Al Jassmi. Uti-
lization of additive manufacturing in evaluating the performance of internally defected
materials. IOP Conf. Ser.: Mater. Sci. Eng., 362:012026, 2018.

[23] M. Al-Refai, M. Syam, and D. Baleanu. Analytical treatment to systems of frac-
tional differential equations with modified atangana-baleanu derivative. Fractals,
31:2340156, 2023.

[24] G. Duffing. Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre
technische Bedeutung. Druck und Verlag von Fridr. Vieweg & Sohn, Braunschweig,
1918.

[25] R. H. Abraham and C. D. Shaw. Dynamics – the geometry of behavior, Part I. Aerial
Press, Santa Cruz, 2000.

[26] D. D. Ganji, A. R. Sahouli, and M. Famouri. A new modification of he’s homotopy
perturbation method for rapid convergence of nonlinear undamped oscillators. J.
Appl. Math. Comput., 30:181–192, 2009.

[27] J. H. He. Homotopy perturbation method: A new nonlinear analytical technique.
Appl. Math. Comput., 135:73–79, 2003.

[28] Muhammad Asif Zahoor Raja, Saleem Abbas, Muhammed Ibrahem Syam, and Ab-
dul Majid Wazwaz. Design of neuro-evolutionary model for solving nonlinear sin-
gularly perturbed boundary value problems. Applied Soft Computing, 62:373–394,
2018.

[29] J. I. Ramos. On the variational iteration method and other iterative techniques for
nonlinear differential equations. Appl. Math. Comput., 199:39–69, 2008.

[30] M. I. Syam and H. I. Siyyam. Numerical differentiation of implicitly defined curves.



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 17 of 18

J. Comput. Appl. Math., 108(1-2):131–144, 1999.
[31] M. Syam. The modified broyden-variational method for solving nonlinear elliptic

differential equations. Chaos, Solitons & Fractals, 32(2):392–404, 2007.
[32] Bothayna S.H. Kashkaria and Muhammed I. Syam. Evolutionary computational intel-

ligence in solving a class of nonlinear volterra–fredholm integro-differential equations.
Journal of Computational and Applied Mathematics, 311:314–323, 2017.

[33] M. F. El-Sayed and M. I. Syam. Electrohydrodynamic instability of a dielectric com-
pressible liquid sheet streaming into an ambient stationary compressible gas. Archive
of Applied Mechanics, 77:613–626, 2007.

[34] H. M. Jaradat, M. Alquran, and M. I. Syam. A reliable study of new nonlinear
equation: Two-mode kuramoto–sivashinsky. International Journal of Applied and
Computational Mathematics, 4:64, 2018.

[35] S. Bourazza, Sami H. Altoum, and et al. Discharging process within a storage con-
tainer considering numerical method. Journal of Energy Storage, 66:107490, 2023.

[36] Sami H. Altoum and et al. Efficacy of magnetic force on nanofluid laminar transporta-
tion and convective flow. Journal of Magnetism and Magnetic Materials, 581:170964,
2023.

[37] Mahmmoud M. Syam, Mohammad Alkhedher, and Muhammed I. Syam. Thermal
and hydrodynamic analysis of mhd nanofluid flow over a permeable stretching surface
in porous media: Comparative study of fe3o4, cu, and ag nanofluids. International
Journal of Thermofluids, 26:101055, 2025.

[38] Mahmmoud M. Syam and Muhammed I. Syam. Impacts of energy transmission
properties on non-newtonian fluid flow in stratified and non-stratified conditions.
International Journal of Thermofluids, 23:100824, 2024.

[39] R. E. Mickens. Nonlinear Oscillations. Cambridge University Press, New York, 1981.
[40] Basem S. Attili, Khalid Furati, and Muhammed I. Syam. An efficient implicit

runge–kutta method for second order systems. Applied Mathematics and Compu-
tation, 178(2):229–238, 2006.

[41] L. Cveticanin. Analytic approach for the solution of the complex-valued strong non-
linear differential equation of duffing type. Physica A, 297:348–360, 2001.

[42] S. M. Syam, Z. Siri, S. H. Altoum, M. A. Aigo, and R. Md. Kasmani. A new method
for solving physical problems with nonlinear phoneme within fractional derivatives
with singular kernel. J. Comput. Nonlinear Dyn., 19, 2024.

[43] S. M. Syam, Z. Siri, S. H. Altoum, and R. Md. Kasmani. An efficient numerical
approach for solving systems of fractional problems and their applications in science.
Mathematics, 11, 2023.

[44] S. M. Syam, Z. Siri, R. M. Kasmani, and K. Yildirim. A new method for solving
sequential fractional wave equations. J. Math., 2023.

[45] M. I. Syam, M. Sharadga, and I. Hashim. A numerical method for solving fractional
delay differential equations based on the operational matrix method. Chaos, Solitons
& Fractals, 147:110977, 2021.

[46] S. M. Syam, Z. Siri, and R. M. Kasmani. Operational matrix method for solving
fractional system of riccati equations. In 2023 International Conference on Frac-



M. I. Syam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5884 18 of 18

tional Differentiation and Its Applications (ICFDA), pages 1–6, Ajman, United Arab
Emirates, 2023.

[47] S. M. Syam, Z. Siri, S. H. Altoum, M. A. Aigo, and R. M. Kasmani. A novel study
for solving systems of nonlinear fractional integral equations. Appl. Math. Sci. Eng.,
31(1), 2023.

[48] M. I. Syam, H. Alahbabi, R. Alomari, S. M. Syam, S. M. Hussein, S. Rabih, and
N. Al Saafeen. A numerical study for fractional problems with nonlinear phenomena
in physics. Progress in Fractional Differentiation and Applications, 10(3):451–461,
2024.

[49] L. Abdelhaq, S. M. Syam, and M. I. Syam. An efficient numerical method for two-
dimensional fractional integro-differential equations with modified atangana–baleanu
fractional derivative using operational matrix approach. Partial Differ. Equ. Appl.
Math., 11, 2024.

[50] M. I. Syam and M. Al-Refai. Fractional differential equations with atangana–baleanu
fractional derivative: Analysis and applications. Chaos Solitons Fractals X, 2, 2019.

[51] M. M. Syam and M. I. Syam. Computational study of magnetohydrodynamic squeeze
flow between infinite parallel disks. Int. J. Thermofluids, 24, 2024.

[52] A.K. Alomari, Muhammed I. Syam, N.R. Anakira, and A.F. Jameel. Homotopy
sumudu transform method for solving applications in physics. Results in Physics,
18:103265, 2020.

[53] M. M. Syam and M. I. Syam. Investigation of slip flow dynamics involving al2o3 and
fe3o4 nanoparticles within a horizontal channel embedded with porous media. Int. J.
Thermofluids, 24, 2024.

[54] Precious Sibanda and A. Khidir. A new modification of the hpm for the duffing
equation with cubic nonlinearity. In Recent Researches in Applied and Computational
Mathematics, pages 139–145, 2011.

[55] A. Verma, W. Sumelka, and P. K. Yadav. The numerical solution of nonlinear frac-
tional lienard and duffing equations using orthogonal perceptron. Symmetry, 15:1753,
2023.

[56] H. Singh and H. M. Srivastava. Numerical investigation of the fractional-order liénard
and duffing equations arising in oscillating circuit theory. Frontiers in Physics, 8:120,
2020.