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

Numerical Solutions of the SIR Mathematical Model of
Computer Viruses Involving Non-linear Fractional

Order Differential Equation

Rajratana M. Kamble1,∗, Pramod R. Kulkarni2

1 Shri Vitthal Rukmini Arts, Commerce and Science College, Sawana, Tq. Mahagoan, Dist.
Yavatmal, Maharashtra, India
2 N.E.S Science College, Nanded, Maharashtra, India

Abstract. In this paper, the non-linear Fractional order susceptible-infected-Recovered(SIR) math-
ematical model of computer viruses is presented. For this Fractional differential transform method
(FDTM) and the Laplace-Adomian decomposition method (LADM) are applied and compared
the obtained results with Runge-Kutta Felhberg Method for γ = 1. Also, graphs of solutions are
plotted up to five iterations from both the method. The fractional derivative used in the model is
Caputo Fractional derivative. Using Lyapunov stability analysis,the stability verified of the given
mathematical model. Solutions are obtained for three different fractional order values of γ. The
validity of the result is confirmed by converting model into integer order. Errors from both the
methods are compared.

2020 Mathematics Subject Classifications: 34A08

Key Words and Phrases: Caputo fractional derivative, fractional differential transform method,
Laplace Adomian decomposition method, Lyapunov stability, Runge-Kutta Felhberg method

1. Introduction

Computer viruses are harmful programs for normal functioning of computers. There
are various types of computer viruses that attack computers. Different types of mathemat-
ical models representing spread of computer viruses have been proposed in the literature[1]
[1],[2],[3],[4],[5],[6],[7],[8],[9],[10],[11],[12],[13],[14],[15],[16][17]. These models can be solved
using mathematical techniques such as the Differential Transform Method, Collection
Method [13], Homotopy Analysis Method [2],[15] Variational Iterational Method[18], and
others[7],[16]. Various mathematical and engineering problems are solved using semi-
analytical techniques [14],[19],[20],[21],[15],[22],[23],[24],[25]. The Differntial Transform
Method and Laplace Adomian Decomposition Method are the tools to solve linear and non-
linear preoblems in engineering, physics, mathematics, etc [25],[26],[27],[28],[29],[30],[31]

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

Email addresses: kamblerajratna2@gmail.com (R. M. Kamble), pramodrkul@gmail.com (P. R. Kulkarni)

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



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 2 of 40

,[32],[33],[34],[35],[36],[37],[38],[39],[40],[41],[42],[43],[44],[45],[46]
In this paper, we have designed and solved the non-linear fractional order system of the
Susceptible-Infected-Recovered (SIR) mathematical model of computer viruses. For this
purpose, the Fractional Differential Transform Method (FDTM) and Laplace Adomian
Decomposition Method (LADM) is applied. The approximate solutions are estimated for
different iterations (up to five iterations), and the results are compared with those ob-
tained by 4/5 Runge-Kutta Frlhberg Algorithm (RKFA). The graphs of the approximate
solutions up to five iterations and the absolute errors are plotted at different values of γ
and t, illustrating the evolution of the system. Comparing the solutions and errors ob-
tained by FDTM and LADM with RKFA, we have concluded that as compared to the
LADM, the FDTM gives more accurate results under certain conditions. The methods
and results presented in this chapter on the non-linear fractional order system for modeling
computer virus propagation through the Susceptible-Infected-Recovered (SIR) framework
have broader implications for mathematical modeling and systems analysis. The use of
non-linear fractional order models is not limited to computer viruses; similar models can be
applied to biological systems such as the spread of infectious diseases in populations. The
methodologies outlined in this chapter can be adapted to various epidemiological contexts,
potentially improving the understanding of disease dynamics and control measures. The
methods employed, particularly the Fractional Differential Transform Method (FDTM)
and Differential Transform Method (DTM), can be integrated with other mathematical
techniques, such as Agent-based Modeling or Network Theory to provide a more com-
prehensive analysis of complex systems. This versatility allows the exploration of various
dynamics and interactions in different contexts. The mathematical framework developed
can be adapted for various applications. The most common and general model is the clas-
sical non-linear fractional order Susceptible-Infected-Recovered (SIR) model for computer
virus propagation is given below.

2. The S-I-R Model and Parameters

The governing equations of the non-linear fractional-order SIR model[47] are given by:

dγs(t)

dtγ
= f1 − λs(t)i(t)− ds(t),

dγi(t)

dtγ
= f2 + λs(t)i(t)− εi(t)− dr(t),

dγr(t)

dtγ
= f3 + εi(t)− dr(t).

(1)

with the initial conditions

s(0) = s0, i(0) = i0, r(0) = r0 (2)

The parameters and initial values of system 1 are given in Table 1.
Here the fractional derivative used is the Caputo fractional derivative and γ ∈ (0, 1].

The generalization of integer order differentiation is known as fractional order differentia-
tion. As the reader can see, calculating derivatives of fractional order is far more difficult



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 3 of 40

Table 1: List of parameters and their values

Parameter Meaning Value

s(t) Susceptible computers to the virus at time t s(0) = 20
i(t) Infected computers to the virus at time t i(0) = 15
r(t) Recovered computers from the virus at time t r(0) = 10
f1, f2, f3 Rate of other computers connecting to the network 0
λ Rate of virus infection for susceptible computers 0.001
ε Rate of recovery from virus for infected computers 0.1
d Rate of removing from the network 0.1

than computing derivatives of traditional integer order. Nevertheless, it has been demon-
strated that fractional order mathematical models many phenomena’ attributes are better
described by integrals and derivatives than by the previously employed integer order mod-
els.This is because systems are typically not flawless and can be influenced by outside
factors, for instance. Consequently, it might not be possible to comprehend the trajecto-
ries of state variables using derivatives of integer order. With fractional derivatives, we
may choose whatever fractional differential equation best captures the model’s dynamics
since we have an infinite number of derivative orders at our disposal. Among other places,
[48]and [49] exhibit experimental data and techniques for certain real-world events demon-
strating that fractional order derivatives offer more effective solution curve modelling.
In [50], Kamble, R., and Kukarni, P. have proved the existence and uniqueness of solutions
for the following equation

DαDβx (τ) = f (t, x (τ) , ϕx (τ) , ψx (τ)) , τ ∈ [0, 1] , x (0) = x (1) = 0

where 0 < α ≤ 1, 0 < β ≤ 1, Dα, Dβ are the Caputo fractional derivatives of order α , β
respectively, f : [0, 1] × R3−→R is a continuous function λ, δ : [0, 1] × [0, 1] → [0,+∞),
and

ϕx (τ) =

∫ τ

0
λ (τ, s)x (s) ds, ψx (τ) =

∫ τ

0
δ (τ, s)x (s) ds

ϕ∗ = sup
t∈[0,1]

∣∣∣∣∫ t

0
λ (τ, s) ds

∣∣∣∣ <∞, ψ∗ = sup
t∈[0,1]

∣∣∣∣∫ t

0
δ (τ, s) ds

∣∣∣∣ <∞

In [51], Kamble, R., and Kukarni, P. have proved the existence and uniqueness of solutions
for the conformable fractional order Volterra-Fredholm type integro-differentail equations
of the form

dαx(t)

dt
= f (t) +

∫ t

0
p (t, s, x (s)) ds+

∫ b

0
q (t, s, x (s)) ds

with the initial condition
x (0) = x0

where the term dαx(t)
dt represents the conformable fractional order derivative of fractional

order α ∈ (0, 1),t ∈ I = [0, b] , f : I → X, p, q : I × I ×X → X are continuous functions



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 4 of 40

and x0 is an element of a real Banach space X, with the norm ∥.∥.
In [52], Kamble, R., and Kukarni, P. have given new definition and Prove there properties

Definition 1. Let f : [0,∞) → R be a real valued function, α ∈ (0, 1], a > 0 ∈ R, a ̸= 1.
We define the (α, a)-Extended fractional derivative(EFD) of f at x, denoted Dα

a f(x), by

Dα
a f(x) = lim

h→0

f(xahx
(−α)

)− f(x)

h
, (3)

provided the limit exists.

.
In [53], Kamble, R., and Kukarni, P. have proved, Existence and uniqueness of solutions
for exponential fractional differential equations.
In [54], Kamble, R., and Kukarni, P. have developed method to solve, (α, 1) Fractional Dif-
ferential Difference Equations With Conditions,Linear And Nonlinear with help of Laplace
Transform And Laplace Decomposition Method.

3. Lyapunov Stability analysis of the given non Linear SIR model of
Computer Viruses.

Using lyapunov stability theory we have investigated the stability of given SIR frac-
tional order non linear mathematical model is stable.
The detail about lyapunov stability analysis of caputo fractional order nonliner system is
given in [55]. Consider the following non linear system

Dγ
t0
Y(t) = g(t,Y(t)) = AY(t) + h(t,Y(t)),Y(t0) = Y0

where γ ∈ (0, 1),Y ∈ Rn called the state vector of the non linear system.
Here A ∈ Rn×n is the constant matrix of the given system, and h : R+ ×Rn → Rn is the
non linear part of the system, where h(t, 0) = 0 for all t ≥ 0. The corresponding linear
system is

Dγ
t0
Y(t) = AY(t),Y(t0) = Y0 (4)

The detailed stability of the above linear system is given in [55].
Now consider the following non-linear system

Dγ
t0
Y(t) = AY(t) +Bu+ h(t,Y(t)),Y(t0) = Y0

where B ∈ Rn×p is the constant matrix of the given system,u ∈ Rp is the control input to
be defined.
Writing the given SIR mathematical model of Computer Viruses in the form

Dγ
t0
Y(t) = AY(t) +Bu+ h(t,Y(t)),Y(t0) = Y0



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 5 of 40

where A =

−d 0 0
0 −ϵ −d
0 ϵ −d

, B =

f1f2
f3

, and Y(t) =
[
s(t) i(t) r(t)

]
and h(t,Y(t)) =[

−λs(t)i(t), λs(t)i(t), 0
]
.

Since B =

f1f2
f3

=
00
0

, the system convert into

Dγ
t0
Y(t) = AY(t) + h(t,Y(t)),Y(t0) = Y0 (5)

Definition 2. [55] The non linear function h(t,Y) is said to be one sided Lipschitz if
there exist ρ ∈ R such that
⟨h(t,Y1)− h(t,Y2),Y1 −Y2⟩ ≤ ρ(∥Y1 −Y2∥)2

Theorem 1. [55] Let Y = 0be an equilibrium point of the system 4. If the state ma-
trix A is Hurwitz then the trivial solution of the fractional linear system 4 is fractional
asymptotically stable.

Theorem 2. [55] Let Y = 0 be an equilibrium point of the system 5. Let that the state
matrix A is Hurwitz and the condition ∥h(t,Y(t))∥ < ρ holds.
If there exist a positive definite matrix P such that the following inequality holds
ρ < λmin(Q)

2λmax(P ) where λ is eigen value of matrix.

where ATP + PA = −Q, then the trivial solution of the fractional linear system 5 is
fractional asymptotically stable.

Theorem 3. [55] If the function h(t,Y(t)) is Lipschitz continuous,L is Lipschitz Con-
stant. Assume that the following assumption is hold : Then There exists a positive sym-
metric matrix P and positive constant ϵ such that the following inequalities hold ATP +
PA+ ϵI < 0 and L < ϵ

2λmaxP
where λ is eigen value of matrix, then the trivial solution of

fractional system 5 is asymptotically stable.

The given SIR mathematical model of Computer Viruses can be written in the form

Dγ
t0
Y(t) = AY(t) +Bu+ h(t,Y(t)),Y(t0) = Y0

where A =

−d 0 0
0 −ϵ −d
0 ϵ −d

, B =

f1f2
f3

, Y(t) =
[
s(t) i(t) r(t)

]
and h(t,Y(t)) =

[
−λs(t)i(t), λs(t)i(t), 0

]
Since B =

f1f2
f3

=
00
0

, the system convert into

Dγ
t0
Y(t) = AY(t) + h(t,Y(t)),Y(t0) = Y0

For the given system A =

−d 0 0
0 −ϵ −d
0 ϵ −d

=
−0.1 0 0

0 −0.1 −0.1
0 0.1 0.1


A is Hurwitz matrix since the eigen values are λ1 = −0.1,



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 6 of 40

λ2 = 0.1 + 0.44721359549996i, λ3 = 0.1− 0.44721359549996i
We can find positive definite matrices P and Q as given below where

P =

10 0 0
0 10 0
0 0 10

, Q =

2 0 0
0 2 0
0 0 2

, ρ = 1
10 .

With these values, all the conditions of theorem 2 are satisfied. Also, with the P given as
above andϵ = 1 L = 1

20 all the condition of theorem 3 is satisfied. Therefore by theorem
3 the given fractional non linear system of SIR model is also asymptotically stable.

4. Numerical Methods

The Fractional Differential Transform Method (FDTM) is used to solve the non-linear
fractional order system of the Susceptible-Infected-Recovered mathematical model of com-
puter viruses. In this section, we apply FDTM to solve the system of equations 1. The
numerical solutions are compared with those obtained using the Differential Transform
Method (DTM), verifying that our proposed results are valid when fractional order deriva-
tive γ = 1. Additionally, we have taken five iterations of the solutions (k = 5) and compute
the values of s(t), i(t), and r(t) for different time points.
The Laplace Adomian Decomposition method (LADM) is also used to solve the non-linear
fractional order system of the Susceptible-Infected-Recovered mathematical model of com-
puter viruses. In this section, we apply LADM to solve the system of equations 1. The
numerical solutions are compared, verifying that our proposed results are valid when frac-
tional order derivative γ = 1. Additionally, we have taken five iterations of the solutions
(k = 5) and compute the values of s(t), i(t), and r(t) for different time points.
The Runge-Kutta Felhberg Algorithm (RKFA) is also used to solve the non-linear frac-
tional order system of the Susceptible-Infected-Recovered mathematical model of computer
viruses for fractional order derivative γ = 1.
The results are visually presented through graphs that illustrate the evolution of the sus-
ceptible, infected, and recovered populations over time for both methods. A comparison of
results obtained from FDTM and LADM is performed to analyze differences and evaluate
the effectiveness of the fractional order model in capturing the dynamics of computer virus
propagation. Future work could explore higher-order fractional derivatives, incorporate
more complex virus spread scenarios, or apply the model to different types of computer
security threats.

5. Fractional Differential Transform Method

The steps in the fractional differential transform method are given in [56]. The frac-
tional differentiation in the Riemann-Liouville sense is defined as:

Dγ
t0
u(t) =

1

(n− γ)

∫ t

t0

dm

dxm
(t− s)n−γ−1u(s)ds , (6)



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 7 of 40

where n − 1 ≤ γ < m, n is a positive integer, and t > t0. The expansion of an analytic
and continuous function f(t) in terms of a fractional power series is given as:

u(t) =
∞∑
k=0

U(k)(t− t0)
k/µ, (7)

where µ is the order of the fraction, and U(k) is the fractional differential transform of
u(t). The relation between the Riemann-Liouville operator and the Caputo operator is
given as follows:

CDγ
x0
u(t) = Dγ

x0

[
u(t)−

n−1∑
k=0

uk(t0)

k!
(t− t0)

k

]
. (8)

Here, the notation C is used for the Caputo operator to distinguish it from the Riemann-
Liouville operator.

Since the initial conditions are implemented to the integer order differential equation,
the transformed initial conditions can be obtained using the following relation:

U(k) =

 1
(k/µ)!

[
d(k/µ)u(t)

dt(k/µ)

]
x=x0

, if k/µ ∈ Z+,

0, if k/µ /∈ Z+.
(9)

where k = 0, 1, 2, . . . , γµ−1, and γ is the order of the fractional differential equation. The
parameter µ is chosen such that γµ is a positive integer.

The following results[56] are applied in the proofs of the results in the next section.

Theorem 4. [56] If w(t) = u(t)±v(t), then its transformed form is W (k) = U(k)±V (k).

Theorem 5. [56] If w(t) = u(t)v(t), then W (k) =
∑k

l=0 U(l)V (k − l).

Theorem 6. [56] If w(t) = (t− t0)
p, then W (k) = δ(k − λp), where:

δ(k) =

{
1, k = 0,

0, k ̸= 0.
(10)

Theorem 7. [56] If w(t) = Dγ
t0
u(t), then:

W (k) =
Γ(γ + 1 + k/µ)

Γ(1 + k/µ)
U(k + γµ). (11)



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 8 of 40

6. Methodology of the Proposed Non-Linear Fractional Mathematical
Model

Using the transformation rules from equations theorem4–theorem7, the transformed
equations are

S(k + γµ) =
Γ(1 + k/µ)

Γ(1 + γ + k/µ)

[
f1δ(k)− λ

k∑
l=0

S(l)I(k − l)− dS(k)

]

I(k + γµ) =
Γ(1 + k/µ)

Γ(1 + γ + k/µ)

[
f2δ(k) + λ

k∑
l=0

S(l)I(k − l)− εI(k)− dR(k)

]

R(k + γµ) =
Γ(1 + k/µ)

Γ(1 + γ + k/µ)
[f3δ(k) + εI(k)− dR(k)]

(12)

The initial conditions transform into:

S(k) = 0, k = 1, 2, ..., γµ− 1

I(k) = 0, k = 1, 2, ..., γµ− 1

R(k) = 0, k = 1, 2, ..., γµ− 1,

S(k) = 20, I(k) = 15, R(k) = 10, k = 0.

(13)

The approximate series solution for system 1 is:

s(t) =
n∑

j=0

S(j)tj/µ, i(t) =
n∑

j=0

I(j)tj/µ, r(t) =
n∑

j=0

R(j)tj/µ (14)

Case I: γ = 1/10
For γ = 1/10, setting µ = 10, using 12, we obtain

S(k + 10γ) =
Γ(1 + k/10)

Γ(1 + γ + k/10)

[
f1δ(k)− λ

k∑
l=0

S(l)I(k − l)− dS(k)

]

I(k + 10γ) =
Γ(1 + k/10)

Γ(1 + γ + k/10)

[
f2δ(k) + λ

k∑
l=0

S(l)I(k − l)− εI(k)− dR(k)

]

R(k + 10γ) =
Γ(1 + k/10)

Γ(1 + γ + k/10)
[f3δ(k) + εI(k)− dR(k)]

Substituting the initial conditions S(0) = 20, I(0) = 15, R(0) = 10, we obtain the values
of S(k), I(k), R(k), 0 ≤ k ≤ 4 as follows. Thus, the approximate solution becomes

s(t) = 20− 2.4177t1/10 + 0.3350t2/10 − 0.04716t3/10 + 0.005841t4/10

i(t) = 15− 2.3126t1/10 + 0.0996t2/10 + 0.03319t3/10 − 0.004836t4/10 (15)

r(t) = 10 + 0.5255t1/10 − 0.2981t2/10 + 0.04037t3/10 − 0.00070805t4/10



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 9 of 40

Table 2: Values of S(k), I(k), R(k), γ = 1/10.

Sr. No. k S(k) I(k) R(k)

1 0 20 15 10
2 1 −2.4177 −2.3126 0.5255
3 2 0.3350 0.0996 −0.2981
4 3 −0.04716 0.03319 0.04037
5 4 0.005841 −0.004836 −0.00070805

Case II: γ = 1/2
For γ = 1/2, setting µ = 2, using (12), we obtain

S(k + 2γ) =
Γ(1 + k/2)

Γ(1 + γ + k/2)

[
f1δ(k)− λ

k∑
l=0

S(l)I(k − l)− dS(k)

]

I(k + 2γ) =
Γ(1 + k/2)

Γ(1 + γ + k/2)

[
f2δ(k) + λ

k∑
l=0

S(l)I(k − l)− εI(k)− dR(k)

]

R(k + 2γ) =
Γ(1 + k/2)

Γ(1 + γ + k/2)
[f3δ(k) + εI(k)− dR(k)]

Using the initial conditions S(0) = 20, I(0) = 15, R(0) = 10, we obtain the values of
S(k), I(k), R(k), 0 ≤ k ≤ 4 as follows.
Thus, the approximate solution takes the form

Table 3: Values of S(k), I(k), R(k), γ = 1/2.

Sr. No. k S(k) I(k) R(k)

1 0 20 15 10
2 1 −2.5953 −2.4825 0.5642
3 2 0.3084 0.0915 −0.2699
4 3 −0.02477 0.01479 0.02718
5 4 0.002359 −0.003505 −0.001411

s(t) = 20− 2.5953t1/2 + 0.3084t− 0.02477t3/2 + 0.002359t2

i(t) = 15− 2.4825t1/2 + 0.0915t+ 0.01479t3/2 − 0.003505t2

r(t) = 10 + 0.5642t1/2 − 0.2699t+ 0.02718t3/2 − 0.001411t2

(16)

Case III: γ = 1
For γ = 1, setting γµ = 1,i.e µ = 1 we obtain (12) as follows

S(k + 1) =
Γ(1 + k)

Γ(2 + k)

[
f1δ(k)− λ

k∑
l=0

S(l)I(k − l)− dS(k)

]



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 10 of 40

I(k + 1) =
Γ(1 + k)

Γ(2 + k)

[
f2δ(k) + λ

k∑
l=0

S(l)I(k − l)− εI(k)− dR(k)

]

R(k + 1) =
Γ(1 + k)

Γ(2 + k)
[f3δ(k) + εI(k)− dR(k)]

The table of values of S(k), I(k), R(k), 0 ≤ k ≤ 4 can be obtained as follows. The

Table 4: Values of S(k), I(k), R(k), γ = 1.

Sr. No. k S(k) I(k) R(k)

1 0 20 15 10
2 1 −2.3 −2.2 0.5
3 2 0.15425 0.04575 −0.135
4 3 −0.007904 0.005737 0.006025
5 4 0.0003097 −0.0004061 −0.0000216

approximate solution up to five iterations in this case is

s(t) = 20− 2.3t+ 0.15425t2 − 0.007904t3 + 0.0003097t4

i(t) = 15− 2.2t+ 0.04575t2 + 0.005737t3 − 0.0004061t4 (17)

r(t) = 10 + 0.5t− 0.135t2 + 0.006025t3 − 0.0000216t4

7. Laplace-Adomian Decomposition Method (LADM)

Now we solve the system of equations (1) using the Laplace-Adomian Decomposition
Method. Denoting the Laplace transform L[ϕ(t)] = Φ(z) and the Laplace transform of
the Caputo fractional derivative of the function ϕ by

L
[dγϕ(t)

dt

]
= zγL[ϕ(t)]−

n−1∑
k=0

zγ−k−1ϕ(k)(0)

Since γ ∈ (0, 1] i.e. n = 1, the above equation become

L
[dγϕ(t)

dt

]
= zγL[ϕ(t)]− zγ−1ϕ(0)

Using this relation and the system of equations (1), we get

L[s(t)] = s(0)

z
− λ

L[s(t)i(t)]
zγ

− d
L[s(t)]
zγ

L[i(t)] = i(0)

z
+ λ

L[s(t)i(t)]
zγ

− ε
L[i(t)]
zγ

− d
L[r(t)]
zγ



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 11 of 40

L[r(t)] = r(0)

z
+ ε

L[i(t)]
zγ

− d
L[r(t)]
zγ

Setting A(t) = s(t)i(t), we have

L[s(t)] = s(0)

z
− λ

L[A(t)]
zγ

− d
L[s(t)]
zγ

L[i(t)] = i(0)

z
+ λ

L[A(t)]
zγ

− ε
L[i(t)]
zγ

− d
L[r(t)]
zγ

L[r(t)] = r(0)

z
+ ε

L[i(t)]
zγ

− d
L[r(t)]
zγ

The non-linear term A is called Adomian polynomial and it is given by

A =
∞∑
j=0

Aj , An =
n∑

j=0

s(j).i(n− j)

and

s =

∞∑
j=0

sj , i =

∞∑
j=0

ij , r =

∞∑
j=0

rj

The terms sj , ij , rj are obtained as follows.

L[s0] =
s(0)

z
=⇒ s0 = s(0) = 20

L[i0] =
i(0)

z
=⇒ s0 = i(0) = 15

L[r0] =
r(0)

z
=⇒ r0 = r(0) = 10

L[s1] = −λL[s0i0]
zγ

− d
L[s0]
zγ

=⇒ s1 = −2.3
tγ

Γ(γ + 1)

L[i1] = λ
L[s0i0]
zγ

− ε
L[i0]
zγ

− d
L[r0]
zγ

=⇒ i1 = −2.2
tγ

Γ(γ + 1)

L[r1] = ε
L[i0]
zγ

− d
L[r0]
zγ

=⇒ r1 = 0.5
tγ

Γ(γ + 1)

L[s2] = −λL[A1]

zγ
− d

L[s1]
zγ

=⇒ s2 = 0.3085
t2γ

Γ(2γ + 1)

L[i2] = λ
L[A1]

zγ
− ε

L[i1]
zγ

− d
L[r1]
zγ

=⇒ i2 = 0.0915
t2γ

Γ(2γ + 1)

L[r2] = ε
L[i1]
zγ

− d
L[r1]
zγ

=⇒ r2 = −0.27
t2γ

Γ(2γ + 1)

L[s3] = −λL[A2]

zγ
− d

L[s2]
zγ

=⇒ s3 = −0.0193325
t3γ

Γ(3γ + 1)



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 12 of 40

L[i3] = λ
L[A2]

zγ
− ε

L[i2]
zγ

− d
L[r2]
zγ

=⇒ i3 = 0.00293675
t3γ

Γ(3γ + 1)

L[r3] = ε
L[i2]
zγ

− d
L[r2]
zγ

=⇒ r3 = 0.03615
t3γ

Γ(3γ + 1)

L[s4] = −λL[A3]

zγ
− d

L[s3]
zγ

=⇒ s4 = 0.0025250375
t4γ

Γ(4γ + 1)

L[i4] = λ
L[A3]

zγ
− ε

L[i3]
zγ

− d
L[r3]
zγ

=⇒ i4 = −0.0045004625
t4γ

Γ(4γ + 1)

L[r4] = ε
L[i3]
zγ

− d
L[r3]
zγ

=⇒ r4 = −0.0033321325
t4γ

Γ(4γ + 1)

It can be verified that for j ≥ 1, we have

L[sj ] = −λL[Aj−1]

zγ
− d

L[sj−1]

zγ

L[ij ] = λ
L[Aj−1]

zγ
− ε

L[ij−1]

zγ
− d

L[rj−1]

zγ

L[rj ] = ε
L[ij−1]

zγ
− d

L[rj−1]

zγ

The solution of the SIR model of computer viruses by Laplace Adomian decomposition
method up to five iteration is given by

s(t) =
4∑

j=0

sj = 20− 2.3
tγ

Γ(γ + 1)
+ 0.3085

t2γ

Γ(2γ + 1)
− 0.0193325

t3γ

Γ(3γ + 1)

+ 0.0025250375
t4γ

Γ(4γ + 1)

i(t) =

4∑
j=0

ij = 15− 2.2
tγ

Γ(γ + 1)
+ 0.0915

t2γ

Γ(2γ + 1)
+ 0.00293675

t3γ

Γ(3γ + 1)

− 0.0045004625
t4γ

Γ(4γ + 1)
(18)

r(t) =

4∑
j=0

rj = 10 + 0.5
tγ

Γ(γ + 1)
− 0.27

t2γ

Γ(2γ + 1)
+ 0.03615

t3γ

Γ(3γ + 1)

− 0.0033321325
t4γ

Γ(4γ + 1)

As in the earlier case, the solutions for different values of γ are obtained as follows.
CaseI: γ = 1/10

s(t) =20− 2.4176145292t1/10 + 0.3359947896t2/10 − 0.021541986t3/10

+ 0.0028458694t4/10



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 13 of 40

i(t) =15− 2.312500854t1/10 + 0.0996548566t2/10 + 0.00322722506t3/10

− 0.0050722925t4/10 (19)

r(t) =10 + 0.5255683759t1/10 − 0.2940635112t2/10 + 0.0402798531t3/10

− 0.0037555141t4/10

Case II: γ = 1/2

s(t) = 20− 2.595271866t1/2 + 0.3084t− 0.0145429311t3/2 + 0.0012625188t2

i(t) = 15− 2.4824339588t1/2 + 0.0915t+ 0.002209179t3/2 − 0.0022502313t2 (20)

r(t) = 10 + 0.5641895361t1/2 − 0.27t+ 0.0271939459t3/2 − 0.0016660663t2

Case III: γ = 1

s(t) = 20− 2.3t+ 0.15425t2 − 0.0032220833t3 + 0.0001052099t4

i(t) = 15− 2.2t+ 0.04575t2 + 0.0004894583t3 − 0.0001875193t4 (21)

r(t) = 10 + 0.5t− 0.135t2 + 0.006025t3 − 0.0001388389t4

8. Numerical Approximations and Estimate of Errors

Using the equations (15), (16), which are obtained using the Fractional Differential
Transform Method and the equations (19), (20), obtained using the Laplace Adomian
Decomposition Method, the approximate values of the solutions s(t), i(t) and r(t) for the
cases γ = 0.1 and γ = 0.5 are as shown by the tables 5 to 10. In these tables, the
values are obtained by using the MATLAB software correct up to fourteen digits after the
decimal. The notations s(t)FDTM, i(t)FDTM and r(t)FDTM are used to represent the
solutions s(t), i(t) and r(t) obtained using the Fractional Differential Transform Method
respectively. Similarly, the notations s(t)LADM, i(t)LADM and r(t)LADM are used to
represent the solutions s(t), i(t) and r(t) obtained using the Laplace Adomian Decompo-
sition Method respectively. The approximate values of the absolute errors in s(t), i(t) and
r(t) are represented by Er[s(t)], Er[i(t)] and Er[r(t)] and are obtained by the formulae

Er[s(t)] = |s(t)FDTM − s(t)LADM |
Er[i(t)] = |i(t)FDTM − i(t)LADM |
Er[r(t)] = |r(t)FDTM − r(t)LADM |



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 14 of 40

Table 5: Estimate of errors in s(t) for γ = 0.1.

t s(t)FDTM s(t)LADM Er[s(t)]

0 20.00000000000000 20.00000000000000 0.00000000000000
5 17.55702229544161 17.59431154300001 0.03728924755840
10 17.40781060925223 17.45308607588309 0.04527546663086
15 17.31712461031913 17.36782451716560 0.05069990684647
20 17.25124112568540 17.30616983119194 0.05492870550654
25 17.19925180577707 17.25769584483516 0.05844403905810
30 17.15619659362801 17.21767482183604 0.06147822820803

Table 6: Estimate of errors in i(t) for γ = 0.1.

t i(t)FDTM i(t)LADM Er[i(t)]

0 15.00000000000000 15.00000000000000 0.00000000000000
5 12.46558082849279 12.41676378653062 0.04881704196217
10 12.30053972810497 12.24037435176146 0.06016537634351
15 12.19982823339885 12.13183803127759 0.06799020212126
20 12.12647775783722 12.05232598332779 0.07415177450943
25 12.06849040276133 11.98917726395060 0.07931313881073
30 12.02039883376837 11.93660290522253 0.08379592854584

Table 7: Estimate of errors in r(t) for γ = 0.1.

t r(t)FDTM r(t)LADM Er[r(t)]

0 10.00000000000000 10.00000000000000 0.00000000000000
5 10.27004286714835 10.26974491443414 0.00029795271420
10 10.26787884169264 10.26652744914304 0.00135139254960
15 10.26545056249648 10.26327199707570 0.00217856542077
20 10.26315827919200 10.26027696818093 0.00288131101106
25 10.26103572422611 10.25753353094654 0.00350219327957
30 10.25906856101875 10.25500465803592 0.00406390298284



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 15 of 40

Table 8: Estimate of errors in s(t) for γ = 0.5.

t s(t)FDTM s(t)LADM Er[s(t)]

0 20.00000000000000 20.00000000000000 0.00000000000000
5 15.52077075898145 15.61380112841204 0.09303036943060
10 14.32954461214130 14.55193192445817 0.22238731231687
15 13.66621435930453 14.02421111980599 0.35799676050146
20 13.28956940388280 13.77791436087322 0.48834495699041
25 13.11162500000000 13.71834853250000 0.60672353250000
30 13.08993014024754 13.79851856567137 0.70858842542383

Table 9: Estimate of errors in i(t) for γ = 0.5.

t i(t)FDTM i(t)LADM Er[i(t)]

0 15.00000000000000 15.00000000000000 0.00000000000000
5 9.98419347279288 9.87505250806161 0.10914096473127
10 8.18184657457090 7.90969179323354 0.27215478133736
15 6.82841519839620 6.38011427914472 0.44830091925149
20 5.64878030720242 5.02572029367717 0.62306001352524
25 4.53312500000000 3.74508301850000 0.78804198150000
30 3.42353249758467 2.48594621404391 0.93758628354076

Table 10: Estimate of errors in r(t) for γ = 0.5.

t r(t)FDTM r(t)LADM Er[r(t)]

0 10.00000000000000 10.00000000000000 0.00000000000000
5 10.18069619104760 10.17395205545793 0.00674413558967
10 9.80356412390077 9.77746541222377 0.02609871167700
15 9.39817751417899 9.34005725374527 0.05812026043372
20 8.99183221094853 8.88900625418141 0.10282595676712
25 8.58912500000000 8.42889948050000 0.16022551950000
30 8.18948040334127 7.95915496540499 0.23032543793628

9. Runge-Kutta Felhberg Algorithm

An important family of predictor-corrector methods to solve linear and non-linear
ordinary differential equations is known to the Runge-Kutta methods, named after the
German mathematicians C. Runge (1856–1927) and M. W. Kutta (1867–1944).

Consider the ordinary differential equation

Y ′(t) = F [t, Y (t)], Y (0) = Y0



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 16 of 40

where Y (t) = (y1(t), y2(t), · · · , yn(t))T , F : [a, b]× Rn → Rn is a continuous function and
a, b ∈ R. In case the function F is non-linear and the analytical solutions do not exist, we
need to obtain an approximate solution by dividing the interval [a, b] into n equal parts,
each of length h = b−a

n , n ∈ N,h is also called as the step size. Considering the mesh
points ti = a + ih, i = 0, 1, · · · , n, the n + 1th approximate solution at the mth stage to
the initial value problem is given by:

Y (tn+1) = Yn+1 = Yn +
m∑
i=1

aiki

where

km = F (tn + cmh+ Yn + h

m−1∑
i=1

dmiki)

The coefficients ci, 2 ≤ i ≤ m, dij , 1 ≤ j < i ≤ m, and ai, i ≤ i ≤ m are to be
determined in the form of a table, known as Butcher Tableau. In order to improve the rate
of convergence, the algorithm based on step-size adjustment was developed by Felhberg
and is known as Runge-Kutta Felhberg Algorithm (RKFA). This algorithm is an adaptive
algorithm where at each step, two different approximate solutions are obtained and are
compared. If the two approximate solutions agree up to a desired level of accuracy, the
approximate solution is accepted. In case the two approximate solutions did not agree,
the step-size is reduced and the same procedure is repeated.

This algorithm evaluates F [t, Y (t)] six times per step using the embedded 4th order
and 5th order Runge-Kutta method. Moreover, the algorithm also evaluates the amount
of errors. The approximate solution to the initial value problem by Runge-Kutta method
of order 4 is given by:

Yk+1 = Yk +
25k1
216

+
1408k2
2565

+
2197k4
4101

− k5
5

where at each step, the following six values are calculated.

k1 = hF (tk, Yk)

k2 = hF

(
tk +

h

4
, Yk +

k1
4

)
k3 = hF

(
tk +

3h

8
, Yk +

3k1
32

+
9k2
32

)
k4 = hF

(
tk +

12h

13
, Yk +

1932k1
2197

− 7200k2
2197

+
7296k3
2197

)
k5 = hF

(
tk + h, Yk +

439k1
216

− 8k2 +
3680k3
513

− 845k4
4104

)
k6 = hF

(
tk +

h

2
, Yk −

8k1
27

+ 2k2 −
3544k3
2565

+
1859k4
4104

− 11k5
40

)



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 17 of 40

A better approximate solution is given by the 5th order Runge-Kutta Felhberg Algorithm,
which is given by

Yk+1 = Yk +
16k1
135

+
6656k3
12825

+
28561k4
56430

− 9k5
50

+
2k6
55

It has been proved that RKFA provides more accurate solutions as compared to the
Runge-Kutta method of order 4. Moreover, RKFA requires less number of iterations
as compared to the R-K method to arrive at the desired level of accuracy. Using the
RKFA, we have obtained the numerical values of s(t), i(t) and r(t) and compared them
with the values obtained by FDTM and LADM for the case γ = 1. Moreover, we
have obtained the estimates of the errors in these solutions with respect to the RKFA.
The approximate values of the absolute errors in s(t), i(t) and r(t) are represented by
ErFDTM [s(t)], ErFDTM [i(t)], ErFDTM [r(t)], etc. and are obtained by the formulae

ErFDTM [s(t)] = |s(t)FDTM − s(t)RKFA|
ErFDTM [i(t)] = |i(t)FDTM − i(t)RKFA|
ErFDTM [r(t)] = |r(t)FDTM − r(t)RKFA|
ErLADM [s(t)] = |s(t)LADM − s(t)RKFA|, etc.

The estimate of errors in s(t), i(t), r(t) for γ = 1 are as shown in the tables 11, 12, 13.

Table 11: Estimate of errors in s(t) for γ = 1.

t s(t)FDTM s(t)LADM s(t)RKFA ErFDTM[s(t)] ErLADM[s(t)]

0.0 20.00000000000000 20.00000000000000 20.00000000000000 0.00000000000000 0.00000000000000
0.2 19.54610726352000 19.54614439166944 19.546107256391014 0.00000000712900 0.00003713527840
0.4 19.10418207232000 19.10447648004224 19.104181955888297 0.00000011643180 0.00029452415400
0.6 18.67386287312000 18.67484766521024 18.673862145750277 0.00000072736980 0.00098551946000
0.8 18.25480000512000 18.25711338732544 18.254797171108123 0.00000283401190 0.00231621621730
1.0 17.84665570000000 17.85113312660000 17.846647378609305 0.00000832139070 0.00448574799070
1.2 17.44910408192000 17.45677040330624 17.449083821756879 0.00002026016320 0.00768658154940
1.4 17.06183116752001 17.07389277777664 17.061787968569263 0.00004319895080 0.01210480920740
1.6 16.68453486592000 16.70237185040384 16.684451412115802 0.00008345380420 0.01792043828800
1.8 16.31692497872000 16.34208326164064 16.316775584419219 0.00014939430080 0.02530767722140
2.0 15.95872320000000 15.99290669200000 15.958471474156529 0.00025172584350 0.03443521784350



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 18 of 40

Table 12: Estimate of errors in i(t) for γ = 1.

t i(t)FDTM i(t)LADM i(t)RKFA ErFDTM[i(t)] ErLADM[i(t)]

0.0 15.00000000000000 15.00000000000000 15.00000000000000 0.00000000000000 0.00000000000000
0.2 14.56187524624000 14.56183361563552 14.561875256612106 0.00000001037210 0.00004164097660
0.4 14.12767677184000 14.12734652483712 14.127676929730777 0.00000015789070 0.00033040489360
0.6 13.69765656144000 13.69655142049152 13.697657518861158 0.00000095742110 0.00110609836960
0.8 13.27205100544000 13.26945379474432 13.272054690053119 0.00000368461310 0.00260089530880
1.0 12.85108090000000 12.84605193900000 12.851091662011378 0.00001076201130 0.00503972301130
1.2 12.43495144704000 12.42633694392192 12.538540600282305 0.10358915324230 0.11220365636040
1.4 12.02385225424000 12.01029269943232 12.023907974631463 0.00005572039140 0.01361527519910
1.6 11.61795733504000 11.59789589471232 11.618065010548600 0.00010767550860 0.02016911583630
1.8 11.21742510864000 11.18911601820192 11.217618014408080 0.00019290576800 0.02850199620610
2.0 10.82239840000000 10.78391535760000 10.822723791238749 0.00032539123870 0.03880843363870

Table 13: Estimate of errors in r(t) for γ = 1.

t r(t)FDTM r(t)LADM r(t)RKFA ErFDTM[r(t)] ErLADM[r(t)]

0.0 10.00000000000000 10.00000000000000 10.00000000000000 0.00000000000000 0.00000000000000
0.2 10.09464816544000 10.09464817778578 10.094648185984170 0.00000002054410 0.00000000819840
0.4 10.17878504704000 10.17878524457242 10.178785335767941 0.00000028872790 0.00000009119550
0.6 10.25269860064000 10.25269960064786 10.252699863108273 0.00000126246820 0.00000026246040
0.8 10.31667595264000 10.31667911315866 10.316679322626605 0.00000336998660 0.00000020946800
1.0 10.37100340000000 10.37101111611000 10.371010136930472 0.00000673693040 0.00000097917960
1.2 10.41596641024000 10.41598241036570 10.415977336803497 0.00001092656340 0.00000507356230
1.4 10.45184962144000 10.45187926364818 10.451864314231031 0.00001469279100 0.00001494941710
1.6 10.47893684224000 10.47898741053850 10.478952588015369 0.00001574577530 0.00003482252320
1.8 10.49751105184000 10.49759205247634 10.497521581722806 0.00001052988280 0.00007047075350
2.0 10.50785440000000 10.50797785776000 10.507848413694944 0.00000598630510 0.00012944406510

10. Results and Discussion

In the second section of this paper, we have formulated the non-linear fractional order
SIR model of computer viruses with the specific initial conditions as described by the table
1. In the third section, we concluded that the system of equations 1 is asymptotically
stable. In the sixth section, we have obtained the numerical values of the solutions which
are given by the equations (15) and (16) for the cases γ = 0.1, 0.5. Whereas using the
FDTM the solutions are obtained for the cases γ = 0.1, 0.5 as given by equations (19) and
(20). To have a better understanding of the solution curves and errors, we have compared
the numerical values obtained by FDTM and LADM with those obtained by RKFA. As
can be noticed, we don’t get a clear picture of the solutions and the errors just from the
tabulated values. For a better understanding of the system, we plot the graphs of the
approximate solutions and errors for different values of γ.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 19 of 40

Figure 1: Comparison of graphs of s(t) by FDTM and LADM for γ = 0.1

0 5 10 15 20 25 30
17

17.5

18

18.5

19

19.5

20

s(t) by FDTM and LADM for γ=0.1

t

s
(t

)

 

 

s(t)FDTM

s(t)LADM



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 20 of 40

Figure 2: Comparison of graphs of i(t) by FDTM and LADM for γ = 0.1

0 5 10 15 20 25 30
11.5

12

12.5

13

13.5

14

14.5

15

i(t) by FDTM and LADM for γ=0.1

t

i(
t)

 

 

i(t)FDTM

i(t)LADM



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 21 of 40

Figure 3: Comparison of graphs of r(t) by FDTM and LADM for γ = 0.1

0 5 10 15 20 25 30
10

10.05

10.1

10.15

10.2

10.25

10.3

10.35

r(t) by FDTM and LADM for γ=0.1

t

r(
t)

 

 

r(t)FDTM

r(t)LADM

The graphs shown by figures 1, 2 and 3 represent the solution curves for s(t), i(t) and
r(t) by FDTM and LADM for γ = 0.1. From these figures, it can be observed that the
solutions obtained by FDTM and LADM almost agree for first few values of the variable
t, to be specific, for t ∈ [0, 5]. But as t increases, slight differences in the solution curves
are observed. The absolute errors in the solutions by both the methods can be observed
in figure 4.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 22 of 40

Figure 4: Comparison of errors in s(t), i(t), r(t) by FDTM and LADM for γ = 0.1

0 5 10 15 20 25 30
0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

Errors in s(t), i(t), r(t) by FDTM and LADM for γ=0.1

t

E
rr

o
rs

 i
n

 s
(t

),
 i
(t

),
 r

(t
)

 

 

Er[s(t)]

Er[i(t)]

Er[r(t)]

The phase-plane portraits of s(t), i(t), r(t) by both the methods are as plotted in figure
5. From this figure we note that the solution curves for s(t), i(t), r(t) progress almost hand-
in-hand for the initial values of the variable t, but differ later on.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 23 of 40

Figure 5: Phase plane portraits of s(t), i(t), r(t) by FDTM and LADM for γ = 0.1

0 5 10 15 20 25 30
10

11

12

13

14

15

16

17

18

19

20

Phase portratis of s(t), i(t), r(t) by FDTM and LADM for γ=0.1

t

s
(t

),
 i
(t

),
 r

(t
)

 

 

s(t)FDTM

s(t)LADM

i(t)FDTM

i(t)LADM

r(t)FDTM

r(t)LADM

The solution curves, errors and the phase-plain portrait for γ = 0.5 are as shown by
figures 6 to figure 10. In this case also we have a similar type of conclusion as in the case
for γ = 0.1. This type of pattern can be observed for the both the methods for γ ∈ (0, 1).



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 24 of 40

Figure 6: Comparison of graphs of s(t) by FDTM and LADM for γ = 0.5

0 5 10 15 20 25 30
13

14

15

16

17

18

19

20

s(t) by FDTM and LADM for γ=0.5

t

s
(t

)

 

 

s(t)FDTM

s(t)LADM



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 25 of 40

Figure 7: Comparison of graphs of i(t) by FDTM and LADM for γ = 0.5

0 5 10 15 20 25 30
2

4

6

8

10

12

14

16

i(t) by FDTM and LADM for γ=0.5

t

i(
t)

 

 

i(t)FDTM

i(t)LADM



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 26 of 40

Figure 8: Comparison of graphs of r(t) by FDTM and LADM for γ = 0.5

0 5 10 15 20 25 30
7.5

8

8.5

9

9.5

10

10.5

r(t) by FDTM and LADM for γ=0.5

t

r(
t)

 

 

r(t)FDTM

r(t)LADM



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 27 of 40

Figure 9: Comparison of errors in s(t), i(t), r(t) by FDTM and LADM for γ = 0.5

0 5 10 15 20 25 30
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Errors in s(t), i(t), r(t) by FDTM and LADM for γ=0.5

t

E
rr

o
rs

 i
n

 s
(t

),
 i
(t

),
 r

(t
)

 

 

Er[s(t)]

Er[i(t)]

Er[r(t)]



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 28 of 40

Figure 10: Phase plane portraits of s(t), i(t), r(t) by FDTM and LADM for γ = 0.5

0 5 10 15 20 25 30
2

4

6

8

10

12

14

16

18

20

Phase portratis of s(t), i(t), r(t) by FDTM for γ=0.5

t

s
(t

),
 i
(t

),
 r

(t
)

 

 

s(t)FDTM

s(t)LADM

i(t)FDTM

i(t)LADM

r(t)FDTM

r(t)LADM

For γ = 1, figure 11 shows solution curves for s(t) by FDTM, LADM and RKFA. In
the graph, the curves are respectively named as s(t)FDTM, s(t)LADM and s(t)RKFA.
In the figure, we can not find any differences in the curves just by naked eyes and the curves
appear to be coincident and one can not have better understanding of which method among
FDTM and LADM is close enough to RKFD. To overcome this difficulty, we have obtained
a zoom-in figure of figure 11, which is shown in figure 12. From this figure, we observe
that the curve obtained by FDTM is comparatively colse to the one obtained by RKFA.
A similar conclusion for the curves i(t) and r(t) can be drawn form the figure 13 to figure
16.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 29 of 40

Figure 11: Comparison of graphs of s(t) by FDTM, LADM and RKFA for γ = 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
15.5

16

16.5

17

17.5

18

18.5

19

19.5

20

t

s
(t

)

s(t) by FDTM, LADM and RKFA for γ=1

 

 

s(t)FDTM

s(t)LADM

s(t)RKFA

Figure 12: Zoom-in s(t) by FDTM and LADM for γ = 1

1.986 1.988 1.99 1.992 1.994 1.996 1.998 2
15.965

15.97

15.975

15.98

15.985

15.99

15.995

t

s
(t

)

s(t) by FDTM, LADM and RKFA for γ=1

 

 

s(t)FDTM

s(t)LADM

s(t)RKFA



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 30 of 40

Figure 13: Comparison of graphs of i(t) by FDTM and LADM for γ = 1

0 0.5 1 1.5 2
10.5

11

11.5

12

12.5

13

13.5

14

14.5

15

t

i(
t)

i(t) by FDTM, LADM and RKFA for γ=1

 

 

i(t)FDTM

i(t)LADM

i(t)RKFA



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 31 of 40

Figure 14: Zoom-in i(t) by FDTM and LADM for γ = 1

1.974 1.976 1.978 1.98 1.982 1.984 1.986 1.988

10.825

10.83

10.835

10.84

10.845

10.85

10.855

t

i(
t)

i(t) by FDTM, LADM and RKFA for γ=1

 

 

i(t)FDTM

i(t)LADM

i(t)RKFA

Figure 15: r(t) by FDTM and LADM for γ = 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
10

10.1

10.2

10.3

10.4

10.5

10.6

t

r(
t)

r(t) by FDTM, LADM and RKFA for γ=1

 

 

r(t)FDTM

r(t)LADM

r(t)RKFA



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 32 of 40

Figure 16: Zoom-in r(t) by FDTM and LADM for γ = 1

1.9772 1.9773 1.9774 1.9775 1.9776 1.9777 1.9778 1.9779 1.978 1.9781

10.5066

10.5067

10.5067

10.5068

10.5068

t

r(
t)

r(t) by FDTM, LADM and RKFA for γ=1

 

 

r(t)FDTM

r(t)LADM

r(t)RKFA

The figure 17 shows the estimate of absolute errors in s(t) by FDTM and LADM with
respect to RKFA. In the figure, we observe that the absolute error by FDTM is almost
close to zero and the curve is going parallel to the x-axis while the error by LADM is
though in between 10−2 and 10−4, but still growing exponentially. The figure 18 shows
the error in i(t) and the figure 19 shows the error in r(t).



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 33 of 40

Figure 17: Error in s(t) by FDTM and LADM for γ = 1

0 0.5 1 1.5 2
0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

t

E
r[

s
(t

)]

Error in s(t) by FDTM, LADM with respect to RKFA for γ=1

 

 

ErFDTM[s(t)]

ErLADM[s(t)]



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 34 of 40

Figure 18: Error in i(t) by FDTM and LADM for γ = 1

0 0.5 1 1.5 2
0

0.02

0.04

0.06

0.08

0.1

0.12

t

E
r[

i(
t)

]

Error in i(t) by FDTM, LADM with respect to RKFA for γ=1

 

 

ErFDTM[i(t)]

ErLADM[i(t)]



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 35 of 40

Figure 19: Error in r(t) by FDTM and LADM for γ = 1

0 0.5 1 1.5 2
0

0.2

0.4

0.6

0.8

1

1.2

1.4
x 10

−4

t

E
r[

r(
t)

]

Error in r(t) by FDTM, LADM with respect to RKFA for γ=1

 

 

ErFDTM[r(t)]

ErLADM[r(t)]

11. Conclusion

In this paper, the non-linear Fractional order susceptible-infected-Recovered(SIR) math-
ematical model of computer viruses is presented.The fractional derivative used in the model
is Caputo Fractional derivative. Using Lyapunov stability analysis,the stability verified
of the given mathematical model. The FDTM and the LADM are two reliable and use-
ful techniques that were used in this work to solve the non-linear epidemiological model
of computer viruses Involving Non-linear Fractional Order Differential Equation. These
models are among the useful computer engineering models that we may use to monitor,
assess, and forecast computer viruses within a network. To demonstrate the effectiveness
and precision of the approach that was described, the residual errors for five iterations
were displayed using the LADM, FDTM, and RKFA. Additionally, residual error graphs
were shown to illustrate the capabilities of the methodologies that were provided. These
findings suggest that the FDTM is more easy and applicable than the LADM without



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 36 of 40

adomian polynomial. From the graphical analysis of the solution curves and errors by
by FDTM and LADM, we conclude that the FDTM for solving the non-linear system of
equations 1 presents better solutions as compared to the LADM. However, one should note
that due to the complexity of calculations, while finds the solutions by FDTM and LADM,
we have considered only five iterations and a small interval for the values of the variable
t. The conclusions drawn here may vary if one considers a large number of iterations and
find out the long term behaviour of the SIR model. Hence further investigation in this
matter is needed. By including a few more parameters, we may improve the model. In
many other fields, this model may also be used to provide alternative models.

References

[1] Fred Cohen. Computer viruses: theory and experiments. Computers & security,
6(1):22–35, 1987.

[2] Asad A Freihat, M Zurigat, and Ali H Handam. The multi-step homotopy analysis
method for modified epidemiological model for computer viruses. Afrika Matematika,
26(3):585–596, 2015.

[3] Xie Han and Qiulin Tan. Dynamical behavior of computer virus on internet. Applied
mathematics and computation, 217(6):2520–2526, 2010.

[4] Bimal Kumar Mishra and Navnit Jha. Fixed period of temporary immunity after run
of anti-malicious software on computer nodes. Applied Mathematics and Computation,
190(2):1207–1212, 2007.

[5] Jeffrey O Kephart, Tad Hogg, and Bernardo A Huberman. Dynamics of computa-
tional ecosystems. Physical Review A, 40(1):404, 1989.

[6] Jose RC Piqueira, Adolfo A De Vasconcelos, Carlos ECJ Gabriel, and Vanessa O
Araujo. Dynamic models for computer viruses. computers & security, 27(7-8):355–
359, 2008.

[7] José Roberto C Piqueira and Vanessa O Araujo. A modified epidemiological model
for computer viruses. Applied Mathematics and Computation, 213(2):355–360, 2009.

[8] Jianguo Ren, Xiaofan Yang, Lu-Xing Yang, Yonghong Xu, and Fanzhou Yang. A
delayed computer virus propagation model and its dynamics. Chaos, Solitons &
Fractals, 45(1):74–79, 2012.

[9] Jianguo Ren, Xiaofan Yang, Qingyi Zhu, Lu-Xing Yang, and Chunming Zhang. A
novel computer virus model and its dynamics. Nonlinear Analysis: Real World Ap-
plications, 13(1):376–384, 2012.

[10] John C Wierman and David J Marchette. Modeling computer virus prevalence with a
susceptible-infected-susceptible model with reintroduction. Computational statistics
& data analysis, 45(1):3–23, 2004.

[11] Hua Yuan and Guoqing Chen. Network virus-epidemic model with the point-to-group
information propagation. Applied Mathematics and Computation, 206(1):357–367,
2008.

[12] Qingyi Zhu, Xiaofan Yang, and Jianguo Ren. Modeling and analysis of the spread of
computer virus. Applied Mathematics and Computation, 17(12):5117–5124, 2012.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 37 of 40

[13] Yalçın Öztürk and Mustafa Gülsu. Numerical solution of a modified epidemiological
model for computer viruses. Applied Mathematical Modelling, 39(23-24):7600–7610,
2015.

[14] Mohammad Ali Fariborzi Araghi and Samad Noeiaghdam. Fibonacci-regularization
method for solving cauchy integral equations of the first kind. Ain Shams Engineering
Journal, 8(3):363–369, 2017.

[15] Samad Noeiaghdam, Muhammad Suleman, and Hüseyin Budak. Solving a modified
nonlinear epidemiological model of computer viruses by homotopy analysis method.
Mathematical Sciences, 12(3):211–222, 2018.

[16] Jagdev Singh, Devendra Kumar, Zakia Hammouch, and Abdon Atangana. A frac-
tional epidemiological model for computer viruses pertaining to a new fractional
derivative. Applied mathematics and computation, 316:504–515, 2018.

[17] Bimal Kumar Mishra and Samir Kumar Pandey. Fuzzy epidemic model for the trans-
mission of worms in computer network. Nonlinear Analysis: Real World Applications,
11(5):4335–4341, 2010.

[18] Samad Noeiaghdam. A novel technique to solve the modified epidemiological model
of computer viruses. SeMA Journal, 76(1):97–108, 2019.

[19] MA Fariborzi Araghi and Samad Noeiaghdam. A novel technique based on the ho-
motopy analysis method to solve the first kind cauchy integral equations arising in
the theory of airfoils. Journal of Interpolation and Approximation in Scientific Com-
puting, 2016(1):1–13, 2016.

[20] Mohammad Ali Fariborzi Araghi and Samad Noeiaghdam. Homotopy regularization
method to solve the singular volterra integral equations of the first kind. Jordan
Journal of Mathematics and Statistics, 11(1):1–12, 2018.

[21] Mohammad Ali Fariborzi and Samad Noeiaghdam. Valid implementation of the sinc-
collocation method to solve linear integral equations by the cadna library. Journal of
Mathematical Modeling, 7(1):63–84, 2019.

[22] Samad Noeiaghdam and Emran Khoshrouye Ghiasi. An efficient method to
solve the mathematical model of hiv infection for cd8+ t-cells. arXiv preprint
arXiv:1907.01106, 2019.

[23] Samad Noeiaghdam, Mohammad Ali Fariborzi Araghi, and Saeid Abbasbandy. Find-
ing optimal convergence control parameter in the homotopy analysis method to
solve integral equations based on the stochastic arithmetic. Numerical Algorithms,
81(1):237–267, 2019.

[24] Samad Noeiaghdam, Denis Sidorov, and Valery Sizikov. Control of accuracy on taylor-
collocation method to solve the weakly regular volterra integral equations of the first
kind by using the cestac method. arXiv preprint arXiv:1811.09802, 2018.

[25] Eisa Zarei and Samad Noeiaghdam. Solving generalized abel’s integral equa-
tions of the first and second kinds via taylor-collocation method. arXiv preprint
arXiv:1804.08571, 2018.

[26] Mostafa Nourifar, Ahmad Aftabi Sani, and Ali Keyhani. Efficient multi-step differ-
ential transform method: Theory and its application to nonlinear oscillators. Com-
munications in Nonlinear Science and Numerical Simulation, 53:154–183, 2017.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 38 of 40

[27] Mohammad Mehdi Rashidi and Esmaeel Erfani. A new analytical study of mhd
stagnation-point flow in porous media with heat transfer. Computers & Fluids,
40(1):172–178, 2011.

[28] Kunjan Shah, Twinkle Singh, and Adem Kılıçman. Combination of integral and
projected differential transform methods for time-fractional gas dynamics equations.
Ain Shams Engineering Journal, 9(4):1683–1688, 2018.

[29] Abdul-Majid Wazwaz. The combined laplace transform–adomian decomposition
method for handling nonlinear volterra integro–differential equations. Applied Math-
ematics and Computation, 216(4):1304–1309, 2010.

[30] Hassan Eltayeb Gadain. Solving coupled pseudo-parabolic equation using a modified
double laplace decomposition method. Acta Mathematica Scientia, 38(1):333–346,
2018.

[31] Fazal Haq, Kamal Shah, Ghaus ur Rahman, and Muhammad Shahzad. Numerical so-
lution of fractional order smoking model via laplace adomian decomposition method.
Alexandria Engineering Journal, 57(2):1061–1069, 2018.

[32] Susmita Paul, Sankar Prasad Mondal, and Paritosh Bhattacharya. Numerical solution
of lotka volterra prey predator model by using runge–kutta–fehlberg method and
laplace adomian decomposition method. Alexandria Engineering Journal, 55(1):613–
617, 2016.

[33] Saeid Abbasbandy. Extended newton’s method for a system of nonlinear equations
by modified adomian decomposition method. Applied Mathematics and Computation,
170(1):648–656, 2005.

[34] Jun-Sheng Duan, Randolph Rach, and Abdul-Majid Wazwaz. Higher order numeric
solutions of the lane–emden-type equations derived from the multi-stage modified
adomian decomposition method. International Journal of Computer Mathematics,
94(1):197–215, 2017.

[35] Abdelhalim Ebaid, Mona D Aljoufi, and Abdul-Majid Wazwaz. An advanced study on
the solution of nanofluid flow problems via adomian’s method. Applied Mathematics
Letters, 46:117–122, 2015.

[36] S Hosseini, E Babolian, S Abbasbandy, et al. A new algorithm for solving van der pol
equation based on piecewise spectral adomian decomposition method. International
Journal of Industrial Mathematics, 8(3):177–184, 2016.

[37] Hadi Roohani Ghehsareh, Saeid Abbasbandy, and Babak Soltanalizadeh. Analytical
solutions of the slip magnetohydrodynamic viscous flow over a stretching sheet by
using the laplace–adomian decomposition method. Zeitschrift für Naturforschung A,
67(5):248–254, 2012.

[38] Abdul-Majid Wazwaz, Randolph Rach, and Jun-Sheng Duan. Adomian decompo-
sition method for solving the volterra integral form of the lane–emden equations
with initial values and boundary conditions. Applied Mathematics and Computation,
219(10):5004–5019, 2013.

[39] Sunil Kumar, Deepak Kumar, Saeid Abbasbandy, and MM Rashidi. Analytical so-
lution of fractional navier–stokes equation by using modified laplace decomposition
method. Ain Shams Engineering Journal, 5(2):569–574, 2014.



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 39 of 40

[40] Samad Noeiaghdam, Eisa Zarei, and Hasan Barzegar Kelishami. Homotopy analysis
transform method for solving abel’s integral equations of the first kind. Ain Shams
Engineering Journal, 7(1):483–495, 2016.

[41] Samad Noeiaghdam and Emran Khoshrouye Ghiasi. Solving a non-linear model of hiv
infection for cd4+ t cells by combining laplace transformation and homotopy analysis
method. arXiv preprint arXiv:1809.06232, 2018.

[42] Samad Noeiaghdam, Muhammad Suleman, and Hüseyin Budak. Solving a modified
nonlinear epidemiological model of computer viruses by homotopy analysis method.
Mathematical Sciences, 3(18):211–222, 2012.

[43] S. Noeiaghdam M.A. Fariborzi Araghi. Homotopy analysis transform method for
solving generalized abel’s fuzzy integral equations of the first kind. IEEE, 2016.

[44] Sumit Gupta, Devendra Kumar, and Jagdev Singh. Analytical solutions of
convection–diffusion problems by combining laplace transform method and homotopy
perturbation method. Alexandria Engineering Journal, 54(3):645–651, 2015.

[45] Claude Rodrigue Bambe Moutsinga, Edson Pindza, and Eben Mare. Homotopy per-
turbation transform method for pricing under pure diffusion models with affine coef-
ficients. Journal of King Saud University-Science, 30(1):1–13, 2018.

[46] MH Tiwana, K Maqbool, and AB Mann. Homotopy perturbation laplace transform
solution of fractional non-linear reaction diffusion system of lotka-volterra type dif-
ferential equation. Engineering Science and Technology, an International Journal,
20(2):672–678, 2017.

[47] Samad Noeiaghdam. Numerical approximation of modified non-linear sir model of
computer viruses. arXiv preprint arXiv:1901.10804, 2019.

[48] Ricardo Almeida, Nuno RO Bastos, and M Teresa T Monteiro. Modeling some real
phenomena by fractional differential equations. Mathematical Methods in the Applied
Sciences, 39(16):4846–4855, 2016.

[49] Kai Diethelm, Neville J Ford, Alan D Freed, and Yu Luchko. Algorithms for the
fractional calculus: a selection of numerical methods. Computer methods in applied
mechanics and engineering, 194(6-8):743–773, 2005.

[50] Kamble Rajratna M. and Kulkarni Pramod Ramakant. On some existence and unique-
ness results for nonlinear fractional differential equations with boundary conditions.
Econophysics, Sociophysics and Other Multidisciplinary Sciences Journal (ESMSJ),
12(1):11–18, 2023.

[51] Kamble Rajratna M. and Kulkarni Pramod R . Existence and uniqueness of contin-
uous solutions for conformable fractional integro-differential equations in cone metric
spaces. Communications on Applied Nonlinear Analysis, 32(9s):1303–1311, 2025.

[52] Rajratana Maroti Kamble and Pramod Ramakant Kulkarni . Extended fractional
derivative: Some results involving classical properties and applications. Scientific
Temper, 15(4):3026–3033, 2024.

[53] Kamble Rajratna M and Kulkarni Pramod R. Existence and uniqueness of solutions
for exponential fractional differential equations. Scientific Temper, 15(4):3008–3012,
2024.

[54] Kamble Rajratna M and Kulkarni Pramod R. Laplace transform and laplace de-



R. M. Kamble, P. R. Kulkarni / Eur. J. Pure Appl. Math, 18 (4) (2025), 6693 40 of 40

composition method for of order of fractional differential difference equations with
interval conditions linear and nonlinear. IOSR Journal of Mathematics (IOSR-JM),
21(1):04–29, 2025.

[55] Imed Basdouri, Souad Kasmi, Fahmi Mabrouk, and Rim Messaoud. Lyapunov stabil-
ity analysis of caputo fractional-order nonlinear systems. Authorea Preprints, 2020.

[56] Aytac Arikoglu and Ibrahim Ozkol. Solution of fractional differential equations by
using differential transform method. Chaos, Solitons & Fractals, 34(5):1473–1481,
2007.