# Control of Hopf Bifurcation and Chaos in a Delayed Lotka-Volterra Predator-Prey System with Time-Delayed Feedbacks

Authors: Huitao Zhao; Yaowei Sun; Zhen Wang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A delayed Lotka-Volterra predator-prey system with time delayed feedback is studied by using the theory of functional differential equation and Hassard’s method. By choosing appropriate control parameter, we investigate the existence of Hopf bifurcation. An explicit algorithm is given to determine the directions and stabilities of the bifurcating periodic solutions. We find that these control laws can be applied to control Hopf bifurcation and chaotic attractor. Finally, some numerical simulations are given to illustrate the effectiveness of the results found.


---

## Body

## 1. Introduction

Lotka-Volterra system is one of the most classical and important systems in the field of mathematical biology. Since the word of Volterra, there have been extensively detailed investigations on Lotka-Volterra system including stability, attractivity, persistence, periodic oscillation, bifurcation and chaos (see [1–6] and the references therein). In particular, the properties of periodic solutions arising from the Hopf bifurcation are of great interest [7–10]. But the study on chaos control of Lotka-Volterra system is scarce.Reference [3] and the references therein proposed that, for a two-species competition system with delays
(1)u˙i(t)=ui(t)[ri+∑j=1n‍aijuj(t-τij)],i=1,2,
when τij is big enough, the chaotic behavior may occur. For example, Yan and Zhang [9] investigated the following delayed prey-predator system with a single delay:
(2)x˙(t)=x(t)[r1-a11x(t-τ)-a12y(t-τ)],y˙(t)=y(t)[-r2+a21x(t-τ)-a22y(t-τ)],
where r1, r2, a11, a12, a21, and a22 are all positive constants. The delay τ≥0 denotes the gestation period of the predator. Their results show that, taking τ as the bifurcation parameter, when τ passes through a certain critical value, the positive equilibrium loses its stability and Hopf bifurcation takes place. Furthermore, when τ takes a sequence of critical values containing the above critical value, the positive equilibrium of system (2) will undergo a Hopf bifurcation. With the further increase of the delay, the system will show the chaotic phenomenon (see Figure 1).Figure 1
Waveform plot and phase plot of system (2) with τ=1.4.In the sense of biology, chaotic behavior sometimes is to the disadvantage of virtuous cycle and develop of the ecosystem, so we want to control this chaos phenomenon and create periodic orbits. So far, many researchers have proposed chaos control schemes in recent years [11–16]. For example, Song and Wei in [17] investigated the chaos phenomena of Chen’s system using the method of delayed feedback control. Their results show that, when the controlling parameter K to be some value, taking the delay τ as the bifurcation parameter, when τ passes through a certain critical value, the stability of the equilibrium will be changed from unstable to stable, chaos vanishes, and a periodic solution emerges.To the end of controlling chaos in system (2), stimulating by the works of above, we add some delayed feedback terms to system (2), that is, the following delayed feedback control system:
(3)x˙(t)=x(t)[r1-a11x(t-τ)-a12y(t-τ)]+k1(x(t)-x(t-τ)),y˙(t)=y(t)[-r2+a21x(t-τ)-a22y(t-τ)]+k2(y(t)-y(t-τ)),
where ki(i=1,2) denote the capture coefficient when ki<0 (or release coefficient when ki>0). By choosing τ and ki as bifurcation parameter, we get the conditions under which Hopf bifurcation occurs. And then, we derive the explicit algorithm for determining the direction of the Hopf bifurcation and stability of the bifurcating periodic solutions. At last we will give some example showing that when ki is fixed, with τ increasing, the stability of the positive equilibrium will be changed, chaos vanishes, a periodic solution occurs.This paper is organized as follows. In Section2, we first focus on the stability and Hopf bifurcation of the positive equilibrium. In Section 3, we derive the direction and stability of Hopf bifurcation by using normal form and central manifold theory. Finally in Section 4, numerical simulations are performed to support the stability results.


## 2. Stability and Hopf Bifurcation Analysis with Delayed Feedback Control

In this section, by analyzing the characteristic equation of the linearized system of system (3) at the positive equilibrium, we investigate the stability of the positive equilibrium and the existence of the local Hopf bifurcations occurring at the positive equilibrium. To guarantee that system (3) has always a positive equilibrium, throughout this section, we assume that the coefficients of system (3) satisfy the following condition:(H1)
a21r1-a11r2>0.Clearly, under the hypothesis (H1), system (3) has a unique positive equilibrium E*(x*,y*), where
(4)x*=a22r1+a12r2a11a22+a12a21,y*=a21r1-a11r2a11a22+a12a21.
Let u1(t)=x(t)-x*, u2(t)=y(t)-y*; then system (3) can be rewritten as the following equivalent system:
(5)u˙1(t)=(u1(t)+x*)[-a11u1(t-τ)-a12u2(t-τ)]+k1(u1(t)-u1(t-τ)),u˙2(t)=(u2(t)+y*)[a21u1(t-τ)-a22u2(t-τ)]+k2(u2(t)-u2(t-τ)).
Thus, the positive equilibrium E*(x*,y*) of system (3) is transformed into the equilibrium (0,0) of system (5). Linearizing system (5) at the equilibrium (0,0) yields the following linear system:
(6)u˙1(t)=-a11x*u1(t-τ)-a12x*u2(t-τ)+k1(u1(t)-u1(t-τ)),u˙2(t)=a21y*u1(t-τ)-a22y*u2(t-τ)+k2(u2(t)-u2(t-τ)).
The characteristic equation of system (6) is
(7)λ2+a1λ+a2+a3e-λτ+a4λe-λτ+a5e-2λτ=0,
where a1=-(k1+k2), a2=k1k2, a3=-(k1a22y*+k2a11x*+2k1k2), a4=a11x*+a22y*+k1+k2, a5=a11a22x*y*+a12a21x*y*+k1a22y*+k2a11x*+k1k2. Multiplying eλτ by both sides of (8), we have
(8)(λ2+a1λ+a2)eλτ+a3+a4λ+a5e-λτ=0.
Thus, iω(ω>0) is a root of (8) if and only if ω satisfies the following equation:
(9)(-ω2+a1iω+a2)(cos⁡ωτ+isinωτ)+a3+a4iω+a5(cos⁡ωτ-isinωτ)=0.
Separating the real and imaginary parts, we have
(10)(ω2-a2-a5)cos⁡ωτ+a1ωsinωτ=a3,(ω2-a2+a5)sinωτ-a1ωcos⁡ωτ=a4ω.
It follows that
(11)sinωτ=d3ω3+d4ωω4+d1ω2+d2,cos⁡ωτ=d5ω2+d6ω4+d1ω2+d2,
where d1=a12-2a2, d2=a22-a52, d3=a4, d4=a1a3-a2a4-a4a5, d5=a3-a1a4, d6=a3a5-a2a3. From sin2ωτ+cos⁡2ωτ=1, we have
(12)ω8+e1ω6+e2ω4+e3ω2+e4=0,
where e1=2d1-d32, e2=2d2+d12-d52-2d3d4, e3=2d1d2-2d5d6-d42, e4=d22-d62.Denotez=ω2; (12) becomes
(13)z4+e1z3+e2z2+e3z+e4=0.
Let
(14)G(z)=z4+e1z3+e2z2+e3z+e4.
Since lim⁡z→∞⁡G(z)=+∞, we conclude that if e4<0, then (13) has at least one positive root. If all the parameters of system (5) are given, it is easy to get the roots of (13) by using a computer. Suppose the following.(H2)
Equation (13) has at least one positive real root. Without loss of generality, we assume that (14) has four positive roots, defined by z1, z2, z3, z4, respectively. Then (12) has four positive roots
(15)ω1=z1,ω2=z2,ω3=z3,ω4=z4.From (11), we have
(16)cos⁡ωkτ=d5ωk2+d6ωk4+d1ωk2+d2.
If we denote
(17)τk(j)=1ωk{arccos(d5ωk2+d6ωk4+d1ωk2+d2)+2jπ},hhhhhhhhhhhhhhhhhhhhhhhhj=0,1,2,…,
where k=1,…,4;j=0,1,…, then ±iωk is a pair of purely imaginary roots of (8) with τ=τk(j).Define(18)τ0=τk0(0)=min⁡k=1,…,4{τk(0)},ω0=ωk0.
Let λ(τ)=α(τ)+iω(τ) be the root of (8) near τ=τk(j) satisfying α(τk(j))=0, ω(τk(j))=ωk. Substituting λ(τ) into (8) and taking the derivative with respect to τ, we have
(19)(2λ+a1)eλτdλdτ+(λ2+a1λ+a2)eλτ(λ+τdλdτ)+a4dλdτ-a5e-λτ(λ+τdλdτ)=0.
It follows that
(20)dλdτ=a5λe-λτ-(λ2+a1λ+a2)λeλτ(2λ+a1)eλτ+(λ2+a1λ+a2)τeλτ-a5τe-λτ+a4.
Then
(21)[dλdτ]-1=(2λ+a1)eλτ+a4a5λe-λτ-(λ2+a1λ+a2)λeλτ-τλ.Through tedious calculating, we can get(22)[dλdτ]τ=τk(j)-1=E+iFC+iD-τk(j)iωk,
where
(23)C=a5ωksinωkτk(j)-(ωk2-a2)ωksinωkτk(j)+a1ωk2cos⁡ωkτk(j),D=a5ωkcos⁡ωkτk(j)+(ωk2-a2)ωkcos⁡ωkτk(j)+a1ωk2sinωkτk(j),E=a1cos⁡ωkτk(j)-2ωksinωkτk(j)+a4,F=a1sinωkτk(j)+2ωkcos⁡ωkτk(j).
Denote Q=C2+D2; then Q>0, and we have
(24)QRe[dλdτ]τ=τk(j)-1=EC+DF,
and note that
(25)sign⁡{Re[dλdτ]τ=τk(j)}=sign⁡{Re[dλdτ]τ=τk(j)-1}=sign⁡{QRe[dλdτ]τ=τk(j)-1}.In order to get the main results, it is necessary to make the following assumption.(H3)
EC+DF≠0. Then, if (H3) holds, the transversality conditions Re[dλ/dτ]τ=τ0≠0 hold.Note that whenτ=0, (8) becomes
(26)λ2+(a1+a4)λ+a2+a3+a5=0.
From [9], we know that all the roots of (26) have negative real parts; hence, the positive equilibrium E* is locally asymptotically stable for τ=0. Then, we can employ a result of Ruan and Wei [18] to analyze (8). For the convenience of the reader, we state it as follows.Lemma 1.
Consider the exponential polynomial(27)P(λ,e-λτ1,…,e-λτm)=λn+p1(0)λn-1+⋯+pn-1(0)λ+pn(0)+[p1(1)λn-1+⋯+pn-1(1)λ+pn(1)]e-λτ1+⋯+[p1(m)λn-1+⋯+pn-1(m)λ+pn(m)]e-λτm,
where τi⩾0(i=1,2,…,m) and pj(i)(i=0,1,…,m;j=1,2,…,n) are constants. As (τ1,τ2,…,τm) vary, the sum of the order of the zeroes of P(λ,e-λτ1,…,e-λτm) on the open right half plane can change only if a zero appears on or crosses the imaginary axis.From Lemma1 and the above assumption, we can obtain the following theorem.Theorem 2.
Suppose that(H2) and (H3) hold; then the following results hold true.(i)
The positive equilibriumE*(x*,y*) of system (3) is asymptotically stable for τ∈[0,τ0);(ii)
system (3) exhibits Hopf bifurcation at the positive equilibrium for τ=τ0.


## 3. Direction and Stability of the Hopf Bifurcation

In this section, we obtain the conditions under which a family periodic solutions bifurcate form the steady state at the critical value ofτ. As pointed out by Hale and Verduyn Lunel [19] and Hassard et al. [20], it is interesting to determine the direction, stability, and period of these periodic solutions bifurcating from the steady state. Following the ideal of [20], we derive the explicit formulae for determining the properties of the Hopf bifurcation at the critical value of τ using the normal form and the center manifold theory.For the sake of simplicity of notation, we denote the critical valuesτj(k) as τk, and when τ=τk, we denote the pair of purely imaginary roots of (8) as ±iωk. Let μ=τ-τk; then μ=0 is the Hopf bifurcation value of system (5). In the following, we consider the equivalent system (6). Let t=τt; then the system (5) can be rewritten as a functional differential equation in C([-1,0],R2):
(28)x˙(t)=Lμ(xt)+f(μ,xt),
where x=(u1,u2)T, For ϕ=(ϕ1,ϕ2)T∈C([-1,0],R2),
(29)Lμ(ϕ)=(τk+μ)[k100k2][ϕ1(0)ϕ2(0)]+(τk+μ)[-a11x*-k1-a12x*a21y*-a22y*-k2][ϕ1(-1)ϕ2(-1)],(30)f(μ,ϕ)=(τk+μ)[-a11ϕ1(0)ϕ1(-1)-a12ϕ1(0)ϕ2(-1)a21ϕ2(0)ϕ1(-1)-a22ϕ2(0)ϕ2(-1)].Obviously,L(μ) is a continuous linear function mapping C([-1,0],R2) into R2. By the Riesz representation theorem, there exists a 2×2 matrix function η(θ,μ)(-1⩽θ⩽0), whose elements are of bounded variation such that
(31)Lμϕ=∫-10‍dη(θ,0)ϕ(θ),forϕ∈C([-1,0],R2).
In fact, we can choose
(32)η(θ,μ)=(τk+μ)[k100k2]δ(θ)-(τk+μ)[-a11x*-k1-a12x*a21y*-a22y*-k2]δ(θ+1),
where δ denotes Dirac-delta function. For ϕ∈C([-1,0],R2), define
(33)A(μ)ϕ={dϕ(θ)dθ,θ∈[-1,0),∫-10‍dη(s,μ)ϕ(s),θ=0,R(μ)ϕ={0,θ∈[-1,0),f(μ,ϕ),θ=0.
Then when θ=0, the system
(34)x˙t=A(μ)xt+R(μ)xt
is equivalent to the system (29), where xt(θ)=x(t+θ), θ∈[-1,0]. For ψ∈C1([0,1],(R2)*), define
(35)A*ψ(s)={-dψ(s)ds,s∈(0,1],∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(36)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0), and let A=A(0); then A and A* are adjoint operators. By the discussion in Section 2, we know that ±iωkτk are eigenvalues of A. Thus, they are also eigenvalues of A*. We first need to compute the eigenvector of A and A* corresponding to iωkτk and -iωkτk, respectively.Suppose thatq(θ)=(1,α)Teiωkτkθ is the eigenvector of A corresponding to iωkτk. Then Aq(θ)=iωkτkq(θ). It follows from the definition of A, Lμϕ, and η(θ,μ) that(37)τk[iωk-k1+(a11x*+k1)e-iωkτka12x*e-iωkτk-a21y*e-iωkτkiωk-k2+(a22y*+k2)e-iωkτk]q(0)=[00].Thus, we can easily obtain α=-a21y*e-iωkτk/(iωk-k2+(a22y*+k2)e-iωkτk), q(0)=(1,α)T.Similarly, letq*(s)=D(1,α*)eiωkτks be the eigenvector of A* corresponding to -iωkτk. By the definition of A*, we can compute α*=a12x*eiωkτk/(iωk+k2-(a22y*+k2)eiωkτk).In order to assure〈q*(s),q(θ)〉=1, we need to determine the value of D. From (36), we have
(38)〈q*(s),q(θ)〉=D¯(1,α¯*)(1,α)T-∫-10‍∫ξ=0θ‍D¯(1,α¯*)e-iωkτk(ξ-θ)dη(θ)(1,α)Teiωkτkξdξ=D¯{αα¯*)e-iωkτk1+αα¯*+τk(-a11x*-k1+a21y*α¯*-a12x*αvvvvvvvvvvvvvvvvvvvv-(a22y*+k2)αα¯*)e-iωkτk}.
Thus, we can choose
(39)D¯={-(a22y*+k2)αα¯*)e-iωkτk1+αα¯*+τk(-a11x*-k1+a21y*α¯*-a12x*αvvvvvvvvvvvvvv-(a22y*+k2)αα¯*)e-iωkτk}-1,
such that 〈q*(s),q(θ)〉=1, 〈q*(s),q¯(θ)〉=0.In the following, we first compute the coordinates to describe the center manifoldC0 at μ=0. Define
(40)z(t)=〈q*,xt〉,W(t,θ)=xt(θ)-2Re{z(t)q(θ)}.
On the center manifold C0, we have
(41)W(t,θ)=W(z(t),z¯(t),θ)=W20(θ)z22+W11(θ)zz¯+W02(θ)z¯22+W30(θ)z36+⋯,
where z and z¯ are local coordinates for center manifold C0 in the direction of q and q¯. Note that W is real if xt is real. We consider only real solutions. For the solution xt∈C0, since μ=0, we have
(42)z˙=iωkτkz+〈q*(θ),f({z(t)q(θ)}0,W(z(t),z¯(t),θ)vvvvvvvvvvvvvvvvvvvvv+2Re{z(t)q(θ)})〉=iωkτkz+q¯*(0)f({z(t)q(0)}0,W(z(t),z¯(t),0)vvvvvvvvvvvvvvvvvvv+2Re{z(t)q(0)})≜iωkτkz+q¯*(0)f0(z,z¯)=iωkτkz+g(z,z¯),
where
(43)g(z,z¯)=q¯*(0)f0(z,z¯)=g20(θ)z22+g11(θ)zz¯+g02(θ)z¯22+⋯.
From (40) and (41), we have
(44)xt(θ)=(x1t(θ),x2t(θ))T=W(t,θ)+zq(θ)+z¯q¯(θ).
In addition, q(θ)=(1,α)Teiωkτkθ; then
(45)x1t(0)=z+z¯+W20(1)(0)z22+W11(1)(0)zz¯+W02(1)(0)z¯22+O(|(z,z¯)|3),x2t(0)=αz+α¯z¯+W20(2)(0)z22+W11(2)(0)zz¯+W02(2)(0)z¯22+O(|(z,z¯)|3),x1t(-1)=ze-iωkτk+z¯eiωkτk+W20(1)(-1)z22+W11(1)(-1)zz¯+W02(1)(-1)z¯22+O(|(z,z¯)|3),x2t(-1)=αze-iωkτk+α¯z¯eiωkτk+W20(2)(-1)z22+W11(2)(-1)zz¯+W02(2)(-1)z¯22+O(|(z,z¯)|3).
By the definition of f(μ,xt), we have
(46)g(z,z¯)=D¯τk(1,α¯*)[-a11x1t(0)x1t(-1)-a12x1t(0)x2t(-1)a21x2t(0)x1t(-1)-a22x2t(0)x2t(-1)].
Substituting x1t(0), x2t(0), x1t(-1), and x2t(-1) into the above equation and comparing the coefficients with (43), we get
(47)g20=2D¯τk(-a11-a12α+a21αα¯*-a22α2α¯*)e-iωkτk,g11=2D¯τk[(-a11-a12α¯+a21αα¯*-a22αα¯α¯*)eiωkτkvvvvvvvv-(a11+a12α-a21α¯α¯*+a22αα¯α¯*)e-iωkτk],g02=2D¯τk(-a11-a12α+a21α¯α¯*-a22α¯2α¯*)eiωkτk,g21=D¯τk[-a11(2W11(1)(-1)+W20(1)(-1)vvvvvvvvvvvvvvvv+W20(1)(0)eiωkτk+2W11(1)(0)e-iωkτk)vvvvvvv-a12(2W11(2)(-1)+W20(2)(-1)vvvvvvvvvvvvvv+α¯W20(1)(0)eiωkτk+2αW11(1)(0)e-iωkτk)vvvvvvv+a21α¯*(2αW11(1)(-1)+α¯W20(1)(-1)vvvvvvvvvvvvvvv+W20(2)(0)eiωkτk+2W11(2)(0)e-iωkτk)vvvvvv-a22α¯*(2αW11(2)(-1)+α¯W20(2)(-1)vvvvvvvvvvvvvv+α¯W20(2)(0)eiωkτk+2αW11(2)(0)e-iωkτk)].In order to assure the value ofg21, we need to compute W20(θ) and W11(θ). From (34) and (40), we have
(48)W˙=x˙t-z˙q-z¯˙q¯={AW-2Re{q¯*(0)f0q(θ)},θ∈[0,1)AW-2Re{q¯*(0)f0q(θ)}+f0,θ=0≜AW+H(z,z¯,θ),
where
(49)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯.
Notice that near the origin on the center manifold C0, we have
(50)W˙=Wzz˙+Wz¯z¯˙;
thus, we have
(51)(A-2iωkτkI)W20(θ)=-H20(θ),AW11(θ)=-H11(θ).
From (48), for θ∈[-1,0), we have
(52)H(z,z¯,θ)=-q¯*(0)f0q(θ)-q*(0)f¯0q¯(θ)=-gq(θ)-g¯q¯(θ).
Comparing the coefficients with (48) gives that
(53)H20(θ)=-g20q(θ)-g¯02q¯(θ),H11(θ)=-g11q(θ)-g¯11q¯(θ).
From (51), (53), and the definition of A, we can get
(54)W˙20(θ)=2iωkτkW20(θ)+g20q(θ)+g¯02q¯(θ).Notice thatq(θ)=q(0)eiωkτkθ; we have
(55)W20(θ)=ig20ωkτkq(0)eiωkτkθ+ig¯023ωkτkq¯(0)e-iωkτkθ+E1e2iωkτkθ,
where E1=(E1(1),E1(2))∈R2 is a constant vector. In the same way, we can also obtain
(56)W11(θ)=-ig11ωkτkq(0)eiωkτkθ+ig¯11ωkτkq¯(0)e-iωkτkθ+E2,
where E2=(E2(1),E2(2))∈R2 is also a constant vector.In what follows, we will computeE1 and E2. From the definition of A and (51), we have
(57)∫-10‍dη(θ)W20(θ)=2iωkτkW20(0)-H20(0),(58)∫-10‍dη(θ)W11(θ)=-H11(0).From (48) and (49), we have
(59)H20(0)=-g20q(0)-g¯02q¯(0)+2τk[-a11e-iωkτk-a12αe-iωkτka21αe-iωkτk-a22α2e-iωkτk],(60)H11(0)=-g11q(0)-g¯11q¯(0)+2τk[-a11Re{e-iωkτk}-a12Re{αe-iωkτk}a21Re{αeiωkτk}-a22Re{αα¯e-iωkτk}].Substituting (55) and (59) into (57) and noticing that
(61)[iωkτkI-∫-10‍eiωkτkθdη(θ)]q(0)=0,[-iωkτkI-∫-10‍e-iωkτkθdη(θ)]q¯(0)=0,
we obtain
(62)[2iωkτkI-∫-10‍e2iωkτkθdη(θ)]E1=2τk[-a11e-iωkτk-a12αe-iωkτka21αe-iωkτk-a22α2e-iωkτk],
which leads to
(63)BE1=2[-a11e-iωkτk-a12αe-iωkτka21αe-iωkτk-a22α2e-iωkτk],
where(64)B=[2iωk-k1+(a11+k1)x*e-2iωkτka12x*e-2iωkτk-a21y*e-2iωkτk2iωk-k2+(a22y*+k2)e-2iωkτk].It follows that
(65)E1=2B-1[-a11e-iωkτk-a12αe-iωkτka21αe-iωkτk-a22α2e-iωkτk].Similarly, substituting (56) and (60) into (58), we can get the formula of E2, where(66)E2=2[a11x*a12x*-a21y*a22y*]-1×[-a11Re{e-iωkτk}-a12Re{αe-iωkτk}a21Re{αeiωkτk}-a22Re{αα¯e-iωkτk}].Thus, we can determineW20(θ) and W11(θ). Furthermore, we can determine each gij. Therefore, each gij is determined by the parameters and delay in (5). Thus, we can compute the following values:
(67)c1(0)=i2ωkτk(g20g11-2|g11|2-13|g02|2)+12g21,μ2=-Re{c1(0)}Re{λ′(0)},T2=-Im⁡{c1(0)}+μ2Im⁡{λ′(0)}ωkτk,β2=2Re{c1(0)},
which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value τk; that is, μ2 determines the directions of the Hopf bifurcation: if μ2>0(<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcation exists for τ>τ0(<τ0); β2 determines the stability of the bifurcation periodic solutions: the bifurcating periodic solutions are stable (unstable) if β2<0(>0); and T2 determines the period of the bifurcating periodic solutions: the period increases (decrease) if T2>0(<0).


## 4. Numerical Simulations

In this section, we present numerical results to verify the analytical predictions obtained in the previous sections and use the delayed-feedback controller to control the Hopf bifurcation and chaos of system (3).For the convenience of the calculation, we take the parameters of system (3) as r1=r2=1, a11=1, a12=1, a21=2, a22=1. Then system (3) has a positive equilibrium E*=(2/3,1/3).From [9], when k1=k2=0 or τ=0, system (3) becomes system (2), which has a bifurcation point τ0=0.8072, and when τ∈[0,0.8072), E* is asymptotically stable and is unstable when τ>0.8072, and with the increase of the delay τ, chaos occurs via t periodic-doubling bifurcation (Figure 2(a)), when τ=1.4, k1=k2=0, system (3) becomes a chaotic one like illustrated by Figure 1.The bifurcation diagram ofy: (a) k1=k2=0; (b) k1=-0.2, k2=0.5.
(a)
(b)Under the delayed-feedback control, if we choosek1=-0.2, k2=-0.5, from the algorithm of Section 2, we can compute that τ0≈1.4294 and Re[dλ/dτ]τ=τ0=0.1017>0; hence from Theorem 2, when τ∈[0,1.4294), E* is asymptotically stable and is unstable when τ>1.4294, and when τ passes through the critical point τ0, a family of periodic orbits will bifurcate from E*. Furthermore, from the algorithm of Section 3, we can compute that g20=1.4704+2.4314i, g11=6.1433+8.6345i, g02=9.4686+15.7816i, g21=-24.6452+243.1853i; hence, based on the formula of the end of Section 3, we have c1(0)=-27.0460-64.6042i. It is easy to see that μ2>0, β2<0. This means that the Hopf bifurcation is supercritical and the bifurcating periodic solutions are stable (Figure 5(b)).By comparing the two bifurcation plots of Figure2(a) and Figure 2(b), we observe that the bifurcation occurrence is postponed as the system is being controlled. Figure 3 shows that when τ=0.9, if k1=k2=0, E* is unstable, and there are periodic orbits bifurcating from E* (Figure 3(a)); if we take k1=-0.2, k2=-0.5, E* becomes stable (Figure 3(b)). Figure 4 shows that when τ=1.4, if k1=k2=0, E* is unstable, and system (3) becomes chaotic (Figure 4(a)); if we take k1=-0.2, k2=-0.5, then E* becomes stable and chaos vanishes (Figure 4(b)). Figure 5 shows that when τ=1.5, if k1=k2=0, E* is unstable, and system (3) still exhibits chaotic behavior (Figure 5(a)); if we take k1=-0.2, k2=-0.5, then E* becomes unstable and stable periodic orbits bifurcate from E* (Figure 5(b)).Waveform plot and phase plot of system (3) with τ=0.9: (a) k1=k2=0; (b) k1=-0.2, k2=-0.5.
(a)
(b)Waveform plot and phase plot of system (3) with τ=1.4: (a) k1=k2=0; (b) k1=-0.2, k2=-0.5.
(a)
(b)Waveform plot and phase plot of system (3) with τ=1.5: (a) k1=k2=0; (b) k1=-0.2, k2=-0.5.
(a)
(b)


## 5. Conclusion

In this paper, we have studied a delayed Lotka-Volterra predator-prey system with time delayed feedback by using the theory of functional differential equation and Hassard’s method. By analyzing the corresponding characteristic equations, the local stability of the positive equilibrium of system (3) was discussed.We have obtained the estimated length of gestation delay which would not affect the stable coexistence of both prey and predator species at their equilibrium values. The existence of Hopf bifurcation for system (3) at the positive equilibrium was also established. From theoretical analysis it was shown that the larger values of gestation time delay cause fluctuation in individual population density and hence the system becomes unstable.As the estimated length of delay to preserve stability and the critical length of time delay for Hopf-bifurcation are dependent upon the system parameters, it is possible to impose some control, which will prevent the possible abnormal oscillation in population density. Our results show that if we choose some appropriate parameters, the oscillation can be controlled to a stable equilibrium or a stable periodic orbit; that is to say, we can achieve the ecological equilibrium by adjusting the capture (or release) level.



---

*Source: 104156-2014-06-04.xml*