# Dynamic Properties of a Differential-Algebraic Biological Economic System

Authors: Hongyang Zhang; Chunrui Zhang
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We analyze a differential-algebraic biological economic system with time delay. The model has two different Holling functional responses. By considering time delay as bifurcation parameter, we find that there exists stability switches when delay varies, and the Hopf bifurcation occurs when delay passes through a sequence of critical values. Furthermore, we also consider the stability and direction of the Hopf bifurcation by applying the normal form theory and the center manifold theorem. Finally, using Matlab software, we do some numerical simulations to illustrate the effectiveness of our results.


---

## Body

## 1. Introduction

Recently, dynamics of differential-algebraic equations have become a subject of intense research activities because differential-algebraic equation models can describe many kinds of practical systems such as power system, biological economic system, and singular prey-predator economic model [1, 2]. The dynamic relationship between predators and their preys with harvesting has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance. Harvesting has a strong impact on the dynamic evolution of a population. Depending on the nature of applied harvesting strategy, the long-run stationary density of the population may be significantly smaller than the long-run stationary density of a population in the absence of harvesting. Time delays of one type or another have been incorporated into mathematical models of population dynamics due to maturation time, capturing time, or other reasons. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate. Recently, a number of researches have been devoted to the dynamic properties for delayed predator-prey systems with harvesting. For further details, refer to [3–7].Kar and Pahari [8] have investigated the dynamical behavior of an exploited system consisting of a prey and a predator that are harvested:(1.1)ẋ=x(t)g(x(t-τ))-yp(x(t))-q1E1x(t),ẏ=y(t)(-d+αp(x(t)))-q2E2y(t),
where x(t) and y(t) respect prey and predator population densities at time t and g(x) is the growth rate of prey when y(t) is zero. p(x) is the response function, d is the death rate of y(t) and α is the conversion factor. q1E1x and q2E2y represent harvesting of x(t) and y(t).In system (1.1) x(t) and y(t) have the same Holling function. But in real life they may have different adaptation for the environment. So we formulate the aforementioned problem with different Holling functions as follows:(1.2)ẋ=α(1-x(t-τ)k)x(t)-βx(t)y(t)1+bx(t)+cy(t)-E(t)x(t),ẏ=-dy(t)+aβx(t)y(t)1+bx(t).In 1954, Gordon analyzed the effect of harvest effort on ecosystem from an economic perspective and proposed the following economic theory [9]:Net Economic Revenue (NER) = Total Revenue (TR) − Total Cost (TC).Based on the previous economic theory, we letE(t) represent the harvest effort for x(t) at time t. TR=pE(t)x(t) and TC=c0E(t). p represents harvesting reward per unit harvesting effort for unit weight of prey. c0 represents the cost per unit harvest effort for prey. Then a differential-algebraic model which consists of two differential equations and an algebraic equation can be established as follows:(1.3)ẋ=α(1-x(t-τ)k)x(t)-βx(t)y(t)1+bx(t)+cy(t)-E(t)x(t),ẏ=-dy(t)+aβx(t)y(t)1+bx(t),E(t)(px(t)-c0)-m=0.
Define(1.4)g(x,E)=E(t)(px(t)-c0)-m.
From the standpoint of biology, we are only interested in the dynamics of model (1.3) in the region(1.5)R+3={(x,y,E)∣x≥0,y≥0,E≥0}.


## 2. Stability Analysis

Firstly, we consider the location and number of the equilibria of model (1.3) in R+3. The interior equilibrium point is a root of the following equation:(2.1)α(1-x0k)x0-βx0y01+bx0+cy0-E0x0=0,-dy0+aβx0y01+bx0=0,E0(px0-c0)-m=0.
It is obvious that (2.1) has an only real solution Y0=(x0,y0,E0) where(2.2)x0=daβ-bd,y0=(α(1-x0/k)-m/(px0-c0))(1+bx0)β-c(α(1-x0/k)-m/(px0-c0)),E0=mpx0-c0,
and αβ>bd,  0<m<α(px0-c0)(1-x0/k). In order to analyze the local stability of the positive equilibrium point for the system (1.3) we use the linear transformation YT=QNT(2.3)N=(u(t),v(t),E¯(t)),Q=(100010-pE0px0-c001).
Then we have(2.4)u(t)=x(t),v(t)=y(t),E¯(t)=pE0px0-c0x(t)+E(t),Dyg(x)Q=(0,0,px0-c0),
for which system (1.3) yields(2.5)u̇=α(1-u(t-τ)k)u(t)-βu(t)v(t)1+bu(t)+cv(t)-E¯(t)u(t)+pE0u2(t)px0-c0,v̇=-dv(t)+aβu(t)v(t)1+bu(t),(E¯(t)-pE0u(t)px0-c0)(pu(t)-c0)-m=0.Now we derive the formula for determining the properties of the positive equilibrium point of the system (2.5). Firstly, we consider the local parametric [10] ψ, which is defined as follows:(2.6)ψ(z(t))=N0T+U0z(t)+V0h(z(t)),
where U0=(100100), V0=(001), z(t)=(y1(t)y2(t)), and h=(00h3),  h3=m/(p(u0+y1(t))-c0)+pE0(u0+y1(t))/(px0-c0)-E0. Introducing perturbations u(t)=u0+y1andv(t)=v0+y2 and defining  ψ(z(t))=(u(t),v(t),E¯(t))T,  E¯(t)=E0+h3(y(1),y(2)), and  g(ψ(z(t)))=0, then we can obtain the parametric system of system (2.5) as follows:(2.7)ẏ1=α(1-y1(t-τ)+u0k)(y1(t)+u0)-β(y1(t)+u0)(y2(t)+v0)1+b(y1(t)+u0)+c(y2(t)+v0)-(E¯0+h3(y1(t),y2(t)))(y1(t)+u0)+pE0(y1(t)+u0)px0-c0,ẏ2=-d(y2+v0)+aβ(y1(t)+u0)(y2(t)+v0)1+b(y1(t)+u0).
Neglecting 2nd- and higher-order products of y1 and y2, we get(2.8)ẏ1=(α(1-u0k)-βv0(1+cv0)(1+bv0+cv0)2-E0+pE0u0px0-c0)y1(t)-αu0ky1(t-τ)-βu0(1+bu0)(1+bu0+cv0)2y2(t),ẏ2=aβv0(1+bu0)2y2(t).
The characteristic equation of system (2.8) at (0,0) takes the form(2.9)λ2+(A-αu0k+αu0ke(-λτ))λ+B=0.First of all, we know that the roots of (2.9) with τ=0 have always negative real parts when A<0,B>0, where(2.10)A=-α(1-2u0k)-βv0(1+cv0)(1+bu0+cv0)2-E0+pE0u0pu0-c,B=αβv0(1+bu0+cv0)2(1+bu0)2.Next, we will investigate the distribution of the roots of (2.9). If iω(ω>0) is a root of (2.9), then we have(2.11)-ω2+(A-αu0k+αu0keiωτ)iω+B=0.
Separating the real and imaginary parts, we have(2.12)αωu0ksin⁡(ωτ)=ω2-B,αωu0kcos⁡(ωτ)=ω(αu0k-A),
which lead to(2.13)ω4+[(αu0k-A)2-2B-(αu0k)2]ω2+B2=0.
When (αu0/k-A)2<2B+(αu0/k)2 and [(αu0/k-A)2-2B-(αu0/k)2]2>4B2 hold, we can get that (2.13) has two positive roots ω+ and ω-.According to the analysis previous, we have the following lemma 2.1:Lemma 2.1.
(1) IfA<0, B>0 and (αu0/k-A)2>2B+(αu0/k)2,  then all the roots of (2.9) have negative real parts for all τ≥0.
(2) If(αu0/k-A)2<2B+(αu0/k)2 and [(αu0/k-A)2-2B-(αu0/k)2]2>4B2, then (2.13) has two positive roots ω+, ω-. Furthermore, one has
(2.14)τ±=1ω±arccos(kαu0(αu0k-A))+2nπω±,n=0,1,2,….Differentiating the two sides of (2.9) and noticing that λ is a function of τ, we have(2.15)(dλdτ)-1=2λ-(αu0/k)λτe-λτ+(A-αu0/k+(αu0/k)e-λτ)(αu0/k)λ2e-λτ,and then we can obtain(2.16)sign⁡{Re(dλdτ)}τ=τ+,ω=ω+>0,sign⁡{Re(dλdτ)}τ=τ-,ω=ω-<0.Therefore, we can obtain the following results about the stability of the positive equilibrium and the Hopf bifurcation of (2.7).Theorem 2.2.
Assume thatA<0,B>0. (1)
If(αu0/k-A)2<2B+(αu0/k)2 and [(αu0/k-A)2-2B-(αu0/k)2]2>4B2 are satisfied, then there exists an integer N such that the zero solution of system (2.7) is asymptotically stable when τ∈[0,τ0+)⋃(τ0-,τ1+)⋃(τ1-,τ2+)⋃⋯⋃(τN-1-,τN+) and unstable when τ∈(τ0+,τ0-)⋃(τ1+,τ1-)⋃⋯⋃(τN+,+∞).(2)
Whenτ=τn±,n=0,1,2…, there are Hopf bifurcations at the equilibrium point of system (2.7).


## 3. Direction and the Stability of Hopf Bifurcation

In the previous section, we have already obtained the conditions which ensure that model (2.7) undergoes the Hopf bifurcation at the critical values τn(n=0,1,…). In this section, we will study the direction and stability of the Hopf bifurcation based on the normal form approach theory and center manifold theory introduced by Hassard et al. [11].First of all, by the transformationy1=u-u0,  y2=v-v0, t=t/τ, τ=τn+μ, the parametric system (2.7) is equivalent to the following Functional Differential Equation (FDE) system in C=C([-1,0],R2):(3.1)ż(t)=Lμ(zt)+f(μ,zt),
where z(t)=(y1(t),y2(t))T,  zt(θ)=z(t+θ) for θ∈[-1,0], and Lμ:C→R:f:R×C→R.(3.2)Lμ(ϕ)=(τn+μ)(a11a12a210)ϕT(0)+(τn+μ)(b11000)ϕT(-1),a11=(α-αu0k)-βv0(1+cv0)(1+bu0+cv0)2-E0+pE0u0pu0-c0,a12=-βu0(1+bu0)(1+bu0+cv0)2,a21=aβv0(1+bv0)2,b11=-αu0k.
Let ϕ=(ϕ1,ϕ2)T. Then(3.3)f(μ,ϕ)=(τn+μ)(f11,f22)T,
with(3.4)f11=(pE0px0-c0-p2E0u0(px0-c0)2+βbv0(1+bu0+cv0)3)ϕ1(0)2-αkϕ1(0)ϕ1(-1)-β(1+bu0+cv0)2ϕ1(0)ϕ2(0),f22=-βabv0(1+bu0+cv0)3ϕ1(0)2+aβ(1+bu0+cv0)2ϕ1(0)ϕ2(0).
By the Riesz representation theorem, there exists a matrix whose components are bounded variation functions η(θ,μ) such that(3.5)Lμ=∫-10dη(θ,μ)ϕ(θ),ϕ∈C,θ∈[-1,0].
In fact, we can choose(3.6)η(θ,μ)=(τn+μ)(a11a12a210)δ(θ)+(τn+μ)(b11000)δ(θ+1),
where(3.7)δ(θ)={0,θ≠0,1,θ=0,
for ϕ∈C1([-1,0],R2). Define(3.8)A(μ)ϕ(θ)={dϕ(θ)dθ,-1≤θ<0,∫-10dη(θ,μ)ϕ(θ),θ=0,R(μ)ϕ={0,θ∈[-1,0),f(μ,ϕ),θ=0.
Then the system (3.1) is equivalent to(3.9)żt=A(μ)zt+R(μ)zt,
for ψ∈C([0,1],(R2)*), the adjoint operator A* of A is(3.10)A*ψ(s)={-dψ(s)ds,0<s≤1,∫-10dηT(s,0)ψ(-s),s=0,
and a bilinear inner product is given by(3.11)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫θ=10∫ξ=0θψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ.By the discussion in Section2, we know that ±iτnω± are eigenvalues of A(0). Thus, they are also eigenvalues of A*. It is easy to calculate that(3.12)q(θ)=(1,r)Teiω±τnθ,q*(s)=G(r*,1)eiω±τns,
where(3.13)r=iω-a11-b11e-iωτna12,G¯=1r+r*¯-τnr*¯b11e-iωτn,〈q*(s),q(θ)〉=1,〈q*(s),q¯(θ)〉=0.
Let zt be the solution of (3.9) when μ=0. Define(3.14)ż(t)=〈q*,zt〉,w(t,θ)=zt-2Re{ztq(θ)}.
On the center manifold C0, we have(3.15)w(t,θ)=w(z(t),z¯(t),θ)=w20(θ)z22+w11(θ)zz¯+w02(θ)z¯22+⋯.
Note that w is real. When μ=0 we have(3.16)ż=iωτnz+〈q*(θ),f(0,ω+2Rez(t)q(θ))〉.
Then(3.17)ż=iωτnz+g(z,z¯),
where(3.18)g(z,z¯)=g20(θ)z22+g11(θ)zz¯+g02(θ)z¯22+⋯.
From (3.9) and (3.16), we have(3.19)Ẇ=żt-żq-z¯̇q¯={Aw-2Req¯*(0)f(z,z¯)q(θ),-1≤θ≤0,Aw-2Req¯*(0)f(z,z¯)q(θ)+f,θ=0.
Rewrite(3.20)Ẇ=AW+H(z,z¯,θ),
where(3.21)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯.
We can obtain(3.22)(A-2iωτn)w20=-H20(θ),Aw11(θ)=-H11(θ).
Now we have g20, g11, g02, and  g21:(3.23)g20=2G¯τn[Fr*¯-βrr*¯-arβ(1+bu0+cv0)2-αr*¯ke-iωτnθ-βabv0(1+bu0+cv0)3],g11=G¯τn[2Fr*¯+(2aB-2Br*¯)Re(r)-2βabv0(1+bu0+cv0)3-2αr*¯kRe(eiωτnθ)],g02=2G¯τn[Fr*¯-βr¯r*¯-aβr¯(1+bu0+cv0)2-αr*¯keiωτnθ-βabv0(1+bu0+cv0)3],g21=2G¯τn[(2Fr*¯-αr*¯ke-iωτnθ-2βabv0(1+bu0+cv0)3+aβr-βrr*¯(1+bu0+cv0)2)W20(1)(0)+(aβ-βr*¯)(1+bu0+cv0)2W11(2)(0)+(Fr*¯-βabv0(1+bu0+cv0)3-αr*¯2keiωτnθ+aβr¯-βr¯r*¯2(1+bu0+cv0)2)W20(1)(0)-αr*¯kW11(1)(-1)+aβ-βr*¯2(1+bu0+cv0)2W20(2)(0)-αr*¯2kW20(1)(-1)],
where(3.24)F=pE0px0-c0-p2E0u0(px0-c0)2+βbv0(1+bu0+cv0)3,W20(θ)=ig20τnωq(0)eiωτnθ+ig¯023ωτnq¯(0)e-iωτnθ+M1e2iωτnθ,W11(θ)=-ig11τnωq(0)eiωτnθ+ig¯11ωτnq¯(0)e-iωτn+M2,M1(1)=-4A1(1)iω+2a12A1(2)a21a12+2b11iωe-2iωτn+4ω2+2a11iω,M1(2)=-2A1(1)a21+(4iω-2a21-2b11e-2iωτn)A1(2)a21a12+2b11iωe-2iωτn+4ω2+2a11iω,M2(1)=-2A2(2)a21,M2(2)=-2a21A2(1)+2(a11+b11)A2(2)a12a21.A1(1)=F-βr(1+bu0+cv0)2-αke-iωτnθ,A1(2)=aβr(1+bu0+cv0)2-βabv0(1+bu0+cv0)3,A2(1)=F-βRe(r)(1+bu0+cv0)2-αkRe(eiωτnθ),A2(2)=aβRe(γ)(1+bu0+cv0)2-βabv0(1+bu0+cv0)3.At last, we can calculate(3.25)C1(0)=i2ωτn(g11g20-2|g112|-|g02|23)+g212,μ2=-Re{c1(0)}Re{λ′(τn)},β2=2Re{c1(0)},T2=-Im⁡{c1(0)}+μ2Im⁡{λ′(τn)}ωτn.Now we have a conclusion.Conclusion.
(1) The sign ofμ2 determines the direction of Hopf bifurcation. When μ2>0, the Hopf bifurcation is supercritical; when μ2<0, the Hopf bifurcation is subcritical.
(2)β2 determines the stability of bifurcated periodic solutions. When β2<0, the periodic solutions are stable; when β2>0, the periodic solutions is unstable.
(3)T2 determines the period of bifurcated periodic solutions. When T2>0, the period increases; when T2<0, the period decreases.


## 4. Numerical Simulations

In this section, we use some numerical simulations to illustrate the analytical results we obtained in previous sections.Letα=β=2, k=p=2,  a=b=c0=d=1,  c=0, and  m=1/4. The system (1.3) is(4.1)ẋ=2(1-x(t-τ)2)x(t)-2x(t)y(t)1+x(t)-E(t)x(t),ẏ=-y(t)+2x(t)y(t)1+x(t),0=E(t)(2x(t)-1)-14.The only positive equilibrium point of (4.1) is Y0=(1,3/4,1/4). By simple computation we have ω-=0.4164 and  ω+=0.9005. Substituting these parameters into (2.14) gives τ0-=1.2136, τ0+=0.5612,  τ1+=7.5386, and  τ1-=16.2129. Applying Theorem (2.5), we have that the positive equilibrium point of (4.1) is asymptotically stable when τ∈[0,0.5612)∪(1.2136,7.5386). The results can be illustrated in Figures 1 and 3.Whenτ=0.1, the positive equilibrium point Y0 is asymptotically stable.
(a)(b)(c)(d)Whenτ>τ0+=0.5612 and is sufficiently near τ0+, the bifurcating periodic solution from positive equilibrium point of (4.1) occurs. Furthermore, we consider the properties of the Hopf bifurcation at τ=τ0+=0.5612. Based on Theorem (2.5) and conclusion 1, we get c1(0)=0.1673-0.6835i, and  λ′(τ0)=0.9563-0.1856i. This implies that(4.2)μ2=-0.1749,β2=0.3346,T2=1.2883.We conclude that the bifurcating periodic solution is subcritical and unstable. The results are illustrated in Figure2.Whenτ=0.5618, periodic solutions occur from Y0.
(a)(b)(c)(d)Whenτ=1.3, the positive equilibrium point Y0 is stable.
(a)(b)(c)(d)



---

*Source: 205346-2012-01-19.xml*