Three-species predator-prey model with respect to Caputo and Caputo-Fabrizio fractional operators

TL;DR

This paper introduces a new numerical method for analyzing three-species predator-prey models using Caputo and Caputo-Fabrizio fractional derivatives, highlighting the impact of operator type and order on system stability.

Contribution

It develops an enhanced Adams-Bashforth based numerical scheme for Caputo-Fabrizio operators and compares stability conditions with Caputo derivatives.

Findings

Operator type and order significantly affect system stability.

The new numerical method improves stability analysis.

Dynamics vary with fractional operator choice.

Abstract

We study distributed lag effects in three-dimensional Lotka-Volterra systems by applying the concept of fractional calculus. We derive a new numerical method that provides enhanced stability for the Caputo-Fabrizio operator based on Adams-Bashforth method, considering non-singular kernel in the definition of Caputo-Fabrizio operator. We investigate the stability conditions of this system with comparisons to the Caputo fractional derivative. Numerical results show that the type of differential operators and the value of orders significantly influence the stability of the numerical solution, and dynamics of the Lotka-Volterra system.

Tables

Table 1. Table 1: Stability conditions for Caputo and CF operators.

Equilibrium point	Caputo derivative	CF operator
${\epsilon }_0$	Always saddle	$a_1>\frac{1}{1-\alpha }$
${\epsilon }_1$	and	and or $\frac{a_1{a}_4-a_2{a}_3}{a_2}>\frac{\alpha }{1-\alpha }$   and   $\frac{a_{1}a6-a_2{a}_5}{a_2}>\frac{\alpha }{1-\alpha }$
${\epsilon }_2$		or $\left[1-a_3-\frac{a_4}{a_6}\left(1-a_5\right)\right]>\frac{1}{1-\alpha }$  and $\frac{\left[a_2\left(1-a_5\right)\pm \sqrt{a_2^2{\left(1-a_5\right)}^2+4a_6\left(1-a_5\right)\left(a_1{a}_6+a_2\left(1-a_5\right)\right)}\right]}{2a_6}>\frac{1}{1-\alpha }$
${\epsilon }_3$		or $1-a_4-\frac{a_6}{a_4}\left(1-a_3\right)+\frac{a_7}{a_4}\left[a_1{a}_4+a_2\left(1-a_3\right)\right]>\frac{1}{1-\alpha }$    and $\frac{\left[a_2\left(1-a_3\right)\pm \sqrt{a_2^2{\left(1-a_3\right)}^2+4a_4\left(1-a_3\right)\left[a_1{a}_4+a_2\left(1-a_3\right)\right]}\right]}{2a_4}>\frac{1}{1-\alpha }$
${\epsilon }_4$	$a_6>\frac{a_2{a}_4\left(a_3-1\right)\left[w+a_2\left(a_3-1\right)\right]}{w\left(a_2+a_4\right)+a_2{a}_4\left(a_3-1\right)}$ $\left(w=a_4\left(1+a_1{a}_7-a_5\right)+\left(a_6-a_1{a}_7\right)\left(a_3-1\right)\right)$	$a_6>\frac{a_2{a}_4\left(a_3-1\right)\left[w+a_2\left(a_3-1\right)\right]}{w\left(a_2+a_4\right)+a_2{a}_4\left(a_3-1\right)}$ or ${\lambda }_1,{\lambda }_2,{\lambda }_3>\frac{1}{\left(1-\alpha \right)}$ (see equations (7)-(11))

Table 2. Table 2: The summary of examples; C , CF , and U ​ ( 0 ) 𝑈 0 U(0) denote Caputo, Caputo-Fabrizio, and initial values, respectively, and the notation ✓ indicates the system is asymptotically stable, while ✗ implys unstability.

Example 1				C	CF	C	CF
Coefficient	$U\left(0\right)$	Equilibrium	Eigenvalues	$\alpha =0.98$		$\alpha \le 0.66$
$a_1=3$ $a_2=0.5$ $a_3=4$ $a_4=3$ $a_5=4$ $a_6=9$ $a_7=4$	$x_0=0.5$ $y_0=0.9$ $z_0=0.1$	${\epsilon }_0(0,0,0)$	${\lambda }_0(-3,-3,3)$	✗	✗	✗	✓
		${\epsilon }_1(6,0,0)$	${\lambda }_1(-3,15,51)$	✗	✗	✗	✓
		${\epsilon }_2(0.33,0,2.83)$	${\lambda }_2(-0.083-2.914i,-0.083+2.914i,-2+0i)$	✓	✓	✓	✓
		${\epsilon }_3(1,2.5,0)$	${\lambda }_3(-0.25-2.727i,-0.25+2.727i,16+0i)$	✗	✗	✗	✓
		${\epsilon }_4(1,-1.5,4)$ Not Acceptable	${\lambda }_4(-1.239-5.904i,-1.239+5.904i,1.978+0i)$	✗	✗	✗	✗
Example 2				C		CF
Coefficient	$U\left(0\right)$	Equilibrium	Eigenvalues	$\alpha =0.6$
$a_1=3$ $a_2=0.5$ $a_3=4$ $a_4=3$ $a_5=14$ $a_6=9$ $a_7=4$	$x_0=2$ $y_0=2$ $z_0=3$	${\epsilon }_0(0,0,0)$	${\lambda }_0(-13,-3,3)$	✗		✓
		${\epsilon }_1(6,0,0)$	${\lambda }_1(-3,15,41)$	✗		✓
		${\epsilon }_2(1.44,0,2.28)$	${\lambda }_2(-0.361-5.429i,-0.361+5.429i,1.333+0i)$	✗		✗
		${\epsilon }_3(1,2.5,0)$	${\lambda }_3(-0.25-2.727i,-0.25+2.727i,6+0i)$	✗		✓
		${\epsilon }_4(1,1,1.5)$	${\lambda }_4(0.276-4.123i,0.276+4.123i,-1.053+0i)$	✓		✗
Example 3				C		CF
Coefficient	$U\left(0\right)$	Equilibrium	Eigenvalues	$\alpha =0.4$
$a_1=8$ $a_2=0.5$ $a_2=4$ $a_4=1$ $a_5=7$ $a_6=9$ $a_7=4$	$x_0=0.5$ $y_0=0.1$ $z_0=5$	${\epsilon }_0(0,0,0)$	${\lambda }_0(-6,-3,8)$	✗		✓
		${\epsilon }_1(160,0,0)$	${\lambda }_1(-8,157,1434)$	✗		✓
		${\epsilon }_2(0.666,0,7.966)$	${\lambda }_2(-0.016-6.913i,-0.016+6.913i,-2.333+0i)$	✓		✓
		${\epsilon }_3(3,7.85,0)$	${\lambda }_3(-0.075-4.852i,-0.075+4.852i,52.4+0i)$	✗		✓
		${\epsilon }_4(3,-5.25,13.1)$ Not Acceptable	${\lambda }_4(-1.274-18.50i,-1.274+18.50i,2.398+0i)$	✗		✓

Equations

$${}^{\text{C}}{D}_{a^{+}}^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_a^t(t-\tau {)}^{-\alpha }{f}^{\prime }(\tau )d\tau ,0<\alpha <1,t>a.$$

$${}^{\text{CF}}{D}_{a^{+}}^{\alpha }f(t)=\frac{M(\alpha )}{1-\alpha }{\int }_a^t\exp (-\frac{\alpha }{1-\alpha }(t-\tau ))f^{\prime }(\tau )d\tau ,t\ge a,$$

$$D_{a^{+}}^{\alpha }x(t)=Ax(t),$$

$$\{\begin{array}{l}{D}_{a^{+}}^{\alpha }x(t)=x(t)(a_1-a_{2}x(t)-y(t)-z(t))\\ D_{a^{+}}^{\alpha }y(t)=y(t)((1-a_3)+a_{4}x(t))\\ D_{a^{+}}^{\alpha }z(t)=z(t)((1-a_5)+a_{6}x(t)+a_{7}y(t))\end{array}$$

$$x(0)=x_0,y(0)=y_0,z(0)=z_0,$$

$$\{\begin{array}{l}{D}_{a^{+}}^{\alpha }\text{u}(t)=\text{F}(\text{u}(t))0<t<\infty \\ \text{u}(0)={\text{u}}_0,\end{array}$$

$$D_{a^{+}}^{\alpha }\text{u}(t)=\text{A}\text{u}(t)+x(t)\text{B}\text{u}(t)+y(t)\text{C}\text{u}(t)+z(t)\text{D}\text{u}(t)$$

$$\text{A}=\left[\begin{array}{ccc}{a}_1& 0& 0\\ 0& 1-a_3& 0\\ 0& 0& 1-a_5\end{array}\right],\text{B}=\left[\begin{array}{ccc}-a_2& 0& 0\\ 0& a_4& 0\\ 0& 0& a_6\end{array}\right],\text{C}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 0& 0\\ 0& 0& a_7\end{array}\right],\text{D}=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right].$$

$$\Vert \text{F}(\text{u}(t))-\text{F}(\text{v}(t))\Vert$$

$$\displaystyle=\bigg{\|}\textbf{A}\textbf{u}(t)+u_{1}(t)\textbf{B}\textbf{u}(t)+u_{2}(t)\textbf{C}\textbf{u}(t)+u_{3}(t)\textbf{D}\textbf{u}(t)-(\textbf{A}\textbf{v}(t)+v_{1}(t)\textbf{B}\textbf{v}(t)+v_{2}(t)\textbf{C}\textbf{v}(t)+v_{3}(t)\textbf{D}\textbf{v}(t))\bigg{\|}$$

$$⩽\Vert \text{A}(\text{u}(t)-\text{v}(t))\Vert +\Vert u_1(t)\text{B}(\text{u}(t)-\text{v}(t))\Vert +\Vert (u_1(t)-v_1(t))\text{B}\text{v}(t)\Vert +\Vert u_2(t)\text{C}(\text{u}(t)-\text{v}(t))\Vert$$

$$+\Vert (u_2(t)-v_2(t))\text{C}\text{v}(t)\Vert +\Vert u_3(t)\text{D}(\text{u}(t)-\text{v}(t))\Vert +\Vert (u_3(t)-v_3(t))\text{D}\text{v}(t)\Vert$$

$$\displaystyle\leqslant\bigg{[}\|\textbf{A}\|+\|\textbf{B}\|\left(|u_{1}(t)|+\|\textbf{v}(t)\|\right)+\|\textbf{C}\|\left(\left|u_{2}(t)\right|+\left\|\textbf{v}(t)\right\|\right)+\|\textbf{D}\|\left(\left|u_{3}(t)\right|+\left\|\textbf{v}(t)\right\|\right)\bigg{]}\times\left\|\textbf{u}(t)-\textbf{v}(t)\right\|,$$

$$\Vert \text{F}(\text{u}(t))-\text{F}(\text{v}(t))\Vert ⩽\text{L}\Vert \text{u}(t)-\text{v}(t)\Vert ,$$

$$\text{L}=\Vert \text{A}\Vert +(\Vert \text{B}\Vert +\Vert \text{C}\Vert +\Vert \text{D}\Vert )(M_1+M_2)>0,$$

$$D_{a^{+}}^{\alpha }x(t)=0,D_{a^{+}}^{\alpha }y(t)=0,D_{a^{+}}^{\alpha }z(t)=0.$$

$${\epsilon }_0=(0,0,0),$$

$${\epsilon }_1=(\frac{a_1}{a_2},0,0),$$

$${\epsilon }_2=\left(\frac{a_5-1}{a_6},0,\frac{a_1{a}_6-a_2(a_5-1)}{a_6}\right),$$

$${\epsilon }_3=\left(\frac{a_3-1}{a_4},\frac{a_1{a}_4-a_2(a_3-1)}{a_4},0\right),$$

$${\epsilon }_4=\left(\frac{a_3-1}{a_4},\frac{a_4(a_5-1)-a_6(a_3-1)}{a_7{a}_4},\frac{a_4(1+a_1{a}_7-a_5)+(a_6-a_2{a}_7)(a_3-1)}{a_7{a}_4}\right).$$

$$J(x^{\ast },y^{\ast },z^{\ast })=\left[\begin{array}{ccc}{a}_1-2a_2{x}^{\ast }-y^{\ast }-z^{\ast }& -x^{\ast }& -x^{\ast }\\ a_4{y}^{\ast }& 1-a_3+a_4{x}^{\ast }& 0\\ a_6{z}^{\ast }& a_7{z}^{\ast }& a_6{x}^{\ast }-a_5+a_7{y}^{\ast }+1\end{array}\right].$$

$$J({\epsilon }_0)=\left[\begin{array}{ccc}{a}_1& 0& 0\\ 0& 1-a_3& 0\\ 0& 0& 1-a_5\end{array}\right],$$

$$J({\epsilon }_1)=\left[\begin{array}{ccc}-a_1& -\frac{a_1}{a_2}& -\frac{a_1}{a_2}\\ 0& 1-a_3+\frac{a_1{a}_4}{a_2}& 0\\ 0& 0& 1-a_5+\frac{a_1{a}_6}{a_2}\end{array}\right],$$

$$a_1{a}_4<a_2{a}_3-a_2,$$

$$a_1{a}_6<a_2{a}_5-a_2,$$

$$a_2(1-a_3)(1-\alpha )>a_2-a_1{a}_4(1-\alpha ),$$

$$a_2(1-a_5)(1-\alpha )>a_2-a_1{a}_6(1-\alpha ).$$

$$J({\epsilon }_2)=\left[\begin{array}{ccc}\frac{a_2}{a_6}(1-a_5)& \frac{1-a_5}{a_6}& \frac{1-a_5}{a_6}\\ 0& 1-a_3-\frac{a_4}{a_6}(1-a_5)& 0\\ a_1{a}_6+a_2(1-a_5)& \frac{a_7}{a_6}\left(a_1{a}_6+a_2(1-a_5)\right)& 0\end{array}\right],$$

$$J({\epsilon }_2)=\left[\begin{array}{ccc}A& B& C\\ 0& D& 0\\ E& F& 0\end{array}\right].$$

$$(\lambda -D)({\lambda }^2-A\lambda -CE)=0,$$

$$A<0,C<0,E>0,D<0,$$

$${\lambda }_1=D<0,{\lambda }_2+{\lambda }_3=A<0,{\lambda }_2{\lambda }_3=-CE>0.$$

$${\lambda }_1=\left[1-a_3-\frac{a_4}{a_6}(1-a_5)\right],$$

$$\lambda_{2,3}=\frac{\bigg{[}a_{2}(1-a_{5})\pm\sqrt{a_{2}^{2}(1-a_{5})^{2}+4a_{6}(1-a_{5})(a_{1}a_{6}+a_{2}(1-a_{5}))}\bigg{]}}{2a_{6}}.$$

Full text

1 Introduction

The use of fractional calculus has rapidly increased in many fields of science and engineering [1, 2]. Fractional calculus brings more degrees of freedom for differentiation in the modeling of various phenomena, such as complex networks [3, 4], optimal control problems [5, 6], and Viscoelastic systems [7]. Such flexible differential operators do not have unique definitions, however. The Grunwald-Letnikov, Riemann-Liouville, and Caputo definitions are examples of commonly used approaches that have been used in science and engineering.

Here, we study the Caputo-Fabrizio (CF) fractional operator that at first comes from the definition of Caputo fractional derivative [8]. It replaces the singular kernel of the Caputo derivative with an exponential function and has two representations for the temporal and spatial variables [9]. This approach has been used, for instance, to model the behavior of the diffusion-convection equation, fractional Nagumo equation, and to control the wave on a shallow water [10, 11]. The new operator has also been successfully applied in cancer treatment, HIV/AIDS infection, and tumor-obesity model [12, 13, 14]. In [15], a comprehensive overview of the CF operator has been conducted, showing that this operator is applied in economic and physical models.

Another significant application is provided by the classical Lotka-Volterra systems, which are sometimes called predator-prey or parasite-host equations. Such systems play a remarkable role in mathematical biology [16], and in financial systems, for example, biunivoc capital transfer from mother bank to subsiding bank and from subsiding bank to individuals or companies [17]. At first, these models were introduced independently by Alfred J. Lotka and Vito Volterra as a simplified model of two species predator-prey population dynamics [18]. Whereas the classical formulation uses an integer-order differential, this is not always optimal due to the nonlocality of the interactions and the potential existence of memory or lag effects in the real systems. In 2007, Ahmed et al. [19] have introduced the fractional-order Lotka-Volterra system. Recent studies have generalized Lotka-Volterra models to two-predator one-prey dynamics [20] and analysed a Lotka-Volterra fractional-order model using the Caputo fractional derivative [21, 22].

However, due to the appearance of a singularity in the definition of the Caputo fractional derivative, this operator is impractical for the modeling of some nonlocal dynamics [23]. Hence, the new nonsingular CF fractional operator has been proposed to overcome this shortcoming. Recently, scholars have examined the efficiency of this fractional framework through some real practical cases and provided more accurate parameter fitting than the classic integer and noninteger order models [24]. Importantly, Tarasov [15] has shown that a system with the CF operator, in contrast to the Caputo fractional derivatives, cannot describe processes with memory effects but can suitably model processes with continuously distributed lag.

Motivated by the above discussion, we consider three-dimensional Lotka-Volterra differential equations described by the CF operator. We formulate a new corrected numerical method to solve the CF system. We analyse the stability properties of the Lotka-Volterra model under Caputo and CF fractional operators. Consequently, we reveal new insights from differences in the stability region of these frameworks. To investigate the impact of different stability properties on the dynamic, we compare and illustrate some example cases.

2 Definitions

2.1 Fractional calculus

In this section, we outline the key definitions for fractional derivatives. The fractional derivative in the sense of Caputo is defined as follow [25]:

[TABLE]

The kernel $(t-\tau)^{-\alpha}$ in Eq. (1) cause a singularity at $t=\tau$ that can be considered as a drawback in this definition. In 2015, Caputo and Fabrizio defined the following fractional derivative as [9]

[TABLE]

where $f$ is a continuous and differentiable function on $C^{1}[a,b]$ and $M(\alpha)$ is a normalization function such that $M(0)=M(1)=1$.

For more details on the above-mentioned fractional operators, the readers are referred to [8, 9].

2.2 Stability of the fractional-order system

In this part of the paper, we recall some basic theorems to proceed with our goal on the stability of the Lotka-Volterra system with the CF operator.

Definition 2.1**.**

The autonomous system, with $x(t_{0})=x_{0}$ is asymptotically stable if and only if $\displaystyle\lim_{t\to+\infty}\|x(t)\|=0,$ where $\|.\|$ is Euclidean norm.

Consider the linear fractional-order autonomous system as follows

[TABLE]

where $x(t)\in\mathbb{R}^{n}$, $A\in\mathbb{R}^{n\times n}$, $0<\alpha<1$ and $D^{\alpha}_{a^{+}}$ is one of Caputo or CF operators.

Theorem 2.1**.**

The linear autonomous system (3) with Caputo fractional derivative for $0<\alpha<1$ is asymptotically stable if and only if $\left|\mathrm{arg}(\mathrm{spec}(A))\right|>\dfrac{\alpha\pi}{2}$, $\mathrm{spec}(A)$ is the spectrum (set of all eigenvalues) of $A$ [26].

Theorem 2.2**.**

The system (3) with CF operator is asymptotically stable if eigenvalues $\lambda(A)$ of the matrix A satisfy one of the following conditions [27]

$\|(\lambda(A))\|\geqslant\dfrac{1}{1-\alpha},\lambda(A)\neq\dfrac{1}{1-\alpha},$**

• 2)

$\mathrm{Re}(\lambda(A))>\dfrac{1}{1-\alpha},$**

• 3)

$\mathrm{Re}(\lambda(A))<0,$**

• 4)

$\left|\mathrm{Im}(\lambda(A))\right|>\dfrac{1}{2(1-\alpha)}.$**

3 Modeling and its stability

One can find the stability conditions of this system with Caputo derivative in [28]. In this section, we investigate the stability of the Lotka-Volterra system involving the CF operator and Table 1 summarizes all the required conditions for the stability in respect to the both operators.

3.1 Fractional Lotka-Volterra model

In this study, we consider a three-species Lotka-Volterra model as follows

[TABLE]

where $0<\alpha\leq 1$ and $a_{i}>0,i=0,1,\ldots,7$ and ${}_{0}D^{\alpha}_{t}$ is one of the differential operators Caputo or CF with initial conditions

[TABLE]

where $x_{0},y_{0},z_{0}\in\mathbb{R}^{+}.$ In this model $x(t)\geq 0$ represents the population of the prey, $y(t)\geq 0$ and $z(t)\geq 0$ represent the population of predators at the time $t$.

Now we consider system (4) in a compact form as follows

[TABLE]

where $\textbf{u}(t)={(x(t),y(t),z(t))^{T}}\in\mathscr{L}[0,t^{\prime}]$, where $\mathscr{L}[0,t^{\prime}]$ be the set of all continuous vector $\textbf{u}(t)$ defined on the interval $[0,t^{\prime}]$ ($t^{\prime}>0$) and F is a real-valued continuous vector function. Then system (4) can be written in the form

[TABLE]

where

[TABLE]

Theorem 3.3**.**

For $\textbf{u}(t)\in\mathscr{L}[0,t^{\prime}]$, system (6) has a unique solution.

Proof.

Let $\textbf{F}(\textbf{u}(t))=\textbf{A}\textbf{u}(t)+u_{1}(t)\textbf{B}\textbf{u}(t)+u_{2}(t)\textbf{C}\textbf{u}(t)+u_{3}(t)\textbf{D}\textbf{u}(t)$ and $\textbf{F}(\textbf{v}(t))=\textbf{A}\textbf{v}(t)+v_{1}(t)\textbf{B}\textbf{v}(t)+v_{2}(t)\textbf{C}\textbf{v}(t)+v_{3}(t)\textbf{D}\textbf{v}(t)$ where $\textbf{F}(\textbf{u}(t))$, $\textbf{F}(\textbf{v}(t))\in\mathscr{L}[0,t^{\prime}]$. Since $\textbf{u}(t)=(u_{1},u_{2},u_{3})$, $\textbf{v}(t)=(v_{1},v_{2},v_{3})\in\mathscr{L}[0,t^{\prime}]$ such that $\textbf{u}(t)\neq\textbf{v}(t)$. The following inequality holds

[TABLE]

then we have

[TABLE]

where

[TABLE]

and $M_{1}$ and $M_{2}$ are positive constant and satisfy $\|\textbf{u}\|\leq M_{1}$, $\|\textbf{v}\|\leq M_{2}$ as a aresult of $\textbf{u},\textbf{v}\in\mathscr{L}[0,t^{\prime}].$ It means that $\textbf{F}(X(t))$ is continuous and satisfying Lipschitz condition, then the initial value problem (6) has a unique solution. ∎

3.2 Stability of the model

In this subsection, we discuss the stability of non-linear Lotka-Volterra differential equations (4) described by the CF operator. In the case of non-linear systems, we study the local stability of equilibrium points and the following theorems are presented to investigate the local stability of equilibrium points. In order to determine the equilibrium points of system (4), let us consider

[TABLE]

The equilibrium points of system (4) are obtained and denoted as

[TABLE]

To adjust the conditions for the actual situations, the equilibrium points must be nonnegative. In this regard, it is obvious that $\varepsilon_{0}$ and $\varepsilon_{1}$ always exist, and $\varepsilon_{2}$ exists when $a_{3}\geq 1$ and $a_{1}a_{4}\geq a_{2}(a_{3}-1)$, and it happens for $\varepsilon_{3}$ when $a_{5}\geq 1$ and $a_{1}a_{6}\geq a_{2}(a_{5}-1)$. Finally, the conditions $a_{3}\geq 1$, $a_{4}(a_{5}-1)\geq a_{6}(a_{3}-1)$ and $a_{4}\geq\dfrac{(a_{2}a_{7}-a_{6})(a_{3}-1)}{(1+a_{1}a_{7}-a_{5})}$ (or if $(1+a_{1}a_{7}-a_{5})<0$ then $a_{4}\leq\dfrac{(a_{6}-a_{2}a_{7})(a_{3}-1)}{(1+a_{1}a_{7}-a_{5})}$, else if $(1+a_{1}a_{7}-a_{5}=0$ then $a_{6}>a_{2}a_{7}$) are necessary for the existence of $\varepsilon_{4}$.

Theorem 3.4**.**

Let $\varepsilon^{*}$ be an equilibrium point of the nonlinear system (4), with CF operator, then equilibrium point $\varepsilon^{*}$ is asymptotically stable if eigenvalues of Jacobian matrix $\lambda(J(\varepsilon^{*}))$ satisfy one of the following conditions

$\|\lambda(J(\varepsilon^{*}))\|\geqslant\dfrac{1}{1-\alpha},\lambda(J(\varepsilon^{*}))\neq\dfrac{1}{1-\alpha},$**

• 2)

$\mathrm{Re}(\lambda(J(\varepsilon^{*})))>\dfrac{1}{1-\alpha},$**

• 3)

$\mathrm{Re}(\lambda(J(\varepsilon^{*})))<0,$**

• 4)

$\left|\mathrm{Im}(\lambda(J(\varepsilon^{*})))\right|>\dfrac{1}{2(1-\alpha)}.$**

Proof.

The proof is straightforward with Theorem 2.2 and [26]. ∎

To study the local stability of the equilibrium points such as $(x^{\ast},y^{\ast},z^{\ast})$ for system (4) we provide the Jacobian matrix $J(x^{\ast},y^{\ast},z^{\ast})$ as follows

[TABLE]

3.2.1 The first equilibrium

For $\varepsilon_{0}$, the Jacobian can be expressed as

[TABLE]

where eigenvalues are $\lambda_{1}=a_{1}$, $\lambda_{2}=1-a_{3}$, $\lambda_{3}=1-a_{5}$. Since $a_{1}>0$, $1-a_{3}<0$, $1-a_{4}<0$, then $\varepsilon_{0}$ is stable if $a_{1}>\dfrac{1}{1-\alpha}$.

3.2.2 The second equilibrium

For $\varepsilon_{1}$, the Jacobian matrix is

[TABLE]

where $\lambda_{1}=-a_{1}<0$, $\lambda_{2}=1-a_{3}+\dfrac{a_{1}a_{4}}{a_{2}}$, and $\lambda_{3}=1-a_{5}+\dfrac{a_{1}a_{6}}{a_{2}}$. Thus, $\varepsilon_{1}$ is asymptotically stable when

[TABLE]

or

[TABLE]

3.2.3 The third equilibrium

For $\varepsilon_{2}$ the Jacobian matrix is

[TABLE]

we use the below notation

[TABLE]

The characteristic equation is as follows

[TABLE]

by the condition $\dfrac{a_{5}-1}{a_{6}}<\dfrac{a_{1}}{a_{2}}<\dfrac{a_{3}-1}{a_{4}}$ we get

[TABLE]

therefore

[TABLE]

The eigenvalues are

[TABLE]

and

[TABLE]

In this case, we can conclude that $\varepsilon_{2}$ is locally asymptoialy stable. However, when the condition $\dfrac{a_{5}-1}{a_{6}}<\dfrac{a_{1}}{a_{2}}<\dfrac{a_{3}-1}{a_{4}}$ was not available, $\varepsilon_{2}$ could be locally asymtotically stable when $\lambda_{1},\lambda_{2},\lambda_{3}>\dfrac{1}{1-\alpha}$, which leads

[TABLE]

and

[TABLE]

3.2.4 The fourth equilibrium

Jacobian of $\varepsilon_{3}$ is

[TABLE]

where $w=1-a_{5}-\dfrac{a_{6}}{a_{4}}(1-a_{3})+\dfrac{a_{7}}{a_{4}}[a_{1}a_{4}+a_{2}(1-a_{3})]$. Same as $\varepsilon_{2}$, we can provid the stability condition as

[TABLE]

where $\lambda_{1}=1-a_{4}-\dfrac{a_{6}}{a_{4}}(1-a_{3})+\dfrac{a_{7}}{a_{4}}[a_{1}a_{4}+a_{2}(1-a_{3})]$ and

[TABLE]

When the condition $\dfrac{a_{3}-1}{a_{4}}<\dfrac{a_{1}}{a_{2}}<\dfrac{a_{5}-1}{a_{6}}$ is not available, $\varepsilon_{3}$ is locally asymptotically stable when $\lambda_{1}(1-\alpha)>1$, $\lambda_{2}(1-\alpha)>1$, $\lambda_{3}(1-\alpha)>1$.

3.2.5 The fifth equilibrium

For $\varepsilon_{4}$, Jacobian matrix is as follows

[TABLE]

where

[TABLE]

To compute eigenvalues of the above matrix, we consider the characteristic polynomial, $L(\lambda)=\lambda^{3}+a\lambda^{2}+b\lambda+c,$ where $a=-A,b=-B(C+D),c=-BCE.$ It is obvious $a,c>0.$ If $a_{1}a_{2}>a_{3}$ then Routh-Hurwitz criterion shows that the all roots of are negative. The equation $a_{1}a_{2}-a_{3}=B[A(C+D)+CE]$ is positive if

[TABLE]

where

[TABLE]

To investigate the stability of the system in the sense of CF operator we consider the following parameter for $L(\lambda)$ as

[TABLE]

If $\Delta>0$, then we have only one real solution

[TABLE]

If $\Delta=0$, there are repeated roots

[TABLE]

If $\Delta>0$ then roots are same as below

[TABLE]

Consequently, $\varepsilon_{4}$ is locally asymptotically stable if in any case all of the eigenvalues satisfying these conditions

[TABLE]

We end this section by summarizing the stability conditions of all the equilibrium points for Caputo and CF operators in Table 1.

4 Numerical algorithm

The Predictor-Corrector methods are well-known numerical approach so that their extensions can provide accurate numerical solutions of fractional differential equations (FDEs). For instance, in [29], a numerical method based on Adams-Bashforth methods is proposed for solving FDEs with Caputo derivative. In [30], authors have investigated a fractional Adams-Bashforth method for solving FDEs with the CF operator, although their arguments are flawed. For this aim, in this part of the paper, we correct this method to solve the Lotka-Volterra system of Caputo and the CF operator and we compare the solutions using the both operators. Consider the following differential equation:

[TABLE]

which is equivalent to the following equation

[TABLE]

where $T_{n-1}(x)$ is the Taylor expansion of $f(x)$ centered at $x_{0}=0$ and $T_{n-1}(x)=\sum_{i=0}^{n-1}\frac{x^{i}}{i!}f_{0}^{(i)}$.

The corrector formula $f_{k+1}$ can be written as follows

[TABLE]

and by the fractional Adams-Bashforth-multon method [29], $f_{k+1}^{p}$ is determined by

[TABLE]

where

[TABLE]

Now, consider the following fractional-order system involving CF operator

[TABLE]

We consider $0\leqslant\alpha\leqslant 1$ for simplicity and assume that $(x_{0},y_{0},z_{0})$ is the initial point. Applying the above scheme, system (15) can be discretized as follows

[TABLE]

where

[TABLE]

and

[TABLE]

In the following, we apply the proposed numerical technique for simulation the solutions of system 4. It is worth mentioning that to obtain the numerical results in the sense of Caputo derivative we use the equivalent Adams-Bashforth-multon method described in [29].

5 Numerical implementation

In this part, we discuss the numerical results of the fractional-order Lotka-Volterra model (4) with Caputo and CF operators, by using the numerical method described in Sec. 4. It is helpful to classify the numerical results according to the positions of the eigenvalues and discuss the behavior of the system in the sense of Caputo and CF operator. To illustrate such a classification, we provide stability and instability region for both operators in Fig. 1 and set four eigenvalues on the plane for different cases. The unstable domain of the system with Caputo derivative is an unbound region limited by two lines with angles $-\alpha\frac{\pi}{2}$ and $\alpha\frac{\pi}{2}$. On the other hand, the unstable domain of the system with CF operator is a bounded closed circle centered at $(0,\frac{1}{2(1-\alpha)})$ with the radius $\frac{1}{2(1-\alpha)}$. For both cases, it is clear that the stability of the system has an inverse relation with the order derivative; in fact, the smaller $\alpha$ is, the more there is space for stablitity and vice versa. As it is shown in Fig 1, the eigenvalues can be located in four distinct classes: $\lambda_{A}$ is in an area where both systems are stable; the class of $\lambda_{B}$ is where the system with Caputo derivatives is stable, but the system with CF operator is not stable; $\lambda_{C}$ denotes a class of eigenvalues staying at where both systems are unstable; and finally, $\lambda_{D}$ is where the system in the sense of CF operator is stable but the system with Caputo derivatives is not stable.

We collect a summary of the three examples in Table 2 to easily compare the behavior of the model concerning the Caputo and CF operators.

5.1 Example 1

As one can see in Table 2, the parameters of this example give five distinct equilibrium points (and corresponding eigenvalues), while equilibrium $\epsilon_{4}$ is not acceptable since it has a negative value. Thus, we should expect negative-value solutions of the system when we do not impose any constraints on the components. There is a recommended paper [31] to avoid going toward such meaningless solutions and getting a feasible solution.

Furthermore, this example shows the stability of the equilibrium points depends on the value of the fractional-order $\alpha$. As we expect from Fig. 1, the number of stable equilibrium points increases when we reduce the value of $\alpha$. In this case, when $\alpha$ is $0.98$ for both operators, the system is asymptotically stable only at $\epsilon_{2}$ (see Table 2). Indeed, Fig. 2 (left) shows that the system gets steady at $\epsilon_{2}$, with different oscillations which are related to the definition of the operators. Nonetheless, the condition $\alpha\leq 0.66$ provides a larger area for the stability of the system so that three eigenvalues $\lambda_{0}$, $\lambda_{1}$, and $\lambda_{3}$ stay in the class of $\lambda_{D}$ (see Fig. 1 and Table 2). Hence, with appropriate initial values and differential orders, the system could converge to $\epsilon_{3}$ (see Fig. 2, right).

5.2 Example 2

This example confirms the points mentioned in the previous one; by setting $\alpha=0.6$, the equilibrium $\epsilon_{4}$ is the only stable equilibrium point in the sense of Caputo derivative, while the equilibrium points $\epsilon_{0}$, $\epsilon_{1}$, and $\epsilon_{3}$ are stable concerning the CF operator. But, as shown in Fig. 1, it is interesting that we have here an eigenvalue, $\lambda_{4}$, in the class of $\lambda_{B}$ alongside the $\lambda_{D}$, where $\lambda_{0}$, $\lambda_{1}$, and $\lambda_{3}$ are. As a result, Fig. 3 shows that the system can start from a point to converge asymptotically to the only equilibrium that is stable in the sense of Caputo, rather than the CF operator. Therefore, it could make a challenge for one who assumes a system having a more stability region may lead to more potential to achieve a steady-state.

5.3 Example 3

This example can complete the discussion and make clear the substantial role of the initial values on the behavior of the very Lotka-Volterra model. Considering the information in Table 2, determining the location of four eigenvalues in the $\lambda_{D}$ class (Fig. 1), we expect that the system is most stable for the CF operator and noticeably unstable for the Caputo derivative. Although Fig. 4 illustrates this expectation, it is not an absolute scenario when the process is supposed to start from a point leading to equilibrium $\epsilon_{3}$. It depends on the domain of attraction that the initial values stay and the corresponding eigenvalue, which is in the class of $\lambda_{B}$ (see Fig. 5). For more information on finding the domain of attractions to specific equilibrium points, see [31].

6 Conclusion

We have investigated the three-dimensional Lotka-Volterra system for the CF operator. Concerning the existence of a non-singular kernel in the definition of CF operator, we investigated the stability of the system and suggested a new numerical method with improved stability properties based on Adams-Bashforth methods. This numerical scheme improves the efficiency of the simulations for both Caputo and CF operators. Numerical results demonstrate how the behavior of the Lotka-Volterra system can depend on the type of differential operator and the value of fractional order. Moreover, we have shown that the CF operator provides different properties compared to the classical Caputo derivative. Overall, this analysis can enhance our understanding of the exceptional dynamics in complex systems.

Analysing the behavior of Lotka-Volterra models under incommensurate fractional orders, where the different interacting partners may have different degrees of memory or lag effects, is a fascinating line for further research. In real-world complex systems, a time-variable dependency on the past states is likely and may lead to anomalous behaviors that pose challenges for modeling. Besides, fractional calculus provides tools to simulate systems with such incommensurate fractional-order derivatives. Therefore, finding the stability region of the three-dimensional Lotka-Volterra model with incommensurate fractional orders in the sense of Caputo and CF operators as well as examining the domain of attractions are promising directions for future studies.