Next Article in Journal
Eigenvalue −1 of the Vertex Quadrangulation of a 4-Regular Graph
Next Article in Special Issue
Theoretical Investigation of Fractional Estimations in Liouville–Caputo Operators of Mixed Order with Applications
Previous Article in Journal
On the Interpolating Family of Distributions
Previous Article in Special Issue
A Note on the Time-Fractional Navier–Stokes Equation and the Double Sumudu-Generalized Laplace Transform Decomposition Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Note on a Fractional Extension of the Lotka–Volterra Model Using the Rabotnov Exponential Kernel

by
Mohamed M. Khader
1,2,*,
Jorge E. Macías-Díaz
3,4,*,
Alejandro Román-Loera
5,* and
Khaled M. Saad
6,7
1
Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh 11566, Saudi Arabia
2
Department of Mathematics, Faculty of Science, Benha University, Benha 13518, Egypt
3
Department of Mathematics, School of Digital Technologies, Tallinn University, Narva Rd. 25, 10120 Tallinn, Estonia
4
Department of Mathematics and Physics, Autonomous University of Aguascalientes, Ave. Universidad 940, Ciudad Universitaria, Aguascalientes 20100, Mexico
5
Department of Electronics, Autonomous University of Aguascalientes, Ave. Universidad 940, Ciudad Universitaria, Aguascalientes 20100, Mexico
6
Department of Mathematics, College of Sciences and Arts, Najran University, Najran P.O. Box 1988, Saudi Arabia
7
Department of Mathematics, Faculty of Applied Science, Taiz University, Taiz P.O. Box 6803, Yemen
*
Authors to whom correspondence should be addressed.
Axioms 2024, 13(1), 71; https://doi.org/10.3390/axioms13010071
Submission received: 1 December 2023 / Revised: 4 January 2024 / Accepted: 12 January 2024 / Published: 21 January 2024
(This article belongs to the Special Issue Fractional Calculus - Theory and Applications II)

Abstract

:
In this article, we study the fractional form of a well-known dynamical system from mathematical biology, namely, the Lotka–Volterra model. This mathematical model describes the dynamics of a predator and prey, and we consider here the fractional form using the Rabotnov fractional-exponential (RFE) kernel. In this work, we derive an approximate formula of the fractional derivative of a power function ζ p in terms of the RFE kernel. Next, by using the spectral collocation method (SCM) based on the shifted Vieta–Lucas polynomials (VLPs), the fractional differential system is reduced to a set of algebraic equations. We provide a theoretical convergence analysis for the numerical approach, and the accuracy is verified by evaluating the residual error function through some concrete examples. The results are then contrasted with those derived using the fourth-order Runge-Kutta (RK4) method.

1. Introduction

Fractional calculus has been a topic of interest for researchers who have found interesting applications in biology, physics, engineering, fluid mechanics, and viscoelasticity [1,2]. To that end, the usefulness of the different types of fractional derivatives has been examined. Some of those derivatives are the Caputo and the Riemann–Liouville operators, which are fractional operators with a power kernel. In the present work, we will employ a different fractional operator, namely, the so-called Yang–Aty–Cattani (YAC) operator (see, e.g., [3,4]). This fractional derivative is also called the Prabhakar operator or generalized Mittag–Leffler operator. Among those fractional operators, we will employ the YAC derivative based on the Rabotnov function, which is a generalized Mittag–Leffler function. The derivative with an exponential kernel is called the Caputo–Fabrizio derivative, and the Atangana–Baleanu derivative if the kernel is the Mittag–Leffler function. These types of non-singular derivatives have applications to groundwater flow [5], medical sciences [6], chaos theory [7] and other areas [8].
In the context of approximating fractional derivatives numerically, we will employ a spectral collocation technique that relies on specific polynomials. One noteworthy characteristic of these approaches is their capacity to yield precise outcomes with minimal errors. These methods heavily leverage the significance of polynomials, and the orthogonality characteristics of these polynomials are harnessed for the approximation of periodic functions within a closed and bounded interval on the real number line [9,10].
Certainly, this holds true for VLPs. The objective of this work is to assess the precision of the numerical approach when employing the fractional derivative with the RFE kernel. We provide an estimated formula for the fractional-order derivative of the power function ζ p using the RFE kernel. Furthermore, we utilize this approximative formula along with the characteristics of VLPs to calculate the numerical solution for the proposed system.
It is shown that this method is efficient and that there is an excellent agreement with the currently available exact solutions. Finally, we will conclude that the non-singular operator is better suited for numerical simulations of the mathematical model. This conclusion is drawn from a comparison with the RK4 method.
The predator–prey equations, a well-known set of nonlinear first-order differential equations, are employed as a basis for illustrating the dynamics of specific biological systems characterized by the interaction between two species one functioning as a predator and the other as prey. Samardzija [11] extended this model by incorporating the concept of solitary prey being targeted by predators within the Lotka–Volterra system. Our work in this paper centers on the examination of this system.
The paper is organized as follows: In Section 2, we introduce some fundamental concepts related to fractional-order integrals and derivatives, along with information about Vieta–Lucas polynomials. Section 3 contains the numerical results and convergence analysis for the proposed model, with a particular focus on discussing the numerical simulation. The concluding remarks, observations, and considerations for future work are presented in Section 4.

2. Preliminaries

2.1. Definitions

We start this subsection by providing some remarks and recalling definitions of some useful fractional derivatives. For further reading, we refer the reader to [12,13]. Throughout, we let H 1 ( 0 , b ) be the classical Sobolev space. With this convention, the following definition serves to recall the well-known Caputo derivative [14].
Definition 1. 
The Caputo fractional derivative   C D ν of order 0 < ν 1 for φ ( ζ ) H 1 ( 0 , b ) is given by
  C D ν φ ( ζ ) = 1 Γ ( 1 ν ) 0 ζ φ ( τ ) ( ζ τ ) ν d τ , ζ > 0 ,
where Γ ( . ) is the usual Gamma function.
It is evident from this definition that the integral depends on the singular kernel ( ζ τ ) ν . Here, it is worth mentioning that other fractional-order operators have been designed to avoid the singularity of the kernel. For example, the following definition redefines the Caputo derivative using the well-known Rabotnov kernel.
Definition 2. 
The Caputo fractional derivative of order ν on the interval [ 0 , 1 ] for a function Θ ( ζ ) using the Rabotnov kernel is defined by
  R F E D ν Θ ( ζ ) = 0 ζ Θ ( n ) ( τ ) R ν [ Ω ( ζ τ ) ν ] d τ , n 1 < ν n .
Here, Ω R + , and the RFE function R ν is defined as
R ν [ Ω ( ζ ) ν ] = κ = 0 ( Ω ) κ ζ ( κ + 1 ) ( ν + 1 ) 1 Γ [ ( κ + 1 ) ( ν + 1 ) ] .
For more detail on the RFE-operator derivative, we refer to [3,4].
As we mentioned above, the fractional derivative (2) is based on a Rabotnov kernel, and it will be used in this manuscript to study a predator–prey model. However, the derivation of an approximate solution for the fractional dynamical system will require approximating the fractional derivative (2). Obviously, the analytical calculation of that derivative is, in general, a difficult task. For that reason, we concentrate only on the numerical approximation of derivatives of powers. In fact, the approximate formula of the fractional derivative concerning the RFE-fractional is provided in the next theorem. The proof is based on the use of the Simpson integration scheme.
Theorem 1 
(Yang et al. [4]). The RFE derivative of order ν of ζ p with p n , ( n = ν ) is
  R F E D ν ζ p = Γ ( p + 1 ) Γ ( p + 1 ν ) × h 3 [ Ψ ν , p ( ζ , 0 ) + Ψ ν , p ( ζ , m ) + 2 κ = 1 m 2 Ψ ν , p ( ζ , 2 κ ) + 4 κ = 1 m 1 Ψ ν , p ( ζ , 2 κ + 1 ) ] .
The interval [ 0 , 1 ] is divided into m equal subintervals with length h, · is the smallest integer function, κ = κ m , for each κ = 0 , 1 , , m , and
Ψ ν , p ( ζ , ) = p ν R ν [ Ω ( ζ ) ν ] .

2.2. Vieta–Lucas Polynomials

To use a spectral method for solving the fractional dynamical system of interest, we need to express the unknown functions through series in terms of some orthogonal functions. In this work, we will use Vieta–Lucas polynomials [15]. In this stage, we will recall now some important properties of the shifted Vieta–Lucas polynomials (see [15] for more properties).
The VLPs of degree κ N 0 is denoted by V L κ ( z ) , and is defined as
VL κ ( z ) = 2 cos ( κ ψ ) , ψ = arccos ( 0.5 z ) .
It is obvious that ψ [ 0 , π ] and | z | 2 . Moreover, it is easy to prove that the polynomials VL k ( z ) satisfy VL 0 ( z ) = 2 and VL 1 ( z ) = z . Moreover, the following recurrence relation is satisfied, for each κ = 2 , 3 ,
VL k ( z ) = z VL k 1 ( z ) VL k 2 ( z ) .
Notice that if z = 4 ζ 2 , we obtain a new class of orthogonal polynomials on [ 0 , 1 ] . We will denote these polynomials using the notation VL m s ( ζ ) , and they satisfy the following relation:
VL κ + 1 s ( ζ ) = ( 4 ζ 2 ) VL κ 1 s ( ζ ) VL κ 2 s ( ζ ) , κ = 2 , 3 , ,
where VL 0 s ( ζ ) = 2 and VL 1 s ( ζ ) = 4 ζ 2 . Moreover, notice that VL κ s ( 0 ) = 2 ( 1 ) κ and VL κ s ( 1 ) = 2 , for each κ = 0 , 1 , 2 ,
These polynomials can be written also in an analytic form. More precisely, let v κ ( ζ ) represent V L κ s ( ζ ) . Under this convention, the following equation is satisfied:
v κ ( ζ ) = 2 κ j = 0 κ ( 1 ) j 4 κ j Γ ( 2 κ j ) Γ ( j + 1 ) Γ ( 2 κ 2 j + 1 ) ζ κ j , κ = 2 , 3 ,
Notice that v κ ( ζ ) are orthogonal polynomials on [ 0 , 1 ] with respect to the weight function 1 ζ ζ 2 . As a consequence,
v i ( ζ ) , v j ( ζ ) = 0 1 v i ( ζ ) v j ( ζ ) ζ ζ 2 d ζ = 0 , i j 0 , 4 π , i = j = 0 , 2 π , i = j 0 .
These polynomials form a basis for the Euclidean space L 2 ( [ 0 , 1 ] ) . It follows that, for each ψ ( ζ ) L 2 [ 0 , 1 ] , the following identity is satisfied for suitable constants c j :
ψ ( ζ ) = j = 0 c j v j ( ζ ) .
The constants c j must be evaluated to express ψ ( ζ ) in terms of the series on v j ( ζ ) . By considering only the first m + 1 terms in the series of (10), an approximation for the function ψ ( ζ ) is given by
ψ m ( ζ ) = j = 0 m c j v j ( ζ ) ,
where the coefficients c j , for j = 0 , 2 , , m , can be evaluated using the integrals
c j = 1 δ j 0 1 ψ ( ζ ) v j ( ζ ) ζ ζ 2 d ζ ,
and
δ j = 4 π , j = 0 , 2 π , j = 1 , 2 , , m .

3. Numerical Results

3.1. Convergence Analysis

The present subsection is devoted to providing the convergence analysis of our numerical approach. We begin this subsection by recalling a useful result from the literature.
Lemma 1 
(Yang et al. [4]). Assume that ψ ( ζ ) L w ˜ 2 ( [ 0 , 1 ] ) with respect to the weight function w ˜ ( ζ ) = 1 ζ ζ 2 , and that | ψ ( ζ ) | ε , for some constant ε. Then the partial sum (11) of the series (10) is uniformly convergent to ψ ( ζ ) as m . In addition, the following estimates are satisfied:
  • The coefficients in the series (10) are bounded, more precisely,
    c j ε 4 j j 2 1 , j > 2 .
  • An error estimate in the L w ˜ 2 ( [ 0 , 1 ] ) -norm is given by
    ψ ( ζ ) ψ m ( ζ ) w ˜ < ε 12 m 3 .
Theorem 2. 
Let the function ψ ( ζ ) [ 0 , 1 ] have continuous derivatives up to ( m + 1 ) -times and let ψ m ( ζ ) be the best approximation of the function ψ ( ζ ) defined in Equation (10). Then we have the following:
ψ ( ζ ) ψ m ( ζ ) α β m + 1 π ( m + 1 ) ! ,
where
α = max 0 ζ 1 ψ ( m + 1 ) ( ζ ) , β = max { ζ 0 , 1 ζ 0 } .
Proof. 
The Taylor series approximation of the function ψ ( ζ ) in the neighborhood of a point ζ = ζ 0 is given by
ψ ( ζ ) = ψ ( ζ 0 ) + ( ζ ζ 0 ) 1 ! ψ ( 1 ) ( ζ 0 ) + ( ζ ζ 0 ) 2 2 ! ψ ( 2 ) ( ζ 0 ) + + ( ζ ζ 0 ) m m ! ψ ( m ) ( ζ 0 ) + ( ζ ζ 0 ) m + 1 ( m + 1 ) ! ψ ( m + 1 ) ( τ ) ,
where ζ 0 [ 0 , 1 ] , τ ( ζ 0 , ζ ) . Assume
ψ ˜ ( ζ ) = ψ ( ζ 0 ) + ( ζ ζ 0 ) 1 ! ψ ( 1 ) ( ζ 0 ) + ( ζ ζ 0 ) 2 2 ! ψ ( 2 ) ( ζ 0 ) + + ( ζ ζ 0 ) m m ! ψ ( m ) ( ζ 0 ) ,
then
ψ ( ζ ) ψ ˜ m ( ζ ) = ( ζ ζ 0 ) m + 1 ( m + 1 ) ! ψ ( m + 1 ) ( τ ) .
If ψ m ( ζ ) is the best approximation of ψ ( ζ ) , then
ψ ( ζ ) ψ m ( ζ ) 2 ψ ( ζ ) ψ ˜ m ( ζ ) 2 = 0 1 w ( ζ ) ψ ( ζ ) ψ ˜ m ( ζ ) 2 d ζ 0 1 w ( ζ ) ( ζ ζ 0 ) m + 1 ( m + 1 ) ! ψ ( m + 1 ) ( τ ) 2 d ζ .
It is assumed that ψ ( ζ ) has continuous derivatives up to ( m + 1 ) -times. Therefore, there exists a constant α such that
α = max 0 ζ 1 ψ ( m + 1 ) ( ζ ) .
Now we have
ψ ( ζ ) ψ m ( ζ ) 2 0 1 w ( ζ ) α ( m + 1 ) ! ( ζ ζ 0 ) m + 1 2 d ζ .
Considering β = max { ζ 0 , 1 ζ 0 } , then Equation (22) becomes
ψ ( ζ ) ψ m ( ζ ) 2 α 2 β 2 ( m + 1 ) [ ( m + 1 ) ! ] 2 0 1 w ( ζ ) d ζ .
Since w ( ζ ) = 1 ζ ζ 2 , so with some simple computation for the integral, we get:
ψ ( ζ ) ψ m ( ζ ) α β m + 1 π ( m + 1 ) ! .
The proof is completed. □
Readers may refer to [16] for more information on the convergence analysis of the VLPs approximation (11). Now, we are in a position to compute the RFE-fractional derivative of the VLPs approximation (11).
Theorem 3. 
The RFE fractional derivative of order β for the function ψ i ( ζ ) defined in (11) can be calculated as follows:
  R F E D β ψ i ( ζ ) = j = β i χ i , j , β [ Ψ β , p ( ζ , 0 ) + Ψ β , p ( ζ , m ) + 2 κ = 1 m 2 Ψ β , p ( ζ , 2 κ ) + 4 κ = 1 m 1 Ψ β , p ( ζ , 2 κ + 1 ) ] ,
where
χ i , j , β = ( 1 ) j ( 2 i ) 4 i j Γ ( 2 i j ) Γ ( j + 1 ) Γ ( 2 i 2 j + 1 ) × h Γ ( i j + 1 ) 3 Γ ( i j + 1 β ) ,
and Ψ β , p ( ζ , ) = p β R β [ Ω ( ζ ) β ] p = i j .
Proof. 
Using Theorem 1, we obtain that
  R F E D β ζ i j = Γ ( i j + 1 ) Γ ( i j β + 1 ) × h 3 [ G β , p ( ζ , 0 ) + Ψ β , p ( ζ , m ) + 4 κ = 1 m 1 Ψ β , p ( ζ , 2 κ + 1 ) + 2 κ = 1 m 2 Ψ β , p ( ζ , 2 κ ) ] .
The interval [ 0 , 1 ] is divided into m equal subintervals with length h, which means that h = 1 m . Let κ = κ m , for each κ = 0 , 1 , , m , and let
Ψ β , p ( ζ , ) = p β R β [ Ω ( ζ ) β ] p = i j .
Now, using Equations (8) and (27), we can readily evaluate the RFE derivative of the i-th degree ψ i ( ζ ) as follows:
  R F E D β ψ i ( ζ ) = j = 0 i ( 1 ) j ( 2 i ) 4 i j Γ ( 2 i j ) Γ ( j + 1 ) Γ ( 2 i 2 j + 1 )   R F E D β ζ i j = j = β i Γ ( i j + 1 ) Γ ( i j + 1 β ) × ( 1 ) j ( 2 i ) 4 i j Γ ( 2 i j ) Γ ( j + 1 ) Γ ( 2 i 2 j + 1 ) × h 3 [ Ψ β , p ( ζ , 0 ) + Ψ β , p ( ζ , m ) + 4 κ = 1 m 1 Ψ β , p ( ζ , 2 κ + 1 ) + 2 κ = 1 m 2 Ψ β , p ( ζ , 2 k ) ] .
From this result, we can easily get the Equation (25), as desired. □

3.2. Computer Simulations

We start this subsection by describing the computer implementation of our model. The goal is to solve predator–prey equations that are used to describe the dynamics in some biological systems [17]. Our system has two species interacting with each other, one predator and the other prey. Recently, Samardzija [11] proposed an extension of this model, considering two predators and single a prey via the Lotka–Volterra system. Using the fractional operators of this work, the system reads as follows:
  R F E D ν Φ 1 ( ζ ) = σ 1 Φ 1 ( ζ ) σ 2 Φ 1 ( ζ ) Φ 2 ( ζ ) + σ 3 Φ 1 2 ( ζ ) σ 4 Φ 3 ( ζ ) Φ 1 2 ( ζ ) ,
  R F E D ν Φ 2 ( ζ ) = σ 5 Φ 2 ( ζ ) + σ 6 Φ 1 ( ζ ) Φ 2 ( ζ ) ,
  R F E D ν Φ 3 ( ζ ) = σ 7 Φ 3 ( ζ ) + σ 4 Φ 3 ( ζ ) Φ 1 2 ( ζ ) ,
with initial conditions
Φ 1 ( 0 ) = Φ 1 0 , Φ 2 ( 0 ) = Φ 2 0 , Φ 3 ( 0 ) = Φ 3 0 .
Here, population sizes of predators are Φ 1 ( ζ ) and Φ 2 ( ζ ) at the time ζ , and the population size of prey is Φ 3 ( ζ ) . Moreover, σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 , and σ 7 are parameters that describe the interaction between the three species [11].
Let us approximate the unknown functions Φ 1 ( ζ ) , Φ 2 ( ζ ) , and Φ 3 ( ζ ) through the shifted Vieta–Lucas polynomials Φ 1 , N ( ζ ) , Φ 2 , N ( ζ ) and Φ 3 , N ( ζ ) , respectively. To that end, we let
Φ 1 , N ( ζ ) = i = 0 N a i v i ( ζ ) , Φ 2 , N ( ζ ) = i = 0 N b i v i ( ζ ) , Φ 3 , N ( ζ ) = i = 0 N c i v i ( ζ ) .
Then, substituting Equations (25) and (34) into the system (30)–(32), we obtain:
j = ν N a i χ N , j , ν [ Ψ ν , p ( ζ , 0 ) + Ψ ν , p ( ζ , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ , 2 k + 1 ) + 2 k = 1 m 2 Ψ ν , p ( ζ , 2 k ) ] = σ 1 i = 0 N a i v i ( ζ ) σ 2 i = 0 N a i v i ( ζ ) i = 0 N b i v i ( ζ ) + σ 3 i = 0 N a i v i ( ζ ) 2 σ 4 i = 0 N c i v i ( ζ ) i = 0 N a i v i ( ζ ) 2 ,
j = ν N b i χ N , j , ν [ Ψ ν , p ( ζ , 0 ) + Ψ ν , p ( ζ , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ , 2 k + 1 ) + 2 k = 1 m 2 Ψ ν , p ( ζ , 2 k ) ] = σ 5 i = 0 N b i v i ( ζ ) + σ 6 i = 0 N a i v i ( ζ ) i = 0 N b i v i ( ζ ) ,
j = ν N c i χ N , j , ν [ Ψ ν , p ( ζ , 0 ) + Ψ ν , p ( ζ , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ , 2 k + 1 ) + 2 k = 1 m 2 Ψ ν , p ( ζ , 2 k ) ] = σ 7 i = 0 N c i v i ( ζ ) + σ 4 i = 0 N c i v i ( ζ ) i = 0 N a i v i ( ζ ) 2 .
By collocation of these three equations at N points ζ r , all of them being roots of v N ( ζ ) ), the equations will be reduced to
j = ν N a i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ r , 2 k + 1 ) + 2 k = 1 m 2 Ψ ν , p ( ζ r , 2 k ) ] = σ 1 i = 0 N a i v i ( ζ r ) σ 2 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) + σ 3 i = 0 N a i v i ( ζ r ) 2 σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 ,
j = ν N b i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ r , 2 k + 1 ) + 2 k = 2 m 2 Ψ ν , p ( ζ r , 2 k ) ] = σ 5 i = 0 N b i v i ( ζ r ) + σ 6 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) ,
j = ν N c i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 m 1 Ψ ν , p ( ζ r , 2 k + 1 ) + 2 k = 1 m 2 Ψ ν , p ( ζ r , 2 k ) ] = σ 7 i = 0 N c i v i ( ζ r ) + σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 .
Substituting Equation (34) in (33), the initial conditions (33) will be transformed into the following algebraic equations:
i = 0 N 2 ( 1 ) i a i = Φ 1 , 0 , i = 0 N 2 ( 1 ) i b i = Φ 2 , 0 , i = 0 N 2 ( 1 ) i c i = Φ 3 , 0 .
Notice that the system (38)–(41) can be expressed as a constrained optimization problem concerning the cost functions (CFs) defined by
C F 1 = r = 0 N | j = ν N a i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 , k o d d m 1 Ψ ν , p ( ζ r , k ) + 2 k = 2 , k e v e n m 2 Ψ ν , p ( ζ r , k ) ] σ 1 i = 0 N a i v i ( ζ r ) + σ 2 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) σ 3 i = 0 N a i v i ( ζ r ) 2 + σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 | ,
C F 2 = r = 0 N | j = ν N b i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 , k o d d m 1 Ψ ν , p ( ζ r , k ) + 2 k = 2 , k e v e n m 2 Ψ ν , p ( ζ r , k ) ] + σ 5 i = 0 N b i v i ( ζ r ) σ 6 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) | ,
C F 3 = r = 0 N | j = ν N c i χ N , j , ν [ Ψ ν , p ( ζ r , 0 ) + Ψ ν , p ( ζ r , m ) + 4 k = 1 , k o d d m 1 Ψ ν , p ( ζ r , k ) + 2 k = 2 , k e v e n m 2 Ψ ν , p ( ζ r , k ) ] + σ 7 i = 0 N c i v i ( ζ r ) σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 | ,
with the constraints (Cons)
Cons = | i = 0 N 2 ( 1 ) i a i Φ 1 0 | + | i = 0 N 2 ( 1 ) i b i Φ 2 0 | + | i = 0 N 2 ( 1 ) i c i Φ 3 0 | .
We implemented the Penalty Leap-Frog method [18] for solving the constrained optimization problem (42)–(45), for the coefficients a i , b i and c i , for each i = 0 , 1 , , N . In turn, this leads us to formulate the approximate solution by substitution into Equation (34). For more details concerning the Penalty Leap-Frog Method and Constrained Optimization Problem, the reader can see [19,20,21]. Where the Leap-Frog algorithm for solving the constrained optimization problem is given in [19], and the relationship between constrained optimal points and asymptotically stable critical points has been discussed in [20]. Finally, the convergence of the ODE method in constrained optimization is given in [21].
Hence, we define the residual error function (REF) with the help of the approximation Formula (25) as follows:
REF [ Φ 1 ( ζ ) ] =   R F E D ν Φ 1 ( ζ ) σ 1 i = 0 N a i v i ( ζ r ) + σ 2 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) σ 3 i = 0 N a i v i ( ζ r ) 2 + σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 ,
REF [ Φ 2 ( ζ ) ] =   R F E D ν Φ 2 ( ζ ) + σ 5 i = 0 N b i v i ( ζ r ) σ 6 i = 0 N a i v i ( ζ r ) i = 0 N b i v i ( ζ r ) ,
REF [ Φ 3 ( ζ ) ] =   R F E D ν Φ 3 ( ζ ) + σ 7 i = 0 N c i v i ( ζ r ) σ 4 i = 0 N c i v i ( ζ r ) i = 0 N a i v i ( ζ r ) 2 .
The approximation solution is always in near alignment with the precise solution, as indicated by a decreasing absolute relative error ( REF ( ζ ) 0 ). It is noteworthy to emphasize that the original equation model is considered to be well-posed, and the residual is a measure of how far the approximate solution deviates from the precise solution. We use this kind of error measurement in the setting of fractional order with ν , where the exact solution is unknown. It should be noted that REF takes on multiple forms; for a more thorough explanation, see [22].
Finally, we present some numerical results for the model (30)–(33). We study the effect of the parameters ν and N, letting σ 1 = σ 2 = σ 3 = σ 4 = 1 , σ 5 = 2 , σ 6 = 3 and σ 7 = 2.7 . As initial conditions, we choose Φ 1 0 = 1 , Φ 2 0 = 1.4 , Φ 3 0 = 1 .  Figure 1 shows the numerical solutions for different values of ν = 1.0 ,   0.95 ,   0.90 , with N = 6 . Meanwhile, Figure 2 shows the results of our simulations using different values of N = 6 ,   8 ,   10 at ν = 0.93 . To verify the numerical results, we will examine now the residual error function, which is calculated by replacing the approximate solutions in the form (30)–(33). Figure 3 shows the REF of the numerical solution at ν = 0.94 with N = 10 . A comparison of the results obtained using the suggested technique and those obtained using the RK4 method at ( ν = 1 ) with N = 7 is shown in Figure 4.
These findings elucidate that the present approach enhances both the efficiency and outcomes of the method. Additionally, the results affirm that the suggested method is apt for addressing the fractional-order representation of the proposed model. It is noteworthy that the numerical solution’s behavior, generated by the proposed method, is contingent on the values of N and ν , by the inherent characteristics of the model under investigation. In conclusion, we can assert with assurance that all the theoretical investigations conducted in this study have been successfully fulfilled.

4. Conclusions

In this work, we calculated a numerical solution for a predator–prey model of fractional derivatives based on the RFE. The model considers two predators and one prey; mathematically, it is a system of coupled ordinary differential equations of orders in the interval ( 0 , 1 ] . The approach followed in this work is based on the use of the shifted Vieta–Lucas polynomials, and it allowed us to validate the effectiveness of the proposed method. Furthermore, we noted that the precision of the error can be managed and diminished by integrating extra terms from the series of the approximate solution. Moreover, we verified that the RFE operator, devoid of singularity, proves to be more suitable for numerically simulating the model discussed in this article. The graphical results and those obtained through the RK4 approach are comparable. Clearly, various types of polynomials can be employed to examine the same system. Future endeavors aim to expand and generalize the topics this study encompasses, such as
  • Employing the finite element method for addressing the identical issue: Conducting theoretical investigations into the depth that elucidates the suggested model, incorporating optimal control of the results;
  • Adjusting the interpretation of the fractional derivative to alternatives such as the Atangana–Baleanu–Caputo model or a variable-order model.

Author Contributions

Conceptualization, M.M.K. and K.M.S.; methodology, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; software, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; validation, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; formal analysis, M.M.K. and K.M.S.; investigation, M.M.K. and K.M.S.; resources, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; data curation, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; writing—original draft preparation, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; writing—review and editing, M.M.K., J.E.M.-D., A.R.-L. and K.M.S.; visualization, M.M.K. and K.M.S.; supervision, J.E.M.-D., A.R.-L. and K.M.S.; project administration, J.E.M.-D., A.R.-L. and K.M.S.; funding acquisition, J.E.M.-D., A.R.-L. and K.M.S. All authors have read and agreed to the published version of the manuscript.

Funding

One of the authors (J.E.M.-D.) was funded by the National Council of Science and Technology of Mexico (CONACYT) through grant A1-S-45928.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kilbas, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives: Theory and Applications; Gordon & Breach: Philadelphia, PA, USA, 1993. [Google Scholar]
  2. Khader, M.M.; Saad, K.M. On the numerical evaluation for studying the fractional KdV, KdV-Burger’s, and Burger’s equations. Eur. Phys. J. Plus 2018, 133, 335. [Google Scholar] [CrossRef]
  3. Kumar, S.; Nisar, K.S.; Kumar, R.; Cattani, C.; Samet, B. A new Rabotnov fractional-exponential function-based fractional derivative for diffusion equation under external force. Math. Methods Appl. Sci. 2020, 47, 4460–4471. [Google Scholar] [CrossRef]
  4. Yang, X.J.; Abdel-Aty, M.; Cattani, C. A new general fractional-order derivative with Rabotnov fractional-exponential kernel applied to model the anomalous heat transfer. Therm. Sci. 2019, 23, 1677–1681. [Google Scholar] [CrossRef]
  5. Gao, W.; Ghanbari, B.; Baskonus, H.M. New numerical simulations for some real-world problems with Atangana-Baleanu fractional derivative. Chaos Solitons Fractals 2019, 128, 34–43. [Google Scholar] [CrossRef]
  6. Morales-Delgado, V.F.; Gomez-Aguilar, J.F.; Saad, K.M.; Escobar-Jimenez, R.F. Application of the Caputo-Fabrizio and Atangana-Baleanu fractional derivatives to the mathematical model of cancer chemotherapy effect. Math. Methods Appl. Sci. 2019, 42, 1167–1193. [Google Scholar] [CrossRef]
  7. Toufik, M.; Atangana, A. New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models. Eur. Phys. J. Plus 2017, 132, 444. [Google Scholar] [CrossRef]
  8. Kumar, D.; Singh, J.; Baleanu, D. On the analysis of vibration equation involving a derivative with Mittag-Leffler law. Math. Methods Appl. Sci. 2020, 43, 443–457. [Google Scholar] [CrossRef]
  9. Abd-Elhameed, W.M.; Youssri, Y.H. Connection formulae between generalized Lucas polynomials and some Jacobi polynomials: Application to certain types of fourth-order BVPs. Int. J. Appl. Comput. Math. 2020, 6, 45. [Google Scholar] [CrossRef]
  10. Agarwal, F.; El-Sayed, A.A. Vieta–Lucas polynomials for solving a fractional-order mathematical physics model. Adv. Differ. Equ. 2020, 2020, 626. [Google Scholar] [CrossRef]
  11. Samardzija, N.; Greller, L.D. Explosive route to chaos through a fractal torus in a generalized Lotka–Volterra model. Bull. Math. Biol. 1988, 50, 465–491. [Google Scholar] [CrossRef]
  12. Diethelm, K. An algorithm for the numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal. 1997, 5, 1–6. [Google Scholar]
  13. Khader, M.M.; Babatin, M.M. Numerical treatment for solving fractional SIRC model and influenza A. Comput. Appl. Math. 2014, 33, 543–556. [Google Scholar] [CrossRef]
  14. Khader, M.M.; Adel, M. Chebyshev wavelet procedure for solving FLDEs. Acta Appl. Math. 2018, 158, 1–10. [Google Scholar] [CrossRef]
  15. Horadam, A.F. Vieta Polynomials; The University of New England: Armidale, NSW, Australia, 2000. [Google Scholar]
  16. Zakaria, M.; Khader, M.M.; Al-Dayel, I.; Al-Tayeb, W. Solving fractional generalized Fisher-Kolmogorov-Petrovsky-Piskunov’s equation using compact finite difference method together with spectral collocation algorithms. J. Math. 2022, 2022, 1901131. [Google Scholar]
  17. Adel, M.; Sweilam, N.H.; Khader, M.M. Semi-analytical scheme with its stability analysis for solving the fractional-order predator–prey equations by using Laplace-VIM. J. Appl. Math. Comput. Mech. 2023, 22, 5–17. [Google Scholar] [CrossRef]
  18. El-Hawary, H.M.; Salim, M.S.; Hussien, H.S. Ultraspherical integral method for optimal control problems governed by ordinary differential equations. J. Glob. Optim. 2003, 25, 283–303. [Google Scholar] [CrossRef]
  19. Snyman, J.A. The LFOPC Leap-Frog algorithm for constrained optimization. Comput. Math. Appl. 2000, 40, 1085–1096. [Google Scholar] [CrossRef]
  20. Zhou, Z.F.; Shi, Y. An ODE method for solving constrained optimization. Comput. Math. Appl. 1997, 3, 97–102. [Google Scholar] [CrossRef]
  21. Zongfang, Z. A convergence of ODE method in constrained optimization. J. Math. Anal. Appl. 1998, 218, 297–307. [Google Scholar] [CrossRef]
  22. Parand, K.; Delkhosh, M. Operational matrices to solve nonlinear Volterra-Fredholm integrodifferential equations of multi-arbitrary order. Gazi Univ. J. Sci. 2016, 29, 895–907. [Google Scholar]
Figure 1. Approximate solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) for different values of ν .
Figure 1. Approximate solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) for different values of ν .
Axioms 13 00071 g001
Figure 2. Approximate solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) for different values of N.
Figure 2. Approximate solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) for different values of N.
Axioms 13 00071 g002
Figure 3. REF of the solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c).
Figure 3. REF of the solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c).
Axioms 13 00071 g003
Figure 4. REF of the solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) by VLPs and RK4 ( ν = 1 ).
Figure 4. REF of the solutions Φ 1 ( ζ ) (a), Φ 2 ( ζ ) (b), and Φ 3 ( ζ ) (c) by VLPs and RK4 ( ν = 1 ).
Axioms 13 00071 g004
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Khader, M.M.; Macías-Díaz, J.E.; Román-Loera, A.; Saad, K.M. A Note on a Fractional Extension of the Lotka–Volterra Model Using the Rabotnov Exponential Kernel. Axioms 2024, 13, 71. https://doi.org/10.3390/axioms13010071

AMA Style

Khader MM, Macías-Díaz JE, Román-Loera A, Saad KM. A Note on a Fractional Extension of the Lotka–Volterra Model Using the Rabotnov Exponential Kernel. Axioms. 2024; 13(1):71. https://doi.org/10.3390/axioms13010071

Chicago/Turabian Style

Khader, Mohamed M., Jorge E. Macías-Díaz, Alejandro Román-Loera, and Khaled M. Saad. 2024. "A Note on a Fractional Extension of the Lotka–Volterra Model Using the Rabotnov Exponential Kernel" Axioms 13, no. 1: 71. https://doi.org/10.3390/axioms13010071

APA Style

Khader, M. M., Macías-Díaz, J. E., Román-Loera, A., & Saad, K. M. (2024). A Note on a Fractional Extension of the Lotka–Volterra Model Using the Rabotnov Exponential Kernel. Axioms, 13(1), 71. https://doi.org/10.3390/axioms13010071

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop