# Discrete-Time Predator-Prey Model with Bifurcations and Chaos

Authors: K. S. Al-Basyouni; A. Q. Khan
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, local dynamics, bifurcations and chaos control in a discrete-time predator-prey model have been explored inℝ+2. It is proved that the model has a trivial fixed point for all parametric values and the unique positive fixed point under definite parametric conditions. By the existing linear stability theory, we studied the topological classifications at fixed points. It is explored that at trivial fixed point model does not undergo the flip bifurcation, but flip bifurcation occurs at the unique positive fixed point, and no other bifurcations occur at this point. Numerical simulations are performed not only to demonstrate obtained theoretical results but also to tell the complex behaviors in orbits of period-4, period-6, period-8, period-12, period-17, and period-18. We have computed the Maximum Lyapunov exponents as well as fractal dimension numerically to demonstrate the appearance of chaotic behaviors in the considered model. Further feedback control method is employed to stabilize chaos existing in the model. Finally, existence of periodic points at fixed points for the model is also explored.


---

## Body

## 1. Introduction

In the mid-1920s, an Italian Biologist U. D’Ancona studied the population variations in different species of fish that were interacting with one another. In his investigation, he found some percentage data of total catch of different fish species during the World War I which was brought into various Mediterranean ports. In precise, the data determine the percentage (%) of total catch of selachians, which are not adorable as food fish, and during the years 1914–1923, the percentage data for the port of Fiume, Italy, is given in Table1.Table 1
Percentage data for the port of Fiume, Italy, during the years 1914–1923.
Number of years% values191411.9191521.4191622.1191721.2191836.4191927.3192016.0192115.9192214.8192310.7D’Ancona was surprised by a very large increase in the percentage of selachians during the World War I. He reasoned the increase in the percentage of selachians was because of decline in fishing during this period. But how does the intensity of fishing affect the fish populations? The answer was obviously the struggle for existence between competing species, which was of great concern to D’Ancona and also to the fishing industry. Naturally, selachians depend on food fish for their survival as selachians are predators and food fish are prey. D’Ancona thought that this accounted for the large increase in number of selachians during the war period. Since the level of fishing was less, there were more food fish available to selachians, which therefore increased rapidly. D’Ancona just shows the increase in number of selachians when fishing’s level is reduced. D’Ancona did not explain why a declining level of fishing is more helpful to predators than to their prey. Later, paying all possible biological clarifications to this phenomenon, D’Ancona twisted to his coworker, the well-known mathematician Vito Volterra. With hope, Volterra would come up with a mathematical model of the growth of the selachians and food fish, their prey, and his model would provide the answer to D’Ancona’s question. Volterra started his analysis on this problem by separating all the fish into the prey populationxt and the predator population yt. Then, he reasoned that the food fish do not compete very fast among themselves for their food supply since this is very plentiful, and the density of fish population is not very much. Hence, in the absence of predators, their prey would grow according to the Malthusian law of population growth: dx/dt=rx for positive constant r. Next, Volterra reasoned the number of contacts per unit time between selachians (predators) and food fish (prey) is bxy for positive constant b. Therefore, dx/dt=rx−bxy, and so Volterra reasoned that predators have a characteristic pace of lessening −dy relative to their present number and that they likewise increment at a rate cxy relative to their present number y and their nourishment supply x, and hence, one has [1–5](1)dxdt=rx−bxy,dydt=cxy−dy.Now include the effects of fishing in (1). It is observed that fishing decreases the population of food fish at a rate εxt and decreases the population of selachians at a rate εyt, and the parameter ϵ reflects the intensity of fishing. So, by the effect of fishing, the continuous-time model which is depicted in (1) becomes as follows [1]:(2)dxdt=r−εx−bxy,dydt=cxy−d+εy.It is important here to note that generally many biological models are directed by continuous as well as discrete-time systems, and in recent years, many authors gave great contribution towards discrete models [6–15]. The reasons are that, for nonoverlapping generation, discrete models are much convincing than continuous models, and moreover, these models provide more effective computational results for numerical simulations [6, 16–21]. Due to scientific computation, in the present study, our aim is to explore the qualitative behavior of the discrete-time model corresponding to (2). It is predicted that the discrete model is dynamically consistent with the continuous-time model. It is also mentioned in [15] that if population has a nonoverlapping generation, as stated above, it is essential to write a discrete system that corresponds to model (2). So, applying the forward Euler scheme, (2) becomes as follows:(3)xt+1−xth=r−εxt−bxtyt,yt+1−yth=cxtyt−d+εyt.After some manipulations, from (3), one gets(4)xt+1=1+hr−εxt−bhxtyt,yt+1=1−hd+εyt+chxtyt.This paper is organized as follows: Section2 is about the topological classifications at fixed points of the model (4), whereas existence of possible bifurcations at respective fixed points is given in Section 3. The comprehensive bifurcation analysis at a positive fixed point is investigated in Section 4. In Section 5, some simulations are performed to demonstrate obtained theoretical results, and this is also about the study of fractal demission that characterized the strange attractors. In Section 6, we investigated the chaos control by the feedback control method, whereas existence of periodic points is studied in Section 7. The conclusion of the paper is given in Section 8.


## 2. Topological Classifications at Fixed Points

We will study the topological classifications at fixed points of the model (4) in this section. In order to determine fixed points, one need to solve the following system, where x^,y^ is the fixed point of (4):(5)x^=1+hr−εx^−bhx^y^,y^=1−hd+εy^+chx^y^.It is noted thatx^,y^=0,0, satisfying (5) obviously, and so Etrival=0,0 is a trivial fixed point of (4). For a positive fixed point, system (5) reduces into the following form:(6)hr−ε−bhy^=0,−hd+ε+chx^=0.From (6), one gets x=d+ε/c and y=r−ε/b. Therefore, (4) has the unique positive fixed point: Epositive=d+ε/c,r−ε/b if r>ϵ.Hereafter, we will study the topological classifications atEtrival and Epositive of the model (4) by method of linearization. So, the Jacobian matrix Ω|Ex^,y^ evaluated at Ex^,y^ becomes as follows:(7)Ω|Ex^,y^≔1+hr−ε−bhy^−bhx^chy^1−hd+ε+chx^.Now, the auxiliary equation ofΩ|Ex^,y^ at Ex^,y^ is as follows:(8)ζ2−TEx^,y^ζ+DEx^,y^=0,where(9)TEx^,y^=2+hr−2hε−bhy^−hd+chx^,DEx^,y^=1+hr−ε−bhy^1−hd+ε+chx^+bch2x^y^.Now, in the rest of the section, we will give topological classification at equilibria:Etrival and Epositive as follows:


### 2.1. Topological Classifications atEtrival

It is noted that, atEtrival, (7) takes the following form:(10)Ω|Etrival≔1+hr−ε001−hd+ε,whose characteristic roots are ζ1=1+hr−ε and ζ2=1−hd+ε. So based on stability theory, one can conclude the topological classifications at Etrival as the following result.Proposition 1.
ForEtrival, the following classifications hold:(I)
If(11)0<h<min2ε−r,2d+ε,whereε>r,thenEtrival is a sink.(II)
If(12)h>max2ε−r,2d+ε,whereε>r,thenEtrival is a source.(III)
If(13)2ε−r<h<2d+ε,thenEtrival is a saddle.(IV)
If(14)h=2ε−ror(15)h=2d+ε,thenEtrival is nonhyperbolic.


### 2.2. Topological Classifications atEpositive

Now, we will give topological classifications atEpositive for the considered system. After some manipulations, at Epositive, (7) takes the following form:(16)Ω|Epositive≔1−bhcd+εchbr−ε1.Furthermore, atEpositive, (8) becomes(17)ζ2−TEpositiveζ+DEpositive=0,where(18)TEpositive=2,DEpositive=1+h2r−εd+ε.Finally, from (17) one gets(19)ζ1,2=2±Disc2,where(20)Disc≔TEpositive2−4DEpositive=4h2ε−rd+ε.Hereafter, at equilibrium:Epositive, we will summarize the topological classifications by allowing sign of discriminant quantity, i.e., Disc=4h2ε−rd+ε<0respectively≥0, as a following proposition.Proposition 2.
IfDisc=4h2ε−rd+ε<0, then for Epositive, the following classifications hold:(I)
Epositive is never stable focus.(II)
If(21)r>ε,thenEpositive is unstable focus.(III)
If(22)r=ε,thenEpositive is nonhyperbolic.Proposition 3.
IfDisc=4h2ε−rd+ε≥0, then for Epositive, the following classifications hold:(I)
If(23)r>εh2d+ε−4h2d+εwith(24)d>4−ε2h2εh2,thenEpositive is a stable node.(II)
If(25)r<εh2d+ε−4h2d+ε,along with (24) holds, then Epositive is an unstable node.(III)
If(26)r=εh2d+ε−4h2d+ε,thenEpositive is nonhyperbolic.


## 2.1. Topological Classifications atEtrival

It is noted that, atEtrival, (7) takes the following form:(10)Ω|Etrival≔1+hr−ε001−hd+ε,whose characteristic roots are ζ1=1+hr−ε and ζ2=1−hd+ε. So based on stability theory, one can conclude the topological classifications at Etrival as the following result.Proposition 1.
ForEtrival, the following classifications hold:(I)
If(11)0<h<min2ε−r,2d+ε,whereε>r,thenEtrival is a sink.(II)
If(12)h>max2ε−r,2d+ε,whereε>r,thenEtrival is a source.(III)
If(13)2ε−r<h<2d+ε,thenEtrival is a saddle.(IV)
If(14)h=2ε−ror(15)h=2d+ε,thenEtrival is nonhyperbolic.


## 2.2. Topological Classifications atEpositive

Now, we will give topological classifications atEpositive for the considered system. After some manipulations, at Epositive, (7) takes the following form:(16)Ω|Epositive≔1−bhcd+εchbr−ε1.Furthermore, atEpositive, (8) becomes(17)ζ2−TEpositiveζ+DEpositive=0,where(18)TEpositive=2,DEpositive=1+h2r−εd+ε.Finally, from (17) one gets(19)ζ1,2=2±Disc2,where(20)Disc≔TEpositive2−4DEpositive=4h2ε−rd+ε.Hereafter, at equilibrium:Epositive, we will summarize the topological classifications by allowing sign of discriminant quantity, i.e., Disc=4h2ε−rd+ε<0respectively≥0, as a following proposition.Proposition 2.
IfDisc=4h2ε−rd+ε<0, then for Epositive, the following classifications hold:(I)
Epositive is never stable focus.(II)
If(21)r>ε,thenEpositive is unstable focus.(III)
If(22)r=ε,thenEpositive is nonhyperbolic.Proposition 3.
IfDisc=4h2ε−rd+ε≥0, then for Epositive, the following classifications hold:(I)
If(23)r>εh2d+ε−4h2d+εwith(24)d>4−ε2h2εh2,thenEpositive is a stable node.(II)
If(25)r<εh2d+ε−4h2d+ε,along with (24) holds, then Epositive is an unstable node.(III)
If(26)r=εh2d+ε−4h2d+ε,thenEpositive is nonhyperbolic.


## 3. Existence of Possible Bifurcations at Fixed Points:Etrival and Epositive

In view of obtained results in Section2 regarding topological classifications at equilibria Etrival and Epositive, we will study the existence of possible bifurcations in this section, as follows:(I)
Recall that, atEtrival, Ω|Etrival has two characteristic roots in which ζ1|14=−1, but ζ2|14=1−2/ε−rd+ε≠−1or1. This implies that (4) may undergo the flip bifurcation if b,c,d,h,r,ε locate in the set:(27)FB|Etrival=b,c,d,h,r,ε,h=2ε−r,orAgain under the nonhyperbolic condition, which is depicted in (15), one can obtain that ζ1|15=1+2/d+εr−ε≠−1or1, but ζ2|15=−1. This implies that model (4) may undergo the flip bifurcation if b,c,d,h,r,ε locate in the set:(28)FB|Etrival=b,c,d,h,r,ε,h=2d+ε.(II)
Under the hypothesis of Proposition2 and obtained nonhyperbolic condition, which is depicted in (22), one gets ζ1,222=1. Hence, model (4) may undergo a Neimark–Sacker bifurcation if parameters b,c,d,h,r,ε locate in the set:(29)HB|Epositive=b,c,d,h,r,ε,r=ε.(III)
Under the hypothesis of Proposition3 and obtained nonhyperbolic condition, which is depicted in (26), one gets ζ1|26=−1, but ζ2|26=3≠−1or1. So model (4) may undergo the flip bifurcation if b,c,d,h,r,ε locate in the set:(30)FB|Epositive=b,c,d,h,r,ε,r=εh2d+ε−4h2d+ε.Note: for case I, it is easy to see that model does not undergo flip bifurcation ifb,c,d,h,r,ε∈FB|Etrival, and hence, Etrival is degenerated with a higher codimension.In the subsequent section, we will present comprehensive N-S and flip bifurcations analysis when parameters, respectively,b,c,d,h,r,ε∈HB|Epositiveandb,c,d,h,r,ε∈FB|Epositive.


## 4. Comprehensive Bifurcation Analysis atEpositive

### 4.1. N-S Bifurcation atEpositive

In the subsequent section, we will discuss N-S bifurcation of model (4) at Epositive if b,c,d,h,r,ε∈HB|Epositive, and hence, the result can be stated as the following theorem.Theorem 1.
Ifb,c,d,h,r,ε∈HB|Epositive, then at Epositive, model (4) does not undergo N-S bifurcation.Proof.
Sinceb,c,d,h,r,ε∈HB|Epositive and hence if bifurcation parameter r varies in a small neighborhood of r∗, that is, r=r∗+τ, where τ≪1, then (4) becomes as follows:(31)xt+1=1+hr∗+τ−εxt−bhxtyt,yt+1=1−hd+εyt+chxtyt,with Epositiveτ=d+ε/c,r∗+τ−ε/b if r∗+τ>ε. The auxiliary equation of Ω|Epositiveτ at Epositiveτ becomes(32)ζ2−TEpositiveτζ+DEpositiveτ=0,where(33)TEpositiveτ=2,DEpositiveτ=1+h2r∗+τ−εd+ε.
Finally, from (32), one gets(34)ζ1,2τ=1±hr∗+τ−εd+ει.
From (34), we obtain(35)ζ1,2≔DEpositiveτ,dζ1,2dτ|τ=0≔h22d+ε>0.
Hereafter by using the following transformation(36)Φt=xt−d+εc,Ψt=yt−r∗+τ−εb,transform Epositiveτ in origin. In view of (36), from (31), one gets(37)Φt+1=1+hr∗+τ−εΦt+d+εc−bhΦt+d+εcΨt+r∗+τ−εb−d+εc,Ψt+1=1−hd+εΨt+r∗+τ−εb+chΦt+d+εcΨt+r∗+τ−εb−r∗+τ−εb.
Now, we will study normal form of (37) if τ=0. After Taylor series expansion about Φt,Ψt=0,0, one gets(38)Φt+1=Φt−bhcd+εΨt−bhΦtΨt,Ψt+1=chbr−εΦt+Ψt+chΦtΨt.
In order to transform liner part of (38) into canonical form, we use the following transformation that can be easily constructed by computation:(39)ΦtΨt=−bhcd+ε00hr−εd+εXtYt.
In view of (39) and (38),(40)Xt+1Yt+1=1hr−εd+ε−hr−εd+ε1XtYt+FXt,YtGXt,Yt,where(41)FXt,Yt=−bh2r−εd+εXtYt,GXt,Yt=−bh2d+εXtYt.
From (41), one gets(42)FXtXt|Etrival=FYtYt|Etrival=0,FXtYt|Etrival=−bh2r−εd+ε,GXtXt|Etrival=GYtYt|Etrival=0,GXtYt|Etrival=−bh2d+ε,FXtXtXt|Etrival=FXtXtYt|Etrival=FXtYtYt|Etrival=FYtYtYt|Etrival=0,GXtXtXt|Etrival=GXtXtYt|Etrival=GXtYtYt|Etrival=GYtYtYt|Etrival=0.
In order to undergo the said bifurcation, the following quantity should be nonzero [22–29]:(43)Γ=−ℜ1−2ζ¯ζ¯21−ζδ11δ20−12δ112−δ022+ℜζ¯δ21,where(44)δ02=18FXtXt−FYtYt+2GXtYt+ιGXtXt−GYtYt+2FXtYt|Etrival,δ11=14FXtXt+FYtYt+ιGXtXt+GYtYt|Etrival,δ20=18FXtXt−FYtYt+2GXtYt+ιGXtXt−GYtYt−2FXtYt|Etrival,δ21=116FXtXtXt+FXtYtYt+GXtXtYt+GYtYtYt+ιGXtXtXt+GXtYtYt−FXtXtYt−FYtYtYtEtrival.
Utilizing (42) in (44), one gets(45)δ02=−bh24d+ε+r−εd+ει,δ11=δ21=0,δ20=bh24−d−ε+r−εd+ει.
Using (45) in (43), one gets Ω=−bh4/16d+ε2<0. Finally, the model considered undergoing N-S bifurcation requires that ζ1,2m≠1,m=1,2,3,4 if τ=0 which corresponds to TEpositive0≠−2,0,1,2. However,TEpositive0=2 which contradicts to the fact that TEpositive0≠−2,0,1,2, and hence, eigenvalues ζ1,2 of fixed point 0,0 lay in the interaction of the unit circle with the coordinate axes when τ=0. Therefore, model (4) does not undergo N-S bifurcation.


### 4.2. Flip Bifurcation atEpositive

In the subsequent section, we will discuss flip bifurcation of model (4) at Epositive if b,c,d,h,r,ε∈FB|Epositive, and hence, the result can be stated as the following theorem.Theorem 2.
Ifb,c,d,h,r,ε∈FB|Epositive, then at Epositive, model (4) undergoes the flip bifurcation if p varies in a small neighborhood of the origin.Proof.
It is noted that ifr in a small nbhd of r∗, i.e, r=r∗+p, where p≪1, then model (4) becomes(46)xt+1=1+hr∗+p−εxt−bhxtyt,yt+1=1−hd+εyt+chxtyt.
Now, one can transformEpositivep in origin by using the following transformations:(47)Φt˜=xt−d+εc,Ψt˜=yt−r∗+p−εb.
From (47), (46) takes the following form:(48)Φt+1˜=1+hr∗+p−εΦt˜+d+εc−bhΦt˜+d+εcΨt˜+r∗+p−εb−d+εc,Ψt+1˜=1−hd+εΨt˜+r∗+p−εb+chΦt˜+d+εcΨt˜+r∗+p−εb−r∗+p−εb.
By Taylor series aboutΦt˜,Ψt˜=0,0, one gets(49)Φt+1˜=Φt˜−bhcd+εΨt˜−bhΦt˜Ψt˜+hpΦt˜,Ψt+1˜=chbr−εΦt˜+Ψt˜+chΦt˜Ψt˜.
Now by using following transformation(50)Φt˜Ψt˜≔−bcr−εε−rd+εbcr−εε−rd+ε11Xt˜Yt˜,linear part of (49) transforms into canonical form. In view of (50), (49) becomes(51)Xt+1˜Yt+1˜=−1003+F˜Xt˜,Yt˜,pG˜Xt˜,Yt˜,p,where(52)F˜Xt˜,Yt˜,p=−12bhXt˜2+12bhYt˜2−12hpXt˜,G˜Xt˜,Yt˜,p=12chXt˜2−12chYt˜2.
Now, for (51), center manifold McEtrival at Etrival is explored in a small nbhd of p. Therefore, McEtrival can be expressed as follows:(53)McEtrival=Xt˜,Yt˜:Yt˜=c0p+c1Xt˜2+c2Xt˜p+c3p3+OXt˜+p3.
After some manipulations, one gets(54)c0=c2=c3=0,c1=−14ch.
Finally, the map (51) restricting to McEtrival is(55)F˜Xt˜=−Xt˜+h1Xt˜2+h2Xt˜p+h3Xt˜2p+h4Xt˜p2+h5Xt˜3+OXt˜+p4,with(56)h1=−12bh,h2=−h2,h3=h4=h5=0.
In order for the map (51) to undergo flip bifurcation, the following discriminatory quantities are required to be nonzero [28, 29]:(57)∧1=∂2F˜∂Xt˜∂p+12∂F˜∂p∂2F˜∂Xt˜2|Xt˜,p=0,0,∧2=16∂3F˜∂Xt˜3+12∂2F˜∂Xt˜22|Xt˜,p=0,0.
On computation, one gets∧1=−h/2≠0, but ∧2=1/6b2h2>0, which shows that model (4) undergoes flip bifurcation if b,c,d,h,r,ε∈FB|Epositive and, in particular, stable period-2 points bifurcating from Epositive.


## 4.1. N-S Bifurcation atEpositive

In the subsequent section, we will discuss N-S bifurcation of model (4) at Epositive if b,c,d,h,r,ε∈HB|Epositive, and hence, the result can be stated as the following theorem.Theorem 1.
Ifb,c,d,h,r,ε∈HB|Epositive, then at Epositive, model (4) does not undergo N-S bifurcation.Proof.
Sinceb,c,d,h,r,ε∈HB|Epositive and hence if bifurcation parameter r varies in a small neighborhood of r∗, that is, r=r∗+τ, where τ≪1, then (4) becomes as follows:(31)xt+1=1+hr∗+τ−εxt−bhxtyt,yt+1=1−hd+εyt+chxtyt,with Epositiveτ=d+ε/c,r∗+τ−ε/b if r∗+τ>ε. The auxiliary equation of Ω|Epositiveτ at Epositiveτ becomes(32)ζ2−TEpositiveτζ+DEpositiveτ=0,where(33)TEpositiveτ=2,DEpositiveτ=1+h2r∗+τ−εd+ε.
Finally, from (32), one gets(34)ζ1,2τ=1±hr∗+τ−εd+ει.
From (34), we obtain(35)ζ1,2≔DEpositiveτ,dζ1,2dτ|τ=0≔h22d+ε>0.
Hereafter by using the following transformation(36)Φt=xt−d+εc,Ψt=yt−r∗+τ−εb,transform Epositiveτ in origin. In view of (36), from (31), one gets(37)Φt+1=1+hr∗+τ−εΦt+d+εc−bhΦt+d+εcΨt+r∗+τ−εb−d+εc,Ψt+1=1−hd+εΨt+r∗+τ−εb+chΦt+d+εcΨt+r∗+τ−εb−r∗+τ−εb.
Now, we will study normal form of (37) if τ=0. After Taylor series expansion about Φt,Ψt=0,0, one gets(38)Φt+1=Φt−bhcd+εΨt−bhΦtΨt,Ψt+1=chbr−εΦt+Ψt+chΦtΨt.
In order to transform liner part of (38) into canonical form, we use the following transformation that can be easily constructed by computation:(39)ΦtΨt=−bhcd+ε00hr−εd+εXtYt.
In view of (39) and (38),(40)Xt+1Yt+1=1hr−εd+ε−hr−εd+ε1XtYt+FXt,YtGXt,Yt,where(41)FXt,Yt=−bh2r−εd+εXtYt,GXt,Yt=−bh2d+εXtYt.
From (41), one gets(42)FXtXt|Etrival=FYtYt|Etrival=0,FXtYt|Etrival=−bh2r−εd+ε,GXtXt|Etrival=GYtYt|Etrival=0,GXtYt|Etrival=−bh2d+ε,FXtXtXt|Etrival=FXtXtYt|Etrival=FXtYtYt|Etrival=FYtYtYt|Etrival=0,GXtXtXt|Etrival=GXtXtYt|Etrival=GXtYtYt|Etrival=GYtYtYt|Etrival=0.
In order to undergo the said bifurcation, the following quantity should be nonzero [22–29]:(43)Γ=−ℜ1−2ζ¯ζ¯21−ζδ11δ20−12δ112−δ022+ℜζ¯δ21,where(44)δ02=18FXtXt−FYtYt+2GXtYt+ιGXtXt−GYtYt+2FXtYt|Etrival,δ11=14FXtXt+FYtYt+ιGXtXt+GYtYt|Etrival,δ20=18FXtXt−FYtYt+2GXtYt+ιGXtXt−GYtYt−2FXtYt|Etrival,δ21=116FXtXtXt+FXtYtYt+GXtXtYt+GYtYtYt+ιGXtXtXt+GXtYtYt−FXtXtYt−FYtYtYtEtrival.
Utilizing (42) in (44), one gets(45)δ02=−bh24d+ε+r−εd+ει,δ11=δ21=0,δ20=bh24−d−ε+r−εd+ει.
Using (45) in (43), one gets Ω=−bh4/16d+ε2<0. Finally, the model considered undergoing N-S bifurcation requires that ζ1,2m≠1,m=1,2,3,4 if τ=0 which corresponds to TEpositive0≠−2,0,1,2. However,TEpositive0=2 which contradicts to the fact that TEpositive0≠−2,0,1,2, and hence, eigenvalues ζ1,2 of fixed point 0,0 lay in the interaction of the unit circle with the coordinate axes when τ=0. Therefore, model (4) does not undergo N-S bifurcation.


## 4.2. Flip Bifurcation atEpositive

In the subsequent section, we will discuss flip bifurcation of model (4) at Epositive if b,c,d,h,r,ε∈FB|Epositive, and hence, the result can be stated as the following theorem.Theorem 2.
Ifb,c,d,h,r,ε∈FB|Epositive, then at Epositive, model (4) undergoes the flip bifurcation if p varies in a small neighborhood of the origin.Proof.
It is noted that ifr in a small nbhd of r∗, i.e, r=r∗+p, where p≪1, then model (4) becomes(46)xt+1=1+hr∗+p−εxt−bhxtyt,yt+1=1−hd+εyt+chxtyt.
Now, one can transformEpositivep in origin by using the following transformations:(47)Φt˜=xt−d+εc,Ψt˜=yt−r∗+p−εb.
From (47), (46) takes the following form:(48)Φt+1˜=1+hr∗+p−εΦt˜+d+εc−bhΦt˜+d+εcΨt˜+r∗+p−εb−d+εc,Ψt+1˜=1−hd+εΨt˜+r∗+p−εb+chΦt˜+d+εcΨt˜+r∗+p−εb−r∗+p−εb.
By Taylor series aboutΦt˜,Ψt˜=0,0, one gets(49)Φt+1˜=Φt˜−bhcd+εΨt˜−bhΦt˜Ψt˜+hpΦt˜,Ψt+1˜=chbr−εΦt˜+Ψt˜+chΦt˜Ψt˜.
Now by using following transformation(50)Φt˜Ψt˜≔−bcr−εε−rd+εbcr−εε−rd+ε11Xt˜Yt˜,linear part of (49) transforms into canonical form. In view of (50), (49) becomes(51)Xt+1˜Yt+1˜=−1003+F˜Xt˜,Yt˜,pG˜Xt˜,Yt˜,p,where(52)F˜Xt˜,Yt˜,p=−12bhXt˜2+12bhYt˜2−12hpXt˜,G˜Xt˜,Yt˜,p=12chXt˜2−12chYt˜2.
Now, for (51), center manifold McEtrival at Etrival is explored in a small nbhd of p. Therefore, McEtrival can be expressed as follows:(53)McEtrival=Xt˜,Yt˜:Yt˜=c0p+c1Xt˜2+c2Xt˜p+c3p3+OXt˜+p3.
After some manipulations, one gets(54)c0=c2=c3=0,c1=−14ch.
Finally, the map (51) restricting to McEtrival is(55)F˜Xt˜=−Xt˜+h1Xt˜2+h2Xt˜p+h3Xt˜2p+h4Xt˜p2+h5Xt˜3+OXt˜+p4,with(56)h1=−12bh,h2=−h2,h3=h4=h5=0.
In order for the map (51) to undergo flip bifurcation, the following discriminatory quantities are required to be nonzero [28, 29]:(57)∧1=∂2F˜∂Xt˜∂p+12∂F˜∂p∂2F˜∂Xt˜2|Xt˜,p=0,0,∧2=16∂3F˜∂Xt˜3+12∂2F˜∂Xt˜22|Xt˜,p=0,0.
On computation, one gets∧1=−h/2≠0, but ∧2=1/6b2h2>0, which shows that model (4) undergoes flip bifurcation if b,c,d,h,r,ε∈FB|Epositive and, in particular, stable period-2 points bifurcating from Epositive.


## 5. Numerical Simulations

Here, some simulations will be presented in order to verify the obtained results in Sections2 and 4. For instance, choose ε=1.9,h=1, then from (24), one gets d>0.2052631578947369. Furthermore, if d=0.5>0.2052631578947369, then from (23), one gets r>0.23333333333333317. But for the existence of unique positive fixed point, it is required that r>ε=1.9. So, for these numerical values and the condition on parameter r where unique positive fixed point is stable focus is r>max0.23333333333333317,1.9. Hence, if one choose parametric value r=2.1 which satisfies r>max0.23333333333333317,1.9 and b=0.4,c=0.5, then it is clear from Figure 1(a) that Epositive=4.8,0.5 is stable focus. Similarly, Figure 1(b) also shows that Epositive of (4) is stable focus. Now, from (25), one can say that Epositive of (4) is unstable focus if r<0.23333333333333317, and hence, flip bifurcation takes place if r<0.23333333333333317. For instance, if r=0.09<0.23333333333333317, then ζ1=−1.0842264752180841, but ζ2=3.084226475218084≠1or−1. Moreover, ∧1=−0.5≠0 but ∧2=0.026666666666666672>0 which implies that stable period-2 points bifurcating from Epositive. Hence, this simulation agrees with theoretical results in Section 4. The 2D bifurcation diagrams with corresponding maximum Lyapunov exponents for said parametric values are plotted in Figure 2. Finally, the trajectories associated with Figures 2(a) and 2(b) are also plotted in Figures 3(a)–3(f) that indicates (4) exhibits complex dynamics having orbits of period-4, period-6, period-8, period-12, period-17, and period-18.Figure 1
Stable node of model (4). (a) r = 2.1 with (0.24, 0.23). (b) r = 2.3 with (0.4, 0.5).
(a)(b)Figure 2
(a, b) Flip bifurcation diagram withr∈0.2,3.95 and 0.24,0.25. (c) M.L.E corresponding to (a, b).
(a)(b)(c)Figure 3
Complex dynamics of model (4). (a) r = 0.01 with (0.4, 0.5). (b) r = 0.12 with (0.4, 0.5). (c) r = 0.193 with (0.04, 0.02). (d) r = 0.2112 with (0.61, 0.61). (e) r = 0.21 with (0.0023, 0.6). (f) r = 0.221 with (0.24, 0.22).
(a)(b)(c)(d)(e)(f)


### 5.1. Fractal Dimension

It is defined by using Lyapunov exponents as follows [30, 31]:(58)DL=J+∑i=1JζJζJ,with ζ1,…,ζn being Lyapunov exponents, where J is the largest integer s.t. ∑i=1JζJ≥0 and ∑i=1J+1ζJ<0. The fractal dimension for considered model (4) becomes(59)DL=1+ζ1ζ2.Now for values ofd,ε,h,c,b, and r, two Lyapunov exponents are numerically computed. If d=0.5,ε=1.9,h=1,c=0.5,b=0.4, then ζ1=3.031748015872047 (resp., ζ1=3.024648117575002) and ζ2=−1.0317480158720471 (resp., ζ2=−1.0246481175750022 for r=0.18 (resp., r=0.192). So fractal dimension for the model (4) is(60)dL=1+3.031748015872047−1.0317480158720471=3.938457810659876forr=0.18,dL=1+3.024648117575002−1.0246481175750022=3.951889595750518forr=0.192.For above chosen parametric values, strange attractors are also plotted in Figures4(a) and 4(b) that demonstrate (4) has a complex dynamical behavior.Figure 4
Strange attractors ifr=0.18 (resp., r=0.192) with 0.8,0.2.
(a)(b)


## 5.1. Fractal Dimension

It is defined by using Lyapunov exponents as follows [30, 31]:(58)DL=J+∑i=1JζJζJ,with ζ1,…,ζn being Lyapunov exponents, where J is the largest integer s.t. ∑i=1JζJ≥0 and ∑i=1J+1ζJ<0. The fractal dimension for considered model (4) becomes(59)DL=1+ζ1ζ2.Now for values ofd,ε,h,c,b, and r, two Lyapunov exponents are numerically computed. If d=0.5,ε=1.9,h=1,c=0.5,b=0.4, then ζ1=3.031748015872047 (resp., ζ1=3.024648117575002) and ζ2=−1.0317480158720471 (resp., ζ2=−1.0246481175750022 for r=0.18 (resp., r=0.192). So fractal dimension for the model (4) is(60)dL=1+3.031748015872047−1.0317480158720471=3.938457810659876forr=0.18,dL=1+3.024648117575002−1.0246481175750022=3.951889595750518forr=0.192.For above chosen parametric values, strange attractors are also plotted in Figures4(a) and 4(b) that demonstrate (4) has a complex dynamical behavior.Figure 4
Strange attractors ifr=0.18 (resp., r=0.192) with 0.8,0.2.
(a)(b)


## 6. Chaos Control

By a state feedback control method, we will stabilize chaotic orbits at an unstable fixed point motivated from existing literature [32, 33]. By adding control force Ct to model (4), then(61)xt+1=1+hr−εxt−bhxtyt+Ct,yt+1=1−hd+εyt+chxtyt,Ct=−g1xt−d+εc−g2yt−r−εb,with g1 and g2 being feedback gains. Now, ΩC|Epositive at Epositive for the controlled system (61) is(62)ΩC|Epositive≔1−g1−bhcd+ε−g2chbr−ε1,whose auxiliary equation is(63)ζC2−TEpositiveζC+DEpositive=0,with(64)TEpositive=2−g1,DEpositive=1−g1+h2r−εd+ε+chbr−εg2.Ifζ1,2 are roots of (63), then(65)ζ1+ζ2=2−g1,(66)ζ1ζ2=1−g1+h2r−εd+ε+chbr−εg2.Here, it is noted that solution of the following equations(67)ζ1=1,(68)ζ1=−1,(69)ζ1ζ2=1,gives the lines of marginal stability, which further grantee the fact that ζ1,2<1. From (66) and (69), one gets(70)L1:−g1+chbr−εg2+h2r−εd+ε=0.From (65), (66), and (67), one gets(71)L2:chbr−εg2+h2r−εd+ε=0.From (65), (66), and (68), one gets(72)L2:2g1−chbr−εg2−h2r−εd+ε=4.Thus,L1,L2, and L3 in g1,g2-plane determines the triangular region that gives ζ1,2<1 (see Figure 5(a)), whereas Figures 5(b) and 5(c) tell about Epositive that the chaotic trajectories are stabilized. For the qualitative behavior of continuous dynamical systems, we refer the interested readers to [34–36] and references cited therein.Figure 5
Control of chaotic trajectories of (61) for d=0.5,ε=1.9,h=1,c=0.5,b=0.4,andr=0.12 with 0.4,0.5: (a) stability region in g1,g2-plane; (b, c) dynamics for xt and yt, respectively.
(a)(b)(c)


## 7. Periodic Points of Model (4)

Existence of periodic points atEtrival and Epositive is investigated in this section.Theorem 3.
Etrival is a periodic point having prime period-1.Proof.
From (4), we have(73)Fx,y≔fx,y,gx,y,where(74)fx,y=1+hr−εx−bhxy,gx,y=1−hd+εy+chxy.
From (73) along with (74), one gets the required results:(75)F|Etrival=Etrival.Theorem 4.
Epositive is a periodic point having prime period-1.Proof.
From (73) along with (74), one gets the required results:(76)F|Epositive=Epositive.Theorem 5.
Etrival is a periodic point of period-2,3,…,n.Proof.
From (73), one gets(77)F2x,y=1+hr−εfx,y−bhfx,ygx,y,1−hd+εgx,y+chfx,ygx,y⟹F2|Epositive=Epositive,F3x,y=1+hr−εf2x,y−bhf2x,yg2x,y,1−hd+εg2x,y+chf2x,yg2x,y⟹F3|Epositive=Epositive,⋮Fix,y=1+hr−εfi−1x,y−bhfi−1x,ygi−1x,y,1−hd+εgi−1x,y+chfi−1x,ygi−1x,y⟹Fi|Epositive=Epositive.
So, from (77), one gets the required statement. □Theorem 6.
Epositive is a periodic point of period-2,3,…,n.Proof.
From (77), one gets the required statement:(78)F2|Epositive=Epositive,F3|Epositive=Epositive,⋮Fi|Epositive=Epositive.


## 8. Conclusion

In this paper, we have investigated the topological classifications at fixed points, bifurcation analysis, and chaos in a model, which is depicted in (4). It is examined that ∀d,ε,h,c,b,r model has a trivial fixed point: Etrival and the unique positive fixed point: Epositive if r>ε. By existing linear theory of stability, we have studied the topological classifications at Etrival and Epositive, and conclusion is reported in Table 2. Furthermore, we have examined the existence of possible bifurcations at Etrival and Epositive and proved that model (4) does not undergo flip bifurcation if b,c,d,h,r,ε∈FB|Etrival, and hence, Etrival is degenerated with a higher codimension. We have also examined that if b,c,d,h,r,ε∈HB|Epositive, then at Epositive, model (4) does not undergo N-S bifurcation. Moreover, at fixed point Epositive, we have proved that model (4) undergoes the flip bifurcation if b,c,d,h,r,ε∈FB|Epositive. Some numerical simulations are performed not only to demonstrate obtained theoretical results but also to tell the complex behaviors in orbits of period-4, period-6, period-8, period-12, period-17, and period-18. We have also computed maximum Lyapunov exponents numerically. By the feedback control method, we have stabilized chaos existing in the considered model. Finally, existence of periodic points at Etrival and Epositive of model (4) is explored.Table 2
Behavior of model (4) at Etrival and Epositive
Fixed pointsRespective behaviorEtrivalSink if0<h<min2/ε−r;2/d+ε,whereε>r; source if h>max2/ε−r,2/d+ε,whereε>rSaddle if2/ε−r<h<2/d+ε; nonhyperbolic if h=2/ε−rorh=2/d+εEpositiveNever stable focus; unstable focus ifr>ε; nonhyperbolic if r=εStable node ifr>εh2d+ε−4/h2d+ε; unstable node if r<εh2d+ε−4/h2d+εNonhyperbolic ifr=εh2d+ε−4/h2d+ε



---

*Source: 8845926-2020-11-16.xml*