# Bifurcation and Chaos of a Discrete Predator-Prey Model with Crowley–Martin Functional Response Incorporating Proportional Prey Refuge

Authors: P. K. Santra; G. S. Mahapatra; G. R. Phaijoo
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The paper investigates the dynamical behaviors of a two-species discrete predator-prey system with Crowley–Martin functional response incorporating prey refuge proportional to prey density. The existence of equilibrium points, stability of three fixed points, period-doubling bifurcation, Neimark–Sacker bifurcation, Marottos chaos, and Control Chaos are analyzed for the discrete-time domain. The time graphs, phase portraits, and bifurcation diagrams are obtained for different parameters of the model. Numerical simulations and graphics show that the discrete model exhibits rich dynamics, which also present that the system is a chaotic and complex one. This paper attempts to present a feedback control method which can stabilize chaotic orbits at an unstable equilibrium point.


---

## Body

## 1. Introduction

The dynamical behavior of the prey-predator system in ecology has been the subject of study for many researchers and hence they have provided a substantial contribution to the growth of the population models [1–8] in the continuous-time domain. In case of populations that have overlapping generations and the birth processes occuring continuously, the ordinary differential equations are used for modeling of predator-prey interaction [9–30]. Many species, such as monocarpic plants, and semelparous animals have discrete non overlapping generations, and their births occur in regular breeding seasons. Their interactions are characterized by difference equations or formulated as discrete-time mappings. Discrete-time predator-prey models are capable of exhibiting much more complicated dynamics in comparison to the corresponding continuous-time models [31–45].The dynamical study of the prey-predator model with different functional responses is an effective approach for time domain. Upadhyay and Naji [46] presented three-species food chain model with Crowley–Martin type functional response for the dynamical study. Zhou [47] presented qualitative analysis of the predator-prey model for Crowley–Martin functional responses. Tripathi et al. [48] applied Crowley–Martin functional response on a delayed density dependent predator-prey model. Li et al. [49] presented a predator-prey model with Crowley–Martin functional response for positive periodic solutions.A certain degree of protection can be offered by nature to a set number of prey populations by providing refuges. The effects of prey refuges on the population dynamics are intricate and complex in nature; however, it can be considered to be constituted by two components for modeling purposes. The first component of refuge effects, which governs positive growth of prey and negatively that of predators, consists of a reduction in prey mortality owing to decreasing predation success. The second one may arise as the trade-offs and byproducts arising from the prey’s hiding behavior, which could cause advantages or even be detrimental for all the interacting populations. González-Olivares and Ramos-Jiliberto [8] presented prey refuges in a simple prey-predator system. Huang et al. [9] presented stability analysis of a prey-predator system incorporating a prey refuge. Ma et al. [12] discussed the effect of prey refuges on a prey-predator system. Chen et al. [13] studied the prey refuge on a Leslie–Gower predator-prey model.Extensive studies have been carried out by many researchers on refuge concept [9–30], and most of these papers consider refuge as a constant number of prey or proportional to prey. These earlier works have been taken into account that the use of refuges by a fraction of prey, or a constant number prey, is influential in existing a stabilizing effect in the dynamics of the interacting populations. This paper will test the above statement assuming that the quantity of prey in refuges is proportional to the prey. We also conduct an analysis of the dynamic properties of such a system via modification of the well-known Lotka–Volterra predator-prey model with prey self-limitation incorporating Crowley–Martin functional response [46, 47]. In this paper, we consider the refuge term proportional to prey density.The rest of the paper is organized as follows. In Section2, a discrete-time prey-predator model under nonoverlapping generation with refuge is formulated. Analysis of the local stability of the proposed model is dealt with in Section 3. The bifurcation analysis of the proposed model is studied through period-doubling bifurcation and Neimark–Sacker bifurcation in Section 4. Sections 5 and 6 give Marottos chaos and control chaos procedures, respectively. The proposed model is supported with the help of numerical simulations in Section 7. Section 8 makes the conclusion to this paper.


## 2. Description of the Proposed Predator-Prey Model

The proposed model belongs to a generalized predator-prey model which incorporates logistic growth having no stage structure for neither the prey nor the predator, which is represented by the following differential equations:(1)dxdt=rx1−xk−cφx,yy,dydt=dφx,yy−fy,where x and y denote the density of prey and predator populations at any time t, respectively, and r,k,c,d,andf are all positive constants and have their biological meanings accordingly. The intrinsic per capita growth rate of prey population is denoted by r, k is the environmental carrying capacity of prey population, c is the maximal per capita consumption rate of predators, d is the efficiency with which predators convert consumed prey into new predators, and f is the per capita death rate of predators. The term φx,y represents the functional response of the predator population.According to Maynard Smith [4], there is a quantity xr of prey population which incorporates refuges. Thereafter, the functional response is modified to incorporate prey refuges in the above model, and the term φx,y becomes the functional response in the form of φx−xr,y. Based on a laboratory study of dragonfly, Crowley and Martin present functional response φx,y=x/1+αx1+βy, where α>0 describe the handling time and β>0 is the magnitude of interference among predators. In the proposed model, we consider refuge xr=bx, and to derive a discrete-time model let dx/dt=xt+h−xt/h and dy/dt=yt+h−yt/h, where xt and yt are the densities of the prey and predator populations in discrete time t. Let h⟶1 and f=1, then the n+1th generation of the prey and predator populations by replacing t by n is obtained as follows:(2)xn+1=r+1xn1−rkr+1xn−c1−bxnyn1+α1−bxn1+βyn,yn+1=d1−bxnyn1+α1−bxn1+βyn.Let us considerr/kr+1=1 and r+1=a, then from (2) we obtain the following proposed discrete-time predator-prey system as follows:(3)xn+1=axn1−xn−c1−bxnyn1+α1−bxn1+βyn,yn+1=d1−bxnyn1+α1−bxn1+βyn,where a,b,c,d,α, and β are all positive constants. By biological meaning of the model variables, we only consider the system in the region Ω=x,y:x≥0,y≥0 in the x,y−plane.


## 3. Fixed Points and Stability Analysis of the Prey-Predator System

Fixed points of system (3) are determined by solving the following nonlinear system of equations:(4)x=ax1−x−c1−bxy1+α1−bx1+βy,y=d1−bxy1+α1−bx1+βy.By simple calculation, we get three nonnegative fixed points as follows:i P0=0,0, ii P1=a−1/a,0, a>1, and iii P2=x2,y2, where x2 is positive root of the equation Ax3+Bx2+Cx+D=0, where A=αβad1−b, B=αβd1−a1−b, C=βd1−a+c1−bd−α, and D=−c, and y2=1/βd1−bx2/1+α1−bx2−1.Lemma 1.
EquationAx3+Bx2+Cx+D=0 has i a unique positive root if a<1, ii three positive roots if 1<a<1+c1−bd−α/βd,d>α, and iii a unique positive root if a>1+c1−bd−α/βd,d>α.Now, the local behavior of system (3) is studied for each equilibrium point of the prey-predator system. The stability of system (3) is carried out by computing the Jacobian matrix corresponding to each equilibrium point. The Jacobian matrix J for system (3) is(5)J=a11a12a21a22,where(6)a11=a1−2x−c1−by1+α1−bx21+βy=a1−2x−cQ,a12=−c1−bx1+α1−bx1+βy2=−cR,a21=d1−by1+α1−bx21+βy=dQ,a22=d1−bx1+α1−bx1+βy2=dR,where(7)Q=1−by1+α1−bx21+βy,R=1−bx1+α1−bx1+βy2.The characteristic equation of matrixJ is λ2−TrJλ+DetJ=0, whereTrJ= Trace of matrix =a1−2x−cQ+dRDetJ= Determinant of matrix =ad1−2xRHence, model (3) is a dissipative dynamical system, if ad1−2xR<1, is a conservative dynamical one, if and only if ad1−2xR=1, and is an undissipated dynamical system otherwise.Proposition 1.
The equilibrium pointP0=0,0 is Sink if a<1, Saddle if a>1, and Nonhyperbolic if a=1.Proposition 2.
The equilibrium pointP1 is Sink if 2−a<1 and d1−a1−b/a−α1−a1−b<1, Source if 2−a>1 and d1−a1−b/a−α1−a1−b>1, Saddle if 2−a>1 and d1−a1−b/a−α1−a1−b<1 or 2−a<1 and d1−a1−b/a−α1−a1−b>1, and Nonhyperbolic if 2−a=1 or d1−a1−b/a−α1−a1−b=1Proposition 3.
Interior equilibrium pointP2x2,y2 is Sink if one of the following conditions holds:ix2<1/2 and maxdR−cQ−1/1−2x2+a,cQ−dR−1/1−2x2−a<adR<1/1−2x2iix2>1/2 and 1/1−2x2<adR<mindR−cQ−1/1−2x2+a,cQ−dR−1/1−2x2−aProposition 4.
Interior equilibrium pointP2x2,y2 is Source if one of the following conditions holds:ix2<1/2 and adR>maxdR−cQ−1/1−2x2+a,cQ−dR−1/1−2x2−a,1/1−2x2iix2>1/2 and adR<mindR−cQ−1/1−2x2+a,cQ−dR−1/1−2x2−a,1/1−2x2Proposition 5.
Interior equilibrium pointP2x2,y2 is Saddle if one of the following conditions holds:ix2<1/2 and dR−cQ−1/1−2x2+a<adR<cQ−dR−1/1−2x2−aiix2>1/2 and cQ−dR−1/1−2x2−a<adR<dR−cQ−1/1−2x2+aProposition 6.
Interior equilibrium pointP2x2,y2 is Nonhyperbolic if one of the following conditions holds:ix2<1/2 and cQ−dR−1/1−2x2−a=adR≠2/1−2x2>dR−cQ−1/1−2x2+aiix2>1/2 and cQ−dR−1/1−2x2−a=adR≠2/1−2x2<dR−cQ−1/1−2x2+aiiix2<1/2 and adR=1/1−2x2>maxdR−cQ−1/1−2x2+a,a1−2x−cQ+dR2/41−2x2ivx2>1/2 and adR=1/1−2x2<mindR−cQ−1/1−2x2+a,a1−2x−cQ+dR2/41−2x2PD1=a,b,c,d,α,β:x2<1/2andcQ−dR−1/1−2x2−a=adR≠2/1−2x2>dR−cQ−1/1−2x2+aPD2=a,b,c,d,α,β:x2>1/2andcQ−dR−1/1−2x2−a=adR≠2/1−2x2<dR−cQ−1/1−2x2+aNS1=a,b,c,d,α,β:x2<1/2andadR=1/1−2x2>maxdR−cQ−1/1−2x2+a,a1−2x−cQ+dR2/41−2x2NS2=a,b,c,d,α,β:x2>1/2andadR=1/1−2x2<mindR−cQ−1/1−2x2+a,a1−2x−cQ+dR2/41−2x2
The prey-predator model (3) undergoes period-doubling bifurcation at P2x2,y2, when parameters vary in a small neighborhood of PD1 or PD2. Model (3) undergoes Neimark–Sacker bifurcation at P2x2,y2, when parameters vary in a small neighborhood of NS1 or NS2.


## 4. Bifurcation Analysis

This section carries out an investigation of the conditions for the existence of Neimark–Sacker bifurcation (NSB) and period-doubling bifurcation (PDB) at the positive fixed pointP2x2,y2 of proposed system (3). There is an emergence of different kinds of bifurcations from the fixed point in dynamical systems, when a particular parameter passes through its critical value. Many dynamical properties of a system can be discussed owing to emergence of NSB and PDB. We discuss NSB and PDB for the positive fixed point P2x2,y2 of the prey-predator system (3) taking b as a bifurcation parameter.


### 4.1. Neimark–Sacker Bifurcation

This section presents NSB of the proposed prey-predator model (3) at P2x2,y2 for the parameters are located in the following set:(8)NS1=a,b,c,d,α,β:x2<12andadR=11−2x2>maxdR−cQ−11−2x2+a,a1−2x−cQ+dR241−2x2.Similar arguments can be applied to the other case:(9)NS2=a,b,c,d,α,β:x2>12andadR=11−2x2<mindR−cQ−11−2x2+a,a1−2x−cQ+dR241−2x2.In analyzing the NSB,b is used as the bifurcation parameter. Furthermore, b∗ (b∗⋘1) is the perturbation of b, and we consider a perturbation of the model as follows:(10)xn+1=axn1−xn−c1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡fxn,yn,b∗,yn+1=d1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then the equilibrium point P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn=0,0 to the third order, the proposed model (10) becomes(11)un+1=α1un+α2vn+α11un2+α12unvn+α22vn2+α111un3+α112un2vn+α122unvn2+α222vn3+Oun+vn4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β111un3+β112un2vn+β122unvn2+β222vn3+Oun+vn4,where(12)α1=fxx2,y2,0=a1−2x2−c1−by21+α1−bx221+βy2,α2=fyx2,y2,0=−c1−bx21+α1−bx21+βy22,α11=fxxx2,y2,0=−2a+2αc1−b2y21+α1−bx231+βy2,α12=fxyx2,y2,0=−c1−b1+α1−bx221+βy22,α22=fyyx2,y2,0=2βc1−bx21+α1−bx21+βy23,α111=fxxxx2,y2,0=−6α2c1−b3y21+α1−bx241+βy2,α112=fxxyx2,y2,0=2αc1−b21+α1−bx231+βy22,α122=fxyyx2,y2,0=2βc1−b1+α1−bx221+βy23,α222=fyyyx2,y2,0=−6β2c1−bx21+α1−bx21+βy24,β1=gxx2,y2,0=d1−by21+α1−bx221+βy2,β2=gyx2,y2,0=d1−bx21+α1−bx21+βy22,β11=gxxx2,y2,0=−2αd1−b2y21+α1−bx231+βy2,β12=gxyx2,y2,0=d1−b1+α1−bx221+βy22,β22=gyyx2,y2,0=−2βd1−bx21+α1−bx21+βy23,β111=gxxxx2,y2,0=6α2d1−b3y21+α1−bx241+βy2,β112=gxxyx2,y2,0=−2αd1−b21+α1−bx231+βy22,β122=gxyyx2,y2,0=−2βd1−b1+α1−bx221+βy23,β222=gyyyx2,y2,0=6β2d1−bx21+α1−bx21+βy24.The characteristic equation associated with the linearization of model (11) at un,vn=0,0 is given by λ2−TrJ1b∗λ+DetJ1b∗=0. Roots of the characteristic equation are(13)λ1,2b∗=TrJ1b∗±i4DetJ1b∗−TrJ1b∗22.Fromλ1,2b∗=1, when b∗=0, we have λ1,2b∗=DetJ1b∗1/2 and l=dλ1,2b∗/db∗h∗=0≠0.In addition, whenb∗=0, λ1,2i≠1, i=1,2,3,4, this is equivalent to TrJ10≠−2,−1,1,2.Now to study the normal form, letγ=Imλ1,2 and δ=Reλ1,2. We define T=01γδ, and using the transformation unvn=Tx¯ny¯n, model (11) becomes(14)x¯n+1=δx¯n−γy¯n+f1x¯n,y¯n,y¯n+1=γx¯n+δy¯n+g1x¯n,y¯n,where the functions f1 and g1 denote the terms in model (14) in variables x¯n,y¯n with the order at least two.In order to undergo Neimark–Sacker Bifurcation, we require the following discriminatory quantityΩ to be nonzero:(15)Ω=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(16)ξ20=18δ2β22−δα22−α12+4γα22+14γα12+18δi4γα22−2α22−2δα22+18i4γβ22+2γ2α22−2α11+18β12+δα11−2β11+δ3α22−δ2β22−δ2α12+δβ124γ,ξ11=12γβ22−δα22+12iγ2α22+α11+δα12+δ2α22+β11−δα11+δβ12−δ2α122γ−δ2β22−δ3α22γ,ξ02=14γ2δα22+α12+β22+14iβ12+2δβ22−2δα12−α11−β11−δα11+δβ12−δ2α124γ+14α22iγ2−3δ2+δ2β22−δ3α224γ,ξ21=38β222γ2+δ2+18β112+14δα112+14δβ122+α12218γ2+38δ2−14δ+38α111+38α222iγ2+2δ2+38α122γδi−18β122γi−38β222γδi−3β111−3δα111+3δβ112−3δ2α112+3δ2β122−3δ3α122+3δ3β222−3δ4α2228γi.From the above analysis we have the following result.Theorem 1.
IfΩ≠0, then model (3) undergoes NSB at P2x2,y2 when the parameter b∗ varies in a small neighborhood of the origin. Moreover, if Ω<0 Ω>0, then an attracting (repelling) inariant closed curve bifurcates from P2x2,y2 for b∗>0 b∗<0.


### 4.2. Period-Doubling Bifurcation

It is clear that one of the eigenvalues of the positive fixed pointP2x2,y2 is λ1=−1 and the other λ2 is neither 1 nor −1 if parameters of the model are located in the following set:(17)PD1=a,b,c,d,α,β:x2<12andcQ−dR−11−2x2−a=adR≠21−2x2>dR−cQ−11−2x2+aorPD2=a,b,c,d,α,β:x2>12andcQ−dR−11−2x2−a=adR≠21−2x2<dR−cQ−11−2x2+a.Here, we discuss PDB of the proposed model (3) at P2x2,y2, when parameters vary in a small neighborhood of PD1. Similar arguments can be applied to the other cases. In analyzing the PDB, b is used as the bifurcation parameter. Furthermore, b∗ (b∗⋘1) is the perturbation of b, and we consider a perturbation of the model as follows: (18)xn+1=axn1−xn−c1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡fxn,yn,b∗,yn+1=d1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then the equilibrium point P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn,b∗=0,0,0 to the third order, model (18) becomes(19)un+1=α1un+α2vn+α11un2+α12unvn+α13unb∗+α23vnb∗+α111un3+α112un2vn+α113un2b∗+α123unvnb∗+Oun+vn+b∗4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β13unb∗+β23vnb∗+β111un3+β112un2vn+β113un2b∗+β123unvnb∗+β223vn2b∗+Oun+vn+b∗4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α13=fxb∗x2,y2,0, α23=fyb∗x2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α113=fxxb∗x2,y2,0, α123=fxyb∗x2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β13=gxb∗x2,y2,0, β23=gyb∗x2,y2,0, β111=gxxxx2,y2,0, β113=gxxb∗x2,y2,0, β123=gxyb∗x2,y2,0, and β223=gyyb∗x2,y2,0.We defineT=α2α2−1−α1λ2−α1, and it is obvious that T is invertible. Using the transformation unvn=Tx¯ny¯n, then model (19) becomes(20)xn+1=−x¯n+f1un,vn,b∗,y¯n+1=λ2y¯n+g1un,vn,b∗,where the functions f1 and g1 denote the terms in model (20) in variables un,vn,b∗ with the order at least two.From the center manifold theorem, we know that there exists a center manifoldWc0,0,0 of model (20) at 0,0 in a small neighborhood of b∗=0, which can be approximately described as follows:(21)Wc0,0,0=x¯n,y¯n,b∗εR3:y¯n+1=α¯1x¯n2+α¯2x¯nb∗+Ox¯n+b∗3,where(22)α¯1=α21+α1α11+α2β111−λ22+β221+α121−λ22−1+α1α121+α1+α2β121−λ22,α¯2=1+α1α231+α1+α2β23α21+λ22−1+α1α13+α2β131+λ22.Model (20) which is restricted to center manifold Wc0,0,0 has the following form:(23)x¯n+1=−x¯n+h1x¯n2+h2x¯nb∗+h3x¯n2b∗+h4x¯nb∗2+h5x¯n3+Ox¯n+b∗3≡Fx¯n,b∗,where(24)h1=α¯2λ2−α¯1α11−α¯2β111+λ2−β221+α¯121+λ2−1+α¯1λ2−α¯1α12−α¯2β121+λ2,h2=λ2−α¯1α13−α¯2β131+λ2−1+α¯1λ2−α¯1α23−α¯2β23α¯21+λ2,h3=λ2−α1α23−α2β23λ2−α1α¯1α21+λ2+λ2−α1α¯1α13−α2β13+α2λ2−α1α113−α2β113−β2231+α121+λ2+2α2α¯2λ2−α1α11−α2β111+λ2−1+α1λ2−α1α123−α2β1231+λ2−2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h4=α¯2λ2−α1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯2α21+λ2+2α2α¯2λ2−α1α11−α2β111+λ2+2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h5=2α2α¯1λ2−α1α11−α2β11+α¯22λ2−α1α111−α2β1111+λ2+2β22α¯1λ2−α11+α11+λ2+λ2−α1α11−α2β11λ2−1−2α1α¯11+λ2−α¯21+α1λ2−α1α112−α2β1121+λ2.For flip bifurcation, we require the two discriminatory quantitiesξ1 and ξ2 to be nonzero:(25)ξ1=∂2F∂x¯∂b∗+12∂F∂b∗∂2F∂x¯20,0,ξ2=16∂3F∂x¯3+12∂2F∂x¯220,0.Finally, from the above analysis, we have the following result.Theorem 2.
Ifξ1≠0 and ξ2≠0, then model (3) undergoes PDB at P2x2,y2 when the parameter b varies. Moreover, if ξ2>0 ξ2<0, then the period-2 points that bifurcation from P2x2,y2 are stable (unstable).


## 4.1. Neimark–Sacker Bifurcation

This section presents NSB of the proposed prey-predator model (3) at P2x2,y2 for the parameters are located in the following set:(8)NS1=a,b,c,d,α,β:x2<12andadR=11−2x2>maxdR−cQ−11−2x2+a,a1−2x−cQ+dR241−2x2.Similar arguments can be applied to the other case:(9)NS2=a,b,c,d,α,β:x2>12andadR=11−2x2<mindR−cQ−11−2x2+a,a1−2x−cQ+dR241−2x2.In analyzing the NSB,b is used as the bifurcation parameter. Furthermore, b∗ (b∗⋘1) is the perturbation of b, and we consider a perturbation of the model as follows:(10)xn+1=axn1−xn−c1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡fxn,yn,b∗,yn+1=d1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then the equilibrium point P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn=0,0 to the third order, the proposed model (10) becomes(11)un+1=α1un+α2vn+α11un2+α12unvn+α22vn2+α111un3+α112un2vn+α122unvn2+α222vn3+Oun+vn4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β111un3+β112un2vn+β122unvn2+β222vn3+Oun+vn4,where(12)α1=fxx2,y2,0=a1−2x2−c1−by21+α1−bx221+βy2,α2=fyx2,y2,0=−c1−bx21+α1−bx21+βy22,α11=fxxx2,y2,0=−2a+2αc1−b2y21+α1−bx231+βy2,α12=fxyx2,y2,0=−c1−b1+α1−bx221+βy22,α22=fyyx2,y2,0=2βc1−bx21+α1−bx21+βy23,α111=fxxxx2,y2,0=−6α2c1−b3y21+α1−bx241+βy2,α112=fxxyx2,y2,0=2αc1−b21+α1−bx231+βy22,α122=fxyyx2,y2,0=2βc1−b1+α1−bx221+βy23,α222=fyyyx2,y2,0=−6β2c1−bx21+α1−bx21+βy24,β1=gxx2,y2,0=d1−by21+α1−bx221+βy2,β2=gyx2,y2,0=d1−bx21+α1−bx21+βy22,β11=gxxx2,y2,0=−2αd1−b2y21+α1−bx231+βy2,β12=gxyx2,y2,0=d1−b1+α1−bx221+βy22,β22=gyyx2,y2,0=−2βd1−bx21+α1−bx21+βy23,β111=gxxxx2,y2,0=6α2d1−b3y21+α1−bx241+βy2,β112=gxxyx2,y2,0=−2αd1−b21+α1−bx231+βy22,β122=gxyyx2,y2,0=−2βd1−b1+α1−bx221+βy23,β222=gyyyx2,y2,0=6β2d1−bx21+α1−bx21+βy24.The characteristic equation associated with the linearization of model (11) at un,vn=0,0 is given by λ2−TrJ1b∗λ+DetJ1b∗=0. Roots of the characteristic equation are(13)λ1,2b∗=TrJ1b∗±i4DetJ1b∗−TrJ1b∗22.Fromλ1,2b∗=1, when b∗=0, we have λ1,2b∗=DetJ1b∗1/2 and l=dλ1,2b∗/db∗h∗=0≠0.In addition, whenb∗=0, λ1,2i≠1, i=1,2,3,4, this is equivalent to TrJ10≠−2,−1,1,2.Now to study the normal form, letγ=Imλ1,2 and δ=Reλ1,2. We define T=01γδ, and using the transformation unvn=Tx¯ny¯n, model (11) becomes(14)x¯n+1=δx¯n−γy¯n+f1x¯n,y¯n,y¯n+1=γx¯n+δy¯n+g1x¯n,y¯n,where the functions f1 and g1 denote the terms in model (14) in variables x¯n,y¯n with the order at least two.In order to undergo Neimark–Sacker Bifurcation, we require the following discriminatory quantityΩ to be nonzero:(15)Ω=−Re1−2λ¯λ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(16)ξ20=18δ2β22−δα22−α12+4γα22+14γα12+18δi4γα22−2α22−2δα22+18i4γβ22+2γ2α22−2α11+18β12+δα11−2β11+δ3α22−δ2β22−δ2α12+δβ124γ,ξ11=12γβ22−δα22+12iγ2α22+α11+δα12+δ2α22+β11−δα11+δβ12−δ2α122γ−δ2β22−δ3α22γ,ξ02=14γ2δα22+α12+β22+14iβ12+2δβ22−2δα12−α11−β11−δα11+δβ12−δ2α124γ+14α22iγ2−3δ2+δ2β22−δ3α224γ,ξ21=38β222γ2+δ2+18β112+14δα112+14δβ122+α12218γ2+38δ2−14δ+38α111+38α222iγ2+2δ2+38α122γδi−18β122γi−38β222γδi−3β111−3δα111+3δβ112−3δ2α112+3δ2β122−3δ3α122+3δ3β222−3δ4α2228γi.From the above analysis we have the following result.Theorem 1.
IfΩ≠0, then model (3) undergoes NSB at P2x2,y2 when the parameter b∗ varies in a small neighborhood of the origin. Moreover, if Ω<0 Ω>0, then an attracting (repelling) inariant closed curve bifurcates from P2x2,y2 for b∗>0 b∗<0.


## 4.2. Period-Doubling Bifurcation

It is clear that one of the eigenvalues of the positive fixed pointP2x2,y2 is λ1=−1 and the other λ2 is neither 1 nor −1 if parameters of the model are located in the following set:(17)PD1=a,b,c,d,α,β:x2<12andcQ−dR−11−2x2−a=adR≠21−2x2>dR−cQ−11−2x2+aorPD2=a,b,c,d,α,β:x2>12andcQ−dR−11−2x2−a=adR≠21−2x2<dR−cQ−11−2x2+a.Here, we discuss PDB of the proposed model (3) at P2x2,y2, when parameters vary in a small neighborhood of PD1. Similar arguments can be applied to the other cases. In analyzing the PDB, b is used as the bifurcation parameter. Furthermore, b∗ (b∗⋘1) is the perturbation of b, and we consider a perturbation of the model as follows: (18)xn+1=axn1−xn−c1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡fxn,yn,b∗,yn+1=d1−b+b∗xnyn1+α1−b+b∗xn1+βyn≡gxn,yn,b∗.Letun=xn−x2andvn=yn−y2, then the equilibrium point P2x2,y2 is transformed into the origin, and further expanding f and g as a Taylor series at un,vn,b∗=0,0,0 to the third order, model (18) becomes(19)un+1=α1un+α2vn+α11un2+α12unvn+α13unb∗+α23vnb∗+α111un3+α112un2vn+α113un2b∗+α123unvnb∗+Oun+vn+b∗4,vn+1=β1un+β2vn+β11un2+β12unvn+β22vn2+β13unb∗+β23vnb∗+β111un3+β112un2vn+β113un2b∗+β123unvnb∗+β223vn2b∗+Oun+vn+b∗4,where α1=fxx2,y2,0, α2=fyx2,y2,0, α11=fxxx2,y2,0, α12=fxyx2,y2,0, α13=fxb∗x2,y2,0, α23=fyb∗x2,y2,0, α111=fxxxx2,y2,0, α112=fxxyx2,y2,0, α113=fxxb∗x2,y2,0, α123=fxyb∗x2,y2,0, β1=gxx2,y2,0, β2=gyx2,y2,0, β11=gxxx2,y2,0, β12=gxyx2,y2,0, β22=gyyx2,y2,0, β13=gxb∗x2,y2,0, β23=gyb∗x2,y2,0, β111=gxxxx2,y2,0, β113=gxxb∗x2,y2,0, β123=gxyb∗x2,y2,0, and β223=gyyb∗x2,y2,0.We defineT=α2α2−1−α1λ2−α1, and it is obvious that T is invertible. Using the transformation unvn=Tx¯ny¯n, then model (19) becomes(20)xn+1=−x¯n+f1un,vn,b∗,y¯n+1=λ2y¯n+g1un,vn,b∗,where the functions f1 and g1 denote the terms in model (20) in variables un,vn,b∗ with the order at least two.From the center manifold theorem, we know that there exists a center manifoldWc0,0,0 of model (20) at 0,0 in a small neighborhood of b∗=0, which can be approximately described as follows:(21)Wc0,0,0=x¯n,y¯n,b∗εR3:y¯n+1=α¯1x¯n2+α¯2x¯nb∗+Ox¯n+b∗3,where(22)α¯1=α21+α1α11+α2β111−λ22+β221+α121−λ22−1+α1α121+α1+α2β121−λ22,α¯2=1+α1α231+α1+α2β23α21+λ22−1+α1α13+α2β131+λ22.Model (20) which is restricted to center manifold Wc0,0,0 has the following form:(23)x¯n+1=−x¯n+h1x¯n2+h2x¯nb∗+h3x¯n2b∗+h4x¯nb∗2+h5x¯n3+Ox¯n+b∗3≡Fx¯n,b∗,where(24)h1=α¯2λ2−α¯1α11−α¯2β111+λ2−β221+α¯121+λ2−1+α¯1λ2−α¯1α12−α¯2β121+λ2,h2=λ2−α¯1α13−α¯2β131+λ2−1+α¯1λ2−α¯1α23−α¯2β23α¯21+λ2,h3=λ2−α1α23−α2β23λ2−α1α¯1α21+λ2+λ2−α1α¯1α13−α2β13+α2λ2−α1α113−α2β113−β2231+α121+λ2+2α2α¯2λ2−α1α11−α2β111+λ2−1+α1λ2−α1α123−α2β1231+λ2−2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h4=α¯2λ2−α1α13−α2β131+λ2+λ2−α1α23−α2β23λ2−α1α¯2α21+λ2+2α2α¯2λ2−α1α11−α2β111+λ2+2β22α¯21+α1λ2−α11+λ2+α¯2λ2−α1α12−α2β12λ2−1−2α11+λ2,h5=2α2α¯1λ2−α1α11−α2β11+α¯22λ2−α1α111−α2β1111+λ2+2β22α¯1λ2−α11+α11+λ2+λ2−α1α11−α2β11λ2−1−2α1α¯11+λ2−α¯21+α1λ2−α1α112−α2β1121+λ2.For flip bifurcation, we require the two discriminatory quantitiesξ1 and ξ2 to be nonzero:(25)ξ1=∂2F∂x¯∂b∗+12∂F∂b∗∂2F∂x¯20,0,ξ2=16∂3F∂x¯3+12∂2F∂x¯220,0.Finally, from the above analysis, we have the following result.Theorem 2.
Ifξ1≠0 and ξ2≠0, then model (3) undergoes PDB at P2x2,y2 when the parameter b varies. Moreover, if ξ2>0 ξ2<0, then the period-2 points that bifurcation from P2x2,y2 are stable (unstable).


## 5. Existence of Marottos Chaos

This section presents chaotic nature of system (3) in the sense of Marrotto.Definition 1.
Let the functionF:Rn⟶Rn be differentiable in BrZ. The point Z∈Rn is an expanding fixed point of F in BrZ, if FZ=Z and all eigenvalues of DFX exceed 1 in norm for all X∈BrZ.Definition 2.
Assume thatZ is an expanding fixed point of F in BrZ for some r>0. Then, Z is said to be a snap-back repeller of F if there exists a point X0∈BrZ with X0≠Z, FMX0=Z, and DFMX0≠0 for some positive integer M.Theorem 3.
IfP2∈U2=x,y:4ad1−2x2R2>max4,a1−2x2−cQ2+dR22, then P2 is an expanding fixed point of F.Proof.
For map,(26)FXn=axn1−xn−c1−bxnyn1+α1−bxn1+βynd1−bxnyn1+α1−bxn1+βyn,Xn=xnynT.
The eigenvalues corresponding with the fixed pointP2=x2,y2 are given by(27)λ1,2=−px2,y2±p2x2,y2−4qx2,y22,where(28)px2,y2=−a1−2x2−cQ2+dR2,qx2,y2=ad1−2x2R2,R2=1−bx21+α1−bx21+βy22,Q2=1−by21+α1−bx221+βy2.
Therefore, the fixed pointP2x2,y2 has a pair of complex eigenvalues, and the norm of them exceeds unity ifp2x2,y2−4qx2,y2<0 and qx2,y2−1>0i.e.,4ad1−2x2R2>max4,a1−2x2−cQ2+dR22.
Thus, we can state the following theorem.Theorem 4.
P
2
x
2
,
y
2 is a snap-back repeller in U2.Proof.
According to the definition of a snap-back repeller, one needs to find one pointP∗x0,y0∈U2 such that P∗≠P2, FMP∗=P2, and DFMP∗≠0, for some positive integer M, where Map F is defined by (3).
To proceed, notice that(29)x1=ax01−x0−c1−bx0y01+α1−bx01+βy0,y1=d1−bx0y01+α1−bx01+βy0,(30)x2=ax11−x1−c1−bx1y11+α1−bx11+βy1,y2=d1−bx1y11+α1−bx11+βy1.
Now, a mapF2 has been constructed to map the point P∗=x0,y0 to the fixed point P2x2,y2 after two iterations if there are solutions different from P2 for equations (29) and (30). The solutions different from P2 for equation (30) satisfy the following equation:(31)x1=y21+βy11−bdy1−αy21+βy1,y1=1+α1−bx11+βy1ax11−x1−x2c1−bx1.
Substitutingx1 and y1 into equation (29) and solving x0,y0, we have(32)x0=y11+βy01−bdy0−αy11+βy0,y0=1+α1−bx01+βy0ax01−x0−x1c1−bx0.
By simple calculations, we obtain(33)DF2P∗=aB−2AB−1+α1−bA1+βDBD+AEc1−b−c1−bADα1−bB1+βD+1+α1−bAβE1+α1−bA21+βD2×1+α1−bA1+βDd1−bCD+AF−d1−bADα1−bC1+βD+1+α1−bAβF1+α1−bA21+βD2−aC−2AC−1+α1−bA1+βDc1−bCD+AF−c1−bADα1−bC1+βD+1+α1−bAβF1+α1−bA21+βD2×1+α1−bA1+βDBD+AEd1−b−d1−bADα1−bB1+βD+1+α1−bAβE1+α1−bA21+βD2,where(34)A=ax1−x−c1−bxy1+α1−bx1+βy,B=a1−2x−cy1+α1−bx21+βy,C=−c1−bx1+α1−bx1+βy2,D=d1−bxy1+α1−bx1+βy,E=dy1+α1−bx21+βy,F=d1−bx1+α1−bx1+βy2.
The solutions of equations (31) and (32) will farther subject to x0,y0,x1,y1≠x2,y2,x0,y0∈U2 and DF2P∗≠0, then P2x2,y2 is a snap-back repeller in U2.
Thus, the following theorem is established.Theorem 5.
If4ad1−2x2R2>max4,a1−2x2−cQ2+dR22, the solutions x1,y1 and x0,y0 of equations (31) and (32) satisfy in addition x0,y0, x1,y1≠x2,y2, x0,y0∈U2,x0,y0=0,0, and DF2P∗≠0, then P2x2,y2 is a snap-back repeller of Map (3), and hence Map (3) is chaotic in the sense of Marotto.


## 6. Chaos Control of the Proposed System

Several methods can be used for obtaining chaos control in discrete-time models. A few of these methods are the state feedback method, pole-placement technique, and hybrid control method. We give a feedback control method to stabilize chaotic orbits at an unstable positive fixed point of the prey-predator system (3), with the following controlled form:(35)xn+1=axn1−xn−c1−bxnyn1+α1−bxn1+βyn+S,yn+1=d1−bxnyn1+α1−bxn1+βyn,with the following feedback control law as the control force:(36)S=−q1xn−x2−q2yn−y2,where q1 and q2 are the feedback gain and x2,y2 is a positive fixed point of the model.The Jacobian matrixJ for system (35) at x2,y2 is(37)J=a11−q1a12−q2a21a22,where(38)a11=a1−2x−c1−by1+α1−bx21+βy,a12=−c1−bx1+α1−bx1+βy2,a21=d1−by1+α1−bx21+βy,a22=d1−bx1+α1−bx1+βy2.The corresponding characteristic equation of matrixJ is(39)λ2−a11+a22−q1λ+a22a11−q1−a21a12−q2.Letλ1 and λ2 are the eigenvalues of the above characteristic equation:(40)λ1+λ2=a11+a22−q1,λ1λ2=a22a11−q1−a21a12−q2.The lines of marginal stability are determined by solving the equationsλ1=±1 and λ1λ2=1. These conditions guarantee that the eigenvalues λ1 and λ2 have modulus less than 1.Supposeλ1λ2=1, from (40), we have line l1 as follows:(41)a22q1−a21q2=a22a11−a21a12−1.Supposeλ1=±1, from (40), we have line l2 and l3 as follows:(42)1−a22q1+a21q2=a11+a22−1−a22a11+a21a12,1+a22q1−a21q2=a11+a22+1+a22a11−a21a12.The stable eigenvalues lie within a triangular region by linel1,l2, and l3.


## 7. Numerical Simulations

This section presents the qualitative dynamical behaviors of the system by presenting Lyapunov exponent, bifurcation diagram, and phase plane for specific parameter values. We present the numerical simulation to observe rich dynamics of our model, so we consider the hypothetical parameters values of the proposed discrete prey-predator system. The value of the parameters are consider asa=4.2, b=0.12, c=3.0, d=3.5, α=0.1, and β=0.1 with initial population 0.6,0.3 for simulation and graphical presentation of different dynamical observations. Figures 1–3 are iterative plots for different refuge parameters, and we observed that refuge has a stabilizing effect. Figures 4–9 are bifurcation diagrams and corresponding local amplification diagram of the prey and predator for varying parameters a,b,c,d,α, and β, respectively.Figure 1
Phase portraits forb=0.0,0.01,0.02,0.03.
(a)
(b)
(c)
(d)Figure 2
Phase portraits forb=0.05,0.06,0.065,0.10.
(a)
(b)
(c)
(d)Figure 3
Phase portraits forb=0.150,0.152,0.153,0.160.
(a)
(b)
(c)
(d)Figure 4
Orbit diagram of species for parametera and local amplification.
(a)
(b)
(c)
(d)Figure 5
Orbit diagram of species for parameterb and local amplification.
(a)
(b)
(c)
(d)Figure 6
Orbit diagram of species for parameterc and local amplification.
(a)
(b)
(c)
(d)Figure 7
Orbit diagram of species for parameterd and local amplification.
(a)
(b)
(c)
(d)Figure 8
Orbit diagram of species for parameterα and local amplification.
(a)
(b)
(c)
(d)Figure 9
Orbit diagram of species for parameterβ and local amplification.
(a)
(b)
(c)
(d)The phase portraits associated with Figure5 are displayed in Figures 1–3. A chaotic attractor for b=0.00 and 0.01. A five-pieces chaotic attractor for b=0.02 and 0.03. Period-10 orbit for b=0.05. Period-68 orbit for b=0.06. Attracting invariant cycle for b=0.065,0.10,0.15, and 0.152. Fixed point for b=0.153 and 0.160.The orbit diagram of species fora, which is a scaling parameter of the intrinsic per capita growth rate of prey population. Figure 4 shows that the system trajectory evolves from a fixed point to NSB and finally into a chaotic attractor. There is one visible periodic windows in the bifurcation diagram for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 10 is the corresponding Lyapunov exponents of Figure 4.Figure 10
Lyapunov exponents of the system fora.Figure5 shows the orbit diagram of prey and predator population for the refuge parameter b with other fixed parameter values. The system trajectory evolves into a periodic orbit from the chaotic attractor and finally into a fixed point. There is only one obvious periodic window in all the bifurcation process for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 11 is the corresponding Lyapunov exponent of Figure 5.Figure 11
Lyapunov exponents of the system forb.Figure6 shows the orbit diagram of the population for the parameter c of per capita consumption rate of predators. The system trajectory evolves from a fixed point to NSB and finally into a chaotic attractor. We also give local magnification of the corresponding bifurcation figure. Figure 12 is the corresponding Lyapunov exponent of Figure 6.Figure 12
Lyapunov exponents of the system forc.Figure7 shows the orbit diagram of species for the parameter d of efficiency with which predators convert the consumed prey into a new predator. The system trajectory evolves from a fixed point to NSB and finally into a chaotic attractor. There is one big visible periodic window in the bifurcation diagram for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 13 is the corresponding Lyapunov exponent of Figure 7.Figure 13
Lyapunov exponents of the system ford.Figure8 shows the orbit diagram of prey and predator population for the parameter α with other fixed parameter values. The system trajectory evolves into a fixed point from the periodic orbit. We also give local magnification of the corresponding bifurcation figure. Figure 14 is the corresponding Lyapunov exponent of Figure 8.Figure 14
Lyapunov exponents of the system forα.Figure9 shows the orbit diagram of prey and predator population for the parameter β with other fixed parameter values. The system trajectory evolves into a periodic orbit from the chaotic attractor and finally into a fixed point. There is only one obvious periodic window in all the bifurcation process for both prey and predator. We also give local magnification of the corresponding bifurcation figure. Figure 15 is the corresponding Lyapunov exponent of Figure 9.Figure 15
Lyapunov exponents of the system forβ.


### 7.1. Chaos Control

We observe chaotic behavior of prey-predator, as shown in Figure16(a). In this case fixed point 0.3183,0.6919 is unstable. In the feedback control method for feedback gain q1=0.5 and q2=−0.5, we observe the fixed point 0.3183,0.6919 is stable, as shown in Figure 16(b).Figure 16
Chaos control.
(a)
(b)


### 7.2. Fractal Dimension

The fractal dimension is defined by using Lyapunov exponents as follows.Letf be a map on Rm. Consider an orbit with Lyapunov exponents h1≥h2≥….≥hm, and let p denote the largest integer such that ∑i=1i=phi≥0.(43)Fractal dimensiondL=0,if no suchpexist,p+∑i=1i=phihp+1,ifp<m,m,ifp=m.Fractal dimension and Lyapunov exponents for different values ofb for the value of the parameters as a=4.2, c=3.0, d=3.5, α=0.1, and β=0.1 with initial population 0.6,0.3 are presented in Table 1.Table 1
Fractal dimension and Lyapunov exponents for different values ofb.
Value of
b
Lyapunov exponents
Fractal dimension
b
=
0.0
h
1
=
0.1295
h
2
=
0.0251
2
b
=
0.01
h
1
=
0.0956
h
2
=
0.0381
2
b
=
0.02
h
1
=
0.0538
h
2
=
−
0.0447
2.2036
b
=
0.03
h
1
=
0.0383
h
2
=
−
0.0991
1.3865


## 7.1. Chaos Control

We observe chaotic behavior of prey-predator, as shown in Figure16(a). In this case fixed point 0.3183,0.6919 is unstable. In the feedback control method for feedback gain q1=0.5 and q2=−0.5, we observe the fixed point 0.3183,0.6919 is stable, as shown in Figure 16(b).Figure 16
Chaos control.
(a)
(b)


## 7.2. Fractal Dimension

The fractal dimension is defined by using Lyapunov exponents as follows.Letf be a map on Rm. Consider an orbit with Lyapunov exponents h1≥h2≥….≥hm, and let p denote the largest integer such that ∑i=1i=phi≥0.(43)Fractal dimensiondL=0,if no suchpexist,p+∑i=1i=phihp+1,ifp<m,m,ifp=m.Fractal dimension and Lyapunov exponents for different values ofb for the value of the parameters as a=4.2, c=3.0, d=3.5, α=0.1, and β=0.1 with initial population 0.6,0.3 are presented in Table 1.Table 1
Fractal dimension and Lyapunov exponents for different values ofb.
Value of
b
Lyapunov exponents
Fractal dimension
b
=
0.0
h
1
=
0.1295
h
2
=
0.0251
2
b
=
0.01
h
1
=
0.0956
h
2
=
0.0381
2
b
=
0.02
h
1
=
0.0538
h
2
=
−
0.0447
2.2036
b
=
0.03
h
1
=
0.0383
h
2
=
−
0.0991
1.3865


## 8. Conclusion

In this paper, we have studied the effects of refuges on a predator-prey interaction, by using the analytical and graphical approach. The refuges are considered as prey refuge proportional to prey density. We evaluate the effects with regard to the local stability of the interior equilibrium point and the long-term dynamics of the interacting populations. The results show that the effects of refuges can stabilize the interior equilibrium point of the proposed prey-predator discrete-time domain model. We have studied bifurcations in a discrete predator-prey model with refuge. Furthermore, it is also shown that the model exhibits various bifurcations of codimension 1, including period-doubling bifurcation, and Neimark–Sacker bifurcation as the values of parameters vary. Numerical simulations and graphical presentation show the rich dynamics in bifurcation and chaotic nature of the system model.



---

*Source: 5309814-2020-01-22.xml*