Abstract

We analyze the dynamics of a fractional order modified Leslie-Gower model with Beddington-DeAngelis functional response and additive Allee effect by means of local stability. In this respect, all possible equilibria and their existence conditions are determined and their stability properties are established. We also construct nonstandard numerical schemes based on Grünwald-Letnikov approximation. The constructed scheme is explicit and maintains the positivity of solutions. Using this scheme, we perform some numerical simulations to illustrate the dynamical behavior of the model. It is noticed that the nonstandard Grünwald-Letnikov scheme preserves the dynamical properties of the continuous model, while the classical scheme may fail to maintain those dynamical properties.

1. Introduction

The dynamical interaction of predator and prey is one of important subjects in ecological science. In recent years, one of the most important species interactions is predator-prey model [1]. One of well-known mathematical models which describe the dynamics of prey-predator interaction is the modified Leslie-Gower model proposed by Aziz-Alaoui and Okiye [2]. In this model, the growth rate of predator is in the form of logistics-type where its carrying capacity is proportional to the prey number and environment protection for predator. One of important parameters describing the prey-predator interaction is the functional response which describes the predator’s rate of prey consumption per capita. Aziz-Alaoui and Okiye [2] and Yu [3] have considered a modified Leslie-Gower model with Holling type II functional response, while Yu [4] considered the same model but with Beddington-DeAngelis functional response. In the normalized variables, the modified Leslie-Gower equation with Beddington-DeAngelis functional response can be written as where and denote population densities of prey and predator at time , respectively. The parameters , , , , , and are positive constants.

One of other factors that influence the interaction of predator and prey is Allee effect, referring to a decrease in per capita fertility rate at low population densities. Allee effect may occur under several mechanisms, such as difficulties in finding mates when population density is low or social dysfunction at small population sizes. When such a mechanism operates, the per capita fertility rate of the species increases with density; that is, positive interaction among species occurs [58]. Recently Indrajaya et al. [9] investigate a modified Leslie-Gower equation with Beddington-DeAngelis functional response and additive (both weak and strong) Allee effect on preywith initial conditions and . The criteria of the Allee effect are as follows [7, 8]:(i)If , then the Allee effect is weak.(ii)If , then the Allee effect is strong.

It is shown in system (2) that the growth rates of both prey and predator depend only on the current state. In many situations, the growth rate is also dependent on the history of variable or its memory. With the rapid development of fractional calculus, fractional differential equations have been implemented in various fields including biological system. This is due the fact that fractional differential equations are naturally related to the real life phenomena with memory which exists in most of biological system [1016]. To describe such memory effect, we first recall the definition of fractional integral operator as well as fractional differential operator.

Definition 1 (see [17]). The Riemann-Liouville -order fractional integral operator of any function , is defined bywhere is the Euler Gamma function.

Definition 2 (see [17]). Let be an integer which satisfies . The Riemann-Liouville -order fractional derivative of function is defined aswhere is the common -order derivative.

The Riemann-Liouville fractional derivative is historically the first concept of fractional derivative and theoretically well established. However, in the case of Riemann-Liouville fractional differential equation, the initial value is usually given in the form of fractional derivative, which is not practical. Consequently, one applies the Caputo fractional derivative which is defined as follows.

Definition 3 (see [17]). The Caputo fractional differential operator of order , with , , is defined byFor simplicity, the Caputo fractional derivative of function of order is denoted by .

From Definition 3, we see the -order fractional derivative at time t is not defined locally; it relies on the total effects of the commonly used m-order integer derivative on the interval . So it can be used to describe the variation of a system in which the instantaneous change rate depends on the past state, which is called the ‘‘memory effect’’ in a visualized manner [18].

In this paper we reconsider system (2). By assuming that the growth rates of both prey and predator at time do not only depend instantaneously on the current state but also depend on the past state, we replace the first order derivatives in system (2) with the fractional order Caputo type derivatives:where . Hence we have a system of fractional differential equation. In the following we discuss the dynamical properties of system (6). To study the stability of equilibrium points, we apply the following stability theorem.

Theorem 4 (see [19]). Consider the following autonomous nonlinear fractional order system:The equilibrium points of the above system are solutions to the equation . An equilibrium point is locally asymptotically stable if all eigenvalues of the Jacobian matrix evaluated at equilibrium satisfy .

2. Equilibria and Their Stability

Based on Theorem 4, equilibria of model (6) can be determined by solving the following system:System (8) has been solved by Indrajaya et al. [9], and the obtained equilibria are as follows:(1)The equilibria of system (6) for the weak Allee effect case are(a) trivial equilibrium , that is, the extinction of both prey and predator point,(b) two axial equilibria, that is, the prey extinction point and the predator extinction point where (c) positive or coexistence equilibrium , where are all possible real positive solutions of cubic equationwhere(2)The equilibria of system (6) for the strong Allee effect case are(a) trivial equilibrium , that is, the extinction of both prey and predator point,(b) three axial equilibria that are the prey extinction point and two predator extinction points: and , where (c) positive or coexistence equilibrium point , where are also all possible real positive solutions of cubic equation (10).Using transformation , (10) can be reduced towhere and . Implementing Cardan’s method as performed by Cai et al. [7], we obtain the existence of the positive equilibria as follows.

Lemma 5 (existence of positive equilibria). Let be the interior equilibrium of model (6) for both weak and strong Allee effects where is a real positive solution of (10). Then the following statements hold:(a) If , then (13) has a single positive root . As a result, model (6) has a unique positive equilibrium point, that is, , with .(b) Suppose that and , then(b1)If , then (13) has a positive root of multiplicity two. Thus, model (6) has a unique positive equilibrium point, that is, .(b2)If , then (13) has two positive roots and . Thus, model (6) has two positive equilibrium points, namely, and , with .(c) If and , then (13) has a unique positive root . As a result, model (6) has a unique positive equilibrium point, that is, with

Moreover, algebraic computations show that if (13) has two positive roots, then they areIf (13) has a positive root, then it must be

To check the local stability of each equilibrium point, we linearize system (6) around the equilibrium and verify all eigenvalues of the Jacobian matrix evaluated at the equilibrium. The stability properties of trivial and axial equilibrium points for the case of weak and strong Allee effect are, respectively, stated in Theorems 6 and 7.

Theorem 6. Stability of trivial and axial equilibrium for weak Allee effect :(i)the trivial equilibrium and the axial equilibrium are always unstable;(ii)the axial equilibrium is asymptotically stable if .

Proof. (i) The Jacobian matrix at is , and the eigenvalues are and . It is clear that . Hence is unstable. The Jacobian matrix at is where one of its eigenvalues is and therefore is unstable because , .
(ii) The Jacobian matrix at is . The eigenvalues of are and . Thus and whenever . This proves part (ii).

Using the same argument as in the proof of Theorem 6, we obtain the following stability properties of equilibria for the case of strong Allee effect.

Theorem 7. Stability of trivial and axial equilibrium for strong Allee effect : (i)the trivial equilibrium ; the axial equilibria: and are always unstable;(ii)the axial equilibrium is always asymptotically stable.

The stability properties of positive (coexistence) equilibrium for the case of weak and strong Allee effect are stated in Theorems 8 and 9.

Theorem 8. Stability of coexistence equilibrium for weak Allee effect :
suppose is the Jacobian matrix at coexistence equilibrium . Equilibrium is asymptotically stable if one of the following mutually exclusive conditions holds:(i); and .(ii), and

Proof. The characteristics equation of is given by (i)If ; ; and then ; hence and the result follows.(ii)If is an eigenvalue of and , then is also an eigenvalue. Using , we have that . Therefore the stability of follows.

Similarly we have Theorem 9 for the stability of coexistence equilibrium for strong Allee effect case.

Theorem 9. Stability of coexistence equilibrium for strong Allee effect :
suppose is the Jacobian matrix evaluated at the coexistence equilibrium . Equilibrium is asymptotically stable if one of the following mutually exclusive conditions holds:(i); ; and (ii), , and

Based on the above theorems it can be seen that the stability properties of both trivial and axial equilibrium points are not dependent on (order of fractional derivative). But may influence significantly the stability of coexistence equilibrium point. Coexistence point can be asymptotically stable although the eigenvalue of Jacobian matrix has positive real part, provided that conditions of Theorem 8(ii) or Theorem 9(ii) are met. This is in contrast to the coexistence equilibrium point of the integer-order model (2) where coexistence point is asymptotically stable only if all real parts of the eigenvalues of the Jacobian matrix are negative [9].

3. Numerical Simulations

To solve system (6), we implement a nonstandard Grünwald-Letnikov scheme which is a combination of the Grünwald-Letnikov approximation [20] and the nonstandard finite difference (NSFD) method [21, 22]. According to [20], the explicit (or implicit) Grünwald-Letnikov (GL) approximation for a fractional differential equation with initial valueis given bywhere ; ; and Here, represents the time step of numerical integration. The Grünwald-Letnikov approximation is proceeding iteratively but the sum in the scheme becomes longer and longer which represents the memory effects. Scherer et al. [20] have shown that the coefficient is positive and satisfies for . Observe that in the standard Grünwald-Letnikov approximation (21), the right hand side of (20) is approximated locally. We implement a nonstandard method which is adopted from the NSFD method [23]. A numerical scheme for an initial value problem is called a NSFD method if at least one of the following conditions is satisfied [21, 22]:(i) The left hand side is approximated by the generalization of forward difference schemeThe nonnegative denominator function has to satisfy .(ii) The approximation of is nonlocal.

By implementing the Grünwald-Letnikov approximation for the fractional derivative and the nonlocal approximation for the right hand side of system (6) we get the following scheme:Observe that scheme (24) is explicit and hence it is simple and easy to be implemented. Besides that, the nonstandard Grünwald-Letnikov scheme (24) also maintains the positivity solutions.

To verify our stability analysis as well as the effectiveness our numerical scheme, we perform some numerical simulations. First we use hypothetic values of parameters , , , , , , , , and . Model (6) with these parameters has four equilibrium points: , , , and . According to Theorem 7, only axial equilibrium is stable for any order of fractional derivative , where . This stability behavior is confirmed by our numerical solutions; see Figures 1 and 2. This shows that strong Allee effect may lead to an extinction of prey population. It is shown in Figure 2 that our numerical solutions for are convergent to the axial equilibrium , indicating that equilibrium is stable asymptotically for any order of fractional derivative. Detail observation shows that as the order of fractional derivative increases the convergence of solution is faster and the solution of system (6) closes to the integer-order model (2).

Next, we set and and keep the rest of parameters as in Figure 1. This weak Allee case has a trivial equilibrium ; two axial equilibria: and ; and two interior points: and . Theorems 7 and 9 state that axial equilibrium is locally stable for and interior point is locally stable if , and other equilibria are always unstable. Such behavior is in accordance with our numerical results depicted in Figures 3(a) and 3(b). It is clearly seen in Figure 3(a) that produces bistable dynamic where depending on the initial values, solutions may be convergent to the extinction of prey point or to interior point. In other words, the solution of system (6) is highly sensitive to the initial conditions. An initially relatively small prey will converge to the prey extinction point. On the other hand, if the prey is initially relatively large then prey and predator will coexist. If we increase the order of derivative such that , the axial equilibrium is still locally stable but the interior point becomes unstable; see Figure 3(b). In latter case, there exists a stable limit cycle which shows that both prey and predator are fluctuating around the interior point. However, the appearance of limit cycle may be suppressed by increasing the coefficient of predator interference. For example we plot in Figure 4(a) the numerical solution using the same parameters as in Figure 3(b) but with . We see that the interior point is now stable while the axial equilibrium point is unstable. It can be said that relatively large predator interference can stabilize the interior point. On the other hand, strong Allee effect can destabilize or even remove the interior equilibrium and can cause the extinction of prey population. For example, we show numerical simulation using parameters the same as in Figure 3(b) except . This simulation shows that all initial values converge to the axial equilibrium which shows the extinction of prey population but predator species can still survive in the habitat because there is an enough environmental protection; see Figure 4(b).

Finally, we compare our numerical results obtained by the NSGL scheme to those obtained by the standard GL scheme using parameters , , , , , , , and ; see Figure 5. We see that both NSGL and GL schemes using produce solutions where their difference cannot be observed in the scale of Figure 5. Using , the numerical solutions of both NSGL and GL schemes are in excellent agreement with those of both schemes using . If we take time step , both schemes have comparable solutions which are initially distorted from solutions with much smaller time step; see Figure 5(a). In Figure 5(b), we plot solutions using relatively large time step ( and ). Although the NSGL scheme has solutions which are quantitatively different from solution with , nevertheless those solutions still have the same behavior as before; that is, they are always positive and convergent to the correct equilibrium point. However, the GL scheme in this case gives unrealistic negative value for prey population. If the time step is further increased, the GL scheme will be unstable and leads to blowing up solutions.

4. Conclusion

The dynamic of a fractional order modified Leslie-Gower model with Beddington-DeAngelis functional response and additive Allee effect has been analyzed. Our model has four types of equilibria that are the trivial (extinction of both prey and predator) equilibrium, two axial equilibria (the prey extinction point and the predator extinction point), and the interior (coexistence) point. The trivial and the predator extinction for both weak and strong Allee effects are always unstable. For the case of weak Allee effect, the prey extinction is conditionally stable while for that of strong Allee effect, the prey extinction is always stable. Our analysis also shows that the interior point for both weak and Allee effects is conditionally stable. The order of fractional derivative may influence the stability of interior point. Here, when the order is larger than critical order , then the interior point may be destabilized. These dynamical properties are confirmed by our NSGL schemes which shows the effectiveness of NSGL scheme. It is also shown that the NSGL scheme preserves the positivity of numerical solutions. Furthermore, our numerical results show that the NSGL scheme produces numerical solutions which satisfy the dynamical behavior of our model. However, the standard GL scheme may fail to preserve such properties; for example, it can produce nonrealistic negative solutions.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the Directorate of Research and Community Service, The Directorate General of Strengthening Research and Development, and the Ministry of Research, Technology and Higher Education (Brawijaya University), Indonesia, Contract no. 137/SP2H/LT/DRPM/III/2016 dated March 10, 2016, and Contract no. 460.18/UN10.C10/PN/2017 dated April 18, 2017.