Dynamics and simulations of discretized Caputo-conformable fractional-order Lotka–Volterra models

Feras Yousef; Billel Semmar; Kamal Al Nasr

doi:10.1515/nleng-2022-0013

Article Open Access

Dynamics and simulations of discretized Caputo-conformable fractional-order Lotka–Volterra models

Feras Yousef , Billel Semmar and Kamal Al Nasr

Published/Copyright: April 2, 2022

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Nonlinear Engineering Volume 11 Issue 1

Abstract

In this article, a prey–predator system is considered in Caputo-conformable fractional-order derivatives. First, a discretization process, making use of the piecewise-constant approximation, is performed to secure discrete-time versions of the two fractional-order systems. Local dynamic behaviors of the two discretized fractional-order systems are investigated. Numerical simulations are executed to assert the outcome of the current work. Finally, a discussion is conducted to compare the impacts of the Caputo and conformable fractional derivatives on the discretized model.

Keywords: Caputo fractional derivative; conformable fractional derivative; prey–predator system; discretization; stability; bifurcations; chaos

1 Introduction

Dynamic analysis is widely used in engineering and science. Understanding the dynamics of prey–predator interactions is one of the most primary operations that shape the framework and function of ecological societies [1,2]. To describe the dynamic behavior between prey and predator, models can be a very useful and powerful tool. The oldest and most celebrated prey–predator model is the Lotka–Volterra model, independently introduced by Lotka [3] and Volterra [4]. This model is formulated in ordinary differential equations such that any small change of the model will lead to a qualitatively different type of behavior. The most important feature of this model is that it lumps prey birth and death rates into one logistic growth term, and it assumes that the predator birth rate remains a linear function of their per capita consumption [5].

Predation is normally quantified in terms of the functional and numerical responses, which are the effects of predation on the prey and predator growth rates, respectively [6]. The classical model of Lotka–Volterra under the linear functional response of predator reads as follows:

(1) d N d T = N ( a − b N − c R ) , d R d T = R ( − d + c e N ) ,

where N and R are functions of time T , representing the prey and predator population densities, respectively [7]. The parameters a , b , c , d , and e are considered as positive constants, representing the maximum birth rate per capita of prey species, the strength of intraspecific competition of the prey species, the strength of intraspecific between prey and predators, the death rate per capita of predator species, and the efficacy of converting ingested prey to new predators, respectively.

Recently, it has been shown that many mathematical models can be effectively reformulated via noninteger-order differential equations owing to the unsuitability and ability of the integer-order differential equations in formulating many phenomena [8]. The noninteger-order derivatives are oftentimes said to be fractional-order derivatives or artlessly fractional derivatives. It is well-known that the derivatives with integer order are in local nature, unlike the derivatives with fractional order [9]. It has been evinced that the nonlocality merit of the fractional derivatives is employed to make them a suitable tool for describing the dynamic behaviors of many life phenomena and dynamical models that have inherited the properties of memory [10,11]. Due to the memory effect, the fractional models integrate all previous information from the past that makes it easier for the researchers to predict and translate the biological models more accurately [12]. Consequently, many existing differential equations describing various phenomena in engineering and science have been recasted by means of fractional derivatives [13,14, 15,16,17, 18,19,20, 21,22,23, 24,25], and their solutions and dynamic behaviors continue to be of widespread interest today in many other disciplines [26,27,28,29,30,31,32,33,34,35].

In the present work, we nondimensionalize the system (1) by performing the following change of variables.

x ( t ) = e a d N ( T ) , y ( t ) = c d R ( T ) , t = c d a T .

Hence, we are lead to work with the following nondimensionalized system:

(2) d x ( t ) d t = r x ( t ) − p x 2 ( t ) − q x ( t ) y ( t ) , d y ( t ) d t = x ( t ) y ( t ) − q y ( t ) ,

where r = a 2 c d , p = b c e , and q = a c .

Thus, the following fractional-order system is obtained by changing the first-order time derivatives by the derivatives of fractional order.

(3) D α x ( t ) = r x ( t ) − p x 2 ( t ) − q x ( t ) y ( t ) , D α y ( t ) = x ( t ) y ( t ) − q y ( t ) ,

where 0 < α < 1 is the fractional-order derivative parameter and t > 0 . It is easily verified that the equilibrium points of system (3) are the trivial state E 0 = ( 0 , 0 ) , the axial state E 1 = r p , 0 , and the steady state of coexistence E 2 = q , r − p q q for r > p q .

Several studies revealed that the discrete-time system exhibits much fruitful dynamic behaviors, such as bifurcations and chaos, than those of its continuous-time system counterpart. Consequently, in the sequel, we will explore the dynamic behaviors of the discretized fractional-order prey–predator system that includes both Caputo and conformable fractional derivatives.

2 The dynamics of the discretized Caputo fractional-order prey–predator system

In this section, we consider the Caputo fractional-order version of model (3) as follows:

(4) D α x ( t ) = r x ( t ) − p x 2 ( t ) − q x ( t ) y ( t ) , D α y ( t ) = x ( t ) y ( t ) − q y ( t ) ,

where D α is the fractional derivative of Caputo type and defined as follows:

(5) D α f ( t ) = 1 Γ ( 1 − α ) ∫ 0 t ( t − τ ) − α f ( τ ) d τ , 0 < α < 1 .

Our first aim is to discretize the model (4) making use of the piecewise-constant approximation [36] as follows:

(6) D α x ( t ) = x t h h r − p x t h h − q y t h h , D α y ( t ) = y t h h − q + x t h h ,

where [ t h ] denotes the integer part of t ∈ [ n h , ( n + 1 ) h ) , n = 0 , 1 , … , and h > 0 is a discretization parameter.

The nth iterative solution of system (6) is given by:

(7) x n + 1 ( t ) = x n ( n h ) + ( t − n h ) α α Γ ( α ) ( x n ( n h ) ( r − p x n ( n h ) − q y n ( n h ) ) ) , y n + 1 ( t ) = y n ( n h ) + ( t − n h ) α α Γ ( α ) ( y n ( n h ) ( − q + x n ( n h ) ) ) .

If we let t → ( n + 1 ) h in system (7), then we will achieve the discretized version of the model (4):

(8) x n + 1 = x n + h α α Γ ( α ) ( x n ( r − p x n − q y n ) ) , y n + 1 = y n + h α α Γ ( α ) ( y n ( − q + x n ) ) .

We next investigate the local stability and bifurcation to the points of equilibrium for system (8).

Theorem 2.1

The point of equilibrium E 0 = ( 0 , 0 ) is a

Saddle point when 0 < h < 2 α Γ ( α ) q α ;
Source when h > 2 α Γ ( α ) q α ;
Non-hyperbolic when h = 2 α Γ ( α ) q α .

Proof

The Jacobian matrix computed at the point of equilibrium E 0 for the linearization of system (8) is given by

J ( E 0 ) = 1 + h α α Γ ( α ) r 0 0 1 − h α α Γ ( α ) q

and has eigenvalues λ 1 = 1 + h α α Γ ( α ) r and λ 2 = 1 − h α α Γ ( α ) q . Since r > 0 , then ∣ λ 1 ∣ > 1 . Now, since q > 0 , we have the following cases:

If 0 < h < 2 α Γ ( α ) q α , then ∣ λ 2 ∣ < 1 and E 0 is a saddle point.
If h > 2 α Γ ( α ) q α , then ∣ λ 2 ∣ > 1 and E 0 is a source.
If h = 2 α Γ ( α ) q α , then ∣ λ 2 ∣ = 1 and E 0 is a nonhyperbolic.□

Theorem 2.2

The point of equilibrium E 1 = r p , 0 is a

Sink when q − 2 α Γ ( α ) h α p < r < min 2 α Γ ( α ) h α , p q ;
Source when r < q − 2 α Γ ( α ) h α p and r > max 2 α Γ ( α ) h α , p q ;
Non-hyperbolic when r = q − 2 α Γ ( α ) h α p or r = 2 α Γ ( α ) h α or r = p q ;
Saddle point for the other values of parameters except those values in (i)–(iii).

Proof

The Jacobian matrix computed at the point of equilibrium E 1 of the linearization of system (8) is given by

J ( E 1 ) = 1 − h α α Γ ( α ) r − h α α Γ ( α ) q r p 0 1 + h α α Γ ( α ) r p − q

and has eigenvalues λ 1 = 1 − h α α Γ ( α ) r and λ 2 = 1 + h α α Γ ( α ) r p − q . Now, we have the following cases:

If q − 2 α Γ ( α ) h α p < r < min 2 α Γ ( α ) h α , p q , then ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ < 1 . Hence, E 1 is a sink.
If r < q − 2 α Γ ( α ) h α p and r > max 2 α Γ ( α ) h α , p q , then ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ > 1 . Hence, E 1 is a source.
if r = q − 2 α Γ ( α ) h α p or r = 2 α Γ ( α ) h α or r = p q , then ∣ λ 1 ∣ = 1 or ∣ λ 2 ∣ = 1 . Hence, E 1 is a nonhyperbolic.
For other values of parameters except those values in (i)–(iii), we have ∣ λ 1 ∣ < 1 and ∣ λ 2 ∣ > 1 (or ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 ). Hence, E 1 is a saddle point.□

Theorem 2.3

The positive point of equilibrium E 2 = q , r − p q q is local asymptotically stable if and only if

max 2 p q h α α Γ ( α ) − 4 q h α α Γ ( α ) 2 + p q , p q < r < p α Γ ( α ) h α + q .

Proof

The Jacobian matrix computed at the point of equilibrium E 2 for the linearization of system (8) is given by

J ( E 2 ) = 1 − p q h α α Γ ( α ) − q 2 h α α Γ ( α ) h α α Γ ( α ) r − p q q 1 .

Next, the trace and determinant of J ( E 2 ) are computed as follows:

(9) Tr ( J ( E 2 ) ) = 2 − p q h α α Γ ( α ) and Det ( J ( E 2 ) ) = 1 − p q h α α Γ ( α ) + q ( r − p q ) h α α Γ ( α ) 2 .

According to the Jury conditions [37], both eigenvalues of the Jacobian matrix J ( E 2 ) have modulus less than 1 if max 2 p q h α α Γ ( α ) − 4 q h α α Γ ( α ) 2 + p q , p q < r < p α Γ ( α ) h α + q .□

The following result is regarding the bifurcation of system (8). We refer the reader to [38] for more details about the major forms of bifurcations in two-dimensional maps.

Theorem 2.4

The positive point of equilibrium E 2 = q , r − p q q forfeits its stability through a

Transcritical bifurcation if r = p q ;
Flip bifurcation if r = 2 p q h α α Γ ( α ) − 4 q h α α Γ ( α ) 2 + p q ;
Neimark–Sacker bifurcation if r = p α Γ ( α ) h α + q .

Proof

From Eq. (9), we have

Tr ( J ( E 2 ) ) − Det ( J ( E 2 ) ) = 1 when r = p q . Hence, the point of equilibrium E 2 forfeits its stability through a transcritical bifurcation when r = p q .
− Tr ( J ( E 2 ) ) − Det ( J ( E 2 ) ) = 1 when r = 2 p q h α α Γ ( α ) − 4 q h α α Γ ( α ) 2 + p q . Hence, the point of equilibrium E 2 forfeits its stability through a flip bifurcation when r = 2 p q h α α Γ ( α ) − 4 q h α α Γ ( α ) 2 + p q .
Det ( J ( E 2 ) ) = 1 when r = p α Γ ( α ) h α + q . Hence, the point of equilibrium E 2 forfeits its stability through a Neimark–Sacker bifurcation if r = p α Γ ( α ) h α + q .□

3 The dynamics of the discretized conformable fractional-order prey–predator system

Several authors prefer to study the fractional-order models by using Caputo derivatives because the definition of Caputo-fractional derivatives shows the property of linearity and yields zero the output for a fractional derivative to a constant. However, this definition fails to show the usual properties of the derivative such as the chain rule, product rule, and quotient rule. More recently, Khalil et al. [39] proposed a novel local derivative operator recognized as a conformable fractional derivative, which is essentially an extension of the well-known limit-based derivative, and this definition satisfies all the usual properties of the integer-order derivatives [40]. We refer the reader to ref. [41] for a comparative study of all definitions of fractional-order derivatives.

The conformable fractional-order version of system (3) is presented as follows:

(10) D a α x ( t ) = r x ( t ) − p x 2 ( t ) − q x ( t ) y ( t ) , D a α y ( t ) = x ( t ) y ( t ) − q y ( t ) ,

where D a α is the fractional derivative of conformable-type and defined for a function f : [ a , ∞ ) → R , a ≥ 0 :

(11) D a α f ( t ) = lim ε → 0 f ( t + ε ( t − a ) 1 − α ) − f ( t ) ε , 0 < α < 1 .

From the aforementioned definition, it has been evinced in ref. [40] the following necessary fact:

(12) D a α f ( t ) = ( t − a ) 1 − α f ′ ( t ) .

In the following, we will adopt piecewise-constant approximation to discretize the model (10).

(13) D a α x ( t ) = x ( t ) r − p x ( t ) − q y t h h , D a α y ( t ) = y ( t ) − q + x t h h .

Applying the rule in Eq. (12) to the first equation in the system (13), for h > 0 and t ∈ [ n h , ( n + 1 ) h ) , n = 0 , 1 , … , gives the following Bernoulli differential equation:

(14) ( t − n h ) 1 − α d x ( t ) d t + ( q y ( n h ) − r ) x ( t ) = − p x 2 ( t ) .

We obtain by simplifying this equation

(15) − x ′ ( t ) x 2 ( t ) + ( r − q y ( n h ) ) x ( t ) ( t − n h ) 1 − α = p ( t − n h ) 1 − α .

Multiplying Eq. (15) by e ( r − q y ( n h ) ) ( t − n h ) α α , we have

(16) d d t 1 x ( t ) e ( r − q y ( n h ) ) ( t − n h ) α α = p ( t − n h ) 1 − α e ( r − q y ( n h ) ) ( t − n h ) α α , t ∈ [ n h , ( n + 1 ) h ) .

Integrating with respect to t on [ n h , t ) both sides of Eq. (16), we obtain

(17) 1 x ( t ) e ( r − q y ( n h ) ) ( t − n h ) α α − 1 x ( n h ) = p r − q y ( n h ) e ( r − q y ( n h ) ) ( t − n h ) α α − 1 .

Ultimately, let t → ( n + 1 ) h in Eq. (17) and replacing x ( n h ) with x n yields

(18) x n + 1 = x n ( r − q y n ) p x n + ( r − q y n − p x n ) e − ( r − q y n ) h α α .

In a similar way, from the second equation in the model (13), we obtain

(19) d y ( t ) y ( t ) = ( x ( n h ) − q ) ( t − n h ) 1 − α d t .

Integrating with respect to t on [ n h , t ) both sides of Eq. (19), we obtain

(20) ln y ( t ) − ln y ( n h ) = ( x ( n h ) − q ) ( t − n h ) α α , t ∈ [ n h , ( n + 1 ) h ) .

For t → ( n + 1 ) h in Eq. (20) and replacing y ( n h ) with y n provides

(21) y n + 1 = y n e ( x n − q ) h α α .

Consequently, the discretized version of the model (13) is derived as follows:

(22) x n + 1 = x n ( r − q y n ) p x n + ( r − q y n − p x n ) e − ( r − q y n ) h α α , y n + 1 = y n e ( x n − q ) h α α .

We next investigate the local stability and bifurcation to the points of equilibrium for system (22).

Theorem 3.1

The point of equilibrium E 0 = ( 0 , 0 ) is a saddle point.

Proof

The Jacobian matrix computed at the point of equilibrium E 0 for the linearization of system (22) is given by

J ( E 0 ) = e r h α α 0 0 e − q h α α

and has eigenvalues λ 1 = e r h α α and λ 2 = e − q h α α . Thus, E 0 is a saddle point since ∣ λ 1 ∣ > 1 and ∣ λ 2 ∣ < 1 .□

Theorem 3.2

The point of equilibrium E 1 = r p , 0 is a

Saddle point if r > p q ;
Sink if r < p q ;
Non-hyperbolic if r = p q .

Proof

The Jacobian matrix computed at the point of equilibrium E 1 for the linearization of system (22) is given by

J ( E 1 ) = e − r h α α q p − 1 + e − r h α α 0 e r p − q h α α ,

and has eigenvalues λ 1 = e − r h α α , which satisfy ∣ λ 1 ∣ < 1 , and λ 2 = e r p − q h α α . Now, we have the following cases:

If r > p q , then ∣ λ 2 ∣ > 1 and E 1 is a saddle point.
If r < p q , then ∣ λ 2 ∣ < 1 and E 1 is a sink.
If r = p q , then ∣ λ 2 ∣ = 1 and E 1 is a nonhyperbolic.□

Theorem 3.3

The positive point of equilibrium E 2 = q , r − p q q is local asymptotically stable if and only if r < p α h α + q .

Proof

The Jacobian matrix computed at the point of equilibrium E 2 for the linearization of system (22) is given by

J ( E 2 ) = e − p q h α α q p ( − 1 + e − p q h α α ) ( r − p q ) h α q α 1 .

The trace and determinant of J ( E 2 ) are computed as follows:

(23) Tr ( J ( E 2 ) ) = 1 + e − p q h α α and Det ( J ( E 2 ) ) = p α e − p q h α α − h α ( r − p q ) − 1 + e − p q h α α p α .

According to the Jury conditions, both eigenvalues of the Jacobian matrix J ( E 2 ) have modulus less than 1 if and only if r < p α h α + q .□

The following result is regarding the bifurcation for model (22).

Theorem 3.4

The positive point of equilibrium E 2 = q , r − p q q forfeits its stability through a Neimark–Sacker bifurcation if r = p α h α + q .

Proof

From Eq. (23), we have Det ( J ( E 2 ) ) = 1 when r = p α h α + q . Hence, the point of equilibrium E 2 forfeits its stability through a Neimark–Sacker bifurcation if r = p α h α + q .□

4 Numerical simulations

Theoretical studies cannot be verified without numerical investigation of the obtained results. In the present study, numerical computations have been accomplished through the use of MATLAB-R2020a software. Let p = 0.5 and q = 1 be fixed and vary α , h , and r . Suppose that the initial state of systems (8) and (22) is ( 0.25 , 0.2 ) . Figures 1 and 2 demonstrate the local stable dynamic behaviors for the two-dimensional discrete systems (8) and (22), respectively, at the positive point of equilibrium E 2 with lessening the fractional-order parameter α . We note that lessening the fractional-order parameter α and fixing the discretization parameter h lead to destabilize the two-dimensional discrete systems and chaotic behavior occurs. Figures 3 and 4 show the local stable dynamic behaviors for the two-dimensional discrete systems (8) and (22), respectively, at the positive point of equilibrium E 2 with rising the discretization parameter h . We note that rising the discretization parameter h along with fixed fractional-order parameter α leads to destabilize the two-dimensional discrete systems and chaotic behavior occurs.

$Figure 1 Stable dynamical behavior of system (8) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) \left(x\left(0),y\left(0))=\left(0.25,0.2) for the parameters p = 0.5 , q = 1 , h = 0.15 , p=0.5,q=1,h=0.15, and r = 1 r=1 : (a) α = 0.95 \alpha =0.95 , (b) α = 0.75 \alpha =0.75 , (c) α = 0.6 \alpha =0.6 , and (d) α = 0.5 \alpha =0.5 .$

Figure 1

Stable dynamical behavior of system (8) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) for the parameters p = 0.5 , q = 1 , h = 0.15 , and r = 1 : (a) α = 0.95 , (b) α = 0.75 , (c) α = 0.6 , and (d) α = 0.5 .

$Figure 2 Stable dynamical behavior of system (22) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) \left(x\left(0),y\left(0))=\left(0.25,0.2) for the parameters p = 0.5 , q = 1 , h = 0.15 , p=0.5,q=1,h=0.15, and r = 1 r=1 : (a) α = 0.95 \alpha =0.95 , (b) α = 0.75 \alpha =0.75 , (c) α = 0.6 \alpha =0.6 , and (d) α = 0.5 \alpha =0.5 .$

Figure 2

Stable dynamical behavior of system (22) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) for the parameters p = 0.5 , q = 1 , h = 0.15 , and r = 1 : (a) α = 0.95 , (b) α = 0.75 , (c) α = 0.6 , and (d) α = 0.5 .

$Figure 3 Stable dynamical behavior of system (8) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) \left(x\left(0),y\left(0))=\left(0.25,0.2) for the parameters p = 0.5 , q = 1 , α = 0.95 , p=0.5,q=1,\alpha =0.95, and r = 1 r=1 : (a) h = 0.15 h=0.15 , (b) h = 0.35 h=0.35 , (c) h = 0.55 h=0.55 , and (d) h = 0.75 h=0.75 .$

Figure 3

Stable dynamical behavior of system (8) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) for the parameters p = 0.5 , q = 1 , α = 0.95 , and r = 1 : (a) h = 0.15 , (b) h = 0.35 , (c) h = 0.55 , and (d) h = 0.75 .

$Figure 4 Stable dynamical behavior of system (22) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) \left(x\left(0),y\left(0))=\left(0.25,0.2) for the parameters p = 0.5 , q = 1 , α = 0.95 , p=0.5,q=1,\alpha =0.95, and r = 1 r=1 : (a) h = 0.15 h=0.15 , (b) h = 0.35 h=0.35 , (c) h = 0.55 h=0.55 , and (d) h = 0.75 h=0.75 .$

Figure 4

Stable dynamical behavior of system (22) subject to the initial condition ( x ( 0 ) , y ( 0 ) ) = ( 0.25 , 0.2 ) for the parameters p = 0.5 , q = 1 , α = 0.95 , and r = 1 : (a) h = 0.15 , (b) h = 0.35 , (c) h = 0.55 , and (d) h = 0.75 .

The bifurcation diagrams in Figures 5 and 6 demonstrate that rising the values of r may destabilize the point of equilibrium E 2 through a Neimark–Sacker bifurcation. The maximum Lyapunov exponents corresponding to Figures 5(a) and 6(a) are given in Figure 7. The chaotic attractor for system (8) and system (22) is presented in Figure 8. The chaotic behavior exists when lessening the fractional-order parameter α and rising the discretization parameter h .

$Figure 5 Bifurcation diagram of system (8) as a function of r r for the parameters p = 0.5 , q = 1 , p=0.5,q=1, and h = 0.15 h=0.15 : (a) α = 0.95 \alpha =0.95 and (b) α = 0.5 \alpha =0.5 .$

Figure 5

Bifurcation diagram of system (8) as a function of r for the parameters p = 0.5 , q = 1 , and h = 0.15 : (a) α = 0.95 and (b) α = 0.5 .

$Figure 6 Bifurcation diagram of system (22) as a function of r r for the parameters p = 0.5 , q = 1 p=0.5,q=1 , and h = 0.15 h=0.15 : (a) α = 0.95 \alpha =0.95 and (b) α = 0.5 \alpha =0.5 .$

Figure 6

Bifurcation diagram of system (22) as a function of r for the parameters p = 0.5 , q = 1 , and h = 0.15 : (a) α = 0.95 and (b) α = 0.5 .

Figure 7

Maximum Lyapunov exponents corresponding to: (a) Figure 5(a) and (b) Figure 6(a).

$Figure 8 Chaotic attractor for the parameters p = 0.5 p=0.5 , q = 1 q=1 , α = 0.95 \alpha =0.95 , h = 0.15 h=0.15 , and r = 4 r=4 : (a) system (8) and (b) system (22).$

Figure 8

Chaotic attractor for the parameters p = 0.5 , q = 1 , α = 0.95 , h = 0.15 , and r = 4 : (a) system (8) and (b) system (22).

5 Discussion and concluding remarks

In the present work, fractional-order Lotka–Volterra models for two fractional-order derivatives types are considered. Discretization process by means of the piecewise-constant approximation is applied, and discrete versions of these systems are obtained. The local stability of the points of equilibrium of these discrete systems is investigated. Moreover, the necessary and sufficient asymptotically stable condition for the point of equilibrium E 2 of these discrete systems is obtained. Numerical simulations are presented to bolster the analytical results. To further confirm the chaos, we plot the time series of the systems (8) and (22), see Figure 9. It is clear that when the bifurcation parameter r is increasing, it leads to chaotic behavior in the systems. The findings of the current work can be summarized in the following points.

It can be seen that when lessening the fractional-order parameter α and fixing the discretization parameter h (or rising the discretization parameter h and fixing the fractional-order parameter α ), the discrete Caputo-conformable model is destabilized and chaotic behavior occurs (Figure 8). Thus, the discretized system is stabilized only for relatively large fractional order α ( α tends to one) and for relatively small discretization parameter h ( h tends to zero).
It can be deduced that when the Caputo-derivative acts on the fractional-order Lotka–Volterra model, the discrete Caputo-system exhibits richer dynamic behaviors than the discrete conformable-model.
It can be observed that the discrete conformable-system forfeits its stability faster than the discrete Caputo-model (Figures 1–6).
By using bifurcation theory, we showed that the discretized Caputo-system undergoes transcritical, flip, and Neimark–Sacker bifurcation at the positive point of equilibrium E 2 . While, the discrete conformable-system undergoes only a Neimark–Sacker bifurcation.

Figure 9

Time series plot of system (8) with respect to Figure 5(a): (a) asymptotically stable for r = 2 and (b) chaotic for r = 4 , and system (22) with respect to Figure 6(a): (c) asymptotically stable for r = 2 and (d) chaotic for r = 4 .

On conclusion, these fractional derivatives act to some extent the same function in importing some of the inherited properties of the time fractional to the time-integer Lotka–Volterra prey–predator models.

Funding information: The authors state no funding involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors declare that they have no known conflict of interest that could appear to influence the work presented in this article.

References

[1] Xu C, Yu Y. Stability analysis of time delayed fractional order predator-prey system with Crowley-Martin functional response. J App Anal Comp. 2019;9(3):928–42. 10.11948/2156-907X.20180175Search in Google Scholar

[2] Liu R, Liu G. Dynamics of a stochastic three species prey–predator model with intraguild predation. J App Anal Comp. 2020;10(1):81–103. 10.11948/jaac20190002Search in Google Scholar

[3] Lotka AJ. Elements of physical biology. Baltimore: Williams and Wilkins; 1925. Search in Google Scholar

[4] Volterra V. Fluctuations in the abundance of a species considered mathematically. Nature. 1926;118:558–60. 10.1038/119012b0Search in Google Scholar

[5] De Boer RJ. Modeling Population Dynamics: a Graphical Approach. Utrecht, Netherlands: Utrecht University; 2006. Search in Google Scholar

[6] Leeuwen EV, Jansen VAA, Bright PW. How population dynamics shape the functional response in a one-predator-two-prey system. Ecology. 2007;88(6):1571–81. 10.1890/06-1335Search in Google Scholar PubMed

[7] Elsadany AA, Matouk AE. Dynamical behaviors of fractional-order Lotka–Volterra predator-prey model and its discretization. J App Math Comp. 2015;49(1):269–83. 10.1007/s12190-014-0838-6Search in Google Scholar

[8] Yousef F, Alquran M, Jaradat I, Momani S, Baleanu D. Ternary-fractional differential transform schema: theory and application. Adv Differ Equ. 2019;2019:197. 10.1186/s13662-019-2137-xSearch in Google Scholar

[9] Yousef F, Alquran M, Jaradat I, Momani S, Baleanu D. New fractional analytical study of three-dimensional evolution equation equipped with three memory indices. J Comput Nonlinear Dynam. 2019;14(11):111008. 10.1115/1.4044585Search in Google Scholar

[10] Du M, Wang Z, Hu H. Measuring memory with the order of fractional derivative. Sci Rep. 2013;3(1):1–3. 10.1038/srep03431Search in Google Scholar PubMed PubMed Central

[11] Alquran M, Yousef F, Alquran F, Sulaiman TA, Yusuf A. Dual-wave solutions for the quadratic-cubic conformable-Caputo time-fractional Klein-Fock-Gordon equation. Math Comp Simu. 2021;185:62–76. 10.1016/j.matcom.2020.12.014Search in Google Scholar

[12] Yavuz M, Sene N. Stability analysis and numerical computation of the fractional predator-prey model with the harvesting rate. Frac Fract. 2020;4(3):35. 10.3390/fractalfract4030035Search in Google Scholar

[13] Naik PA, Eskandari Z, Shahraki HE. Flip and generalized flip bifurcations of a two-dimensional discrete-time chemical model. Math Model Numer Simu Appl. 2021;1(2):95–101. 10.53391/mmnsa.2021.01.009Search in Google Scholar

[14] Naik PA, Yavuz M, Qureshi S, Zu J, Townley S. Modeling and analysis of COVID-19 epidemics with treatment in fractional derivatives using real data from Pakistan. Eur Phys J Plus. 2020;135(10):1–42. 10.1140/epjp/s13360-020-00819-5Search in Google Scholar PubMed PubMed Central

[15] Yavuz M, Coşar FÖ, Günay F, Özdemir FN. A new mathematical modeling of the COVID-19 pandemic including the vaccination campaign. Open J Model Simu. 2021;9(3):299–321. 10.4236/ojmsi.2021.93020Search in Google Scholar

[16] Allegretti S, Bulai IM, Marino R, Menandro MA, Parisi K. Vaccination effect conjoint to fraction of avoided contacts for a Sars-Cov-2 mathematical model. Math Model Numer Simu Appl. 2021;1(2):56–66. 10.53391/mmnsa.2021.01.006Search in Google Scholar

[17] Özköse F, Yavuz M. Investigation of interactions between COVID-19 and diabetes with hereditary traits using real data: a case study in Turkey. Comp Bio Med. 2022;141:105044. 10.1016/j.compbiomed.2021.105044Search in Google Scholar PubMed

[18] Joshi H, Jha BK. Chaos of calcium diffusion in Parkinsonas infectious disease model and treatment mechanism via Hilfer fractional derivative. Math Model Numer Simu Appl. 2021;1(2):84–94. 10.53391/mmnsa.2021.01.008Search in Google Scholar

[19] Özköse F, Senel MT, Habbireeh R. Fractional-order mathematical modelling of cancer cells-cancer stem cells-immune system interaction with chemotherapy. Math Model Numer Simu Appl. 2021;1(2):67–83. 10.53391/mmnsa.2021.01.007Search in Google Scholar

[20] Naik PA, Owolabi KM, Yavuz M, Zu J. Chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells. Chaos Solit Fract. 2020;140:110272. 10.1016/j.chaos.2020.110272Search in Google Scholar

[21] Yavuz M, Özdemir N. Analysis of an epidemic spreading model with exponential decay law. Math Scie Appl E-Notes. 2020;8(1):142–54. 10.36753/mathenot.691638Search in Google Scholar

[22] Gurcan F, Kaya G, Kartal S. Conformable fractional order lotka-volterra predator-prey model: discretization, stability and bifurcation. J Comput Nonlinear Dynam. 2019;14(11):111007. 10.1115/1.4044313Search in Google Scholar

[23] Wang Z, Xie Y, Lu J, Li Y. Stability and bifurcation of a delayed generalized fractional-order prey–predator model with interspecific competition. App Math Comp. 2019;347:360–9. 10.1016/j.amc.2018.11.016Search in Google Scholar

[24] Ahmed E, Elgazzar AS. On fractional order differential equations model for nonlocal epidemics. Phys A. 2007;379:607–14. 10.1016/j.physa.2007.01.010Search in Google Scholar PubMed PubMed Central

[25] Ahmed E, El-Sayed AMA, El-Saka HAA. Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models. J Math Anal Appl. 2007;325:542–53. 10.1016/j.jmaa.2006.01.087Search in Google Scholar

[26] Jaradat I, Alquran M, Sulaiman TA, Yusuf A. Analytic simulation of the synergy of spatial-temporal memory indices with proportional time delay. Chaos Solit Fract. 2022;156:111818. 10.1016/j.chaos.2022.111818Search in Google Scholar

[27] Alquran M, Alsukhour M, Ali M, Jaradat I. Combination of Laplace transform and residual power series techniques to solve autonomous n-dimensional fractional nonlinear systems. Nonlinear Eng. 2021;10(1):282–92. 10.1515/nleng-2021-0022Search in Google Scholar

[28] Hammouch Z, Yavuz M, Özdemir N. Numerical solutions and synchronization of a variable-order fractional chaotic system. Math Model Numer Simu Appl. 2021;1(1):11–23. 10.53391/mmnsa.2021.01.002Search in Google Scholar

[29] Kumar P, Erturk VS. Dynamics of cholera disease by using two recent fractional numerical methods. Math Model Numer Simu Appl. 2021;1(2):102–11. 10.53391/mmnsa.2021.01.010Search in Google Scholar

[30] Maayah B, Yousef F, Arqub OA, Momani S, Alsaedi A. Computing bifurcations behavior of mixed type singular time-fractional partial integrodifferential equations of Dirichlet functions types in Hilbert space with error analysis. Filomat. 2019;33(12):3845–53. 10.2298/FIL1912845MSearch in Google Scholar

[31] Jaradat I, Alquran M, Katatbeh Q, Yousef F, Momani S, Baleanu D. An avant-garde handling of temporal-spatial fractional physical models. Int J Nonlinear SciNumer Simu. 2020;21(2):183–94. 10.1515/ijnsns-2018-0363Search in Google Scholar

[32] Jaradat I, Alquran M, Yousef F, Momani S, Baleanu D. On (2+1)-dimensional physical models endowed with decoupled spatial and temporal memory indices. Eur Phys J Plus. 2019;134(7):360. 10.1140/epjp/i2019-12769-8Search in Google Scholar

[33] Momani S, Arqub OA, Maayah B, Yousef F, Alsaedi A. A reliable algorithm for solving linear and nonlinear Schrödinger equations. Appl Comput Math. 2018;17(2):151–60. Search in Google Scholar

[34] Yousef F, Alkam O, Saker I. The dynamics of new motion styles in the time-dependent four-body problem: weaving periodic solutions. Eur Phys J Plus. 2020;135(9):742. 10.1140/epjp/s13360-020-00774-1Search in Google Scholar

[35] Yousef F, Momani S, Abdalmohsen R. Analytic solution of spatial-temporal fractional Klein-Gordon equation arising in physical models. Proc Int Conf Frac Differ Appl. 2018:1–4. 10.2139/ssrn.3277393.Search in Google Scholar

[36] Kartal S, Gurcan F. Discretization of conformable fractional differential equations by a piecewise constant approximation. Int J Comput Math. 2019;96(9):1849–60. 10.1080/00207160.2018.1536782Search in Google Scholar

[37] Kot M. Elements of mathematical ecology. Cambridge: Cambridge University Press; 2001. 10.1017/CBO9780511608520Search in Google Scholar

[38] Elaydi S. Discrete Chaos: with Applications in Science and Engineering. Boca Raton: Chapman and Hall/CRC; 2008. 10.1201/9781420011043Search in Google Scholar

[39] Khalil R, AlHorani M, Yousef A, Sababheh M. A new definition of fractional derivative. J Compu Appl Math. 2014;264:65–70. 10.1016/j.cam.2014.01.002Search in Google Scholar

[40] Abdeljawad T. On conformable fractional calculus. J Comput Appl Math. 2015;279:57–66. 10.1016/j.cam.2014.10.016Search in Google Scholar

[41] Teodoro GS, TenreiroMachado JA, De Oliveira EC. A review of definitions of fractional derivatives and other operators. J Comput Phys. 2019;338:195–208. 10.1016/j.jcp.2019.03.008Search in Google Scholar

Received: 2022-01-22

Revised: 2022-03-06

Accepted: 2022-03-14

Published Online: 2022-04-02

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/nleng-2022-0013

Keywords for this article

Caputo fractional derivative; conformable fractional derivative; prey–predator system; discretization; stability; bifurcations; chaos

Creative Commons

BY 4.0