# Bifurcation Analysis and Chaos Control in a Discrete Epidemic System

Authors: Wei Tan; Jianguo Gao; Wenjun Fan
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

The dynamics of discreteSI epidemic model, which has been obtained by the forward Euler scheme, is investigated in detail. By using the center manifold theorem and bifurcation theorem in the interior R+2, the specific conditions for the existence of flip bifurcation and Neimark-Sacker bifurcation have been derived. Numerical simulation not only presents our theoretical analysis but also exhibits rich and complex dynamical behavior existing in the case of the windows of period-1, period-3, period-5, period-6, period-7, period-9, period-11, period-15, period-19, period-23, period-34, period-42, and period-53 orbits. Meanwhile, there appears the cascade of period-doubling 2, 4, 8 bifurcation and chaos sets from the fixed point. These results show the discrete model has more richer dynamics compared with the continuous model. The computations of the largest Lyapunov exponents more than 0 confirm the chaotic behaviors of the system x→x+δ[rN(1-N/K)-βxy/N-(μ+m)x], y→y+δ[βxy/N-(μ+d)y]. Specifically, the chaotic orbits at an unstable fixed point are stabilized by using the feedback control method.


---

## Body

## 1. Introduction

Within the field of epidemic theory, both the continuous-time model and the discrete-time model were described by the different equations. The continuous-time epidemic models have been widely investigated by many researchers (e.g., [1–10]). A Leslie-Gower predator-prey model with disease in prey incorporating a prey refuge was investigated in [7]. A ratio-dependent predator-prey model with Allee effect and disease in prey was investigated in [8]. Meanwhile, plenty of scholars used some discretization scheme to the continuous-time systems to study the consistency, convergence, permanence, and stability of the discrete system in [11–16]. In the recent years, there appear to be a number of articles about the flip bifurcation, fold bifurcation, and pitchfork bifurcation and in the sense of Marottos chaos of the discrete system which were presented in [17–28]. Both a discrete-time predator-prey system without Allee effect and a discrete-time predator-prey system with Allee effect were studied in [22, 23]. Dynamics of a system of three interacting populations with Allee effects and stocking was researched in [26]. Stability analysis and rich oscillation patterns in discrete-time FitzHugh-Nagumo excitable system with delayed coupling were studied in [24].In this paper by using the center manifold theorem and bifurcation theorem, the local stability and bifurcation behaviors of a simpleSI epidemic model with surprising dynamic is investigated in detail. The paper is organized as follows. Section 2 gives the specific conditions for existence of flip bifurcation and Neimark-Sacker bifurcation. Section 3 analyzes the two bifurcations by choosing the time step size δ as the parameter of the bifurcation. Section 4 verifies the theoretical analysis and displays the complex and surprising dynamics. In Section 5, chaos approaching to the unable fixed point is controlled by using the feedback control method.A simpleSI epidemic model with surprising dynamic was showed in [1] as(1)S˙t=rN1-NK-βSIN-μ+mS,I˙t=βSIN-μ+dI,where N=S+I is the total population size; the birth process incorporates density dependent effects through a logistic equation with intrinsic growth rate r and the carrying capacity K, and other parameters are positive constant.It is obvious that system (1) is a two-dimensional nonlinear ODE, where it is an impossible task to obtain its analytical solutions to study its qualitative properties including the global stability, so it is necessary to solve numerical solutions or approximate solutions of system (1) according to different discrete difference schemes. Therefore, in order to study the dynamics behavior of system (1), we will focus on the complex dynamical behaviors of a simple discrete SI epidemic model approached by the forward Euler scheme.In this paper, we apply the forward Euler scheme to system (1) to get the following form:(2)x⟶x+δrN1-NK-βxyN-μ+mx,y⟶y+δβxyN-μ+dy,where δ is the step size.


## 2. Stability of Fixed Points

The fixed points of system (2) are A(rN(1-N/K)/μ+m,0) and B(x∗,y∗), where x∗,y∗ satisfy(3)rN1-NK-βxyN-μ+mx=0,βxyN-μ+dy=0.The Jacobian matrix of (2) at any point (x,y) is written as(4)J=1-δβNy+μ+m-δβNxδβNy1+δβNx-μ-d.The characteristic equation of Jacobian matrix (4) is given by(5)Fλ=λ2+px,yλ+qx,y=0,where(6)px,y=-2+δβNy-x+2μ+d+m,qx,y=1-δβNy+μ+m1+δβNx-μ-d+δ2β2N2xy.Let us give a lemma as similar as in [21].Lemma 1.
Supposeλ1 and λ2 are two roots of F(λ)=0; the local stability of the positive point can be determined by the modules of |λ1| and |λ2|.(i)
Ifλ1⩽1 and λ2⩽1, then not only is B(X∗,Y∗) a sink, but also it is locally asymptotically stable.(ii)
Ifλ1⩾1 and λ2⩾1, then not only is B(X∗,Y∗) a source, but also it is unstable.(iii)
Ifλ1⩽1 and λ2⩾1(λ1⩾1 and λ2⩽1), then B(X∗,Y∗) is a saddle.(iv)
Ifλ1=1 or λ2=1, then B(X∗,Y∗) is nonhyperbolic.In order to discuss the stability of the fixed point of system (2), the following lemma is also needed which can be easily proved by the relation between roots and coefficient of a quadratic equation [21].Lemma 2.
LetF(λ)=λ2+pλ+q. Suppose that F(1)>0(F(1)<0); λ1 and λ2 are root of F(λ)=0. Then(i)
λ1⩽1 and λ2⩽1 if F(-1)>0 and q<1;(ii)
λ1⩽1 and λ2⩾1(λ1⩾1 and |λ2|⩽1) if and only if F(-1)<0;(iii)
λ1⩾1 and λ2⩾1 if and only if F(-1)>0 and q>1;(iv)
λ1=1 and λ2≠1 if and only if F(-1)=0 and p≠0,2;(v)
λ1=1 or λ2=1 and λ1 and λ2 are complex conjugate if and only if p2-4q<0 and q=1.Proposition 3.
The eigenvalues of the fixed pointA(rN(1-N/K)/μ+m,0) are λ1=1-δ(μ+m), λ2=1+δβr1-N/K/μ+m-μ-m.(i)
A(rN(1-N/K)/μ+m,0) is sink if (μ+m)(μ+d)>βr(1-N/k) and 0<δ<min⁡{2/μ+m,2(μ+d)(μ+m)/βr(1-N/k)}.(ii)
A(rN(1-N/K)/μ+m,0) is source if (μ+m)(μ+d)>βr(1-N/k) and δ>max⁡{2/μ+m,2(μ+d)(μ+m)/βr(1-N/k)}.(iii)
A(rN(1-N/K)/μ+m,0) is nonhyperbolic if (μ+m)(μ+d)>βr(1-N/k) and δ=2/μ+m or δ=2(μ+d)(μ+m)/βr(1-N/k).(iv)
A(rN(1-N/K)/μ+m,0) is a saddle, except for that values of parameters which (i) to (iii).Let(7)FA=r,N,K,β,μ,m,δ:μ+mμ+d>βr1-Nk,δ=2μ+m,δ≠2μ+dμ+mβr1-N/k,r,N,K,β,μ,m,δ>0.If term (iii) of Proposition 3 holds, then one of the eigenvalues of the fixed point A(rN(1-N/K)/μ+m,0) is −1 and the magnitude of the other is not equal to 1. The point A(rN(1-N/K)/μ+m,0) undergoes flip bifurcation when the parameter changes in small neighborhood of FA.The characteristic equation of the Jacobian matrixJ of system (2) at the positive B(X∗,Y∗) is written as(8)λ2+px∗,y∗λ+qx∗,y∗=0,where(9)px∗,y∗=-2-δE,qx∗,y∗=1+δE+δ2F,E=βNx∗-y∗-2μ-d-m,F=-βNy∗+μ+mβNx∗-μ-d+β2N2x∗y∗.Let(10)Fλ=λ2-2+δEλ+1+δE+δ2F.Therefore,(11)F1=δ2F,F-1=4+2δE+δ2F.According to Lemma 2, we get the following proposition.Proposition 4.
There exists the fixed pointB(x∗,y∗) for every possible parameter:(i)
B(x∗,y∗) is a sink if satisfying
(i.1)
E2-4F⩾0 and 0⩽δ<-E-E2-4F/F;(i.2)
E2-4F<0 and 0⩽δ<-E/F;(ii)
B(x∗,y∗) is a source if satisfying
(ii.1)
E2-4F⩾0 and δ>-E+E2-4F/F;(ii.2)
E2-4F<0 and δ>-E/F;(iii)
B(x∗,y∗) is nonhyperbolic if satisfying
(iii.1)
E2-4F⩾0 and δ=-E∓E2-4F/F;(iii.2)
E2-4F<0 and δ=-E/F;(iv)
B(x∗,y∗) is a saddle, except for that values of parameters which (i) to (iii). If (iii.1) holds, then λ1=-1 and λ2≠1.Let(12)FB1=r,N,K,β,μ,m,δ:E2-4F⩾0,δ=-E-E2-4FF,r,N,K,β,μ,m,δ>0,FB2=r,N,K,β,μ,m,δ:E2-4F⩾0,δ=-E+E2-4FF,r,N,K,β,μ,m,δ>0.If (iii.2) holds, then λ1 and λ2 are a pair of complex conjugate numbers and λ1=-1 and λ2=1.
Let(13)EB1=r,N,K,β,μ,m,δ:E2-4F<0,δ=-EF,r,N,K,β,μ,m,δ>0.


## 3. Bifurcation Analysis on the PositiveB(x∗,y∗)

In this section, we will investigate the flip bifurcation and Neimark-Sacker bifurcation corresponding to the positive fixed point.


### 3.1. Flip Bifurcation

Consider system (2) with arbitrary parameter r,N,K,β,μ,m,δ∈FB1, then(14)x⟶x+δ1rN1-NK-βxyN-μ+mx,y⟶y+δ1βxyN-μ+dy,where(15)δ1=-E-E2-4FF,λ1=-1,λ2≠1.Now, give a perturbation of model (14) as follows:(16)x⟶x+δ1+δ∗rN1-NK-βxyN-μ+mx,y⟶y+δ1+δ∗βxyN-μ+dy,where δ∗≪1 is a limited perturbation parameter.LetX=x∗-x, Y=y-y∗. After the transformation of the fixed B(x∗,y∗) to (0,0), as follows:(17)X⟶a11X+a12Y+a13XY+a14δ∗X+a15δ∗Y+a16δ∗XY,Y⟶a21X+a22Y+a23XY+a24δ∗X+a25δ∗Y+a26δ∗XY,where(18)a11=1-δ1βNy∗+μ+m,a12=-δ1βNx∗,a13=-δ1βN,a14=-βNy∗+μ+m,a15=-βNx∗,a16=-βN,a21=δ1βNy∗,a22=1+δ1βNx∗-μ-d,a23=βNδ1,a24=βNy∗,a25=βNx∗-μ-d,a26=βN,and δ=δ1.Let(19)T=a12a12-1-a11λ2-a11.Consider the following translation:(20)XY=Tuv.TakingT-1 on both sides of (17), obtain(21)uv=-100λ2uv+fX,Y,δ∗gX,Y,δ∗,where(22)fX,Y,δ∗=1a12λ2+1a13λ2-a11-a12a23XY+λ2-a11a14-a12a24δ∗X+λ2-a11a15-a12a25δ∗Y+λ2-a11a16-a12a26δ∗XY,gX,Y,δ∗=1a12λ2+1a131+a11+a12a23XY+1+a11a14+a12a24δ∗X+1+a11a15+a12a25δ∗Y+1+a11a16+a12a26δ∗XY,X=a12u+v,Y=-1+a11u+λ2-a11v.Apply the center manifold theorem to (17) at the origin in the limited neighborhood of δ∗=0. The center manifold Wc(0,0,0) can be approximately presented as(23)Wc0,0,0=u,v,δ∗∈R3∣v=hu,δ∗,h0,0=0,Dh0,0=0,for u, δ∗ sufficiently small.We assume a center manifold of the form(24)hu,δ∗=a1u2+a2δ∗u+0u+δ∗3,which must satisfy(25)Nhu,δ∗=h-u+fu,hu,δ∗,δ∗-λ2hu,δ∗-gu,hu,δ∗,δ∗=0.By the simple calculations, obtain(26)a1=a131+a11+a12a231+a11λ22-1,a2=-1+a11a14+a12a24a12+1+a11a15+a12a251+a11a121+λ22.Consider the map restricted to the center manifoldWc(0,0,0):(27)f:⟶-u+h1a12λ2+1u2+h2a12λ2+1uδ∗+h3a12λ2+1u3+h4a12λ2+1δ∗u2+h5a12λ2+1δ∗2u,where(28)h1=-a13λ2-a11-a12a23a121+a11,h2=λ2-a11a14-a12a24a12-λ2-a11a15-a12a251+a11,h3=a13λ2-a11-a12a23a12a1λ2-2a11-1,h4=a13λ2-a11-a12a23a12a2λ2-2a11-1+λ2-a11a14-a12a24a12a1+λ2-a11a15-a12a25a1λ2-a11-λ2-a11a16-a12a26a121+a11,h5=λ2-a11a14-a12a24a12a2+λ2-a11a15-a12a25λ2-a11a2,L=a12λ2+1.According to flip bifurcation, we get(29)l1=∂2f∂u∂δ∗+12∂f∂δ∗∂2f∂u2∣0,0,l2=16∂3f∂u3+12∂2f∂u22∣0,0.After simple calculations, we obtain(30)l1=h2L+h12L2,l2=h3L+h12L2.Theorem 5.
Ifα2≠0, system (2) will undergo flip bifurcation at the positive point B(x∗,y∗). Also the period-2 orbits that bifurcate from fixed B(x∗,y∗) are stable (resp., unstable) if α2>0 (resp., α2<0).


### 3.2. Neimark-Sacker Bifurcation

Consider model (2) with arbitrary parameter (r,N,K,β,μ,m,δ)∈EB1, which is described by(31)x⟶x+δ2rN1-NK-βxyN-μ+mx,y⟶y+δ2βxyN-μ+dy.B(x∗,y∗) is the fixed point of model (2), which is given by (3) and δ2=-E/F.We consider the perturbation of model (31) as follows:(32)x⟶x+δ2+δ⋆rN1-NK-βxyN-μ+mx,y⟶y+δ2+δ⋆βxyN-μ+dy,where δ⋆≪1 is a limited perturbation parameter.LetX=x-x∗, Y=y-y∗. After the transformation of the fixed B(x∗,y∗) to (0,0), as follows:(33)X⟶b11X+b12Y+b13XY,Y⟶b21X+b22Y+b23XY,where(34)b11=1-δ2+δ⋆βNy∗+μ+m,b12=-δ2+δ⋆βNx∗,b13=-δ2+δ⋆βN,b21=δ2+δ⋆βNy∗,b22=1+δ2+δ⋆βNx∗-μ-d,b23=δ2+δ⋆βN,and δ=δ2.The characteristic equation of model (31) at B(x∗,y∗) is as follows:(35)λ2+pδ⋆λ+qδ⋆=0.Since the parameter(r,N,K,β,μ,m,δ2)∈EB1, the characteristic values of B(x∗,y∗) are a pair of complex conjugate numbers λ3 and λ4 with modulus 1:(36)λ3,4=1+Eδ2+δ⋆2∓iδ2+δ⋆4F-E22.Now, we have(37)λ3,4=qδ⋆,l=dλdδ⋆=-E2>0.Whenδ⋆ varies in limited neighborhood of δ⋆=0, let θ=Re⁡(λ)=1+δ2E/2, ω=Im⁡(λ)=δ24F-E2/2. Neimark-Sacker bifurcation requires that when δ⋆=0, then λ¯n,λn≠1 (n=1,2,3,4), which is equivalent to p(0)≠0,2,-2,1. Since p(0)2-4q(0)<0 and q(0)=1, therefore p(0)2<4; then p(0)≠±2, which follows that(38)E2≠2F,3F.Let(39)T=b120θ-b11-ω.Consider the translation below:(40)XY=Tuv.TakingT-1 on both sides of (19), obtain(41)uv=θ-ωωθuv+fu,vgu,v,where(42)fu,v=b13θ-b11u2-ωuv,gu,v=-1ωb13b11-θ+b12b23uθ-b11u-ωv,X=b12u,Y=θ-b11u-ωv.Now(43)fuu=2b13θ-b11,fuv=-b13ω,fvv=0,fuuu=fuuv=fuvv=fvvv=0,guu=-2ωb13b11-θ+b12b23θ-b11,guv=b12b11-θ+b12b23,gvv=0,guuu=guuv=guvv=gvvv=0.According to Neimark-Sacker bifurcation, the discriminatory quantities are given by(44)μ2=-Re⁡1-2λλ¯21-λΦ11Φ20-12Φ112-Φ022+Re⁡λ¯Φ21,where(45)Φ20=18fuu-fvv+2guv+iguu-gvv-2fuv,Φ11=14fuu+fvv+iguu+gvv,Φ02=18fuu-fvv-2guv+iguu-gvv+2fuv,Φ21=116fuuu+fuvv+guuv+gvvv+iguuu+guvv-fuuv-fvvv.Then we can conclude the following theorem.Theorem 6.
If condition (39) holds, μ2≠0, when the parameter δ⋆ alters around the limit region of point (0,0), system (2) undergoes a Neimark-Sacker bifurcation at the point B(X∗,Y∗). If μ2<0, an attracting invariant closed curve bifurcates from the point B(X∗,Y∗) for δ⋆>0.


## 3.1. Flip Bifurcation

Consider system (2) with arbitrary parameter r,N,K,β,μ,m,δ∈FB1, then(14)x⟶x+δ1rN1-NK-βxyN-μ+mx,y⟶y+δ1βxyN-μ+dy,where(15)δ1=-E-E2-4FF,λ1=-1,λ2≠1.Now, give a perturbation of model (14) as follows:(16)x⟶x+δ1+δ∗rN1-NK-βxyN-μ+mx,y⟶y+δ1+δ∗βxyN-μ+dy,where δ∗≪1 is a limited perturbation parameter.LetX=x∗-x, Y=y-y∗. After the transformation of the fixed B(x∗,y∗) to (0,0), as follows:(17)X⟶a11X+a12Y+a13XY+a14δ∗X+a15δ∗Y+a16δ∗XY,Y⟶a21X+a22Y+a23XY+a24δ∗X+a25δ∗Y+a26δ∗XY,where(18)a11=1-δ1βNy∗+μ+m,a12=-δ1βNx∗,a13=-δ1βN,a14=-βNy∗+μ+m,a15=-βNx∗,a16=-βN,a21=δ1βNy∗,a22=1+δ1βNx∗-μ-d,a23=βNδ1,a24=βNy∗,a25=βNx∗-μ-d,a26=βN,and δ=δ1.Let(19)T=a12a12-1-a11λ2-a11.Consider the following translation:(20)XY=Tuv.TakingT-1 on both sides of (17), obtain(21)uv=-100λ2uv+fX,Y,δ∗gX,Y,δ∗,where(22)fX,Y,δ∗=1a12λ2+1a13λ2-a11-a12a23XY+λ2-a11a14-a12a24δ∗X+λ2-a11a15-a12a25δ∗Y+λ2-a11a16-a12a26δ∗XY,gX,Y,δ∗=1a12λ2+1a131+a11+a12a23XY+1+a11a14+a12a24δ∗X+1+a11a15+a12a25δ∗Y+1+a11a16+a12a26δ∗XY,X=a12u+v,Y=-1+a11u+λ2-a11v.Apply the center manifold theorem to (17) at the origin in the limited neighborhood of δ∗=0. The center manifold Wc(0,0,0) can be approximately presented as(23)Wc0,0,0=u,v,δ∗∈R3∣v=hu,δ∗,h0,0=0,Dh0,0=0,for u, δ∗ sufficiently small.We assume a center manifold of the form(24)hu,δ∗=a1u2+a2δ∗u+0u+δ∗3,which must satisfy(25)Nhu,δ∗=h-u+fu,hu,δ∗,δ∗-λ2hu,δ∗-gu,hu,δ∗,δ∗=0.By the simple calculations, obtain(26)a1=a131+a11+a12a231+a11λ22-1,a2=-1+a11a14+a12a24a12+1+a11a15+a12a251+a11a121+λ22.Consider the map restricted to the center manifoldWc(0,0,0):(27)f:⟶-u+h1a12λ2+1u2+h2a12λ2+1uδ∗+h3a12λ2+1u3+h4a12λ2+1δ∗u2+h5a12λ2+1δ∗2u,where(28)h1=-a13λ2-a11-a12a23a121+a11,h2=λ2-a11a14-a12a24a12-λ2-a11a15-a12a251+a11,h3=a13λ2-a11-a12a23a12a1λ2-2a11-1,h4=a13λ2-a11-a12a23a12a2λ2-2a11-1+λ2-a11a14-a12a24a12a1+λ2-a11a15-a12a25a1λ2-a11-λ2-a11a16-a12a26a121+a11,h5=λ2-a11a14-a12a24a12a2+λ2-a11a15-a12a25λ2-a11a2,L=a12λ2+1.According to flip bifurcation, we get(29)l1=∂2f∂u∂δ∗+12∂f∂δ∗∂2f∂u2∣0,0,l2=16∂3f∂u3+12∂2f∂u22∣0,0.After simple calculations, we obtain(30)l1=h2L+h12L2,l2=h3L+h12L2.Theorem 5.
Ifα2≠0, system (2) will undergo flip bifurcation at the positive point B(x∗,y∗). Also the period-2 orbits that bifurcate from fixed B(x∗,y∗) are stable (resp., unstable) if α2>0 (resp., α2<0).


## 3.2. Neimark-Sacker Bifurcation

Consider model (2) with arbitrary parameter (r,N,K,β,μ,m,δ)∈EB1, which is described by(31)x⟶x+δ2rN1-NK-βxyN-μ+mx,y⟶y+δ2βxyN-μ+dy.B(x∗,y∗) is the fixed point of model (2), which is given by (3) and δ2=-E/F.We consider the perturbation of model (31) as follows:(32)x⟶x+δ2+δ⋆rN1-NK-βxyN-μ+mx,y⟶y+δ2+δ⋆βxyN-μ+dy,where δ⋆≪1 is a limited perturbation parameter.LetX=x-x∗, Y=y-y∗. After the transformation of the fixed B(x∗,y∗) to (0,0), as follows:(33)X⟶b11X+b12Y+b13XY,Y⟶b21X+b22Y+b23XY,where(34)b11=1-δ2+δ⋆βNy∗+μ+m,b12=-δ2+δ⋆βNx∗,b13=-δ2+δ⋆βN,b21=δ2+δ⋆βNy∗,b22=1+δ2+δ⋆βNx∗-μ-d,b23=δ2+δ⋆βN,and δ=δ2.The characteristic equation of model (31) at B(x∗,y∗) is as follows:(35)λ2+pδ⋆λ+qδ⋆=0.Since the parameter(r,N,K,β,μ,m,δ2)∈EB1, the characteristic values of B(x∗,y∗) are a pair of complex conjugate numbers λ3 and λ4 with modulus 1:(36)λ3,4=1+Eδ2+δ⋆2∓iδ2+δ⋆4F-E22.Now, we have(37)λ3,4=qδ⋆,l=dλdδ⋆=-E2>0.Whenδ⋆ varies in limited neighborhood of δ⋆=0, let θ=Re⁡(λ)=1+δ2E/2, ω=Im⁡(λ)=δ24F-E2/2. Neimark-Sacker bifurcation requires that when δ⋆=0, then λ¯n,λn≠1 (n=1,2,3,4), which is equivalent to p(0)≠0,2,-2,1. Since p(0)2-4q(0)<0 and q(0)=1, therefore p(0)2<4; then p(0)≠±2, which follows that(38)E2≠2F,3F.Let(39)T=b120θ-b11-ω.Consider the translation below:(40)XY=Tuv.TakingT-1 on both sides of (19), obtain(41)uv=θ-ωωθuv+fu,vgu,v,where(42)fu,v=b13θ-b11u2-ωuv,gu,v=-1ωb13b11-θ+b12b23uθ-b11u-ωv,X=b12u,Y=θ-b11u-ωv.Now(43)fuu=2b13θ-b11,fuv=-b13ω,fvv=0,fuuu=fuuv=fuvv=fvvv=0,guu=-2ωb13b11-θ+b12b23θ-b11,guv=b12b11-θ+b12b23,gvv=0,guuu=guuv=guvv=gvvv=0.According to Neimark-Sacker bifurcation, the discriminatory quantities are given by(44)μ2=-Re⁡1-2λλ¯21-λΦ11Φ20-12Φ112-Φ022+Re⁡λ¯Φ21,where(45)Φ20=18fuu-fvv+2guv+iguu-gvv-2fuv,Φ11=14fuu+fvv+iguu+gvv,Φ02=18fuu-fvv-2guv+iguu-gvv+2fuv,Φ21=116fuuu+fuvv+guuv+gvvv+iguuu+guvv-fuuv-fvvv.Then we can conclude the following theorem.Theorem 6.
If condition (39) holds, μ2≠0, when the parameter δ⋆ alters around the limit region of point (0,0), system (2) undergoes a Neimark-Sacker bifurcation at the point B(X∗,Y∗). If μ2<0, an attracting invariant closed curve bifurcates from the point B(X∗,Y∗) for δ⋆>0.


## 4.  Numerical Simulation

In this section, we will give the numerical simulations for model (2), which not only illustrate our theoretical analysis but also present the complex, rich, and surprising dynamical behavior. The bifurcation parameters are considered for the following two cases:(1)
varyingδ in the range 4<δ<7 and fixing β=0.1, r=0.03, d=0.0004, μ=0.006, m=0.0004, N=10, and K=500;(2)
varyingδ in the range 2.8<δ<3.6 and fixing β=0.45, r=6, d=0.6, μ=0.1, m=0.3, N=1, and K=1.25.Case 1.
For Case1  β=0.1, r=0.03, d=0.0004, μ=0.006, m=0.0004, N=10, and K=500 with fixed point (0.64,45.2975) at the initial value point (1,1). By a simple calculation, we get δ1=4.438, α1≠0, α2=1.01298931>0, which established Theorem 5. From bifurcation diagram (see Figure 1(c)) and the largest Lyapunov exponent of model (2) (see Figure 1(a)), the fixed point of system (2) is stable when δ1<4.438 and unstable with δ1=4.438 and there exists flip bifurcation when δ1>4.438. Largest Lyapunov exponents (Figure 1(a)) correspond to Figure 1(c).Figure 1
(a) Largest Lyapunov exponents corresponding to (c). (c) Flip bifurcation diagrams of system (2) in (δ,x) at the fixed point B(0.64,45.2975), and not only the initial value point (1,1), but also β=0.1, r=0.03, d=0.0004, μ=0.006, m=0.0004, N=10, and K=500; meanwhile, α2=1.0129893>0. (b) Largest Lyapunov exponents corresponding to (d). (d) Local amplification of (b) for δ∈(5.8,6.3).
(a)
(b)
(c)
(d)The phase portraits which are associated with Figures1(a) and 1(c) are displayed in Figure 3. There are the orbits of periods 2, 4, and 8 with δ1∈(4.438,5.675). Moreover, when δ1∈5.675,7, stable period windows and chaotic sets alternate in appearing. For example, there emerge the stable window of period-9, period-7, period-3, period-5, period-6, period-1 orbits. The phase portraits which are associated with Figures 1(b) and 1(d) are displayed in Figure 3. With δ1 increasing, chaos and period-doubling bifurcation take place by turn. Finally, the system enters into chaos.Case 2.
For Case2  β=0.45, r=6, d=0.6, μ=0.1, m=0.3, N=1, and K=1.25. The initial value is (1,1), after a simple calculation for the fixed point (1.55555556,0.82539683), we get δ2=1.43694715, λ3,4=0.44574894∓i0.89515802, θ=0.44574894, ω=0.89515802, -E/2=0.3857143>0, Δ2=E2-4F=-1.55230772<0, E2≠2F,3F, μ2<0 (see Figure 2(b)); it established Theorem 6. From Figures 2(a) and 2(c), we observe that when δ2<2.96, the fixed point of map (2) is stable, it loses its stability at δ2=2.96, and there appears an invariant circle when δ2 exceeds 2.96. Figures 2(b) and 2(d) are local amplifications of Figures 2(a) and 2(c) respectively for δ2∈(3.25,3.4).Figure 2
(a) Largest Lyapunov exponents corresponding to (c). (c) Neimark-Sacker bifurcation diagrams of system (2) in (δ-x) plane at the point B(1.55555556,0.82539683), and not only the initial value point (1,1), but also β=0.45, r=6, d=0.6, μ=0.1, m=0.3, N=1, and K=1.25. (b) Largest Lyapunov exponents corresponding to (d). (d) Local amplification of (b) for δ∈(3.25,3.37).
(a)
(b)
(c)
(d)Figure 3
Phase portrait of system (2) at the initial value (1,1); β=0.1, r=0.03, d=0.0004, μ=0.006, m=0.0004, δ=4.3, δ=5, δ=5.5, δ=5.66, δ=5.75, δ=5.8, δ=5.9, δ=6.125, δ=6.15, δ=6.25, δ=6.285, and δ=6.35.The phase portraits which are associated with Figures1(a) and 1(c) are displayed in Figure 4. Lyapunov exponents are smaller than 0 when δ∈(0.25,0.296); that is to say, the nonchaotic region is smaller than the chaotic region (0.296,0.34). For δ2∈(0.296,0.34), some Lyapunov exponents are larger than 0 and some Lyapunov exponents are smaller than 0, so there exists stable fixed point or stable period windows in the chaotic region. As we all know, Lyapunov exponents are bigger than 0 which implies the existence of chaos. With δ2 increases, circle disappears and period-53, period-34, period-19, period-23, period-19, period-42, period-15, and period-11 orbits appear, attracting chaotic sets (see Figure 4). From Figure 4, we can see that the attractive cycle is smooth in the vicinity of the bifurcation point. However, when δ2 decreases, the cycle inflates gradually and loses its smoothness, and the system enters chaos. This shows that not only the stability changes of system but also the limit cycle occur when system crosses the bifurcation point.Figure 4
Neimark-Sacker bifurcation phase portrait of system (2) with the initial value point (1,1) and δ=2.955, δ=2.967, δ=2.97, δ=3.15, δ=3.2, δ=3.261, δ=3.2685, δ=3.28, δ=3.33, δ=3.379, δ=3.391, δ=3.429, δ=3.54, δ=3.5485, and δ=3.56, respectively.


## 5. Chaos Control

In this section, a feedback control law as the control forceun is added to (2): the following controlled form of model (2) is(46)xn+1=xn+δrN1-NK-βxnynN-μ+mxn+un,yn+1=yn+δβxnynN-μ+dyn,(47)un=-k1xn-x∗-k2yn-y∗,where k1 and k2 are the feedback gain and (x∗,y∗) is the positive point of model (2). The Jacobian matrix J of controlled system (46) is given by(48)J=a11-k1a12-k2a21a22,where a11,a12,a21,a22 are given in (17). The characteristic equation of (48) is(49)λ2-a11+a22-k1λ+a22a11-k1-a21a12-k2=0.Let λ1 and λ2 be the two roots of (49); then(50)λ1+λ2=a11+a22-k1,(51)λ1λ2=a22a11-k1-a21a12-k2.The lines of marginal stability are determined by solving the equation of λ1=±1 and λ1λ2=1.These conditions guarantee that the eigenvaluesλ1 and λ2 have modulus less than 1.Assume thatλ1λ2=1, and from (51) we havel1:k1a22-k2a21=a11a22-a12a21-1.Assume thatλ1=1, and from (50) and (51) we havel2:k1(1-a22)-k2a21=a11+a22-1-a11a22+a12a21.Assume thatλ1=-1, and from (50) and (51) we havel3:k1(1+a22)-k2a21=a11+a22+1+a11a22-a12a21.The stable eigenvalues lie within a triangular region by lines l1,l2, and l3 (see Figure 5(a)).Figure 5
(a) The bounded region for the eigenvalues of control model (48) and (49) in the (k1,k2) plane for β=0.045, r=6, d=0.6, μ=0.1, m=0.3, N=1, K=1.25, and δ=3.55. (b) The time responses for the state x of control model (48) and (49) in the (n,x) plane for β=0.045, r=6, d=0.6, μ=0.1, m=0.3, N=1, K=1.25, δ=3.556, k1=0.1, and k2=-0.3415. The initial value point (1,1). (c) The time responses for the state y of control model (48) and (49) in the (n,y) plane.
(a)
(b)
(c)We have performed some numerical simulations to see how the state feedback method controls the unstable fixed point. Parameter values are fixed asβ=0.45, r=6, d=0.6, μ=0.1, m=0.3, N=1, K=1.25, and δ=3.55. The initial value is (1,1), and the feedback gain k1=0.2, k2=-0.34. It is shown in Figures 5(b) and 5(c) that a chaotic trajectory is stabilized at the fixed point (1.55555556,0.82539683).


## 6. Discussion

The dynamics of discreteSI epidemic model, which has been obtained by the forward Euler scheme, is investigated in detail. By using the center manifold theorem and bifurcation theorem in the closed first quadrant R+2, the specific conditions for the existence of flip bifurcation and Neimark-Sacker bifurcation have been derived. Moreover, numerical simulation not only presents our theoretical analysis but also exhibits rich and complex dynamical behavior existing in the case of the windows of period-1, period-3, period-5, period-6, period-7, period-9, period-11, period-15, period-19, period-23, period-34, period-42, and period-53 orbits. Meanwhile, there appears the cascade of period-doubling 2, 4, 8 bifurcation and chaos sets from the fixed point. These results show that the discrete model has more richer dynamics compared with the continuous model. The computations of the largest Lyapunov exponents more than 0 confirm the chaotic behaviors of system (2). Specifically, the chaotic orbits at an unstable fixed point are stabilized by using the feedback control method. However, in this paper, respect to the interpretation of these results with sound meaning of general disease dynamics and implications for disease control still need further researches. In our future work, we expect to obtain some more results based on real data from known epidemic disease to illustrate the validity of our theoretic results, such as how to predict the occurrence of disease in which way study global stability, limit cycle, bifurcation, chaos, strange attractors, hidden attractor impact on dynamics of disease, and so on.



---

*Source: 974868-2015-08-31.xml*