Home Mathematics Bifurcation and chaos in a discrete predator-prey system of Leslie type with Michaelis-Menten prey harvesting
Article Open Access

Bifurcation and chaos in a discrete predator-prey system of Leslie type with Michaelis-Menten prey harvesting

  • Jialin Chen , Zhenliang Zhu , Xiaqing He and Fengde Chen EMAIL logo
Published/Copyright: August 18, 2022

Abstract

In this paper, a discrete Leslie-Gower predator-prey system with Michaelis-Menten type harvesting is studied. Conditions on the existence and stability of fixed points are obtained. It is shown that the system can undergo fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation by using the center manifold theorem and bifurcation theory. Numerical simulations are presented to illustrate the main theoretical results. Compared to the continuous analog, the discrete system here possesses much richer dynamical behaviors including orbits of period-16, 21, 35, 49, 54, invariant cycles, cascades of period-doubling bifurcation in orbits of period-2, 4, 8, and chaotic sets.

MSC 2010: 34C25; 92D25; 34D20; 34D40

1 Introduction

Due to its universal existence and importance, the predation relationship between predator and prey is one of the dominant themes in ecology. Since the pioneering works of Lotka and Volterra, their classical models (called Lotka-Volterra predator-prey models) have been reasonably modified to incorporate real biological backgrounds. Among these modifications are Leslie-Gower models, which have been extensively studied (see, for example, [1,2, 3,4] and the references cited therein). In such models, the environmental carrying capacity of the predator species is determined by the abundance of the prey species. A general Leslie-Gower predator-prey model can be written as follows:

d x d t = r x 1 x K ϕ ( x ) y , d y d t = s y 1 h y x ,

where x and y are the densities of the prey and predator, respectively; ϕ is a functional response, and all the parameters are positive constants.

As we know, humans need to exploit natural resources to satisfy their own needs. To ensure the sustainable development of the ecosystem and to maximize economic benefits, harvesting models have therefore attracted the attention of many scholars [1,5,6, 7,8,9, 10,11,12, 13,14,15, 16,17,18, 19,20]. The study of harvesting models can provide suggestions to relevant government agencies on the development of policies to regulate and discipline. Among different forms of harvesting are the constant harvesting and linear harvesting proposed by May et al. [5]. In constant harvesting, the harvest rate is independent of the number of the harvested population, which is not practical. In linear harvesting [7,9,10, 11,12,13, 14,15,19], the harvest rate is proportional to the size of the harvested species, and it has the expression h = q E x . Note that when E or x tends to infinity and the other one is fixed, h ( x ) would tend to infinity. Obviously, this is against the fact that the harvesting capacity and the size of species are limited in reality. To overcome the drawbacks of these two harvesting forms, Clark and Mangel [6] proposed nonlinear harvesting, i.e., Michaelis-Menten type harvesting h ( x ) = q E x m E + n x . If E , then h q x m , the linear harvesting and while if x , then h q E n , the constant harvesting. Therefore, nonlinear harvesting is much more realistic, and it has attracted the attention of many researchers (to name a few, see [1,8,11,12,14,16, 17,18,20]).

The aforementioned works mainly focus on continuous predator-prey systems and obtained results include stability, bifurcation, limit cycles, and so on [21,22,23,24,25]. However, when species have nonoverlapping generations or their sizes are too small, discrete models described by difference equations are more appropriate than continuous-time ones. Over the past few decades, dynamic behaviors of discrete predator-prey systems have been widely studied, some of which can be found in [26,27,28, 29,30,31, 32,33,34, 35,36,37, 38,39,40] and the references cited therein. All these investigations have demonstrated that discrete systems tend to have more complex dynamic behaviors than continuous ones. In particular, He and Lai [28] studied a discrete Lotka-Volterra predator-prey system of Leslie-Gower with Holling type II functional response, whereas Ren et al. [41] proposed a system with Crowley-Martin functional response. In [36], Zhao and Yan investigated the aforementioned system with a modified Holling-Tanner functional response. In addition, many scholars have considered adding harvesting terms to these systems [42,43, 44,45]. For example, Hu and Cao [42] studied a discrete system of Holling type II functional response and Leslie type with constant-yield prey harvesting; Zhu et al. [44] proposed a discrete May type cooperative model incorporating Michaelis-Menten type harvesting. These papers mainly investigated the dynamical behaviors of bifurcation and chaos of the systems.

The aforementioned discussion motivates us to study a discrete Leslie-Gower predator-prey model with nonlinear harvesting. Precisely, the model is based on the following continuous Leslie-Gower predator-prey model with Michaelis-Menten type prey-harvesting and linear functional response proposed by Gupta et al. [46],

(1.1) d x d t = r x 1 x K α x y q E x m E + n x , d y d t = s y 1 h y x ,

where x and y denote the densities of prey and predator at time t , respectively, and r , K , α , q , E , m , n , s , and h are all positive constants. Here, r is the intrinsic growth rate of the prey and K represents the carrying capacity of the prey in the absence of predators; α is the maximum predation rate of the predator and s is the intrinsic grow rate of the predator; q stands for the harvesting coefficient and E is the harvesting effort of the prey species; h represents the number of prey required to provide one predator at equilibrium when y equals to x h .

For the sake of simplicity, we first rescale system (1.1) by introducing

x ¯ = x K , y ¯ = h K y , t ¯ = r t .

Dropping the bars, system (1.1) becomes

(1.2) d x d t = x 1 x b y e c + x , d y d t = p y 1 y x ,

where b = α K h r , c = m E n K , e = q E r n K , and p = s r are positive constants. Then applying the forward Euler scheme to system (1.2) yields the model to be studied in this paper,

(1.3) x x + δ x 1 x b y e c + x , y y + δ p y 1 y x ,

where δ is the step size. The aim of this paper is to study the dynamical behaviors of system (1.3) including the stability of fixed points and bifurcation phenomena.

The rest of the paper is arranged as follows. We first study the existence and stability of fixed points in Section 2. Then, in Section 3, we show that system (1.3) can undergo fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation under appropriate conditions on parameter values. The feasibility of the main theoretic results is illustrated with numerical simulations in Section 4. The paper ends with a brief conclusion.

2 The existence and stability of fixed points

2.1 The existence of fixed points

We start with the existence of fixed points of system (1.3). Note that fixed points of (1.3) are solutions of

(2.1) x = x + δ x 1 x b y e c + x , y = y + δ p y 1 y x .

First, we consider boundary fixed points, that is, x 0 and y = 0 . It follows easily from (2.1) that x satisfies

(2.2) x 2 + ( c 1 ) x + e c = 0 .

Let Δ 1 denotes the discriminant, namely,

(2.3) Δ 1 = ( c + 1 ) 2 4 e .

A necessary condition on the existence of boundary fixed points is Δ 1 0 , i.e., ( c + 1 ) 2 4 e . In this case, the roots of (2.2) can be written as follows:

(2.4) x 1 = 1 c Δ 1 2 , x 2 = 1 c + Δ 1 2 .

Specifically, when Δ 1 = 0 , x 1 , and x 2 collide, and we denote it as

(2.5) x 3 = 1 c 2 .

Taking into consideration of seeking positive solutions of (2.2), we easily obtain the following result on the existence of boundary fixed points of (1.3).

Theorem 2.1

The following statements on the boundary fixed points of (1.3) hold.

  1. If e > c , c < 1 , and ( c + 1 ) 2 > 4 e , then there are two boundary fixed points E 1 ( x 1 , 0 ) and E 2 ( x 2 , 0 ) ;

  2. There is only the boundary fixed point E 3 ( x 3 , 0 ) if c < 1 and ( c + 1 ) 2 = 4 e ;

  3. There is only the boundary fixed point E 2 ( x 2 , 0 ) if either e < c or e = c < 1 .

Here, x 1 and x 2 are given in (2.4), and x 3 is given in (2.5).

Now, we consider positive fixed points, which satisfy

(2.6) 1 x b y e c + x = 0 , y = x

by (2.1). Substituting y = x into the first equation of system (2.6) produces

(2.7) ( b + 1 ) x 2 + ( b c + c 1 ) x + e c = 0 .

Similarly as for the discussion of boundary fixed points, let

Δ 2 = [ c ( b + 1 ) + 1 ] 2 4 e ( b + 1 ) .

Then we require Δ 2 0 or [ c ( b + 1 ) + 1 ] 2 4 e ( b + 1 ) . Under this condition, when Δ 2 > 0 , equation (2.7) has two solutions,

(2.8) x 1 = 1 2 1 b + 1 c Δ 2 b + 1 , x 2 = 1 2 1 b + 1 c + Δ 2 b + 1 ,

and when Δ 2 = 0 , it has a unique root,

(2.9) x 3 = 1 2 1 b + 1 c .

Simple calculations give the existence of positive fixed points of (1.3), which is summarized below.

Theorem 2.2

The following statements are valid for the positive fixed points of (1.3).

  1. If e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , there are two positive fixed points E 1 ( x 1 , y 1 ) and E 2 ( x 2 , y 2 ) ;

  2. If 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) , there is only one positive fixed point E 3 ( x 3 , y 3 ) ;

  3. If either ( e = c and 1 b + 1 > c ) or e < c , there is only one positive fixed point E 2 ( x 2 , y 2 ) .

Here, x i are given in (2.8) or (2.9), and y i = x i , i = 1 , 2 , 3 .

2.2 The stability of fixed points

In this section, we study the (local) stability of the fixed points obtained earlier. Denote the Jacobian matrix of system (1.3) evaluated at a fixed point E ( x , y ) by J ( E ) , where

J ( E ) = J 11 J 12 J 21 J 22

with

J 11 = 1 + δ 1 x b y e c + x + δ x 1 + e ( c + x ) 2 , J 12 = δ b x , J 21 = δ p y 2 x 2 , J 22 = 1 + δ p 1 y x δ p y x .

We write the characteristic equation of J ( E ) as F ( λ ) = λ 2 + B λ + C = 0 . E can be classified according to the two roots λ 1 and λ 2 of F ( λ ) as follows (see, for example, [47, pp. 27]).

Definition 2.1

The fixed point E is called

  1. a sink if λ 1 < 1 and λ 2 < 1 , and it is locally asymptotically stable;

  2. a source if λ 1 > 1 and λ 2 > 1 , and it is unstable;

  3. a saddle if either ( λ 1 > 1 and λ 2 < 1 ) or ( λ 1 < 1 and λ 2 > 1 );

  4. nonhyperbolic if either λ 1 = 1 or λ 2 = 1 .

The following result gives sufficient and necessary conditions on the moduli of λ 1 and λ 2 with respect to 1 in the case of F ( 1 ) > 0 .

Lemma 2.1

[47, Lemma 2] Assume that F ( 1 ) > 0 . Then

  1. λ 1 < 1 and λ 2 < 1 if and only if F ( 1 ) > 0 and C < 1 ;

  2. λ 1 > 1 and λ 2 > 1 if and only if F ( 1 ) > 0 and C > 1 ;

  3. ( λ 1 > 1 and λ 2 < 1 ) or ( λ 1 < 1 and λ 2 > 1 ) if and only if F ( 1 ) < 0 ;

  4. λ 1 = 1 and λ 2 1 if and only if F ( 1 ) = 0 and B 0 , 2;

  5. λ 1 and λ 2 are conjugate complex roots and λ 1 = λ 2 = 1 if and only if B 2 4 C < 0 and C = 1 .

We mention that F ( 1 ) > 0 and F ( 1 ) = 0 imply that B 0 . Therefore, B 0 is redundant in statement 4 of the aforementioned lemma, which will be dropped in the coming discussion.

Similarly, we have the following result for the case of F ( 1 ) < 0 .

Lemma 2.2

Assume that F ( 1 ) < 0 . Then

  1. λ 1 > 1 and λ 2 > 1 if and only if F ( 1 ) < 0 ;

  2. ( λ 1 > 1 and λ 2 < 1 ) or ( λ 1 < 1 and λ 2 > 1 ) if and only if F ( 1 ) > 0 ;

  3. λ 1 = 1 and λ 2 1 if and only if F ( 1 ) = 0 .

Proof

Note that F ( λ ) as a quadratic with the leading coefficient being positive. When F ( 1 ) < 0 , we know that F ( λ ) = 0 has two real distinct roots λ 1 and λ 2 on the two sides of 1. Without loss of generality, say λ 1 < 1 < λ 2 . Then the results follow immediately from the intermediate value theorem.□

Theorem 2.3

Under the conditions on the existence of boundary fixed points of (1.3) in Theorem 2.1,

  1. E 1 ( x 1 , 0 ) is always a source and it is unstable;

  2. E 2 ( x 2 , 0 ) is

    1. a saddle if δ < 2 ( c + x 2 ) Δ 1 x 2 ;

    2. a source if δ > 2 ( c + x 2 ) Δ 1 x 2 ;

    3. nonhyperbolic if δ = 2 ( c + x 2 ) Δ 1 x 2 ;

  3. E 3 ( x 3 , 0 ) is always nonhyperbolic.

Proof

At a boundary fixed point E i ( i = 1 , 2 , 3 ), the Jacobian matrix reduces to

J ( E i ) = 1 + δ x i 1 + e ( c + x i ) 2 δ b x i 0 1 + δ p

as 1 x i = e c + x i . Thus, the eigenvalues of J ( E i ) are λ 1 = 1 + δ x i 1 + e ( c + x i ) 2 and λ 2 = 1 + δ p > 1 .

  1. Since x 1 = 1 2 ( 1 c Δ 1 ) , we obtain

    ( c + x 1 ) 2 e = 1 4 ( 1 + c Δ 1 ) 2 e = 1 2 ( Δ 1 Δ 1 c Δ 1 ) = Δ 1 2 [ Δ 1 ( 1 + c ) ] < 0

    or ( c + x 1 ) 2 < e . It follows that λ 1 > 1 . According to Definition 2.1, E 1 is a source and it is unstable when it exists.

  2. It follows from x 2 = 1 c + Δ 1 2 , we can obtain ( c + x 2 ) 2 e = Δ 1 2 ( Δ 1 + 1 + c ) > 0 , and hence, λ 1 < 1 . Furthermore, since x 2 satisfies 1 x 2 e c + x 2 = 0 , we have

    δ x 2 1 + e ( c + x 2 ) 2 = δ x 2 1 + 1 x 2 c + x 2 = δ Δ 1 x 2 c + x 2 .

    Therefore, if δ < 2 ( c + x 2 ) Δ 1 x 2 holds, then λ 1 < 1 , and hence, E 2 is a saddle. This proves (ii)(a). In a similar way, we can prove (ii)(b) and (ii)(c).

  3. It is easy to obtain that ( c + x 3 ) 2 = e since Δ 1 = 0 . Hence, the eigenvalues of J ( E 3 ) are λ 1 = 1 and λ 2 = 1 + δ p > 1 . According to Definition 2.1, E 3 ( x 3 , 0 ) is always nonhyperbolic.

This completes the proof.□

Now, the characteristic equation of the Jacobian matrix J of system (1.3) at each positive fixed point E i ( i = 1 , 2 , 3 ) is

F ( λ ) λ 2 ( 2 + G δ ) λ + ( 1 + G δ + H p δ 2 ) = 0 ,

where

G = x i e ( c + x i ) 2 1 p , H ( x i ) = x i ( b + 1 ) e ( c + x i ) 2 .

Thus,

F ( 1 ) = H p δ 2 , F ( 1 ) = 4 + 2 G δ + H p δ 2 .

We use Lemmas 2.1 and 2.2 to investigate the stability of the positive fixed points one by one.

Theorem 2.4

Assume that e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) . Then the positive fixed point E 1 in Theorem 2.2 is

  1. a saddle if 0 < δ < G + G 2 4 H p H p ;

  2. a source if δ > G + G 2 4 H p H p ;

  3. nonhyperbolic if δ = G + G 2 4 H p H p .

Proof

We first determine the sign of F ( 1 ) , or equivalently the sign of H . In fact, since E 1 ( x 1 , y 1 ) satisfies (2.7), we obtain

H ( x 1 ) = x 1 ( b + 1 ) e ( c + x 1 ) 2 = x 1 ( b + 1 ) 1 ( b + 1 ) x 1 c + x 1 = x 1 c + x 1 [ 2 ( b + 1 ) x 1 + c ( b + 1 ) 1 ] .

Noting x 1 = y 1 = 1 2 1 b + 1 c Δ 2 b + 1 , we have H ( x 1 ) = Δ 2 x 1 c + x 1 < 0 , and hence, F ( 1 ) < 0 .

As H ( x 1 ) < 0 , F ( 1 ) = 0 as a quadratic equation of δ has a negative solution G G 2 4 H p H p and a positive solution G + G 2 4 H p H p . If 0 < δ < G + G 2 4 H p H p then F ( 1 ) > 0 . By Lemma 2.2(ii), E 1 is a saddle. This proves (i). (ii) and (iii) can be proved similarly by using Lemma 2.2(i) and (iii), respectively.□

For the positive fixed point E 2 ( x 2 , y 2 ) , a similar calculation as for H ( x 1 ) will produce H ( x 2 ) = Δ 2 x 2 c + x 2 > 0 , and hence, F ( 1 ) > 0 . With the assistance of Lemma 2.1, we easily obtain

Theorem 2.5

Suppose one of the following three conditions holds,

  1. e > c , 1 b + 1 > c , and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) ;

  2. e = c and 1 b + 1 > c ;

  3. e < c .

Then the positive fixed point E 2 of system (1.3) is
  1. a sink if one of the following conditions holds,

    1. 2 H p G < 0 and 0 < δ < G H p ;

    2. G < 2 H p and 0 < δ < G G 2 4 H p H p ;

  2. a source if one of the following conditions holds,

    1. 2 H p G < 0 and δ > G H p ;

    2. G < 2 H p and δ > G + G 2 4 H p H p ;

    3. G 0 ;

  3. a saddle if G < 2 H p and G G 2 4 H p H p < δ < G + G 2 4 H p H p ;

  4. nonhyperbolic if one of the following conditions holds:

    1. G < 2 H p , δ = G ± G 2 4 H p H p , and δ 4 G ;

    2. 2 H p < G < 0 and δ = G H p .

Finally, from Theorem 2.2, the positive fixed point E 3 exists when [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) and 1 b + 1 > c . It follows that F ( 1 ) = 0 , and hence,

Theorem 2.6

When [ c ( b + 1 ) + 1 ] 2 = 4 e ( b + 1 ) and 1 b + 1 > c , the positive fixed point E 3 of (1.3) is always nonhyperbolic.

3 Bifurcation analysis

In this section, we analyze different bifurcation types at fixed points of system (1.3) by using the center manifold theorem [48] and bifurcation theory [49,50]. We begin with the fold bifurcation.

3.1 Fold bifurcation

Recall from Theorem 2.2.2 that if

(3.1) e = [ c ( b + 1 ) + 1 ] 2 4 ( b + 1 ) , 1 b + 1 > c ,

then system (1.3) only has one positive fixed point E 3 ( x 3 , y 3 ) and the eigenvalues of the Jacobian matrix J ( E 3 ) are λ 1 = 1 and λ 2 = 1 + G δ . Suppose that

(3.2) e ( x 3 + p ) ( x 3 + c ) 2 x 3 , e ( δ p + δ x 3 2 ) ( x 3 + c ) 2 δ x 3 .

Then, λ 2 1 .

Let u = x x 3 , v = y y 3 , and η = e e . Then system (1.3) can be rewritten as follows:

(3.3) u η v a 1 u + a 2 η + a 3 v + a 11 u 2 + a 12 u η + a 13 u v + O ( ( u + v + η ) 3 ) η b 1 u + b 3 v + b 11 u 2 + b 13 u v + b 33 v 2 + O ( ( u + v + η ) 3 ) ,

where

a 1 = 1 + δ x 3 e ( c + x 3 ) 2 1 , a 2 = δ x 3 c + x 3 , a 3 = δ b x 3 , a 11 = δ e c ( c + x 3 ) 3 1 , a 12 = δ c ( c + x 3 ) 2 , a 13 = δ b , b 1 = δ p , b 3 = 1 δ p , b 11 = δ p x 3 , b 13 = 2 δ p x 3 , b 33 = δ p x 3 .

We choose

T 1 = a 2 ( a 1 λ 2 ) a 2 λ 2 b 3 0 1 λ 2 0 a 2 b 1 0 b 1 ,

which is invertible. Then with the transformation

u η v = T 1 x ˜ η 1 y ˜ ,

we transform (3.3) into

(3.4) x ˜ η 1 y ˜ 1 1 0 0 1 0 0 0 λ 2 x ˜ η 1 y ˜ + ϕ ( x ˜ , y ˜ , η 1 ) 0 ψ ( x ˜ , y ˜ , η 1 ) ,

where

ϕ ( x ˜ , y ˜ , η 1 ) = a 11 b 1 + b 11 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 u 2 + a 12 a 2 ( 1 λ 2 ) u η + a 13 b 1 + b 13 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 u v + b 33 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 v 2 + O ( ( u + v + η ) 3 ) , ψ ( x ˜ , y ˜ , η 1 ) = b 11 ( a 1 λ 2 ) a 11 b 1 ( 1 λ 2 ) b 1 u 2 a 12 1 λ 2 u η + b 13 ( a 1 λ 2 ) a 13 b 1 ( 1 λ 2 ) b 1 u v + b 33 ( a 1 λ 2 ) ( 1 λ 2 ) b 1 v 2 + O ( ( u + v + η ) 3 ) , u = a 2 ( a 1 λ 2 ) x ˜ + a 2 η 1 + ( λ 2 b 3 ) y ˜ , η = ( 1 λ 2 ) η 1 , v = a 2 b 1 x ˜ + b 1 y ˜ .

By the center manifold theory, in a small neighborhood of η 1 = 0 , there exists a center manifold W c ( 0 ) of (3.4) at the fixed point ( x ˜ , y ˜ ) = ( 0 , 0 ) , which can be represented as follows:

W c ( 0 ) = { ( x ˜ , y ˜ , η 1 ) R 3 y ˜ = h ( x ˜ , η 1 ) , h ( 0 , 0 ) = 0 , D h ( 0 , 0 ) = 0 } ,

where x ˜ and η 1 are sufficiently small. Assume that the expression of h is

(3.5) h ( x ˜ , η 1 ) = m 1 x ˜ 2 + m 2 x ˜ η 1 + m 3 η 1 2 + O ( ( x ˜ + η 1 ) 3 ) ,

which must satisfy

(3.6) h ( x ˜ + η 1 + ϕ ( x ˜ , h ( x ˜ , η 1 ) , η 1 ) , η 1 ) = λ 2 h ( x ˜ , η 1 ) + ψ ( x ˜ , h ( x ˜ , η 1 ) , η 1 ) .

Substituting (3.5) into (3.6) and comparing the coefficients of terms like x ˜ k η 1 l in the resultant, we obtain

m 1 = a 2 2 ( a 1 λ 2 ) b 1 ( 1 λ 2 ) 2 { ( a 1 λ 2 ) [ ( a 1 λ 2 ) b 11 a 11 b 1 + b 13 b 1 ] a 13 b 1 2 + b 33 b 1 2 } , m 2 = 2 a 2 2 ( a 1 λ 2 ) ( a 1 λ 2 ) b 11 a 11 b 1 b 1 ( 1 λ 2 ) 2 a 12 a 2 ( a 1 λ 2 ) 1 λ 2 + a 2 2 [ ( a 1 λ 2 ) b 13 a 13 b 1 ] ( 1 λ 2 ) 2 2 m 1 1 λ 2 , m 3 = a 2 2 [ ( a 1 λ 2 ) b 11 a 11 b 1 ] b 1 ( 1 λ 2 ) 2 a 12 a 2 1 λ 2 m 1 + m 2 1 λ 2 .

Therefore, the map (3.4) restricted to the center manifold can be written as follows:

F 1 : x ˜ x ˜ + η 1 + k 1 x ˜ 2 + k 2 x ˜ η 1 + k 3 η 1 2 + O ( ( x ˜ + η 1 ) 3 ) ,

where

(3.7) k 1 = a 2 2 ( a 1 λ 2 ) ( a 1 λ 2 ) a 11 b 1 + b 11 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 + a 13 b 1 + b 13 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) + a 2 b 1 b 33 ( b 3 λ 2 ) 1 λ 2 , k 2 = ( a 1 λ 2 ) 2 a 2 2 a 11 b 1 + 2 a 2 2 b 11 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 + a 12 + a 2 2 a 13 b 1 + a 2 2 b 13 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) , k 3 = a 2 2 a 11 b 1 + a 2 2 b 11 ( b 3 λ 2 ) a 2 ( 1 λ 2 ) b 1 + a 12 .

Since F 1 ( 0 , 0 ) = 0 , F 1 x ˜ ( 0 , 0 ) = 1 , F 1 η 1 ( 0 , 0 ) = 1 , 2 F 1 x ˜ 2 ( 0 , 0 ) = 2 k 1 , and 2 F 1 x ˜ η 1 ( 0 , 0 ) = k 2 , we obtain the following result.

Theorem 3.1

In addition to (3.1) and (3.2), suppose that k 1 0 or k 2 0 . Then system (1.3) undergoes a fold bifurcation at E 3 , where k 1 and k 2 are given by (3.7). Moreover, the fixed points E 1 and E 2 bifurcate from E 3 for e < e , coalesce at E 3 for e = e , and disappear for e > e .

3.2 Flip bifurcation

Next we discuss the flip bifurcations of system (1.3).

Let

F A = ( c , e , δ ) δ = 2 ( c + x 2 ) Δ 1 x 2 , c , e > 0 ,

where x 2 and Δ 1 are given by (2.3) and (2.4), respectively. System (1.3) can undergo flip bifurcation at the boundary fixed point E 2 ( x 2 , 0 ) when parameters belong to F A . Since a center manifold of system (1.3) at E 2 is y = 0 and system (1.3) restricted to it is the logistic model:

x f ( x ) = x + δ x 1 x e c + x .

Its nontrivial fixed point is x 2 = 1 c + Δ 1 2 . It is easy to see that

f x 2 = 1 + δ x 2 1 + e ( c + x 2 ) 2 < 0

when parameters vary in a small neighborhood of F A . Thus, flip bifurcation can occur (see Figure 2). In this case, the predator species becomes extinct, and the prey species undergoes the flip bifurcation to chaos by choosing the appropriate value of the bifurcation parameter δ .

Since E 1 is always unstable, we now focus on flip bifurcation at E 2 due to the biological significance. Here, we choose δ as the bifurcation parameter.

From Theorem 2.5(iv(a)), we can easily obtain that one of the eigenvalues of J ( E 2 ) is 1 and the other one is neither 1 nor 1 . We rewrite the conditions in Theorem 2.5(iv(a)) as the following two sets:

F B 1 = ( b , c , e , p , δ ) (i) b , c , e , p > 0 , δ = G G 2 4 H p H p , G < 2 H p , δ 4 G , (ii) e > c > 0 , 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , or e = c and 1 b + 1 > c , or e < c

and

F B 2 = ( b , c , e , p , δ ) (i) b , c , e , p > 0 , δ = G + G 2 4 H p H p , G < 2 H p , δ 4 G , (ii) e > c > 0 , 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) , or e = c and 1 b + 1 > c , or e < c .

We shall see that flip bifurcation may undergo when parameters belong to F B 1 or F B 2 . We only provide the detail on discussing parameters belonging to F B 1 and omit that for F B 2 as it is similar.

Take parameter values ( b , c , e , p , δ 1 ) arbitrarily from F B 1 . Then system (1.3) with these parameter values becomes

(3.8) x x + δ 1 x 1 x b y e c + x , y y + δ 1 p y 1 y x .

The map (3.8) has a positive fixed point E 2 , whose eigenvalues are λ 1 = 1 and λ 2 = 3 + G δ 1 with λ 2 1 by Theorem 2.5. We consider a perturbation of (3.8) as follows:

(3.9) x x + ( δ 1 + δ ) x 1 x b y e c + x , y y + ( δ 1 + δ ) p y 1 y x ,

where δ 1 , which is a small perturbation parameter.

Let u = x x 2 and v = y y 2 , which transform the fixed point E 2 of map (3.9) into the origin. Moreover, (3.9) is transformed into

(3.10) u v c 1 u + c 2 v + c 11 u 2 + c 12 u v + c 111 u 3 + c 13 u δ + c 23 v δ + c 113 u 2 δ + c 123 u v δ + O ( ( u + v + δ ) 4 ) d 1 u + d 2 v + d 11 u 2 + d 12 u v + d 22 v 2 + d 111 u 3 + d 112 u 2 v + d 122 u v 2 + d 13 u δ + d 23 v δ + d 113 u 2 δ + d 123 u v δ + d 223 v 2 δ + O ( ( u + v + δ ) 4 ) ,

where

(3.11) c 1 = 1 + δ 1 x 2 e ( c + x 2 ) 2 1 , c 2 = δ 1 b x 2 , c 11 = δ 1 e c ( c + x 2 ) 3 1 , c 12 = δ 1 b , c 111 = δ 1 e c ( c + x 2 ) 4 , c 13 = x 2 e ( c + x 2 ) 2 1 , c 23 = b x 2 , c 113 = e c ( c + x 2 ) 3 1 , c 123 = b , d 1 = δ 1 p , d 2 = 1 δ 1 p , d 11 = d 22 = δ 1 p x 2 , d 12 = 2 δ 1 p x 2 , d 111 = d 122 = δ 1 p ( x 2 ) 2 , d 112 = 2 δ 1 p ( x 2 ) 2 , d 13 = d 23 = p , d 113 = d 223 = p x 2 , d 123 = 2 p x 2 .

With the transformation

u v = T 2 x ˜ y ˜ , where T 2 = c 2 c 2 1 c 1 λ 2 c 1 ,

the map (3.10) becomes

(3.12) x ˜ y ˜ 1 0 0 λ 2 x ˜ y ˜ + f ( x ˜ , y ˜ , δ ) g ( x ˜ , y ˜ , δ ) ,

where

f ( x ˜ , y ˜ , δ ) = [ c 11 ( λ 2 c 1 ) c 2 d 11 ] c 2 ( λ 2 + 1 ) u 2 + [ c 12 ( λ 2 c 1 ) c 2 d 12 ] c 2 ( λ 2 + 1 ) u v d 22 λ 2 + 1 v 2 + [ c 13 ( λ 2 c 1 ) c 2 d 13 ] c 2 ( λ 2 + 1 ) u δ + [ c 23 ( λ 2 c 1 ) c 2 d 23 ] c 2 ( λ 2 + 1 ) v δ + [ c 111 ( λ 2 c 1 ) c 2 d 111 ] c 2 ( λ 2 + 1 ) u 3 d 112 λ 2 + 1 u 2 v d 122 λ 2 + 1 u v 2 + [ c 113 ( λ 2 c 1 ) c 2 d 113 ] c 2 ( λ 2 + 1 ) u 2 δ + [ c 123 ( λ 2 c 1 ) c 2 d 123 ] c 2 ( λ 2 + 1 ) u v δ d 223 λ 2 + 1 v 2 δ + O ( ( u + v + δ ) 4 ) , g ( x ˜ , y ˜ , δ ) = [ c 11 ( 1 + c 1 ) + c 2 d 11 ] c 2 ( λ 2 + 1 ) u 2 + [ c 12 ( 1 + c 1 ) + c 2 d 12 ] c 2 ( λ 2 + 1 ) u v + d 22 λ 2 + 1 v 2 + [ c 13 ( 1 + c 1 ) + c 2 d 13 ] c 2 ( λ 2 + 1 ) u δ + [ c 23 ( 1 + c 1 ) + c 2 d 23 ] c 2 ( λ 2 + 1 ) v δ + [ c 111 ( 1 + c 1 ) + c 2 d 111 ] c 2 ( λ 2 + 1 ) u 3 + d 112 λ 2 + 1 u 2 v + d 122 λ 2 + 1 u v 2 + [ c 113 ( 1 + c 1 ) + c 2 d 113 ] c 2 ( λ 2 + 1 ) u 2 δ + [ c 123 ( 1 + c 1 ) + c 2 d 123 ] c 2 ( λ 2 + 1 ) u v δ + d 223 λ 2 + 1 v 2 δ + O ( ( u + v + δ ) 4 ) , u = c 2 ( x ˜ + y ˜ ) , v = ( 1 c 1 ) x ˜ + ( λ 2 + c 1 ) y ˜ .

Now we arrive at determining the center manifold W c ( 0 ) of (3.12) at the fixed point ( x ˜ , y ˜ ) = ( 0 , 0 ) in a small neighborhood of δ = 0 , which is expressed as follows:

W c ( 0 ) = { ( x ˜ , y ˜ , δ ) R 3 y ˜ = h ( x ˜ , δ ) , h ( 0 , 0 ) = 0 , D h ( 0 , 0 ) = 0 }

for x ˜ and δ sufficiently small. The function h must satisfy

(3.13) h ( x ˜ + f ( x ˜ , h ( x ˜ , δ ) , δ ) , δ ) = λ 2 h ( x ˜ , δ ) + g ( x ˜ , h ( x ˜ , δ ) , δ ) .

We write h as follows:

(3.14) h ( x ˜ , δ ) = n 1 x ˜ 2 + n 2 x ˜ δ + n 3 δ * 2 + O ( ( x ˜ + δ ) 3 ) .

Substituting (3.14) into (3.13) and comparing the corresponding coefficients of the left and right sides of the resultant, we can obtain

n 1 = c 2 [ c 11 ( 1 + c 1 ) + c 2 d 11 ] + ( 1 + c 1 ) [ ( d 22 c 12 ) ( 1 + c 1 ) c 2 d 12 ] 1 λ 2 2 , n 2 = ( 1 + c 1 ) [ c 23 ( 1 + c 1 ) + c 2 d 23 ] c 2 [ c 13 ( 1 + c 1 ) + c 2 d 13 ] c 2 ( λ 2 + 1 ) 2 , n 3 = 0 .

Therefore, the restricted map of (3.12) on the center manifold W c ( 0 ) is

(3.15) F 2 : x ˜ x ˜ + h 1 x ˜ 2 + h 2 x ˜ δ + h 3 x ˜ 2 δ + h 4 x ˜ δ * 2 + h 5 x ˜ 3 + O ( ( x ˜ + δ ) 4 ) ,

where

h 1 = 1 λ 2 + 1 { c 2 [ c 11 ( λ 2 c 1 ) c 2 d 11 ] ( 1 + c 1 ) [ c 12 ( λ 2 c 1 ) c 2 d 12 ] d 22 ( 1 + c 1 ) 2 } , h 2 = 1 c 2 ( λ 2 + 1 ) { c 2 [ c 13 ( λ 2 c 1 ) c 2 d 13 ] ( 1 + c 1 ) [ c 23 ( λ 2 c 1 ) c 2 d 23 ] } , h 3 = n 2 λ 2 + 1 { 2 c 2 [ c 11 ( λ 2 c 1 ) c 2 d 11 ] + ( λ 2 1 2 c 1 ) [ c 12 ( λ 2 c 1 ) c 2 d 12 ] } + n 1 c 2 ( λ 2 + 1 ) { c 2 [ c 13 ( λ 2 c 1 ) c 2 d 13 ] + ( λ 2 c 1 ) [ c 23 ( λ 2 c 1 ) c 2 d 23 ] } + 1 λ 2 + 1 { 2 n 2 d 22 ( 1 + c 1 ) ( λ 2 c 1 ) + c 2 [ c 113 ( λ 2 c 1 ) c 2 d 113 ] ( 1 + c 1 ) [ c 123 ( λ 2 c 1 ) c 2 d 123 ] d 223 ( 1 + c 1 ) 2 } , h 4 = n 2 c 2 ( λ 2 + 1 ) { c 2 [ c 13 ( λ 2 c 1 ) c 2 d 13 ] + ( λ 2 c 1 ) [ c 23 ( λ 2 c 1 ) c 2 d 23 ] } , h 5 = 1 λ 2 + 1 { 2 c 2 n 1 [ c 11 ( λ 2 c 1 ) c 2 d 11 ] + n 1 ( λ 2 1 2 c 1 ) [ c 12 ( λ 2 c 1 ) c 2 d 12 ] + c 2 2 [ c 111 ( λ 2 c 1 ) c 2 d 111 ] + c 2 2 d 112 ( 1 + c 1 ) + 2 n 1 d 22 ( 1 + c 1 ) ( λ 2 c 1 ) d 122 c 2 ( 1 + c 1 ) 2 } .

In order for map (3.15) to undergo a flip bifurcation, we require that the two discriminatory quantities α 1 and α 2 are not zero, where

α 1 = 2 F 2 x ˜ δ + 1 2 F 2 δ 2 F 2 x ˜ 2 ( 0 , 0 ) = h 2

and

α 2 = 1 6 3 F 2 x ˜ 3 + 1 2 2 F 2 x ˜ 2 2 ( 0 , 0 ) = h 5 + h 1 2 .

In summary, from the aforementioned discussion and theory in [49,50], we have arrived at the following result.

Theorem 3.2

If α 1 0 and α 2 0 , then system (1.3) undergoes a flip bifurcation at the fixed point E 2 ( x 2 , y 2 ) when the parameter δ varies in a small neighborhood of δ 1 . Moreover, if α 2 > 0 (resp., α 2 < 0 ), then the period-2 orbits that bifurcate from E 2 ( x 2 , y 2 ) are stable (resp., unstable).

3.3 Neimark-Sacker bifurcation

In the remaining of this section, we consider Neimark-Sacker bifurcation.

Since H ( x 1 ) < 0 , the characteristic equation of J ( E 1 ) cannot have two conjugate complex roots on the unit circle. From Theorem 2.5(iv(b)), we can obtain that the eigenvalues of J ( E 2 ) are a pair of conjugate complex numbers with moduli 1. The condition there can be written as follows:

H B = ( b , c , e , p , δ ) (i) b , c , e , p > 0 , δ = G H p , 2 H p < G < 0 , (ii) e > c > 0 , 1 b + 1 > c and [ c ( b + 1 ) + 1 ] 2 > 4 e ( b + 1 ) or e = c and 1 b + 1 > c , or e < c .

It is possible that Neimark-Sacker bifurcation can occur at the positive fixed point E 2 when parameter values belong to H B .

As before, take ( b , c , e , p , δ 2 ) arbitrarily from H B . We have δ = G H p . Choosing δ as a bifurcation parameter, we consider the following perturbation of (1.3),

(3.16) x x + ( δ 2 + δ ¯ ) x 1 x b y e c + x , y y + ( δ 2 + δ ¯ ) p y 1 y x ,

where δ ¯ 1 , which is a small perturbation parameter. Letting u = x x 2 and v = y y 2 , we rewrite (3.16) as follows:

(3.17) u v c 1 u + c 2 v + c 11 u 2 + c 12 u v + c 111 u 3 + O ( ( u + v ) 4 ) d 1 u + d 2 v + d 11 u 2 + d 12 u v + d 22 v 2 + d 111 u 3 + d 112 u 2 v + d 122 u v 2 + O ( ( u + v ) 4 ) ,

where c 1 , c 2 , c 11 , c 12 , c 111 , d 1 , d 2 , d 11 , d 12 , d 22 , d 111 , d 112 , and d 122 are given in (3.11) with replacing δ 1 with δ 1 for δ 2 + δ ¯ .

Note that the characteristic equation associated with the linearization of system (3.17) at ( u , v ) = ( 0 , 0 ) is given as follows:

(3.18) λ 2 + P ( δ ¯ ) λ + Q ( δ ¯ ) = 0 ,

where

P ( δ ¯ ) = 2 G ( δ 2 + δ ¯ ) , Q ( δ ¯ ) = 1 + G ( δ 2 + δ ¯ ) + H p ( δ 2 + δ ¯ ) 2 .

Since ( b , c , e , p , δ 2 ) H B , the eigenvalues are a pair of complex conjugate numbers λ and λ ¯ with moduli 1 by Theorem 2.5(iv(b)), where

λ , λ ¯ = P ( δ ¯ ) 2 ± i 2 4 Q ( δ ¯ ) P 2 ( δ ¯ ) = 1 + G ( δ 2 + δ ¯ ) 2 ± i ( δ 2 + δ ¯ ) 2 4 H p G 2 .

Then

λ = Q ( δ ¯ ) , l = d λ d δ ¯ δ ¯ = 0 = G 2 > 0 .

In addition, it is required that when δ ¯ = 0 , λ m 1 , and λ ¯ m 1 ( m = 1 , 2 , 3 , 4 ), which is equivalent to P ( 0 ) 2 , 0 , 1 , 2 . Note that ( b , c , e , p , δ 2 ) H B . Thus, P ( 0 ) 2 , 2. It suffices to require that P ( 0 ) 0 , 1 , which lead to

(3.19) G 2 2 H p , 3 H p .

Therefore, the eigenvalues λ and λ ¯ of (3.18) do not lie in the intersection of the unit circle with the coordinate axes when δ ¯ = 0 and (3.19) holds.

Next, we derive the normal form of (3.17) at δ ¯ = 0 .

Let δ ¯ = 0 , μ = 1 + G δ 2 2 , ω = δ 2 2 4 H p G 2 , and

T 3 = c 2 0 μ c 1 ω .

With the transformation

u v = T 3 x ˜ y ˜ ,

and we transform the map (3.17) into

(3.20) x ˜ y ˜ μ ω ω μ x ˜ y ˜ + f ˜ ( x ˜ , y ˜ ) g ˜ ( x ˜ , y ˜ ) ,

where

f ˜ ( x ˜ , y ˜ ) = c 11 c 2 u 2 + c 12 c 2 u v + c 111 c 2 u 3 + O ( ( x ˜ + y ˜ ) 4 ) , g ˜ ( x ˜ , y ˜ ) = c 11 ( μ c 1 ) c 2 d 11 c 2 ω u 2 + c 12 ( μ c 1 ) c 2 d 12 c 2 ω u v d 22 ω v 2 + [ c 111 ( μ c 1 ) c 2 d 111 ] c 2 ω u 3 d 112 ω u 2 v d 122 ω u v 2 + O ( ( x ˜ + y ˜ ) 4 ) , u = c 2 x ˜ , v = ( μ c 1 ) x ˜ ω y ˜ .

Straightforward computations give

f ˜ x ˜ x ˜ = 2 c 2 c 11 + 2 c 12 ( μ c 1 ) , f ˜ x ˜ y ˜ = c 12 ω , f ˜ y ˜ y ˜ = 0 , f ˜ x ˜ x ˜ x ˜ = 6 c 2 2 c 111 , f ˜ x ˜ x ˜ y ˜ = f ˜ x ˜ y ˜ y ˜ = f ˜ y ˜ y ˜ y ˜ = 0 , g ˜ x ˜ x ˜ = 2 ω { c 2 [ c 11 ( μ c 1 ) c 2 d 11 ] + ( μ c 1 ) [ c 12 ( μ c 1 ) c 2 d 12 ] d 22 ( μ c 1 ) 2 } , g ˜ x ˜ y ˜ = ( μ c 1 ) ( 2 d 22 c 12 ) + c 2 d 12 , g ˜ y ˜ y ˜ = 2 d 22 ω , g ˜ x ˜ x ˜ x ˜ = 6 c 2 ω { c 2 [ c 111 ( μ c 1 ) c 2 d 111 ] c 2 d 112 ( μ c 1 ) d 122 ( μ c 1 ) 2 } , g ˜ x ˜ x ˜ y ˜ = 2 c 2 [ c 2 d 112 + 2 d 122 ( μ c 1 ) ] , g ˜ x ˜ y ˜ y ˜ = 2 d 122 c 2 ω , g ˜ y ˜ y ˜ y ˜ = 0 .

In order for map (3.20) to undergo Neimark-Sacker bifurcation, we require that the following discriminatory quantity is not zero [49,50],

α = Re ( 1 2 λ ) λ ¯ 2 1 λ ξ 20 ξ 11 1 2 ξ 11 2 ξ 02 2 + Re ( λ ¯ ξ 21 ) δ ¯ = 0 ,

where

ξ 20 = 1 8 [ ( f ˜ x ˜ x ˜ f ˜ y ˜ y ˜ + 2 g ˜ x ˜ y ˜ ) + i ( g ˜ x ˜ x ˜ g ˜ y ˜ y ˜ 2 f ˜ x ˜ y ˜ ) ] , ξ 11 = 1 4 [ ( f ˜ x ˜ x ˜ + f ˜ y ˜ y ˜ ) + i ( g ˜ x ˜ x ˜ + g ˜ y ˜ y ˜ ) ] , ξ 02 = 1 8 [ ( f ˜ x ˜ x ˜ f ˜ y ˜ y ˜ 2 g ˜ x ˜ y ˜ ) + i ( g ˜ x ˜ x ˜ g ˜ y ˜ y ˜ + 2 f ˜ x ˜ y ˜ ) ] , ξ 21 = 1 16 [ ( f ˜ x ˜ x ˜ x ˜ + f ˜ x ˜ y ˜ y ˜ + g ˜ x ˜ x ˜ y ˜ + g ˜ y ˜ y ˜ y ˜ ) + i ( g ˜ x ˜ x ˜ x ˜ + g ˜ x ˜ y ˜ y ˜ f ˜ x ˜ x ˜ y ˜ f ˜ y ˜ y ˜ y ˜ ) ] .

The following result follows directly from the aformentioned analysis and the theory in [49,50].

Theorem 3.3

Suppose (3.19) holds and α 0 . Then system (1.3) undergoes Neimark-Sacker bifurcation at the positive fixed point E 2 when the parameter δ varies in a small neighborhood of δ 2 . Moreover, if α < 0 (resp., α > 0 ), then an attracting (resp., repelling) invariant closed curve bifurcates from the fixed point for δ > δ 2 (resp., δ < δ 2 ).

4 Numerical simulations

The purpose of this section is to present bifurcation diagrams and phase portraits of system (1.3) to confirm the feasibility of the main theoretical results and to show the complex dynamical behaviors through numerical simulations.

Example 4.1

(Fold bifurcation at the positive fixed point E 3 ). We choose e as the bifurcation parameter. With

(4.1) e [ 0 , 0.3 ] , b = 0.1 , c = 0.05 , p = 1.5 , δ = 1 ,

one obtain the bifurcation value e 0.253 and system (1.3) has only one positive fixed point E 3 ( 0.43 , 0.43 ) . It is easy to verify (3.1) and (3.2). In addition, the eigenvalues of J ( E 3 ) are λ 1 = 1 and λ 2 = 0.457 . Note that k 1 = 1.363 0 . Thus, all the conditions of Theorem 3.1 hold, and hence, fold bifurcation occurs at E 3 . Figure 1 agrees with Theorem 3.1 very well. Moreover, we see that when e < e , E 1 is unstable, while E 2 is stable.

Figure 1 
               Fold bifurcation diagram of system (1.3) in the 
                     
                        
                        
                           
                              (
                              
                                 e
                                 ,
                                 x
                              
                              )
                           
                        
                        \left(e,x)
                     
                   plane where parameters are given in (4.1) and the initial value is 
                     
                        
                        
                           
                              (
                              
                                 0.5
                                 ,
                                 0.2
                              
                              )
                           
                        
                        \left(0.5,0.2)
                     
                  . The dash curve corresponds to the unstable fixed point 
                     
                        
                        
                           
                              
                                 E
                              
                              
                                 1
                              
                              
                                 ∗
                              
                           
                        
                        {E}_{1}^{\ast }
                     
                  , and the solid curve corresponds to the stable fixed point 
                     
                        
                        
                           
                              
                                 E
                              
                              
                                 2
                              
                              
                                 ∗
                              
                           
                        
                        {E}_{2}^{\ast }
                     
                  . The fold bifurcation value is 
                     
                        
                        
                           
                              
                                 e
                              
                              
                                 ∗
                              
                           
                           ≈
                           0.253
                        
                        {e}^{\ast }\approx 0.253
                     
                  .
Figure 1

Fold bifurcation diagram of system (1.3) in the ( e , x ) plane where parameters are given in (4.1) and the initial value is ( 0.5 , 0.2 ) . The dash curve corresponds to the unstable fixed point E 1 , and the solid curve corresponds to the stable fixed point E 2 . The fold bifurcation value is e 0.253 .

Example 4.2

(Flip bifurcation at a boundary fixed point E 2 ). Take c = 0.045 and e = 0.048 . Then ( c + 1 ) 2 = 1.09 > 4 e = 0.192 . On the basis of Theorem 2.1 and (2.4), we know that system (1.3) has the boundary fixed points E 1 ( 0.003 , 0 ) and E 2 ( 0.95 , 0 ) . One can check that flip bifurcation emerges from the fixed point E 2 ( 0.95 , 0 ) at δ = 2.21 and ( c , e , δ ) = ( 0.045 , 0.048 , 2.21 ) F A (see Figure 2).

Figure 2 
               Bifurcation diagram of the logistic model 
                     
                        
                        
                           x
                           →
                           x
                           +
                           δ
                           x
                           
                              
                                 1
                                 −
                                 x
                                 −
                                 
                                    
                                       e
                                    
                                    
                                       c
                                       +
                                       x
                                    
                                 
                              
                           
                        
                        x\to x+\delta x\left(1-x-\frac{e}{c+x}\right)
                     
                   with 
                     
                        
                        
                           δ
                           ∈
                           
                              [
                              
                                 1
                                 ,
                                 3.5
                              
                              ]
                           
                        
                        \delta \in \left[1,3.5]
                     
                  .
Figure 2

Bifurcation diagram of the logistic model x x + δ x 1 x e c + x with δ [ 1 , 3.5 ] .

Example 4.3

(Flip bifurcation at the positive fixed point E 2 ). This time, we choose

(4.2) δ [ 0.8 , 1.5 ] , b = 0.1 , c = 0.45 , e = 0.4 , p = 2 .

According to Theorem 2.2(iv) and (2.8), system (1.3) has only one positive fixed point E 2 . After simple calculations, we can verify that flip bifurcation emerges from E 2 ( 0.54 , 0.54 ) at δ = 1.03 with α 1 = 1.93 , α 2 = 37.49 , and ( b , c , e , p , δ ) = ( 0.1 , 0.45 , 0.4 , 2 , 1.03 ) F B 2 . Figure 3 shows the feasibility of Theorem 3.2.

We can observe that the fixed point is stable for δ < 1.03 . When δ reaches 1.03, with the increase of δ , first, two points with a period-2 cycle bifurcated, then from there points with period-4 and period-8 bifurcated in sequence. The phase portraits associated with Figure 3 are displayed in Figure 4. For δ ( 1.03 , 1.312 ) , there are orbits of periods 2, 4, and 8. When δ = 1.35 , we can see chaotic sets in Figure 4(f).

Figure 3 
               Flip bifurcation diagram of system (1.3) with parameters given in (4.2) and the initial value 
                     
                        
                        
                           
                              (
                              
                                 0.3
                                 ,
                                 0.4
                              
                              )
                           
                        
                        \left(0.3,0.4)
                     
                  . (a) In the 
                     
                        
                        
                           
                              (
                              
                                 δ
                                 ,
                                 x
                              
                              )
                           
                        
                        \left(\delta ,x)
                     
                   plane and (b) in the 
                     
                        
                        
                           
                              (
                              
                                 δ
                                 ,
                                 y
                              
                              )
                           
                        
                        \left(\delta ,y)
                     
                   plane.
Figure 3

Flip bifurcation diagram of system (1.3) with parameters given in (4.2) and the initial value ( 0.3 , 0.4 ) . (a) In the ( δ , x ) plane and (b) in the ( δ , y ) plane.

Figure 4 
               Phase portraits for various values of 
                     
                        
                        
                           δ
                        
                        \delta 
                     
                   corresponding to Figure 3. (a) 
                     
                        
                        
                           δ
                           =
                           0.9
                        
                        \delta =0.9
                     
                  , (b) 
                     
                        
                        
                           δ
                           =
                           1.03
                        
                        \delta =1.03
                     
                  , (c) 
                     
                        
                        
                           δ
                           =
                           1.037
                        
                        \delta =1.037
                     
                  , (d) 
                     
                        
                        
                           δ
                           =
                           1.26
                        
                        \delta =1.26
                     
                  , (e) 
                     
                        
                        
                           δ
                           =
                           1.31
                        
                        \delta =1.31
                     
                  , and (f) 
                     
                        
                        
                           δ
                           =
                           1.35
                        
                        \delta =1.35
                     
                  .
Figure 4

Phase portraits for various values of δ corresponding to Figure 3. (a) δ = 0.9 , (b) δ = 1.03 , (c) δ = 1.037 , (d) δ = 1.26 , (e) δ = 1.31 , and (f) δ = 1.35 .

Example 4.4

(Neimark-Sacker bifurcation at E 2 ) We choose the bifurcation parameters and the other parameters as follows:

(4.3) δ [ 0 , 4.5 ] , b = 2 , c = 0.2 , e = 0.18 , p = 0.0859 .

Obviously, e < c . It follows from Theorem 2.2(iv) and (2.8) that system (1.3) has only one positive fixed point E 2 ( 0.17 , 0.17 ) . By simple calculation, we can obtain that the bifurcation parameter is δ = 1.36 , and the eigenvalues of J ( E 2 ) are 0.98 ± 0.21 i . Moreover, we have l = 0.017 > 0 , α = 0.21 , and ( b , c , e , p , δ ) = ( 2 , 0.2 , 0.18 , 2 , 0.0859 ) H B . By Theorem 3.3, Neimark-Sacker bifurcation can undergo at E 2 with δ = 1.36 , which is demonstrated in Figure 5. We can observe from Figure 5(a) and (c) that the fixed point E 2 of system (1.3) goes from stable to chaotic as δ increases. The phase portraits associated with Figure 5(a) and (c) are displayed in Figure 6, which clearly shows how a smooth invariant circle bifurcates from the stable fixed point E 2 ( 0.17 , 0.17 ) . From Figure 6, we see that there are orbits of period-16, 21, 35, 49, 54, invariant cycles, and chaotic sets.

Figure 5 
               Neimark-Sacker bifurcation of (1.3) at 
                     
                        
                        
                           
                              
                                 E
                              
                              
                                 2
                              
                              
                                 ∗
                              
                           
                        
                        {E}_{2}^{\ast }
                     
                  , where the values of parameters are give in (4.3), and the initial value is 
                     
                        
                        
                           
                              (
                              
                                 0.15
                                 ,
                                 0.1
                              
                              )
                           
                        
                        \left(0.15,0.1)
                     
                  . Note that the local amplifications are for 
                     
                        
                        
                           δ
                           ∈
                           
                              [
                              
                                 4
                                 ,
                                 4.38
                              
                              ]
                           
                        
                        \delta \in \left[4,4.38]
                     
                  . (a) In 
                     
                        
                        
                           δ
                           ,
                           x
                        
                        \delta ,x
                     
                   plane, (b) local amplification of (a), (c) in 
                     
                        
                        
                           
                              (
                              
                                 δ
                                 ,
                                 y
                              
                              )
                           
                        
                        \left(\delta ,y)
                     
                   plane, and (d) local amplification of (c).
Figure 5

Neimark-Sacker bifurcation of (1.3) at E 2 , where the values of parameters are give in (4.3), and the initial value is ( 0.15 , 0.1 ) . Note that the local amplifications are for δ [ 4 , 4.38 ] . (a) In δ , x plane, (b) local amplification of (a), (c) in ( δ , y ) plane, and (d) local amplification of (c).

Figure 6 
               Phase portraits for various values of 
                     
                        
                        
                           δ
                        
                        \delta 
                     
                   corresponding to Figure 5. (a) 
                     
                        
                        
                           δ
                           =
                           1.24
                        
                        \delta =1.24
                     
                  , (b) 
                     
                        
                        
                           δ
                           =
                           1.34
                        
                        \delta =1.34
                     
                  , (c) 
                     
                        
                        
                           δ
                           =
                           1.48
                        
                        \delta =1.48
                     
                  , (d) 
                     
                        
                        
                           δ
                           =
                           1.9436
                        
                        \delta =1.9436
                     
                  , (e) 
                     
                        
                        
                           δ
                           =
                           2.30
                        
                        \delta =2.30
                     
                  , (f) 
                     
                        
                        
                           δ
                           =
                           2.4121
                        
                        \delta =2.4121
                     
                  , (g) 
                     
                        
                        
                           δ
                           =
                           2.6339
                        
                        \delta =2.6339
                     
                  , (h) 
                     
                        
                        
                           δ
                           =
                           2.7073
                        
                        \delta =2.7073
                     
                  , (i) 
                     
                        
                        
                           δ
                           =
                           3.32
                        
                        \delta =3.32
                     
                  , (j) 
                     
                        
                        
                           δ
                           =
                           3.40
                        
                        \delta =3.40
                     
                  , (k) 
                     
                        
                        
                           δ
                           =
                           3.818
                        
                        \delta =3.818
                     
                  , (l) 
                     
                        
                        
                           δ
                           =
                           3.82
                        
                        \delta =3.82
                     
                  , (m) 
                     
                        
                        
                           δ
                           =
                           4.12
                        
                        \delta =4.12
                     
                  , (n) 
                     
                        
                        
                           δ
                           =
                           4.24
                        
                        \delta =4.24
                     
                  , (o) 
                     
                        
                        
                           δ
                           =
                           4.35
                        
                        \delta =4.35
                     
                  .
Figure 6

Phase portraits for various values of δ corresponding to Figure 5. (a) δ = 1.24 , (b) δ = 1.34 , (c) δ = 1.48 , (d) δ = 1.9436 , (e) δ = 2.30 , (f) δ = 2.4121 , (g) δ = 2.6339 , (h) δ = 2.7073 , (i) δ = 3.32 , (j) δ = 3.40 , (k) δ = 3.818 , (l) δ = 3.82 , (m) δ = 4.12 , (n) δ = 4.24 , (o) δ = 4.35 .

5 Conclusion

In this paper, a discrete Leslie-Gower predator-prey system with Michaelis-Menten type harvesting has been proposed and analyzed. Conditions on the existence and stability of fixed points are given. We obtained that all the boundary fixed points are unstable when they exist, and while for positive fixed points, only E 2 is locally stable within the appropriate parameter range, and the others are all unstable. It is proven that the fixed points of system (1.3) possess fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation by using the center manifold theorem and bifurcation theory. Numerical simulations are performed to support the main theoretic results. Example 4.1 illustrates that the number of fixed points changes from two to one to none as the bifurcation parameter e changes through fold bifurcation (see Figure 1). Example 4.2 indicates that flip bifurcation occurs at the boundary fixed point E 2 , and in other words, the predator becomes extinct while the prey experiences the flip bifurcation from stable to chaotic (see Figure 2). Examples 4.3 and 4.4 from numerical simulations show that choosing the integral step size δ as the bifurcation parameter, system (1.3) undergoes a flip bifurcation at the positive fixed point E 2 , which includes orbits of period-2, 4, 8 (see Figure 4), and a Neimark-Sacker bifurcation, which includes orbits of period-16, 21, 35, 49, 54, invariant cycles, and chaotic sets (see Figure 6). These results reveal that the integral step size δ plays a vital role in the stability of the system (1.3). It reminds us that it is necessary to clarify the integral step size assumed in advance when dealing with numerical solutions or approximate solutions of the original continuous predator-prey system with the Holling and Leslie type functional responses.

The continuous version of this system has been investigated by Gupta et al. [46]. They explored the local stability, saddle-node bifurcation, limit cycle, and Hopf bifurcation. In our work, we exhibited various bifurcations specific to discrete systems, including fold bifurcation, flip bifurcation, and Neimark-Sacker bifurcation. This tells us that the discrete system has much richer dynamical behaviors than the continuous analog has.

At the end of the paper, we would like to mention that a referee pointed out it may be a good idea to study the dynamic behaviors of the discrete fast-slow system, and recently, we read the papers [51,52, 53,54] about the continuous slow-fast system; however, it seems that we could not do some research immediately on this direction, and we will leave this for future study.

Acknowledgements

The authors would like to thank anonymous referees for carefully reading the paper, and for bring our attention to the slow-fast system.

  1. Funding information: This work was supported by the Natural Science Foundation of Fujian Province (2020J01499).

  2. Author contributions: All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

  3. Conflict of interest: The authors declare that there is no conflict of interests.

  4. Data availability statement: Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

References

[1] D. P. Hu and H. J. Cao, Stability and bifurcation analysis in a predator-prey system with Michaelis-Menten type predator harvesting, Nonlinear Anal. Real World Appl. 33 (2017), 58–82. 10.1016/j.nonrwa.2016.05.010Search in Google Scholar

[2] C. R. Zhu and K. Lei, Bifurcations analysis of Leslie-Gower predator-prey models with nonlinear predator-harvesting, Discrete Contin. Dyn. Syst. Ser. S 10 (2017), 1187–1206. 10.3934/dcdss.2017065Search in Google Scholar

[3] S. B. Yu, Effect of predator mutual interference on an autonomous Leslie-Gower predator-prey model, IAENG Int. J. Appl. Math. 49 (2019), 229–233. Search in Google Scholar

[4] Q. Yue, Dynamics of a modified Leslie-Gower predator-prey model with Holling-type II schemes and a prey refuge, SpringerPlus. 5 (2016), no. 1, 1–12. 10.1186/s40064-016-2087-7Search in Google Scholar PubMed PubMed Central

[5] R. M. May, J. R. Beddington, C. W. Clark, S. J. Holt, and R. M. Laws, Management of multispecies fisheries, Science 205 (1979), no. 4403, 267–277. 10.1126/science.205.4403.267Search in Google Scholar PubMed

[6] C. W. Clark and M. Mangel, Aggregation and fishery dynamics: a theoretical study of schooling and the purse seine tuna fisheries, Fish. Bull. 77 (1979), no. 2, 317–337. Search in Google Scholar

[7] N. Zhang, Y. Kao, F. Chen, B. Xie, and S. Li, On a predator-prey system interaction under fluctuating water level with nonselective harvesting, Open Math. 18 (2020), no. 1, 458–475. 10.1515/math-2020-0145Search in Google Scholar

[8] Z. Zhu, F. Chen, L. Lai, and Z. Li, Dynamic behaviors of a discrete May type cooperative system incorporating Michaelis-Menten type harvesting, IAENG Int. J. Appl. Math. 50 (2020), 1–10. Search in Google Scholar

[9] X. D. Xie, F. D. Chen, and Y. L. Xue, Note on the stability property of a cooperative system incorporating harvesting, Discrete Dyn. Nat. Soc. 2014 (2014), Article ID 327823. 10.1155/2014/327823Search in Google Scholar

[10] F. D. Chen, H. L. Wu, and X. D. Xie, Global attractivity of a discrete cooperative system incorporating harvesting, Adv. Differential Equations 2016 (2016), no. 1, 268. 10.1186/s13662-016-0996-ySearch in Google Scholar

[11] F. Chen, Q. Zhou, and S. Lin, Global stability of symbiotic medel of commensalism and parasitism with harvesting in commensal populations, WSEAS Trans. Math. 21 (2022), 424–432. 10.37394/23206.2022.21.50Search in Google Scholar

[12] Q. Zhou, S. Lin, F. Chen, and R. U. Wu, Positive periodic solution of a discrete Lotka-Volterra commensal symbiosis model with Michaelis-Menten type harvesting, WSEAS Trans. Math. 21 (2022), 515–523. 10.37394/23206.2022.21.57Search in Google Scholar

[13] J. Chen, Y. Chen, Z. Zhu, and F. Chen, Stability and bifurcation of a discrete predator-prey system with Allee effect and other food resource for the predators, J. Appl. Math. Comput. 2022 (2022), 1–10. 10.1007/s12190-022-01764-5Search in Google Scholar

[14] X. Yu, Z. Zhu, and F. Chen, Dynamic behaviors of a single species stage structure model with Michaelis-Menten-type juvenile population harvesting, Math. 8 (2020), 1281. 10.3390/math8081281Search in Google Scholar

[15] Q. F. Lin, Dynamic behaviors of a commensal symbiosis model with non-monotonic functional response and non-selective harvesting in a partial closure, Commun. Math. Biol. Neurosci. 2018 (2018), Article ID 4. Search in Google Scholar

[16] Y. Liu, L. Zhao, X. Y. Huang, and H. Deng, Stability and bifurcation analysis of two species amensalism model with Michaelis-Menten type harvesting and a cover for the first species, Adv. Differential Equations 2018 (2018), no. 1, 295. 10.1186/s13662-018-1752-2Search in Google Scholar

[17] X. Yu, Z. Zhu, L. Lai, and F. Chen, Stability and bifurcation analysis in a single-species stage structure system with Michaelis-Menten-type harvesting, Adv. Differential Equations 2020 (2020), 238. 10.1186/s13662-020-02652-7Search in Google Scholar

[18] B. G. Chen, The influence of commensalism on a Lotka-Volterra commensal symbiosis model with Michaelis-Menten type harvesting, Adv. Differential Equations 2019 (2019), no. 1, 43. 10.1186/s13662-019-1989-4Search in Google Scholar

[19] X. Y. Huang, F. D. Chen, X. D. Xie, and L. Zhao, Extinction of a two species competitive stage-structured system with the effect of toxic substance and harvesting, Open Math. 17 (2019), no. 1, 856–873. 10.1515/math-2019-0067Search in Google Scholar

[20] Y. Liu, X. Y. Guan, X. D. Xie, and Q. F. Lin, On the existence and stability of positive periodic solution of a nonautonomous commensal symbiosis model with Michaelis-Menten type harvesting, Commun. Math. Biol. Neurosci. 2019 (2019), Article ID 2. Search in Google Scholar

[21] C. Liu, Q. L. Zhang, Y. Zhang, and X. D. Duan, Bifurcation and control in a differential-algebraic harvested prey-predator model with stage structure for predator, Int. J. Bifurcation and Chaos 18 (2008), no. 10, 3159–3168. 10.1142/S0218127408022329Search in Google Scholar

[22] G. J. Peng, Y. L. Jiang, and C. P. Li, Bifurcations of a Holling-type II predator-prey system with constant rate harvesting, Int. J. Bifurcation and Chaos 19 (2009), no. 8, 2499–2514. 10.1142/S021812740902427XSearch in Google Scholar

[23] R. K. Roy and B. Roy, Analysis of prey-predator three species fishery model with harvesting including prey refuge and migration, Int. J. Bifurcation and Chaos, 26 (2016), no. 2, 1650022. 10.1142/S021812741650022XSearch in Google Scholar

[24] M. M. Haque and S. Sarwardi, Dynamics of a harvested prey-predator model with prey refuge dependent on both species, Int. J. Bifurcation and Chaos, 28 (2018), no. 12, 1830040. 10.1142/S0218127418300409Search in Google Scholar

[25] C. X. Huang, H. Zhang, J. D. Cao, and H. J. Hu, Stability and hopf bifurcation of a delayed prey-predator model with disease in the predator, Int. J. Bifurcation and Chaos 29 (2019), no. 7, 1950091. 10.1142/S0218127419500913Search in Google Scholar

[26] X. L. Liu and D. M. Xiao, Bifurcations in a discrete time Lotka-Volterra predator-prey system, Discrete Contin. Dyn. Syst. Ser. B 6 (2006), no. 3, 559. 10.3934/dcdsb.2006.6.559Search in Google Scholar

[27] X. L. Liu and D. M. Xiao, Complex dynamic behaviors of a discrete-time predator-prey system, Chaos Solitons Fractals 32 (2007), no. 1, 80–94. 10.1016/j.chaos.2005.10.081Search in Google Scholar

[28] Z. M. He and X. Lai, Bifurcation and chaotic behavior of a discrete-time predator-prey system, Nonlinear Anal. Real World Appl. 12 (2011), no. 1, 403–417. 10.1016/j.nonrwa.2010.06.026Search in Google Scholar

[29] B. G. Chen, Global attractivity of a discrete competition model, Adv. Differential Equations 2016 (2016), 273. 10.1186/s13662-016-1000-6Search in Google Scholar

[30] L. F. Cheng and H. J. Cao, Bifurcation analysis of a discrete-time ratio-dependent predator-prey model with allee effect, Commun. Nonlinear Sci. Numer. Simul. 38 (2016), 288–302. 10.1016/j.cnsns.2016.02.038Search in Google Scholar

[31] Q. Q. Cui, Q. Zhang, Z. P. Qiu, and Z. Y. Hu, Complex dynamics of a discrete-time predator-prey system with Holling IV functional response, Chaos Solitons Fractals 87 (2016), 158–171. 10.1016/j.chaos.2016.04.002Search in Google Scholar

[32] T. S. Huang and H. Y. Zhang, Bifurcation, chaos and pattern formation in a space-and time-discrete predator-prey system, Chaos Solitons Fractals 91 (2016), 92–107. 10.1016/j.chaos.2016.05.009Search in Google Scholar

[33] S. M. Salman, A. M. Yousef, and A. A. Elsadany, Stability, bifurcation analysis and chaos control of a discrete predator-prey system with square root functional response, Chaos Solitons Fractals, 93 (2016), 20–31. 10.1016/j.chaos.2016.09.020Search in Google Scholar

[34] R. Banerjee, P. Das, and D. Mukherjee, Stability and permanence of a discrete-time two-prey one-predator system with Holling type-III functional response, Chaos Solitons Fractals 117 (2018), 240–248. 10.1016/j.chaos.2018.10.032Search in Google Scholar

[35] J. C. Huang, S. H. Liu, S. G. Ruan, and D. M. Xiao, Bifurcations in a discrete predator-prey model with nonmonotonic functional response, J. Math. Anal. Appl. 464 (2018), no. 1, 201–230. 10.1016/j.jmaa.2018.03.074Search in Google Scholar

[36] J. L. Zhao and Y. Yan, Stability and bifurcation analysis of a discrete predator-prey system with modified Holling-Tanner functional response, Adv. Differential Equations 2018 (2018), no. 1, 1–18. 10.1186/s13662-018-1819-0Search in Google Scholar

[37] S. S. Rana, Bifurcations and chaos control in a discrete-time predator-prey system of Leslie type, J. Appl. Anal. Comput. 9 (2019), no. 1, 31–44. 10.11948/2019.31Search in Google Scholar

[38] P. K. Santra, G. S. Mahapatra, and G. R. Phaijoo, Bifurcation and chaos of a discrete predator-prey model with Crowley-Martin functional response incorporating proportional prey refuge, Math. Probl. Eng. 2020 (2020), Article ID 5309814. 10.1155/2020/5309814Search in Google Scholar

[39] A. Singh and P. Deolia, Dynamical analysis and chaos control in discrete-time prey-predator model, Commun. Nonlinear Sci. Numer. Simul. 90 (2020), 105313. 10.1016/j.cnsns.2020.105313Search in Google Scholar

[40] J. L. Chen, X. Q. He, and F. D. Chen, The influence of fear effect to a discrete-time predator-prey system with predator has other food resource, Mathematics, 9 (2021), no. 8, 865. 10.3390/math9080865Search in Google Scholar

[41] J. L. Ren, L. P. Yu, and S. Siegmund, Bifurcations and chaos in a discrete predator-prey model with Crowley-Martin functional response, Nonlinear Dynam. 90 (2017), no. 1, 19–41. 10.1007/s11071-017-3643-6Search in Google Scholar

[42] D. P. Hu and H. J. Cao, Bifurcation and chaos in a discrete-time predator-prey system of Holling and Leslie type, Commun. Nonlinear Sci. Numer. Simul. 22 (2015), 702–715. 10.1016/j.cnsns.2014.09.010Search in Google Scholar

[43] M. B. Ajaz, U. Saeed, Q. Din, I. Ali, and M. I. Siddiqui, Bifurcation analysis and chaos control in discrete-time modified Leslie-Gower prey harvesting model, Adv. Differential Equations 2020 (2020), 45. 10.1186/s13662-020-2498-1Search in Google Scholar

[44] Z. L. Zhu, F. D. Chen, L. Y. Lai, and Z. Li, Dynamic behaviors of a discrete may type cooperative system incorporating Michaelis-Menten type harvesting, Int. J. Appl. Math. 50 (2020), no. 3, 1–10. Search in Google Scholar

[45] A. Singh and P. Malik, Bifurcations in a modified Leslie-Gower predator-prey discrete model with Michaelis-Menten prey harvesting, J. Appl. Math. Comput. 67 (2021), no. 1, 1–32. 10.1007/s12190-020-01491-9Search in Google Scholar

[46] R. P. Gupta, M. Banerjee, and P. Chandra, Bifurcation analysis and control of Leslie-Gower predator-prey model with Michaelis-Menten type prey-harvesting, Differ. Equ. Dyn. Sys. 20 (2012), no. 3, 339–366. 10.1007/s12591-012-0142-6Search in Google Scholar

[47] G. Y. Chen and Z. D. Teng, On the stability in a discrete two-species competition system, J. Appl. Math. Comput. 38 (2012), no. 1, 25–39. 10.1007/s12190-010-0460-1Search in Google Scholar

[48] D. C. Liaw, Application of center manifold reduction to nonlinear system stabilization, Appl. Math. Comput. 91 (1998), 243–258. 10.1016/S0096-3003(97)10021-2Search in Google Scholar

[49] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer Science and Business Media, Berlin, Germany, 2003. Search in Google Scholar

[50] C. Robinson, Dynamical Systems: Stability, Symbolic Dynamics, and Chaos, CRC Press, 1998. 10.1201/9781482227871Search in Google Scholar

[51] L. Zhao and J. Shen, Relaxation oscillations in a slow-fast predator-prey model with weak Allee effect and Holling-IV functional response, Commun. Nonlinear. Sci. Numer. Simul. 112 (2022), 106517. 10.1016/j.cnsns.2022.106517Search in Google Scholar

[52] C. Li and H. Zhu, Canard cycles for predator-prey systems with Holling types of functional response, J. Differential Equations 254 (2013), no. 2, 879–910. 10.1016/j.jde.2012.10.003Search in Google Scholar

[53] R. Huzak, Predator-prey systems with small predatoras death rate, Electronic J. Qual. Theory Differ. Equ. 86 (2018), 1–16. 10.14232/ejqtde.2018.1.86Search in Google Scholar

[54] M. Krupa and P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points-fold and canard points in two dimensions, SIAM J. Math. Anal. 33 (2001), no. 2, 286–314. 10.1137/S0036141099360919Search in Google Scholar

Received: 2021-08-21
Revised: 2022-06-30
Accepted: 2022-07-01
Published Online: 2022-08-18

© 2022 Jialin Chen et al., published by De Gruyter

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

Articles in the same Issue

  1. Regular Articles
  2. A random von Neumann theorem for uniformly distributed sequences of partitions
  3. Note on structural properties of graphs
  4. Mean-field formulation for mean-variance asset-liability management with cash flow under an uncertain exit time
  5. The family of random attractors for nonautonomous stochastic higher-order Kirchhoff equations with variable coefficients
  6. The intersection graph of graded submodules of a graded module
  7. Isoperimetric and Brunn-Minkowski inequalities for the (p, q)-mixed geominimal surface areas
  8. On second-order fuzzy discrete population model
  9. On certain functional equation in prime rings
  10. General complex Lp projection bodies and complex Lp mixed projection bodies
  11. Some results on the total proper k-connection number
  12. The stability with general decay rate of hybrid stochastic fractional differential equations driven by Lévy noise with impulsive effects
  13. Well posedness of magnetohydrodynamic equations in 3D mixed-norm Lebesgue space
  14. Strong convergence of a self-adaptive inertial Tseng's extragradient method for pseudomonotone variational inequalities and fixed point problems
  15. Generic uniqueness of saddle point for two-person zero-sum differential games
  16. Relational representations of algebraic lattices and their applications
  17. Explicit construction of mock modular forms from weakly holomorphic Hecke eigenforms
  18. The equivalent condition of G-asymptotic tracking property and G-Lipschitz tracking property
  19. Arithmetic convolution sums derived from eta quotients related to divisors of 6
  20. Dynamical behaviors of a k-order fuzzy difference equation
  21. The transfer ideal under the action of orthogonal group in modular case
  22. The multinomial convolution sum of a generalized divisor function
  23. Extensions of Gronwall-Bellman type integral inequalities with two independent variables
  24. Unicity of meromorphic functions concerning differences and small functions
  25. Solutions to problems about potentially Ks,t-bigraphic pair
  26. Monotonicity of solutions for fractional p-equations with a gradient term
  27. Data smoothing with applications to edge detection
  28. An ℋ-tensor-based criteria for testing the positive definiteness of multivariate homogeneous forms
  29. Characterizations of *-antiderivable mappings on operator algebras
  30. Initial-boundary value problem of fifth-order Korteweg-de Vries equation posed on half line with nonlinear boundary values
  31. On a more accurate half-discrete Hilbert-type inequality involving hyperbolic functions
  32. On split twisted inner derivation triple systems with no restrictions on their 0-root spaces
  33. Geometry of conformal η-Ricci solitons and conformal η-Ricci almost solitons on paracontact geometry
  34. Bifurcation and chaos in a discrete predator-prey system of Leslie type with Michaelis-Menten prey harvesting
  35. A posteriori error estimates of characteristic mixed finite elements for convection-diffusion control problems
  36. Dynamical analysis of a Lotka Volterra commensalism model with additive Allee effect
  37. An efficient finite element method based on dimension reduction scheme for a fourth-order Steklov eigenvalue problem
  38. Connectivity with respect to α-discrete closure operators
  39. Khasminskii-type theorem for a class of stochastic functional differential equations
  40. On some new Hermite-Hadamard and Ostrowski type inequalities for s-convex functions in (p, q)-calculus with applications
  41. New properties for the Ramanujan R-function
  42. Shooting method in the application of boundary value problems for differential equations with sign-changing weight function
  43. Ground state solution for some new Kirchhoff-type equations with Hartree-type nonlinearities and critical or supercritical growth
  44. Existence and uniqueness of solutions for the stochastic Volterra-Levin equation with variable delays
  45. Ambrosetti-Prodi-type results for a class of difference equations with nonlinearities indefinite in sign
  46. Research of cooperation strategy of government-enterprise digital transformation based on differential game
  47. Malmquist-type theorems on some complex differential-difference equations
  48. Disjoint diskcyclicity of weighted shifts
  49. Construction of special soliton solutions to the stochastic Riccati equation
  50. Remarks on the generalized interpolative contractions and some fixed-point theorems with application
  51. Analysis of a deteriorating system with delayed repair and unreliable repair equipment
  52. On the critical fractional Schrödinger-Kirchhoff-Poisson equations with electromagnetic fields
  53. The exact solutions of generalized Davey-Stewartson equations with arbitrary power nonlinearities using the dynamical system and the first integral methods
  54. Regularity of models associated with Markov jump processes
  55. Multiplicity solutions for a class of p-Laplacian fractional differential equations via variational methods
  56. Minimal period problem for second-order Hamiltonian systems with asymptotically linear nonlinearities
  57. Convergence rate of the modified Levenberg-Marquardt method under Hölderian local error bound
  58. Non-binary quantum codes from constacyclic codes over 𝔽q[u1, u2,…,uk]/⟨ui3 = ui, uiuj = ujui
  59. On the general position number of two classes of graphs
  60. A posteriori regularization method for the two-dimensional inverse heat conduction problem
  61. Orbital stability and Zhukovskiǐ quasi-stability in impulsive dynamical systems
  62. Approximations related to the complete p-elliptic integrals
  63. A note on commutators of strongly singular Calderón-Zygmund operators
  64. Generalized Munn rings
  65. Double domination in maximal outerplanar graphs
  66. Existence and uniqueness of solutions to the norm minimum problem on digraphs
  67. On the p-integrable trajectories of the nonlinear control system described by the Urysohn-type integral equation
  68. Robust estimation for varying coefficient partially functional linear regression models based on exponential squared loss function
  69. Hessian equations of Krylov type on compact Hermitian manifolds
  70. Class fields generated by coordinates of elliptic curves
  71. The lattice of (2, 1)-congruences on a left restriction semigroup
  72. A numerical solution of problem for essentially loaded differential equations with an integro-multipoint condition
  73. On stochastic accelerated gradient with convergence rate
  74. Displacement structure of the DMP inverse
  75. Dependence of eigenvalues of Sturm-Liouville problems on time scales with eigenparameter-dependent boundary conditions
  76. Existence of positive solutions of discrete third-order three-point BVP with sign-changing Green's function
  77. Some new fixed point theorems for nonexpansive-type mappings in geodesic spaces
  78. Generalized 4-connectivity of hierarchical star networks
  79. Spectra and reticulation of semihoops
  80. Stein-Weiss inequality for local mixed radial-angular Morrey spaces
  81. Eigenvalues of transition weight matrix for a family of weighted networks
  82. A modified Tikhonov regularization for unknown source in space fractional diffusion equation
  83. Modular forms of half-integral weight on Γ0(4) with few nonvanishing coefficients modulo
  84. Some estimates for commutators of bilinear pseudo-differential operators
  85. Extension of isometries in real Hilbert spaces
  86. Existence of positive periodic solutions for first-order nonlinear differential equations with multiple time-varying delays
  87. B-Fredholm elements in primitive C*-algebras
  88. Unique solvability for an inverse problem of a nonlinear parabolic PDE with nonlocal integral overdetermination condition
  89. An algebraic semigroup method for discovering maximal frequent itemsets
  90. Class-preserving Coleman automorphisms of some classes of finite groups
  91. Exponential stability of traveling waves for a nonlocal dispersal SIR model with delay
  92. Existence and multiplicity of solutions for second-order Dirichlet problems with nonlinear impulses
  93. The transitivity of primary conjugacy in regular ω-semigroups
  94. Stability estimation of some Markov controlled processes
  95. On nonnil-coherent modules and nonnil-Noetherian modules
  96. N-Tuples of weighted noncommutative Orlicz space and some geometrical properties
  97. The dimension-free estimate for the truncated maximal operator
  98. A human error risk priority number calculation methodology using fuzzy and TOPSIS grey
  99. Compact mappings and s-mappings at subsets
  100. The structural properties of the Gompertz-two-parameter-Lindley distribution and associated inference
  101. A monotone iteration for a nonlinear Euler-Bernoulli beam equation with indefinite weight and Neumann boundary conditions
  102. Delta waves of the isentropic relativistic Euler system coupled with an advection equation for Chaplygin gas
  103. Multiplicity and minimality of periodic solutions to fourth-order super-quadratic difference systems
  104. On the reciprocal sum of the fourth power of Fibonacci numbers
  105. Averaging principle for two-time-scale stochastic differential equations with correlated noise
  106. Phragmén-Lindelöf alternative results and structural stability for Brinkman fluid in porous media in a semi-infinite cylinder
  107. Study on r-truncated degenerate Stirling numbers of the second kind
  108. On 7-valent symmetric graphs of order 2pq and 11-valent symmetric graphs of order 4pq
  109. Some new characterizations of finite p-nilpotent groups
  110. A Billingsley type theorem for Bowen topological entropy of nonautonomous dynamical systems
  111. F4 and PSp (8, ℂ)-Higgs pairs understood as fixed points of the moduli space of E6-Higgs bundles over a compact Riemann surface
  112. On modules related to McCoy modules
  113. On generalized extragradient implicit method for systems of variational inequalities with constraints of variational inclusion and fixed point problems
  114. Solvability for a nonlocal dispersal model governed by time and space integrals
  115. Finite groups whose maximal subgroups of even order are MSN-groups
  116. Symmetric results of a Hénon-type elliptic system with coupled linear part
  117. On the connection between Sp-almost periodic functions defined on time scales and ℝ
  118. On a class of Harada rings
  119. On regular subgroup functors of finite groups
  120. Fast iterative solutions of Riccati and Lyapunov equations
  121. Weak measure expansivity of C2 dynamics
  122. Admissible congruences on type B semigroups
  123. Generalized fractional Hermite-Hadamard type inclusions for co-ordinated convex interval-valued functions
  124. Inverse eigenvalue problems for rank one perturbations of the Sturm-Liouville operator
  125. Data transmission mechanism of vehicle networking based on fuzzy comprehensive evaluation
  126. Dual uniformities in function spaces over uniform continuity
  127. Review Article
  128. On Hahn-Banach theorem and some of its applications
  129. Rapid Communication
  130. Discussion of foundation of mathematics and quantum theory
  131. Special Issue on Boundary Value Problems and their Applications on Biosciences and Engineering (Part II)
  132. A study of minimax shrinkage estimators dominating the James-Stein estimator under the balanced loss function
  133. Representations by degenerate Daehee polynomials
  134. Multilevel MC method for weak approximation of stochastic differential equation with the exact coupling scheme
  135. Multiple periodic solutions for discrete boundary value problem involving the mean curvature operator
  136. Special Issue on Evolution Equations, Theory and Applications (Part II)
  137. Coupled measure of noncompactness and functional integral equations
  138. Existence results for neutral evolution equations with nonlocal conditions and delay via fractional operator
  139. Global weak solution of 3D-NSE with exponential damping
  140. Special Issue on Fractional Problems with Variable-Order or Variable Exponents (Part I)
  141. Ground state solutions of nonlinear Schrödinger equations involving the fractional p-Laplacian and potential wells
  142. A class of p1(x, ⋅) & p2(x, ⋅)-fractional Kirchhoff-type problem with variable s(x, ⋅)-order and without the Ambrosetti-Rabinowitz condition in ℝN
  143. Jensen-type inequalities for m-convex functions
  144. Special Issue on Problems, Methods and Applications of Nonlinear Analysis (Part III)
  145. The influence of the noise on the exact solutions of a Kuramoto-Sivashinsky equation
  146. Basic inequalities for statistical submanifolds in Golden-like statistical manifolds
  147. Global existence and blow up of the solution for nonlinear Klein-Gordon equation with variable coefficient nonlinear source term
  148. Hopf bifurcation and Turing instability in a diffusive predator-prey model with hunting cooperation
  149. Efficient fixed-point iteration for generalized nonexpansive mappings and its stability in Banach spaces
Downloaded on 7.12.2025 from https://www.degruyterbrill.com/document/doi/10.1515/math-2022-0054/html
Scroll to top button