Dynamical Behavior of Discrete Ratio-Dependent Predator-Prey System

Mingxia Duan¹, Jiying Ma²

¹Department of Mathematics, University of Shanghai for Science and Technology, Shanghai, China.
²College of Science, University of Shanghai for Science and Technology, Shanghai, China.

DOI: 10.4236/ojapps.2023.133032                

Abstract

In this paper, we propose a discrete ratio-dependent predator-prey system. The stability of the fixed points of this model is studied. At the same time, it is shown that the discrete model undergoes fold bifurcation and flip bifurcation by using bifurcation theory and the method of approximation by a flow. Numerical simulations are presented not only to demonstrate the consistence with our theoretical analyses, but also to exhibit the complex dynamical behaviors, such as the cascade of period-doubling bifurcation in period-2 and the chaotic sets. The Maximum Lyapunov exponents are numerically computed to confirm further the complexity of the dynamical behaviors. These results show that the direct discrete method has more rich dynamic behaviors than the discrete model obtained by Euler method.

Keywords

Discrete-Time Predator-Prey System, Ratio-Dependent Functional Response, Stability and Bifurcation Analysis

1. Introduction and Preliminaries

The interaction between predators and their prey has long been the studied topics of ecology and mathematical ecology. In recent years, discrete dynamic models have been widely studied. Compared with continuous dynamic systems, discrete dynamic systems have the following three advantages: Firstly, when populations have non-overlapping generations or the population number is small, discrete dynamic systems are more appropriate than continuous dynamic systems. Secondly, the numerical simulation results obtained from the discrete dynamic system are more accurate than the continuous dynamic system. Thirdly, discrete models exhibit more complex dynamic characteristics than corresponding continuous models [1]. For example, the single-species discrete Logistic model has more complex dynamic behaviors such as bifurcation and chaos.

There are two ways to get discrete system. On the one hand, we start from the continuous model and transform the continuous model into a discrete model by some method, such as Euler iteration [2] and semi-discrete methods [3]. On the other hand, we start from the discrete model directly. Elettreby [4] [5] studied the discrete prey-predator model with mixed functional response from the following discrete model.

$\{\begin{array}{l}x\left(n+1\right)=x\left(n\right)\left(a-bx\left(n\right)\right)-p\left(x\left(n\right)\right)y\left(n\right),\\ y\left(n+1\right)=q\left(x\left(n\right)\right)y\left(n\right)-dy\left(n\right),\end{array}$ (1.1)

where $x\left(n\right)$ and $y\left(n\right)$ represent population densities of prey and predator at a discrete time step n, respectively; d is the natural death rate of the predator; $p\left(x\left(n\right)\right)y\left(n\right)$ and $q\left(x\left(n\right)\right)y\left(n\right)$ represent the predator-prey interaction, $p\left(x\left(n\right)\right)$ is the predator functional responses, $q\left(x\left(n\right)\right)$ describes how predator transforms the consumed prey into the growth of predator. Because the function $p\left(x\left(n\right)\right)$ only depends on the density of prey, it is often called a prey-dependent response function. However, more and more evidence shows that in many situations, especially when predators have to search for food, a more appropriate functional response is a function which is dependent on both prey and predator densities.

Differing from the prey-dependent predator-prey systems, the ratio-dependent predator-prey models have two main predictions: a) the equilibrium abundance is positively correlated with the enrichment gradient (see Arditi [6] ) and b) the “paradox of enrichment” (see Rosenzweig [7] ) either completely disappears or enrichment is linked to stability in a more complex way. In this paper, we will study some mathematical characteristic rather than discuss the ecological significance of the model. In 2001, Xiao [8] [9] [10] further studied the ratio-dependent predator-prey model based on Kuang [11]. Xiao made a global qualitative analysis of the model depending on all parameters, and gave the conditions for the existence and nonexistence of limit cycles in the model. Later, the discrete-time dynamical models described by difference equations are widely investigated. Chen [3] showed that positive equilibrium is globally asymptotically stable and established a new sufficient condition. Sohel Rana [12] discretized the continuous model using the forward Euler scheme. It showed that under some parametric conditions, the system undergoes flip and Neimark-Sacker bifurcation.

Arditi [6] proposed the ratio-dependent function in the following form

$p\left(\frac{x}{y}\right)=\frac{c\frac{x}{y}}{m+\frac{x}{y}}=\frac{cx}{my+x}.$

Then we get the following discrete predator-prey model with the ratio-dependent:

$\{\begin{array}{l}x\left(n+1\right)=x\left(n\right)\left(a-bx\left(n\right)\right)-\frac{cx\left(n\right)y\left(n\right)}{x\left(n\right)+my\left(n\right)},\\ y\left(n+1\right)=y\left(n\right)\left(-d+\frac{fx\left(n\right)}{x\left(n\right)+my\left(n\right)}\right),\end{array}$ (1.2)

where $\frac{a}{b}>0$ is the carrying capacity of the prey, d is the natural death rate of the predator, a is prey intrinsic growth rate, c is capturing rate, m is half saturation constant and f is conversion rate. So a, b, c, d, m and f are positive constants.

The paper is organized as follows. In Section 1, as preliminaries, an important lemma is given. In Section 2, we devote to the existence and local stability of the fixed points of system (1.2). In Sections 3 - 4, we analyze the bifurcation of system (1.2) at the boundary equilibrium point by using the bifurcation theory. Moreover, we proved that the system (1.2) undergoes fold and flip bifurcations for some parameters, respectively. In Section 5, numerical simulations are performed to illustrate the obtained theoretical results and reveal some new dynamical properties of the system (1.2). In Section 6, some conclusions close the paper.

Before analyzing the fixed points of the system (1.2), we introduce an important lemma [2], which will be useful later.

Lemma 1.1 Let $F\left(\lambda \right)={\lambda }^2+B\lambda +C$ , where B and C are two real constants. Suppose ${\lambda }_1$ and ${\lambda }_2$ are two roots of $F\left(\lambda \right)=0$ . Then the following statements hold. If $F\left(1\right)>0$ , then

1) $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|<1$ if and only if $F\left(-1\right)>0$ and $C<1$ ;

2) $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|>1$ (or $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$ ) if and only if $F\left(-1\right)<0$ ;

3) $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|>1$ if and only if $F\left(-1\right)>0$ and $C>1$ ;

4) ${\lambda }_1=-1$ and $\left|{\lambda }_2\right|\ne 1$ if and only if $F\left(-1\right)=0$ and $B\ne \mathrm{0,2}$ ;

5) ${\lambda }_1=-1$ and ${\lambda }_2=-1$ if and only if $F\left(-1\right)=0$ and $B=2$ ;

6) $\left|{\lambda }_1\right|$ and $\left|{\lambda }_2\right|$ are a pair of conjugate complex roots and $\left|{\lambda }_1\right|=\left|{\lambda }_2\right|=1$ if and only if $-2<B<2$ and $C=1$ .

2. The Existence and Stability of Fixed Points

In this section, we give the conditions for existence of fixed points of system (1.2) in $R_{+}^2$ and discuss the property of these fixed points.

It is clear that the fixed points of the system (1.2) satisfy the following equations:

$\{\begin{array}{l}x=x\left(a-bx\right)-\frac{cxy}{x+my},\\ y=y\left(-d+\frac{fx}{x+my}\right).\end{array}$ (2.1)

By a simple computation, we have

Lemma 2.1 System (1.2) has two fixed points:

(i) the boundary equilibrium point $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ exists if $a\ne 1$ ;

(ii) the unique positive fixed point $E_2\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ , where $x^{\mathrm{<mo>\; *\; </mo>}}=\frac{a-1}{b}-\frac{c\left(f-d-1\right)}{mbf}$ and $y^{\mathrm{<mo>\; *\; </mo>}}=\frac{f-d-1}{\left(d+1\right)m}{x}^{\mathrm{<mo>\; *\; </mo>}}$ . The fixed point exists if and only if any one of the following conditions holds:

(ii.1) $d+1<f$ and $c\le m\left(a-1\right)$ ;

(ii.2) $d+1<f<\frac{c\left(d+1\right)}{c-m\left(a-1\right)}$ and $c>m\left(a-1\right)$ .

Next, we study the stability of the system (1.2) at fixed points. The Jacobian matrix of the system (1.2) evaluated at $E\left(x\mathrm{,}y\right)$ is as follows:

$J\left(x\mathrm{,}y\right)=\left(\begin{array}{cc}a-2bx-\frac{cm{y}^2}{{\left(x+my\right)}^2}& -\frac{c{x}^2}{{\left(x+my\right)}^2}\\ \frac{fm{y}^2}{{\left(x+my\right)}^2}& -d+\frac{f{x}^2}{{\left(x+my\right)}^2}\end{array}\right)\mathrm{.}$ (2.2)

The characteristic equation of the matrix can be written as

${\lambda }^2+p\left(x,y\right)\lambda +q\left(x,y\right)=0,$

where

$p\left(x\mathrm{,}y\right)=d-a+2bx+\frac{cm{y}^2-f{x}^2}{{\left(x+my\right)}^2}\mathrm{,}$

$q\left(x\mathrm{,}y\right)=-ad+2bdx+\frac{af{x}^2+cmd{y}^2-2bf{x}^3}{{\left(x+my\right)}^2}\mathrm{.}$

Proposition 2.2 When $a\ne 1$ , system (1.2) has a boundary equilibrium point $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ and

(i) it is a sink if $1<a<3$ and $d-1<f<d+1$ ;

(ii) it is a source if one of the following conditions holds:

(ii.1) $0<a<1$ and $0<f<d-1$ ;

(ii.2) $0<a<1$ and $f>d+1$ ;

(ii.3) $a>3$ and $0<f<d-1$ ;

(ii.4) $a>3$ and $f>d+1$ ;

(iii) it is non-hyperbolic if $a=3$ or $f=d-1$ or $f=d+1$ ;

(iv) it is a saddle except if conditions (1)-(3) are satisfied.

Proof. The Jacobian matrix of the fixed point $E_1$ is

${J|}_{E_1}=\left(\begin{array}{cc}2-a& -c\\ 0& f-d\end{array}\right)\mathrm{.}$ (2.3)

The two eigenvalues of the characteristic equation are:

${\lambda }_1=2-a;{\lambda }_2=f-d.$

Then, the fixed point $\left(\frac{a-1}{b}\mathrm{,0}\right)$ is a sink if $\left|{\lambda }_i\right|<1$ , where $i=1,2$ .

$\left|{\lambda }_1\right|=\left|2-a\right|<1\text{then}1<a<3,$

$\left|{\lambda }_2\right|=\left|f-d\right|<1\text{then}d-1<f<d+1.$

So, the fixed point $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ is stable if the following condition are satisfied:

$1<a<3,d-1<f<d+1.$

The other three items of the propositions are proved by using the above method.

Condition (iii) can be written as

$F{B}_{E_1}^1=\left\{\left(a,b,c,d,m,f\right)\in \left(0,+\infty \right):f=d+1,a\ne 3,a\ne 1\right\},$ (2.4)

$F{B}_{E_1}^2=\left\{\left(a,b,c,d,m,f\right)\in \left(0,+\infty \right):f=d-1,a\ne 3,a\ne 1\right\},$ (2.5)

$F{B}_{E_1}^3=\left\{\left(a,b,c,d,m,f\right)\in \left(0,+\infty \right):a=3,a\ne 1,f\ne d\pm 1\right\}.$ (2.6)

For simple calculation, we denote $R_1=2+\frac{c}{m}+\frac{m{f}^2-c{\left(d+1\right)}^2\left(f-d-2\right)}{mf\mathrm{<mo\; stretchy="false">\; [\; </mo>}df-{\mathrm{<mo\; stretchy="false">\; (\; </mo>}d+\mathrm{1\; <mo\; stretchy="false">\; )\; </mo>}}^2\mathrm{<mo\; stretchy="false">\; ]\; </mo>}}$ and $R_2=3+\frac{c}{m}-\frac{c{\left(d+1\right)}^2\left(f-d-3\right)}{mf\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]}$ .

Proposition 2.3. When $f>d+1$ and $a>1+\frac{c\left(f-d-1\right)}{mf}$ , system (1.2) has a unique positive fixed point $E_2\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ and

(i) it is a sink if one of the following conditions holds:

(i.1) $d\ge 1$ , $f<\frac{{\left(d+1\right)}^2}{d}$ and $R_1<a<R_2$ ;

(i.2) $d\ge 1$ , $f=\frac{{\left(d+1\right)}^2}{d}$ , $m>\frac{cd}{{\left(d+1\right)}^2}$ and $a<R_2$ ;

(i.3) $0<d<1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a<\min \left\{R_1,R_2\right\}$ ;

(i.4) $d=1$ , $f>\frac{{\left(d+1\right)}^2}{d}$ and $a<\min \left\{R_1,R_2\right\}$ ;

(i.5) $d>1$ , $\frac{{\left(d+1\right)}^2}{d}<f<\frac{{\left(d+1\right)}^2}{d-1}$ and $a<\min \left\{R_1,R_2\right\}$ ;

(i.6) $d>1$ , $f=\frac{{\left(d+1\right)}^2}{d-1}$ and $a<R_1$ ;

(i.7) $d>1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $R_2<a<R_1$ ;

(ii) it is a source if one of the following conditions holds:

(ii.1) $d\ge 1$ , $f<\frac{{\left(d+1\right)}^2}{d}$ and $a<\min \left\{R_1,R_2\right\}$ ;

(ii.2) $d\ge 1$ , $f=\frac{{\left(d+1\right)}^2}{d}$ , $m<\frac{cd}{{\left(d+1\right)}^2}$ and $a<R_2$ ;

(ii.3) $0<d<1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $R_1<a<R_2$ ;

(ii.4) $d=1$ , $f>\frac{{\left(d+1\right)}^2}{d}$ and $R_1<a<R_2$ ;

(ii.5) $d>1$ , $\frac{{\left(d+1\right)}^2}{d}<f<\frac{{\left(d+1\right)}^2}{d-1}$ and $R_1<a<R_2$ ;

(ii.6) $d>1$ , $f=\frac{{\left(d+1\right)}^2}{d-1}$ and $a>R_1$ ;

(ii.7) $d>1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a>\max \left\{R_1,R_2\right\}$ ;

(iii) it is a saddle if one of the following conditions holds:

(iii.1) $0<d<1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a>R_2$ ;

(iii.2) $d=1$ and $a>R_2$ ;

(iii.3) $d>1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a<R_2$ ;

(iii.4) $d>1$ , $f<\frac{{\left(d+1\right)}^2}{d-1}$ and $a>R_2$ ;

(iv) it is a non-hyperbolic if one of the following conditions holds:

(iv.1) $f=\frac{{\left(d+1\right)}^2}{d}$ , $c=mf$ and $2+\frac{1}{d}<a<6+\frac{1}{d}$ ;

(iv.2) $f>\frac{{\left(d+1\right)}^2}{d}$ and $mf-\frac{4md{f}^3-4m{\left(d+1\right)}^2{f}^2}{{\left(d+1\right)}^2{\left(f-d-1\right)}^2}<c<mf$ ;

(iv.3) $f\ne \frac{{\left(d+1\right)}^2}{d-1}$ and $a=R_2$ .

Proof. The Jacobian matrix of the fixed point $E_2$ is

${J|}_{E_2}=\left(\begin{array}{cc}1-b{x}^{\mathrm{<mo>\; *\; </mo>}}+\frac{c{x}^{\mathrm{<mo>\; *\; </mo>}}{y}^{\mathrm{<mo>\; *\; </mo>}}}{{\left(x^{\mathrm{<mo>\; *\; </mo>}}+m{y}^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}& -\frac{c{\left(x^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}{{\left(x^{\mathrm{<mo>\; *\; </mo>}}+m{y}^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}\\ \frac{fm{\left(y^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}{{\left(x^{\mathrm{<mo>\; *\; </mo>}}+m{y}^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}& 1-\frac{mf{x}^{\mathrm{<mo>\; *\; </mo>}}{y}^{\mathrm{<mo>\; *\; </mo>}}}{{\left(x^{\mathrm{<mo>\; *\; </mo>}}+m{y}^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}\end{array}\right)\mathrm{.}$ (2.7)

By computation, we have the characteristic equation of $J\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ are follows:

${\lambda }^2+B\left(x^{*},y^{*}\right)\lambda +C\left(x^{*},y^{*}\right)=0,$

where $B\left(x^{*},y^{*}\right)=-2+b{x}^{*}-\frac{\left(c-mf\right)x^{*}{y}^{*}}{{\left(x^{*}+m{y}^{*}\right)}^2}$ , $C\left(x^{*},y^{*}\right)=1-b{x}^{*}+\frac{\left(c-mf\right)x^{*}{y}^{*}}{{\left(x^{*}+m{y}^{*}\right)}^2}+\frac{bmf{\left(x^{*}\right)}^2{y}^{*}}{{\left(x^{*}+m{y}^{*}\right)}^2}$ , $x^{\mathrm{<mo>\; *\; </mo>}}=\frac{c\left(d+1-f\right)+mf\left(a-1\right)}{mbf}$ and $y^{*}=\frac{f-d-1}{\left(d+1\right)m}{x}^{*}$ .

Let

$F\left(\lambda \right)={\lambda }^2+B\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)\lambda +C\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{.}$

Obviously, $F\left(1\right)=\frac{bmf{\left(x^{*}\right)}^2{y}^{*}}{{\left(x^{*}+m{y}^{*}\right)}^2}>0$ ,

$\begin{array}{c}F\left(-1\right)=\frac{\left(a-3\right)\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]}{f}-\frac{c\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]}{mf}\\ +\frac{c{\left(d+1\right)}^2\left(f-d-3\right)}{m{f}^2},\end{array}$ (2.8)

$\begin{array}{c}C=\frac{a\left[df-{\left(d+1\right)}^2\right]}{f}-\frac{2\left[df-{\left(d+1\right)}^2\right]}{f}-\frac{c\left[df-{\left(d+1\right)}^2\right]}{f}\\ +\frac{c{\left(d+1\right)}^2\left(f-d-2\right)}{m{f}^2},\end{array}$ (2.9)

$B=a+d-2-\frac{c}{m}-\frac{{\left(d+1\right)}^2}{f}+\frac{c{\left(d+1\right)}^2}{m{f}^2}.$ (2.10)

From Lemma 1.1, we can see that both $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|<1$ if and only if $F\left(-1\right)>0$ and $C<1$ . Then

$\begin{array}{l}F\left(-1\right)>0\Rightarrow \left[\left(d-1\right)f-{\left(d+1\right)}^2\right]mfa>3\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]mf\\ +\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]cf-c{\left(d+1\right)}^2\left(f-d-3\right).\end{array}$

If $\left(d-1\right)f-{\left(d+1\right)}^2>0$ , then $d>1$ , $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a>R_2$ ;

If $\left(d-1\right)f-{\left(d+1\right)}^2=0$ , then $f=\frac{{\left(d+1\right)}^2}{d-1}$ and $d>1$ ;

If $\left(d-1\right)f-{\left(d+1\right)}^2<0$ , then we can divided d into there cases: $d>1$ , $d=1$ and $0<d<1$ .

Case I: if $d>1$ , then $f<\frac{{\left(d+1\right)}^2}{d-1}$ and $a<R_2$ ;

Case II: if $d=1$ , then $a<R_2$ ;

Case III: if $0<d<1$ , then $f>\frac{{\left(d+1\right)}^2}{d-1}$ and $a<R_2$ .

$\begin{array}{l}C<1\Rightarrow \left[df-{\left(d+1\right)}^2\right]mfa<2\left[df-{\left(d+1\right)}^2\right]mf+\left[df-{\left(d+1\right)}^2\right]cf\\ -c{\left(d+1\right)}^2\left(f-d-2\right)+m{f}^2.\end{array}$

If $df-{\left(d+1\right)}^2>0$ , then $f>\frac{{\left(d+1\right)}^2}{d}$ and $a<R_1$ ;

If $df-{\left(d+1\right)}^2=0$ , then $f=\frac{{\left(d+1\right)}^2}{d}$ and $m>\frac{cd}{{\left(d+1\right)}^2}$ ;

If $df-{\left(d+1\right)}^2<0$ , then $f<\frac{{\left(d+1\right)}^2}{d}$ and $a>R_1$ .

From $F\left(-1\right)>0$ and $C<1$ , we get conclusion (i). The propositions (ii) and (iii) are proved by using the above method.

Note that the fixed point $E_2\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is non-hyperbolic if either $\left|{\lambda }_1\right|=1$ or $\left|{\lambda }_2\right|=1$ . There also are two cases: $\Delta <0$ and $\Delta \ge 0$ . When $\Delta <0$ , $E_2\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is non-hyperbolic if and only if $C=1$ and $-2<B<2$ . Then

$-2<B<2\Rightarrow -d+\frac{c}{m}+\frac{{\left(d+1\right)}^2}{f}-\frac{c{\left(d+1\right)}^2}{m{f}^2}<a<4-d+\frac{c}{m}+\frac{{\left(d+1\right)}^2}{f}-\frac{c{\left(d+1\right)}^2}{m{f}^2}$

$\begin{array}{l}C=1\Rightarrow \left[df-{\left(d+1\right)}^2\right]mfa-2\left[df-{\left(d+1\right)}^2\right]mf-\left[df-{\left(d+1\right)}^2\right]cf\\ +c{\left(d+1\right)}^2\left(f-d-2\right)=m{f}^2.\end{array}$

If $f=\frac{{\left(d+1\right)}^2}{d}$ , then $c=mf$ ; If $f\ne \frac{{\left(d+1\right)}^2}{d}$ , then $a=R_1$ . From $C=1$ and $-2<B<2$ , we get conclusion (iv.1, 2).

When $\Delta \ge 0$ , $E_2\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is non-hyperbolic if and only if $F\left(-1\right)=0$ . Then

$\begin{array}{l}F\left(-1\right)=0\Rightarrow \left[\left(d-1\right)f-{\left(d+1\right)}^2\right]mfa-3\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]mf\\ -\left[\left(d-1\right)f-{\left(d+1\right)}^2\right]cf+c{\left(d+1\right)}^2\left(f-d-3\right)=0.\end{array}$

From $\Delta \ge 0$ and $F\left(-1\right)=0$ , we get conclusion (iv.3). Thus, we finish the proof.

3. Fold Bifurcation

From Proposition 2.2 (iii), we know that the eigenvalues at $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ are ${\lambda }_1=2-a$ and ${\lambda }_2=1$ if both $f=d+1$ and $a\ne \mathrm{1,3}$ .

Let $u\left(n\right)=x\left(n\right)-\frac{a-1}{b}$ , $v\left(n\right)=y\left(n\right)$ and $f^{*}=f-d-1$ . We transform the fixed pointed $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ into the origin. After Taylor expansion, the system (1.2) becomes

$\left(\begin{array}{c}u\left(n+1\right)\\ v\left(n+1\right)\end{array}\right)=\left(\begin{array}{cc}2-a& -c\\ 0& 1\end{array}\right)\left(\begin{array}{c}u\left(n\right)\\ v\left(n\right)\end{array}\right)+\left(\begin{array}{c}{F}_1\left(u\left(n\right)\mathrm{,}v\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ F_2\left(u\left(n\right)\mathrm{,}v\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\end{array}\right)\mathrm{,}$ (3.1)

where $U={\left(u\left(n\right)\mathrm{,}v\left(n\right)\right)}^{\text{T}}$ and

$\begin{array}{l}{F}_1\left(u\left(n\right)\mathrm{,}v\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ =-b{u}^2\left(n\right)+\frac{bcm}{a-1}{v}^2\left(n\right)-\frac{b^{2}cm}{{\left(a-1\right)}^2}u\left(n\right)v^2\left(n\right)-\frac{b^{2}c{m}^2}{{\left(a-1\right)}^2}{v}^3\left(n\right)+O\left({\Vert U\Vert }^4\right),\end{array}$

$\begin{array}{l}{F}_2\left(u\left(n\right)\mathrm{,}v\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ =f^{\mathrm{<mo>\; *\; </mo>}}v\left(n\right)-\frac{bm\left(f^{\mathrm{<mo>\; *\; </mo>}}+d+1\right)}{a-1}{v}^2\left(n\right)+\frac{b^{2}m\left(f^{\mathrm{<mo>\; *\; </mo>}}+d+1\right)}{{\left(a-1\right)}^2}u\left(n\right)v^2\left(n\right)\\ +\frac{b^2{m}^2\left(f^{\mathrm{<mo>\; *\; </mo>}}+d+1\right)}{{\left(a-1\right)}^2}{v}^3\mathrm{<mo\; stretchy="false">\; (\; </mo>}n\mathrm{<mo\; stretchy="false">\; )\; </mo>}+O\left({\Vert U\Vert }^4\right).\end{array}$

When $a\ne \mathrm{1,3}$ , we construct an invertible matrix

$T_1=\left(\begin{array}{cc}1& c\\ 0& 1-a\end{array}\right).$

Let

$\left(\begin{array}{c}u\left(n\right)\\ v\left(n\right)\end{array}\right)=T_1\left(\begin{array}{c}\tilde{x}\left(n\right)\\ \tilde{y}\left(n\right)\end{array}\right)\mathrm{.}$

Then system (3.1) can be transformed into

$\left(\begin{array}{c}\tilde{x}\left(n+1\right)\\ \tilde{y}\left(n+1\right)\end{array}\right)=\left(\begin{array}{cc}2-a& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\tilde{x}\left(n\right)\\ \tilde{y}\left(n\right)\end{array}\right)+\left(\begin{array}{c}{g}_1\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),f^{*}\right)\\ g_2\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),f^{*}\right)\end{array}\right),$ (3.2)

where

$\begin{array}{l}{g}_1\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),f^{*}\right)\\ =-c{f}^{\mathrm{<mo>\; *\; </mo>}}\tilde{y}\left(n\right)+b^{2}cm\frac{f^{\mathrm{<mo>\; *\; </mo>}}+d-a+2}{a-1}\tilde{x}\left(n\right){\tilde{y}}^2\left(n\right)\\ -bc\left[c+m\left(f^{\mathrm{<mo>\; *\; </mo>}}+d-a+2\right)\right]{\tilde{y}}^2\left(n\right)-b{\tilde{x}}^2\left(n\right)-2bc\tilde{x}\left(n\right)\tilde{y}\left(n\right)\\ +\left[b^2{c}^{2}m\frac{f^{\mathrm{<mo>\; *\; </mo>}}+d-a+2}{a-1}-b^{2}c{m}^2\left(f^{\mathrm{<mo>\; *\; </mo>}}+d-a+2\right)\right]{\tilde{y}}^3\left(n\right)\mathrm{,}\end{array}$

$\begin{array}{l}{g}_2\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),f^{*}\right)\\ =f^{*}\tilde{y}\left(n\right)+bm\left(f^{*}+d+1\right){\tilde{y}}^2\left(n\right)+\frac{b^{2}m\left(f^{*}+d+1\right)}{1-a}\tilde{x}{\tilde{y}}^2\left(n\right)\\ +\left[b^{2}cm\frac{f^{*}+d+1}{1-a}+b^2{m}^2\left(f^{*}+d+1\right)\right]{\tilde{y}}^3\left(n\right).\end{array}$

By the center manifold theorem and the method of the references [13], we know that there exists a center manifold $M^c\left(0\right)$ of system (3.2) at the origin in a small neighborhood of $f^{*}=0$ , which can be approximately represented as follows:

$M^c\left(0\right)=\left\{\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),f^{*}\right)\in R^3|\tilde{x}\left(n\right)=h\left(\tilde{y}\left(n\right),f^{*}\right),h\left(0,0\right)=0,Dh\left(0,0\right)=0\right\}.$

For $\tilde{y}\left(n\right)$ and $f^{\mathrm{<mo>\; *\; </mo>}}$ sufficiently small, we assume that a center manifold of the form

$h\left(\tilde{y}\left(n\right),f^{*}\right)=b_0{f}^{*}+b_1{f}^{*}\tilde{y}\left(n\right)+b_2{\left(f^{*}\right)}^2+b_3{\tilde{y}}^2\left(n\right)+O\left({\left(\left|\tilde{y}\left(n\right)\right|+\left|f^{*}\right|\right)}^3\right),$

where $O\left({\left(\left|\tilde{y}\left(n\right)\right|+\left|f^{*}\right|\right)}^3\right)$ is a function with order at least three of its variables.

$\begin{array}{c}\tilde{x}\left(n+1\right)=h\left(\tilde{y}\left(n+1\right),f^{*}\right)\\ =b_0{f}^{*}+b_1{f}^{*}\tilde{y}\left(n\right)+b_2{\left(f^{*}\right)}^2+b_3{\tilde{y}}^2\left(n\right)+O\left({\left(\left|\tilde{y}\left(n\right)\right|+\left|f^{*}\right|\right)}^3\right).\end{array}$ (3.3)

Also,

$\begin{array}{c}\tilde{x}\left(n+1\right)=\left(2-a\right)h\mathrm{<mo\; stretchy="false">\; (\; </mo>}\tilde{y}\mathrm{<mo\; stretchy="false">\; (\; </mo>}n\mathrm{<mo\; stretchy="false">\; )\; </mo>,}{f}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo\; stretchy="false">\; )\; </mo>}+g_1\left(h\left(\tilde{y}\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{,}\tilde{y}\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)\\ =\left(2-a\right)b_0{f}^{\mathrm{<mo>\; *\; </mo>}}+\left[\left(2-a\right)b_1-c-2b_{0}bc\right]f^{\mathrm{<mo>\; *\; </mo>}}\tilde{y}\left(n\right)\\ +\left[\left(2-a\right)b_2-b{b}_0^2\right]{\left(f^{\mathrm{<mo>\; *\; </mo>}}\right)}^2+[\left(2-a\right)b_3-b{c}^2-mbcd\\ {}_{}{}^{}+mabc-2mbc]{\tilde{y}}^2\left(n\right)+O\left({\left(\left|\tilde{y}\left(n\right)\right|+\left|f^{*}\right|\right)}^3\right)\mathrm{.}\end{array}$ (3.4)

Compare the coefficients for (3.3) and (3.4), we have

$b_0=0,b_1=\frac{c}{1-a},b_2=0,b_3=\frac{b{c}^2+mbcd-mabc+2mbc}{1-a}.$

System (3.2) restricted to the center manifold is given by

$\tilde{y}\left(n+1\right)=\left(1+f^{\mathrm{<mo>\; *\; </mo>}}\right)\tilde{y}\left(n\right)+bm\left(f^{\mathrm{<mo>\; *\; </mo>}}+d+1\right){\tilde{y}}^2\left(n\right)+O\left({\left(\left|\tilde{y}\left(n\right)\right|+\left|f^{*}\right|\right)}^3\right)\mathrm{.}$ (3.5)

Let

$F\left(\tilde{y}\left(n\right)\mathrm{,}{f}^{\mathrm{<mo>\; *\; </mo>}}\right)=\tilde{y}\left(n+1\right)\mathrm{,}$

then

${\frac{\partial F}{\partial \tilde{y}\left(n\right)}|}_{\left(0,0\right)}=1+f^{*}+2bm\tilde{y}\left(n\right)\left(f^{*}+d+1\right)=1,$

${\frac{\partial F}{\partial f^{*}}|}_{\left(0,0\right)}=\tilde{y}\left(n\right)+bm{\tilde{y}}^2\left(n\right)=0,$

${\frac{{\partial }^{2}F}{\partial \tilde{y}{\left(n\right)}^2}|}_{\left(0,0\right)}=2bm\left(f^{*}+d+1\right)=2bm\left(d+1\right)>0,$

$\frac{{\partial }^{2}F}{\partial \tilde{y}\left(n\right)\partial f^{*}}{|}_{\left(0,0\right)}=1+2bm\tilde{y}\left(n\right)=1,{\frac{{\partial }^{2}F}{\partial {\left(f^{*}\right)}^2}|}_{\left(0,0\right)}=0.$

Theorem 3.1 System (1.2) undergoes a fold bifurcation at $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ when parameter f varies in a small neighbourhood of $F{B}_{E_1}^1$ .

4. Flip Bifurcation

From Proposition 2.2 (iii), we know that the eigenvalues at $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ are ${\lambda }_1=-1$ and ${\lambda }_2=f-d$ if both $a=3$ and $f-d\ne \pm 1$ (similarly for the case of $F{B}_{E_1}^2$ ).

Let $u\left(n\right)=x\left(n\right)-\frac{a-1}{b}$ , $v\left(n\right)=y\left(n\right)$ and $a^{*}=a-3$ . We transform the fixed pointed $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ into the origin. After Taylor expansion, the system (1.2) becomes

$\left(\begin{array}{c}u\left(n+1\right)\\ v\left(n+1\right)\end{array}\right)=\left(\begin{array}{cc}-1& -c\\ 0& f-d\end{array}\right)\left(\begin{array}{c}u\left(n\right)\\ v\left(n\right)\end{array}\right)+\left(\begin{array}{c}{F}_1\left(u\left(n\right),v\left(n\right),a^{*}\right)\\ F_2\left(u\left(n\right),v\left(n\right),a^{*}\right)\end{array}\right),$ (4.1)

where $U={\left(u\left(n\right)\mathrm{,}v\left(n\right)\right)}^{\text{T}}$ and

$\begin{array}{l}{F}_1\left(u\left(n\right),v\left(n\right),a^{*}\right)\\ =-a^{\mathrm{<mo>\; *\; </mo>}}u\left(n\right)-b{u}^2\left(n\right)+\frac{bcm}{a^{\mathrm{<mo>\; *\; </mo>}}+2}{v}^2\left(n\right)-\frac{b^{2}cm}{{\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)}^2}u\left(n\right)v^2\left(n\right)\\ -\frac{b^{2}c{m}^2}{{\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)}^2}{v}^3\left(n\right)+O\left({\Vert U\Vert }^4\right),\end{array}$

$\begin{array}{l}{F}_2\left(u\left(n\right),v\left(n\right),a^{*}\right)\\ =-\frac{bmf}{a^{\mathrm{<mo>\; *\; </mo>}}+2}{v}^2\left(n\right)+\frac{b^{2}mf}{{\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)}^2}u\left(n\right)v^2\left(n\right)+\frac{b^2{m}^{2}f}{{\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)}^2}{v}^3\left(n\right)+O\left({\Vert U\Vert }^4\right).\end{array}$

When $f-d\ne 1$ , we construct an invertible matrix

$T_2=\left(\begin{array}{cc}1& -\left(a^{*}+2\right)c\\ 0& \left(a^{*}+2\right)\left(f-d+1\right)\end{array}\right).$

Let

$\left(\begin{array}{c}u\left(n\right)\\ v\left(n\right)\end{array}\right)=T_2\left(\begin{array}{c}\tilde{x}\left(n\right)\\ \tilde{y}\left(n\right)\end{array}\right).$

Then system (4.1) can be transformed into

$\left(\begin{array}{c}\tilde{x}\left(n+1\right)\\ \tilde{y}\left(n+1\right)\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& f-d\end{array}\right)\left(\begin{array}{c}\tilde{x}\left(n\right)\\ \tilde{y}\left(n\right)\end{array}\right)+\left(\begin{array}{c}{g}_1\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),a^{*}\right)\\ g_2\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),a^{*}\right)\end{array}\right),$ (4.2)

where

$\begin{array}{l}{g}_1\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),a^{*}\right)\\ =-a^{\mathrm{<mo>\; *\; </mo>}}\tilde{x}\left(n\right)-b{\tilde{x}}^2\left(n\right)+\left[bcm\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)\left(1-d\right)\left(f-d+1\right)-b{c}^2{\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)}^2\right]{\tilde{y}}^2\left(n\right)\\ +c{a}^{\mathrm{<mo>\; *\; </mo>}}\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)\tilde{y}\left(n\right)+2bc\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)\tilde{x}\left(n\right)\tilde{y}\left(n\right)\\ +b^{2}cm\mathrm{<mo\; stretchy="false">\; (\; </mo>}d-\mathrm{1\; <mo\; stretchy="false">\; )\; </mo>}\left(f-d+1\right)\tilde{x}\left(n\right){\tilde{y}}^2\left(n\right)\\ +b^{2}cm\left(d-1\right)\left(a^{\mathrm{<mo>\; *\; </mo>}}+2\right)\left(f-d+1\right)\left[m\left(f-d+1\right)-c\right]{\tilde{y}}^3\left(n\right)\mathrm{,}\end{array}$

$\begin{array}{l}{g}_2\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),a^{*}\right)\\ =-bmf\left(f-d+1\right){\tilde{y}}^2\left(n\right)+\frac{b^{2}mf\left(f-d+1\right)}{a^{\mathrm{<mo>\; *\; </mo>}}+2}\tilde{x}\left(n\right){\tilde{y}}^2\left(n\right)\\ +b^{2}mf\left(f-d+1\right)\left[m\left(f-d+1\right)-c\right]{\tilde{y}}^3\left(n\right)\mathrm{.}\end{array}$

By the center manifold theorem and the method of the references [13], we know that there exists a center manifold $M^c\left(0\right)$ of system (4.2) at the origin in a small neighborhood of $a^{*}=0$ , which can be approximately represented as follows:

$M^c\left(0\right)=\left\{\left(\tilde{x}\left(n\right),\tilde{y}\left(n\right),a^{*}\right)\in R^3|\tilde{y}\left(n\right)=h\left(\tilde{x}\left(n\right),a^{*}\right),h\left(0,0\right)=0,Dh\left(0,0\right)=0\right\}.$

For $\tilde{x}\left(n\right)$ and $a^{\mathrm{<mo>\; *\; </mo>}}$ sufficiently small, we assume that a center manifold of the form

$h\left(\tilde{x}\left(n\right),a^{*}\right)=b_0{a}^{*}+b_1{a}^{*}\tilde{x}\left(n\right)+b_2{\left(a^{*}\right)}^2+b_3{\tilde{x}}^2\left(n\right)+O\left({\left(\left|\tilde{x}\left(n\right)\right|+\left|a^{*}\right|\right)}^3\right),$

where $O\left({\left(\left|\tilde{x}\left(n\right)\right|+\left|a^{*}\right|\right)}^3\right)$ is a function with order at least three of its variables.

$\begin{array}{c}\tilde{y}\left(n+1\right)=h\left(\tilde{x}\left(n+1\right),a^{*}\right)\\ =b_0{a}^{*}-b_1{a}^{*}\tilde{x}\left(n\right)+b_2{\left(a^{*}\right)}^2+b_3{\tilde{x}}^2\left(n\right)+O\left({\left(\left|\tilde{x}\left(n\right)\right|+\left|a^{*}\right|\right)}^3\right).\end{array}$ (4.3)

Also,

$\begin{array}{c}\tilde{y}\left(n+1\right)=\left(f-d\right)h\left(\tilde{x}\left(n\right),a^{*}\right)+g_1\left(\tilde{x}\left(n\right),h\left(\tilde{x}\left(n\right),a^{*}\right),a^{*}\right)\\ =\left(f-d\right)b_0{a}^{*}+\left(f-d\right)b_1{a}^{*}\tilde{x}\left(n\right)\\ +\left[\left(f-d\right)b_2-b{b}_0^{2}mf\left(f-d+1\right)\right]{\left(a^{*}\right)}^2\\ +b_3\left(f-d\right){\tilde{x}}^2\left(n\right)+O\left({\left(\left|\tilde{x}\left(n\right)\right|+\left|a^{*}\right|\right)}^3\right).\end{array}$ (4.4)

Compare the coefficients for (4.3) and (4.4), we have

$b_0=0,b_1=0,b_2=0,b_3=0.$

System (4.2) restricted to the center manifold is given by

$\tilde{x}\left(n+1\right)=\left(-1-a^{*}\right)\tilde{x}\left(n\right)-b{\tilde{x}}^2\left(n\right)+O\left({\left(\left|\tilde{x}\left(n\right)\right|+\left|a^{*}\right|\right)}^3\right).$ (4.5)

Let

$F\left(\tilde{x}\left(n\right),a^{*}\right)=\tilde{x}\left(n+1\right).$

then

$F\left(0,0\right)=0,{\frac{\partial F}{\partial \tilde{x}\left(n\right)}|}_{\left(0,0\right)}=-1,$

$\alpha ={\left(\frac{{\partial }^{2}F}{\partial \tilde{x}\left(n\right)\partial a^{*}}+\frac{1}{2}\frac{\partial F}{\partial a^{*}}\frac{{\partial }^{2}F}{\partial {\tilde{x}}^2\left(n\right)}\right)|}_{\left(0,0\right)}=1\ne 0,$

$\beta ={\left[\frac{1}{6}\frac{{\partial }^{3}F}{\partial {\tilde{x}}^3\left(n\right)}+{\left(\frac{1}{2}\frac{{\partial }^{2}F}{\partial {\tilde{x}}^2\left(n\right)}\right)}^2\right]|}_{\left(0,0\right)}=b^2>0.$

Based on the above analysis, we have the following results.

Theorem 4.1 System (1.2) undergoes a flip bifurcation at $E_1\left(\frac{a-1}{b}\mathrm{,0}\right)$ when parameter a varies in a small neighbourhood of $F{B}_{E_1}^3$ .

5. Numerical Simulations

In this section, we use the bifurcation diagrams, phase portraits diagrams and the Maximum Lyapunov exponents diagrams of the system (1.2) to validate the above analytical result. Then, we consider bifurcation parameters in the following two cases.

Case I: Suppose that parameters $a=2$ , $b=1$ , $c=1$ , $d=3$ , $m=1$ in system (1.2). We see that a fold bifurcation occurs at the fixed point $\left(\frac{a-1}{b}\mathrm{,0}\right)$ when $f=d+1=4$ . This confirms Theorem 3.1.

Figure 1 shows the bifurcation diagram of system (1.2). We can see that the fixed point $E_1\left(\mathrm{1,0}\right)$ is stable on the center manifold for $3.5<f\le 4$ . When $f>4$ the fixed point $E_1\left(\mathrm{1,0}\right)$ on the center manifold becomes unstable. At the same time, there occurs a new fixed point $\left(\frac{4}{f}\mathrm{,1}-\frac{4}{f}\right)$ when $f>4$ . Hence, the transcritical type of fold bifurcation takes place. The phase diagrams of system (1.2) which are associated with Figure 1 are displayed in Figures 2(a)-(h) with different f.

Case II: Suppose that parameters $f=2.9$ , $b=1$ , $c=1$ , $d=2$ , $m=1$ in system (1.2). We see that a flip bifurcation occurs at the fixed point $\left(\frac{a-1}{b}\mathrm{,0}\right)$ when $a=3$ . This confirms Theorem 4.1.

Figure 3(a) shows the bifurcation diagram of system (1.2). We can see that the fixed point $E_1\left(a-\mathrm{1,0}\right)$ is stable for $a<3$ . When $a=3$ , system (1.2) occurs a subcritical flip bifurcation at the fixed point $E_1\left(a-\mathrm{1,0}\right)$ . When a increases to 3.6, system (1.2) undergoes a cascade of inverse period-doubling bifurcations. Subsequently, the system (1.2) near the fixed point $E_1\left(a-\mathrm{1,0}\right)$ is chaotic with narrow periodic windows. The Maximum Lyapunov exponents corresponding to Figure 3(a) are calculated in Figure 3(b), which is consistent with the bifurcation diagram. The phase diagrams of system (1.2) which are associated with Figure 3 are displayed in Figure 4(a)-(h) with different f.

(a) (b)

Figure 1. (a) Bifurcation diagram of the system (1.2) in (f, x)-plane when a = 2, b = 1, c = 1, d = 3, m = 1 and f varies in [3.5, 4.5], the initial value is (0.9, 0.1). (b) Bifurcation diagram of the system (1.2) in (f, y)-plane with same parameters in (a).

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 3. Phase portraits for various values of f, a = 2, b = 1, c = 1, d = 3, m = 1, the initial value is (0.9, 0.1), (a) f = 1.496; (b) f = 1.54; (c) f = 1.56; (d) f = 1.57; (e) f = 1.59; (f) f = 1.6; (g) f = 1.62; (h) f = 2.

(a) (b)

Figure 3. (a) Supercritical flip bifurcation diagram of the system (1.2) on the (a, x)-plane with f = 2.9, b = 1, c = 1, d = 2, m = 1, x₀ = 1.9, y₀ = 0.1; (b) Maximum Lyapunov exponents corresponding to (a).

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 4. Phase portraits for various values of f, a = 3, b = 1, c = 1, d = 2, m = 1, the initial value is (1.9, 0.1), (a) f = 0.79; (b) f = 0.8; (c) f = 0.9; (d) f = 1.1; (e) f = 2.9; (f) f = 3.01; (g) f = 3.153; (h) f = 3.154.

6. Conclusion

In this paper, a discrete prey-predator model with ratio-dependent functional response in the closed first quadrant is established. We obtain the equilibrium points of system (1.2) and the conditions for its existence. Moreover, we show that the boundary equilibrium point of system (1.2) can undergo flip bifurcation and Hopf bifurcation by using center manifold theorem and bifurcation theory. The system (1.2) displays interesting dynamical behaviors, including cascade of period-doubling, chaotic with narrow periodic windows and period-2 orbits, which implie that the boundary equilibrium points of system (1.2) lose stability via bifurcation. Nevertheless, the bifurcation of positive equilibrium point has not been discussed. It is a further research topic on how to set reasonable parameters to analyze the bifurcation at the positive equilibrium point.

Conflicts of Interest

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

References

[1]	May, R. (1976) Simple Mathematical Models with Very Complicated Dynamics, Nature, 261, 459-467. https://doi.org/10.1038/261459a0
[2]	Liu, X.L. and Xiao, D.M. (2007) Complex Dynamic Behaviors of a Discrete-Time Predator-Prey System, Chaos, Solitons & Fractals, 32, 80-94. https://doi.org/10.1016/j.chaos.2005.10.081
[3]	Chen, G., Teng, Z. and Hu, Z. (2011) Analysis of Stability for a Discrete Ratio-Dependent Predator-Prey System. Indian Journal of Pure and Applied Mathematics, 42, 1-26. https://doi.org/10.1007/s13226-011-0001-0
[4]	Elettreby, M.F., Nabil, T. and Khawagi, A. (2020) Stability and Bifurcation Analysis of a Discrete Predator-Prey Model with Mixed Holling Interaction. Computer Modeling in Engineering and Sciences, 122, 907-921. https://doi.org/10.32604/cmes.2020.08664
[5]	Elettreby, M.F., Khawagi, A. and Nabil, T. (2019) Dynamics of a Discrete Prey-Predator Model with Mixed Functional Response. International Journal of Bifurcation and Chaos, 29, 1950199. https://doi.org/10.1142/S0218127419501992
[6]	Arditi, R. and Ginzburg, L.R. (1989) Coupling in Predator-Prey Dynamics: Ratio-Dependence. Journal of Theoretical Biology, 139, 311-326. https://doi.org/10.1016/S0022-5193(89)80211-5
[7]	Rosenzweig, M.L. (1971) Paradox of Enrichment: Destabilization of Exploitation Ecosystems in Ecological Time. Science, 171, 385-387. https://doi.org/10.1126/science.171.3969.385
[8]	Xiao, D.M. and Ruan, S.G. (2001) Global Dynamics of a Ratio-Dependent Predator-Prey System. Journal of Mathematical Biology, 43, 268-290. https://doi.org/10.1007/s002850100097
[9]	Xiao, D.M. and Jennings, L.S. (2005) Bifurcations of a Ratio-Dependent Predator-Prey System with Constant Rate Harvesting. SIAM Journal on Applied Mathematics, 65, 737-753. https://doi.org/10.1137/S0036139903428719
[10]	Xiao, D.M., Li, W.X. and Han, M.A. (2006) Dynamics in a Ratio-Dependent Predator-Prey Model with Predator Harvesting. Journal of Mathematical Analysis and Applications, 324, 14-29. https://doi.org/10.1016/j.jmaa.2005.11.048
[11]	Kuang, Y. and Beretta, E. (1998) Global Qualitative Analysis of a Ratio-Dependent Predator-Prey System. Journal of Mathematical Biology, 36, 389-406. https://doi.org/10.1007/s002850050105
[12]	Sohel Rana, S.M. (2017) Chaotic Dynamics and Control of Discrete Ratio-Dependent Predator-Prey System, Discrete Dynamics in Nature and Society, 2017, Article ID: 4537450. https://doi.org/10.1155/2017/4537450
[13]	Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. Applied Mathematical Sciences, 2nd Edition, Springer-Verlag, New York, 112.