# The Modeling and Control of a Singular Biological Economic System with Time Delay in a Polluted Environment

Authors: Yi Zhang; Huiying Zhao; Qingling Zhang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

This paper brings up the idea of a biological economic system with time delay in a polluted environment. Firstly, by proper linear transformation and parametric method, the singular time-delay systems are transformed to differential time-delay systems. Then, using center manifold theory and Poincare normal form method, the direction of Hopf bifurcation and the stability and period of its periodic orbits are analysed. At last, we have performed numerical simulation to support the analytical results.


---

## Body

## 1. Introduction

Environmental pollution has been increasingly influencing the biological systems. In order to investigate the development and dynamics of population of the biological systems, it is necessary to consider the factor of pollution when establishing a mathematical model. In addition, delay is also a kind of common phenomenon in reality and it has great influence on the dynamic behavior of system. Therefore, the delay differential equations are needed to describe the system when the influence of time delay is considered. Time delay can lead to the imbalance of the system and the emergence of a variety of bifurcations, among which Hopf bifurcation is the most common. The properties of Hopf bifurcation consist of the stability of the periodic solutions, the direction of bifurcations, the period, and so forth. In recent years, the theory of delay system has gradually been generalised to many important fields by domestic and foreign scholars, including the applications in circuit communiment-system [1], electrodynamics [2], optical [3], ecological-system [4], and economics [5]. Many research findings on biological applications also emerged, such as the analysis of the stability of a class of stochastic system with time delay [6], investigation on nonautonomous competitive Lotka-Volterra systems with infinite delay [7], researching on dynamic behavior of a class of prey-predator model with time delay in a polluted environment [8], the analysis and control of a class of singular prey-predator model with discrete delay [9] which studies the prey-predator system with commercial harvesting and double time delays, and the dynamic behavior analysis and optimal control of a class of economic model with stage structure and pregnancy delay [10].There are many kinds of research methods for delay differential systems. Among them, the most commonly used ones [11] are the center manifold method and the Poincare normal form method. Being always an important mathematical means to investigate the bifurcation problems with parameter and the qualitative theory of differential equations, more attention has been paid to the Poincare normal form method for a long time, home and abroad. In [12], the author lays the foundation of the center manifold standard method by combining the normal form theory and the center manifold theorem, and the method was used on the investigation of Hopf bifurcation. When it comes to related properties of the Hopf bifurcation, the center manifold standard method is usually used to reduce the dimension of high-dimension system, which isolates the asymptotic behaviors of complex systems, so that we can investigate the original system in a center manifold of low dimension, which is much simpler. This paper takes a singular biological economy system with time-delay in a polluted environment and analyses it using the stability theory of singular system, the theory of economic system, the theory of the Hopf bifurcation of delay differential system, and so forth.


## 2. Model Formulation

A single-creature model with stage structure is investigated in [5, 13]:(1)
x
˙
t
=
a
y
t
-
b
x
t
-
r
1
x
t
,
y
˙
t
=
b
x
t
-
r
2
y
t
-
β
y
2
t
,where x
(
t
) and y
(
t
) are the densities of the immature number of creatures and mature number of creatures at time t, a denotes the birth rate of immature creatures, r
1, r
2 are the death rates of immature creatures and mature creatures, b denotes the conversion rate from immature creatures into mature creatures, and β denotes the intraspecific effect coefficient. All coefficients are positive.A single-creature model in the polluted environment is investigated in [5]:(2)
x
˙
t
=
r
x
t
1
-
x
t
K
-
r
1
u
t
x
t
,
u
˙
t
=
θ
-
h
u
t
,where x
(
t
) is the creatures density, u
(
t
) is the concentration of environment pollutants, r denotes the intrinsic growth rate when there is no pollution, K denotes the capacity of the environment, r
1
u
(
t
) can be interpreted as the measuring response function of the reduction of creatures because of the pollution factor, θ denotes the amount of pollutants that are inputted by the outside, and h
u
(
t
) can be interpreted as the reduction of pollutant concentration because of other factors. Assume that endotoxin excretion rate and purification rate are relatively small in an organism body, and thus it can be neglected.Considering the need of a period of time when the immature creatures change into mature creatures, based on system (1) and system (2), the following system is proposed:(3)
x
˙
t
=
a
y
t
-
b
x
t
-
τ
-
r
1
x
t
-
η
1
u
t
x
t
,
y
˙
t
=
b
x
t
-
τ
-
r
2
y
t
-
β
y
2
t
-
E
t
y
t
-
η
2
u
t
y
t
,
u
˙
t
=
θ
-
h
u
t
,
0
=
E
t
p
y
t
-
c
-
m
,where E
(
t
) is the capture capability of mature creatures at time t, p denotes the unit price, c denotes the unit cost, and m denotes the economic profit. p
E
(
t
)
y
(
t
) is the total revenue, and c
E
(
t
) is the total cost. All the parameters are positive [14, 15].


## 3. Stability Analysis

Theorem 1.
The positive equilibrium of system (3) is locally asymptotically stable when τ
∈
0
,
τ
0; the positive equilibrium of system (3) is unstable when τ
∈
(
τ
0
,
+
∞
); system (3) undergoes a Hopf bifurcation at the positive equilibrium when τ
=
τ
0.Proof.
System (3) will be investigated in this chapter [16].
Considering system (3), let(4)
a
y
-
b
x
-
r
1
x
-
η
1
u
x
=
0
,
b
x
-
r
2
y
-
β
y
2
-
E
y
-
η
2
u
y
=
0
,
θ
-
h
u
=
0
,
E
p
y
-
c
-
m
=
0
.From (4), the following root can be obtained: (5)
u
∗
=
θ
h
,
E
∗
=
m
p
y
∗
-
c
,
x
∗
=
a
b
+
r
1
+
η
1
θ
/
h
y
∗
,where y
∗ is decided by the following equation: (6)
β
p
y
∗
2
-
p
a
b
b
+
r
1
+
θ
/
h
η
1
-
p
r
2
-
p
θ
h
η
2
+
c
β
y
∗
+
a
b
c
b
+
r
1
+
θ
/
h
η
1
-
c
r
2
-
c
θ
h
η
2
+
m
=
0
.The discriminant of the equation is(7)
Δ
=
p
a
b
b
+
r
1
+
θ
/
h
η
1
-
p
r
2
-
p
θ
h
η
2
+
c
β
2
-
4
β
p
a
b
c
b
+
r
1
+
θ
/
h
η
1
-
c
r
2
-
c
θ
h
η
2
+
m
.
WhenΔ
=
0, the parameters of system (3) meet conditions P
a
b
-
(
p
r
2
+
p
θ
/
h
η
1
-
c
β
)
(
b
+
r
1
+
θ
/
h
η
1
)
>
0, and then P
(
x
∗
,
y
∗
,
u
∗
,
E
∗
) is the unique positive equilibrium point of system (1), where(8)
u
∗
=
θ
h
,
E
∗
=
m
p
y
∗
-
c
,
x
∗
=
a
b
+
r
1
+
θ
/
h
η
1
y
∗
,
y
∗
=
a
b
2
β
b
+
r
1
+
θ
/
h
η
1
-
r
2
2
β
-
θ
/
h
η
2
2
β
+
c
2
p
.
WhenΔ
>
0, system (3) has two equilibrium points P
x
i
∗
,
y
i
∗
,
u
i
∗
,
E
i
∗, i
=
1,2, where(9)
u
i
∗
=
θ
h
,
E
i
∗
=
m
p
y
i
∗
-
c
,
x
i
∗
=
a
b
+
r
1
+
θ
/
h
η
1
y
i
∗
,
y
1
∗
=
a
b
2
β
b
+
r
1
+
θ
/
h
η
1
-
r
2
2
β
-
θ
/
h
η
2
2
β
+
c
2
p
+
Δ
2
β
p
>
y
∗
,
y
2
∗
=
a
b
2
β
b
+
r
1
+
θ
/
h
η
1
-
r
2
2
β
-
θ
/
h
η
2
2
β
+
c
2
p
-
Δ
2
β
p
<
y
∗
.
The stability of the unique positive equilibrium point of systemP
(
x
∗
,
y
∗
,
u
∗
,
E
∗
) is investigated as an example; the stability of the other equilibrium points can be investigated in the same way.
In order to make the research more convenient, system (3) can be written as [17](10)
f
X
=
f
1
x
,
y
,
u
,
E
f
2
x
,
y
,
u
,
E
f
3
x
,
y
,
u
,
E
=
a
y
-
b
x
t
-
τ
-
r
1
x
-
η
1
u
x
b
x
t
-
τ
-
r
2
y
-
β
y
2
-
E
y
-
η
2
u
y
θ
-
h
u
,
g
X
=
E
p
y
-
c
-
m
,where X
=
(
x
,
y
,
u
,
E
)
T.
In order to investigate the local stability of the positive equilibrium point, make the following transformation on system (10):(11)
N
t
=
x
1
t
y
1
t
u
1
t
E
1
t
,
Q
=
1
0
0
0
0
1
0
0
0
0
1
0
0
-
p
E
∗
p
y
∗
-
c
0
1
,
X
t
=
Q
N
t
.Then,(12)
x
1
t
=
x
t
,
y
1
t
=
y
t
,
E
1
t
=
E
t
+
p
E
∗
p
y
∗
-
c
y
t
.By generating system (3), the following system can be obtained:(13)
x
˙
1
t
=
a
y
1
t
-
b
x
1
t
-
τ
-
r
1
x
1
t
-
η
1
u
1
t
x
1
t
,
y
˙
1
t
=
b
x
1
t
-
τ
-
r
2
y
1
t
-
β
y
1
2
t
-
E
1
t
-
p
E
∗
p
y
∗
-
c
y
1
t
y
1
t
-
η
2
u
1
t
y
1
t
,
u
˙
1
t
=
θ
-
h
u
1
t
0
=
E
1
t
-
p
E
∗
p
y
∗
-
c
y
1
t
p
y
1
t
-
c
-
m
.
In order to derive the formula determining the properties of the positive equilibrium of system (13), we consider local parametric Ψ of the fourth equation of system (13) as literatures [18], which is given as follows:(14)
N
t
=
Ψ
Z
t
=
N
0
+
u
0
Z
t
+
h
Z
t
,
g
ψ
Z
t
=
0
,where(15)
u
0
=
1
0
0
0
1
0
0
0
1
0
0
0
,
Z
t
=
x
2
t
y
2
t
u
2
t
,
N
0
=
x
∗
y
∗
u
∗
E
∗
,
h
(
Z
(
t
)
)
=
0,0
,
0
,
h
4
x
1
t
,
y
1
t
,
u
1
t
T, and R
3
→
R
Z is a continuous map.(16)
x
1
t
=
x
∗
+
x
2
t
,
y
1
t
=
y
∗
+
y
2
t
,
u
1
t
=
u
∗
+
u
2
t
,
E
1
t
=
E
∗
+
h
4
x
1
t
,
y
1
t
,
u
1
t
.By generating system (13), the following system can be obtained:(17)
x
˙
2
t
=
a
y
∗
+
y
2
t
-
b
x
∗
+
x
2
t
-
τ
-
r
1
x
∗
+
x
2
t
-
η
1
x
∗
+
x
2
t
u
∗
+
u
2
t
,
y
˙
2
t
=
b
x
∗
+
x
2
t
-
τ
-
r
2
y
∗
+
y
2
t
-
β
y
∗
+
y
2
t
2
-
E
∗
+
h
4
x
1
t
,
y
1
t
,
u
1
t
-
p
E
∗
p
y
∗
-
c
y
∗
+
y
2
t
y
∗
+
y
2
t
-
η
2
y
∗
+
y
2
t
u
∗
+
u
2
t
,
u
˙
2
t
=
θ
-
h
u
∗
+
u
2
t
.The linearized system of parametric system (17) at 0,0
,
0 can be given as follows:(18)
x
˙
2
t
=
a
y
2
t
-
b
x
2
t
-
τ
-
r
1
+
η
1
u
∗
+
η
1
u
2
x
2
t
-
η
1
x
∗
u
2
t
+
a
y
∗
-
b
x
∗
-
r
1
x
∗
-
η
1
u
∗
x
∗
,
y
˙
2
t
=
b
x
2
t
-
τ
-
r
2
-
2
p
E
∗
y
∗
p
y
∗
-
c
+
E
∗
+
2
β
y
∗
+
η
2
u
∗
+
β
y
2
t
-
2
p
E
∗
p
y
∗
-
c
+
η
2
u
2
t
y
2
t
-
η
2
u
∗
+
u
2
t
y
∗
+
b
x
∗
-
r
2
y
∗
-
β
y
∗
2
-
E
∗
y
∗
+
p
E
∗
y
∗
2
p
y
∗
-
c
,
u
˙
2
t
=
θ
-
h
u
∗
-
h
u
2
t
.
The characteristic equation of system (18) is as follows:(19)
λ
2
+
A
11
λ
+
A
12
+
b
e
-
λ
τ
λ
+
B
11
e
-
λ
τ
=
0
,where A
=
r
2
-
2
p
E
∗
y
∗
/
p
y
∗
-
c
+
E
∗
+
2
β
y
∗
+
η
2
u
∗
+
η
2
+
β
y
2
(
t
)
-
2
p
E
∗
/
p
y
∗
-
c, A
11
=
A
+
r
1
+
η
1
u
∗
+
η
1
u
2, A
12
=
A
(
r
1
+
η
1
u
∗
+
η
1
u
2
), B
11
=
A
-
a
b, and let λ
=
i
ω to be a root of (19), ω
>
0, λ in (19); then(20)
-
ω
2
+
A
11
i
ω
+
A
12
+
i
ω
b
cos
⁡
ω
τ
+
b
ω
sin
⁡
ω
τ
+
B
11
cos
⁡
ω
τ
-
i
B
11
sin
⁡
ω
τ
=
0
.
Separate the imaginary part and real part from (20) and the following equations can be obtained:(21)
-
ω
2
+
A
12
+
b
ω
sin
⁡
ω
τ
+
B
11
cos
⁡
ω
τ
=
0
,
A
11
ω
+
ω
b
cos
⁡
ω
τ
-
B
11
sin
⁡
ω
τ
=
0
.Then, the following equation can be obtained by calculating(22)
ω
4
-
2
A
12
-
A
11
2
+
b
2
ω
2
+
A
12
2
-
B
11
2
=
0
.Let v
=
ω
2; then (22) can be changed as the following form:(23)
v
2
-
2
A
12
-
A
11
2
+
b
2
v
+
A
12
2
-
B
11
2
=
0
.When -
(
2
A
12
-
A
11
2
+
b
2
)
>
0,  A
12
2
-
B
11
2
>
0, (22) characteristic has no real root; when -
2
A
12
-
A
11
2
+
b
2
>
0, A
12
2
-
B
11
2
≤
0, (22) characteristic has one real root ω
1
+, in which(24)
ω
1
+
=
2
A
12
-
A
11
2
+
b
2
+
2
A
12
-
A
11
2
+
b
2
2
-
4
A
12
2
-
B
11
2
2
.When -
(
2
A
12
-
A
11
2
+
b
2
)
<
0, A
12
2
-
B
11
2
≤
0, (22) characteristic has one real root ω
+, in which(25)
ω
+
=
2
A
12
-
A
11
2
+
b
2
+
2
A
12
-
A
11
2
+
b
2
2
-
4
A
12
2
-
B
11
2
2
.When -
(
2
A
12
-
A
11
2
+
b
2
)
<
0, A
12
2
-
B
11
2
>
0, (22) characteristic have two real roots ω
2
±, in which(26)
ω
2
±
=
2
A
12
-
A
11
2
+
b
2
+
2
A
12
-
A
11
2
+
b
2
2
-
4
A
12
2
-
B
11
2
2
.
According to (21), if -
(
2
A
12
-
A
11
2
+
b
2
)
<
0, A
12
2
-
B
11
2
>
0, τ can be written as the following form: (27)
τ
j
=
1
ω
2
±
arccos
⁡
B
11
2
ω
2
±
2
-
b
A
11
ω
2
±
2
-
A
12
B
11
B
11
2
+
b
2
ω
2
±
2
+
2
j
π
ω
2
±
.Then, (19) has a pair of pure imaginary roots ±
i
ω
2
± at τ
j; assuming that (19) has a solution λ
(
τ
)
=
α
(
τ
)
+
i
ω
(
τ
), when α
(
τ
j
)
=
0, ω
(
τ
j
)
=
±
i
ω
2
±, λ
(
τ
) in (19), taking the derivative of λ with respect to τ in (19),(28)
2
λ
d
λ
d
τ
+
A
11
+
b
e
-
λ
τ
d
λ
d
τ
-
b
λ
τ
e
-
λ
τ
d
λ
d
τ
-
b
λ
2
e
-
λ
τ
-
B
11
λ
e
-
λ
τ
-
B
11
τ
e
-
λ
τ
d
λ
d
τ
=
0
,when(29)
d
λ
d
τ
=
b
λ
2
e
-
λ
τ
+
B
11
λ
e
-
λ
τ
-
A
11
2
λ
+
b
e
-
λ
τ
-
b
λ
τ
e
-
λ
τ
-
B
11
τ
e
-
λ
τ
.When λ
=
i
ω
2
±, τ
j
=
τ
0, (τ
0 is the minimum of τ
j)(30)
d
λ
d
τ
λ
=
i
ω
2
±
τ
=
τ
0
=
b
ω
2
2
+
B
11
ω
2
±
cos
⁡
ω
2
±
τ
-
A
11
-
i
b
ω
2
2
+
B
11
ω
2
±
sin
⁡
ω
2
±
τ
2
ω
2
±
+
b
τ
sin
⁡
ω
2
±
τ
+
B
11
τ
sin
⁡
ω
2
±
τ
-
b
sin
⁡
ω
2
±
τ
i
-
b
τ
+
B
11
τ
-
b
cos
⁡
ω
2
±
τ
.Denote(31)
C
1
=
b
ω
2
2
+
B
11
ω
2
±
cos
⁡
ω
2
±
τ
-
A
11
,
C
2
=
b
ω
2
2
+
B
11
ω
2
±
sin
⁡
ω
2
±
τ
,
C
3
=
2
ω
2
±
+
b
τ
sin
⁡
ω
2
±
τ
+
B
11
τ
sin
⁡
ω
2
±
τ
-
b
sin
⁡
ω
2
±
τ
,
C
4
=
b
τ
+
B
11
τ
-
b
cos
⁡
ω
2
±
τ
.Then, we can get(32)
R
e
d
λ
d
τ
λ
=
i
ω
2
±
τ
=
τ
0
=
R
e
b
ω
2
2
+
B
11
ω
2
±
cos
⁡
ω
2
±
τ
-
A
11
-
i
b
ω
2
2
+
B
11
ω
2
±
sin
⁡
ω
2
±
τ
2
ω
2
±
+
b
τ
sin
⁡
ω
2
±
τ
+
B
11
τ
sin
⁡
ω
2
±
τ
-
b
sin
⁡
ω
2
±
τ
i
-
b
τ
+
B
11
τ
-
b
cos
⁡
ω
2
±
τ
=
R
e
C
1
-
C
2
i
C
3
i
-
C
4
=
-
C
1
C
4
+
C
2
C
3
C
3
2
+
C
4
2
<
0
.According to the above analysis, we can prove the theorem.


## 4. Direction and Stability of the Hopf Bifurcation

Theorem 2.
The properties of Hopf bifurcation are determined by (70), the detailed contents are as follows:(1)
The direction of Hopf bifurcation: the Hopf bifurcation is supercritical (resp., subcritical) whenμ
2
>
0
(resp., μ
2
<
0
) and the bifurcation periodic solutions exist when τ
>
τ
0  (
τ
<
τ
0
)
.
(2)
The stability of the bifurcating periodic solutions: the bifurcation periodic solutions are stable (resp., unstable) ifβ
2
<
0 (resp., β
2
>
0 ).
(3)
The period of the bifurcating periodic solutions: the period increases (resp., decreases) ifτ
2
>
0 (resp., τ
2
<
0).Proof.
Considering system (19), let z
1
=
x
2
-
x
∗,  z
2
=
y
2
-
y
∗,   z
3
=
u
2
-
u
∗; then system (3) can be written as(33)
z
˙
t
=
L
μ
z
t
+
f
μ
,
z
t
,in which z
(
t
)
=
z
1
t
,
z
2
t
,
z
3
t
T
∈
R
3,  L
μ
:
C
→
R
3,  f : R
×
C
→
R
3.(34)
L
μ
ϕ
θ
=
-
r
1
a
0
0
-
r
2
-
E
∗
0
0
0
-
h
ϕ
1
0
ϕ
2
0
ϕ
3
0
+
-
b
0
0
b
0
0
0
0
0
ϕ
1
-
τ
0
ϕ
2
-
τ
0
ϕ
3
-
τ
0
,
(35)
f
μ
,
z
t
=
-
η
1
ϕ
1
0
ϕ
3
0
-
β
ϕ
2
2
0
-
η
2
ϕ
2
0
ϕ
3
0
0
.
By Riesz representation theorem, there existρ
(
θ
,
μ
) in interval θ
∈
[
-
τ
0
,
0
] and(36)
L
μ
ϕ
=
∫
-
τ
0
0
d
ρ
θ
,
μ
ϕ
θ
ϕ
∈
C
,in which ρ
(
θ
,
μ
) is a bounded variation function. Let(37)
ρ
θ
,
μ
=
-
r
1
a
0
0
-
r
2
-
E
∗
0
0
0
-
h
δ
θ
+
-
b
0
0
b
0
0
0
0
0
δ
θ
+
τ
0for ϕ
∈
C
1
-
τ
0
,
0
,
R
3, defining(38)
B
μ
=
d
ϕ
θ
d
θ
,
θ
∈
-
τ
0
,
0
,
∫
-
τ
0
0
d
ρ
μ
,
t
ϕ
t
,
θ
=
0
,
ψ
s
,
ϕ
θ
=
ψ
¯
T
0
ϕ
0
-
∫
-
τ
0
0
∫
ξ
=
0
θ
ψ
¯
T
ξ
-
θ
d
ρ
θ
ϕ
ξ
d
ξ
.
Then, (33) is equivalent to(39)
z
˙
t
=
B
μ
z
t
+
R
μ
z
t
,in which z
t
(
θ
)
=
z
(
t
+
θ
),  θ
∈
[
-
τ
0
,
0
].
Forψ
∈
C
1
(
[
0
,
τ
0
]
,
R
3
),(40)
B
∗
ψ
s
=
d
ψ
s
d
s
,
s
∈
-
τ
0
,
0
,
∫
-
τ
0
0
d
ρ
T
μ
,
t
ϕ
-
t
,
s
=
0
,
(41)
ψ
s
,
ϕ
θ
=
ψ
¯
T
0
ϕ
0
-
∫
-
τ
0
0
∫
ξ
=
0
θ
ψ
¯
T
ξ
-
θ
d
ρ
θ
ϕ
ξ
d
ξ
.
It is easy to calculate thati
ω
2 and -
i
ω
2 are the eigenvalues of B
(
0
) and B
∗
(
0
); then we can obtain the feature vectors of i
ω
2 and -
i
ω
2 of q
(
θ
) and q
∗
(
s
); let q
(
θ
)
=
q
(
0
)
e
i
ω
2
θ
=
(
1
,
α
1
,
β
1
)
T
e
i
ω
2
θ; according to (34) and (37), we can get(42)
i
ω
2
+
r
1
+
b
e
i
ω
2
θ
-
a
0
-
b
e
i
ω
2
θ
i
ω
2
+
r
2
0
0
0
i
ω
2
+
h
1
α
1
β
1
=
0
0
0
.Then,(43)
q
0
=
1
,
α
1
,
β
1
T
=
1
,
i
ω
2
+
r
1
a
-
i
ω
2
-
r
2
,
0
T
.By the same way, let q
∗
(
s
)
=
q
∗
(
0
)
e
i
ω
2
s
=
(
1
,
α
1
∗
,
β
1
∗
)
T
e
i
ω
2
s; we can get(44)
q
∗
0
=
1
,
α
1
∗
,
β
1
∗
T
=
1
,
a
-
i
ω
2
+
r
2
,
0
T
.By (37), we can get(45)
q
∗
s
q
θ
=
D
¯
1
,
α
¯
1
∗
,
β
¯
1
∗
1
,
α
1
,
β
1
T
-
∫
-
τ
0
0
∫
ξ
=
0
θ
D
¯
1
,
α
¯
1
∗
,
β
¯
1
∗
e
-
i
ω
2
ξ
-
θ
d
ρ
θ
1
,
α
1
,
β
1
T
e
i
ω
2
ξ
d
ξ
=
D
¯
1
+
α
¯
1
∗
α
1
+
β
¯
1
∗
β
1
-
∫
-
τ
0
0
1
,
α
¯
1
∗
,
β
¯
1
∗
θ
e
i
θ
ω
2
d
ρ
θ
1
,
α
1
,
β
1
T
=
D
¯
1
+
α
¯
1
∗
α
1
+
β
¯
1
∗
β
1
-
τ
0
α
¯
1
∗
α
1
θ
e
-
i
ω
2
τ
0due to q
∗
(
s
)
q
(
θ
)
=
1; then,(46)
D
=
1
+
α
1
∗
α
¯
1
+
β
1
∗
β
¯
1
-
τ
0
α
1
∗
α
¯
1
θ
e
-
i
ω
2
τ
0
-
1
.
Next, the coordinates to describe the center manifoldC
0 at μ
=
0 will be calculated. Assuming that z
t is the root of (33) at μ
=
0, define(47)
s
t
=
q
∗
,
z
t
,
W
t
,
θ
=
z
t
θ
-
2
R
e
s
t
q
θ
.Equation (47) satisfied the following function in center manifold C
0:(48)
W
t
,
θ
=
W
s
t
,
s
¯
t
,
θin which W
(
s
,
s
¯
,
θ
)
=
W
20
(
θ
)
s
2
/
2
+
W
11
(
θ
)
s
s
¯
+
W
02
s
¯
2
/
2
+
W
30
(
θ
)
s
3
/
6
+
⋯, s
a
n
d
s
¯ are the local coordinates of center manifold C
0 that, in the directions of q
∗
a
n
d
q
¯
∗, if z
t is a real root of (39), when z
t
∈
C
0, μ
=
0, the following equation can be obtained:(49)
s
˙
t
=
i
ω
2
s
+
q
¯
∗
θ
f
0
,
W
s
,
s
¯
,
0
+
2
R
e
s
q
θ
=
i
ω
2
s
+
q
¯
∗
0
f
0
,
W
s
,
s
¯
,
0
+
2
R
e
s
q
0
=
d
e
f
i
ω
2
s
+
q
¯
∗
0
f
s
,
s
¯
.The equation can be written as(50)
s
˙
t
=
i
ω
2
s
t
+
g
s
,
s
¯
,in which(51)
g
s
,
s
¯
=
q
¯
∗
0
f
0
s
,
s
¯
=
g
20
s
2
2
+
g
11
s
s
¯
+
g
02
s
¯
2
2
+
g
21
θ
z
2
z
¯
2
+
⋯due to z
t
(
θ
)
=
(
z
1
t
(
θ
)
,
z
2
t
(
θ
)
,
z
3
t
(
θ
)
)
T
=
W
(
t
,
θ
)
+
s
q
(
θ
)
+
s
¯
q
(
θ
), q
(
θ
)
=
(
1
,
α
1
,
β
1
)
T
e
i
θ
ω
2; then(52)
z
1
t
0
=
s
+
s
¯
+
W
20
1
0
s
2
2
+
W
11
1
s
s
¯
+
W
02
1
s
¯
2
2
+
O
s
,
s
¯
3
,
z
2
t
0
=
α
s
+
α
¯
1
s
¯
+
W
20
2
0
s
2
2
+
W
11
2
s
s
¯
+
W
02
2
s
¯
2
2
+
O
s
,
s
¯
3
,
z
3
t
0
=
β
s
+
β
¯
1
s
¯
+
W
20
3
0
s
2
2
+
W
11
3
s
s
¯
+
W
02
3
s
¯
2
2
+
O
s
,
s
¯
3
.According to (51), we can get that(53)
g
s
,
s
¯
=
q
¯
∗
0
f
0
s
,
s
¯
=
D
¯
1
,
α
¯
1
∗
,
β
¯
1
∗
-
η
1
z
1
t
0
z
3
t
0
-
β
z
2
t
0
2
-
η
2
z
2
t
0
z
3
t
0
0
=
-
D
¯
η
1
s
+
s
¯
+
W
20
1
0
s
2
2
+
W
11
1
0
s
s
¯
+
W
02
1
s
¯
2
2
+
O
s
,
s
¯
3
β
1
s
+
β
¯
1
s
¯
+
W
20
3
0
s
2
2
+
W
11
3
0
s
s
¯
+
W
02
3
0
s
¯
2
2
+
O
s
,
s
¯
3
-
D
¯
α
¯
1
∗
β
α
1
s
+
α
¯
1
s
¯
+
W
20
2
0
s
2
2
+
W
11
2
s
s
¯
0
+
W
02
2
0
s
¯
2
2
+
O
s
,
s
¯
3
α
1
s
+
α
¯
1
s
¯
+
W
20
2
0
s
2
2
+
W
11
2
s
s
¯
0
+
W
02
2
0
s
¯
2
2
+
O
s
,
s
¯
3
-
D
¯
α
¯
1
∗
η
2
α
1
s
+
α
¯
1
s
¯
+
W
20
2
0
s
2
2
+
W
11
2
0
s
s
¯
+
W
02
2
0
s
¯
2
2
+
O
s
,
s
¯
3
β
1
s
+
β
¯
1
s
¯
+
W
20
3
0
s
2
2
+
W
11
3
0
s
s
¯
+
W
02
3
0
s
¯
2
2
+
O
s
,
s
¯
3compared with the parameter of (51); then we get(54)
g
20
=
-
2
D
¯
η
1
β
1
+
β
α
¯
1
∗
α
1
2
+
η
2
α
¯
1
∗
α
1
β
1
,
g
11
=
-
D
¯
η
1
β
1
+
β
¯
1
-
2
D
¯
α
1
∗
β
α
1
α
¯
1
-
D
¯
α
¯
1
∗
η
2
α
1
β
¯
1
+
α
¯
1
β
1
,
g
02
=
-
2
D
¯
η
1
β
¯
1
+
β
α
¯
1
∗
α
¯
1
2
+
η
2
α
¯
1
∗
α
¯
1
β
¯
1
,
g
21
=
-
D
¯
η
1
W
11
3
0
+
W
20
3
0
+
β
¯
1
W
20
1
0
+
β
1
W
11
1
0
-
2
D
¯
α
1
∗
β
1
α
1
W
11
2
0
+
α
¯
1
W
20
2
0
-
D
¯
α
¯
1
∗
η
2
α
1
W
11
3
0
+
α
¯
1
W
20
3
0
+
β
¯
1
W
20
2
0
+
β
1
W
11
2
0 In order to calculate g
21, we need to get W
20
(
θ
) and W
11
(
θ
) by (39) and (41):(55)
W
˙
=
s
˙
t
-
s
˙
q
-
s
¯
˙
q
˙
=
B
W
-
2
R
e
q
¯
∗
0
f
0
q
θ
,
θ
∈
-
τ
0
,
0
,
B
W
-
2
R
e
q
¯
∗
0
f
0
q
0
+
f
0
,
θ
=
0
=
d
e
f
B
W
+
H
s
,
s
¯
,
θin which(56)
H
s
,
s
¯
,
θ
=
H
20
θ
s
2
2
+
H
11
θ
z
z
¯
+
H
02
θ
s
2
2
+
⋯
.By comparing the parameters, then(57)
B
-
2
i
ω
2
W
20
θ
=
H
-
20
θ
,
B
W
11
θ
=
-
H
11
θ
⋯
.According to (51) and (55), when θ
∈
-
τ
0
,
0,(58)
H
s
,
s
¯
,
θ
=
-
q
¯
∗
0
f
0
q
θ
-
q
∗
0
f
¯
0
q
¯
θ
=
-
g
q
θ
-
g
¯
q
¯
θ
.Considering the parameters of (56), then(59)
H
20
θ
=
-
g
20
q
θ
-
g
¯
02
q
¯
θ
,
H
11
θ
=
-
g
11
q
θ
-
g
¯
11
q
¯
θ
.In the same way,(60)
W
11
θ
=
i
g
11
ω
2
q
0
e
i
θ
ω
2
+
i
g
¯
11
3
ω
2
q
¯
0
e
-
i
θ
ω
2
+
E
2
. Calculate vectors E
1 and E
2 and let θ
=
0, by (51), (56), and (57) and we can get(61)
H
20
0
=
-
g
20
q
0
-
g
¯
02
q
¯
0
+
2
0
-
β
1
α
1
,
(62)
H
11
0
=
-
g
11
q
0
-
g
¯
11
q
¯
0
+
2
0
-
R
e
β
1
R
e
α
1
.By the definition of (57), then(63)
-
r
1
a
0
0
-
r
2
0
0
0
-
h
W
20
0
+
-
b
0
0
b
0
0
0
0
0
W
20
-
τ
0
=
2
i
W
20
0
-
H
20
0
,
(64)
-
r
1
a
0
0
-
r
2
0
0
0
-
h
W
11
0
+
-
b
0
0
b
0
0
0
0
0
W
11
-
τ
0
=
-
H
11
0
. Put (60) and (61) to (65), and we can get(65)
2
i
ω
2
+
r
1
+
b
e
-
2
i
ω
2
-
a
0
-
b
e
-
2
i
ω
2
2
i
ω
2
+
r
2
0
0
0
2
i
ω
2
+
h
E
1
=
2
0
-
β
1
α
1 Let E
1
=
E
1
(
1
)
,
E
1
(
2
)
,
E
1
(
3
)
T; according to (66), we can get(66)
E
1
1
=
2
a
β
1
a
e
-
2
i
ω
2
-
r
2
+
2
i
ω
2
r
1
+
b
e
-
2
i
ω
2
+
2
i
ω
2
,
E
1
2
=
2
β
1
r
1
+
b
e
-
2
i
ω
2
+
2
i
ω
2
r
2
+
2
i
w
2
r
1
+
b
e
-
2
i
ω
2
+
2
i
w
2
-
a
e
-
2
i
ω
2
,
E
1
3
=
2
α
1
2
i
ω
2
+
h
. In the same way, put (61) and (63) to (65), we can get(67)
r
1
+
b
-
a
0
-
b
r
2
0
0
0
h
E
2
=
2
0
-
R
e
β
1
R
e
α
1
.Let E
2
=
E
2
(
1
)
,
E
2
(
2
)
,
E
2
(
3
)
T; the root of (67) can be obtained (68)
E
2
1
=
2
a
R
e
β
1
a
b
-
b
r
2
-
r
1
r
2
,
E
2
2
=
2
b
+
r
1
R
e
β
1
a
b
-
b
r
2
-
r
1
r
2
,
E
2
3
=
2
R
e
α
1
h
. According to E
1, E
2, (60) and (61), W
20
(
θ
)
a
n
d
W
11
(
θ
) can be easily obtained. Furthermore, we can see that each g
i
j in (54) is determined by parameters and delays. Thus, we can compute the following quantities.
We have the following results:(69)
C
1
0
=
i
2
ω
2
g
11
g
20
-
2
g
11
2
-
g
02
2
2
+
g
21
2
,
μ
2
=
-
R
e
C
1
0
R
e
λ
′
τ
0
,
β
2
=
2
R
e
C
1
0
,
τ
2
=
-
I
m
C
1
0
+
μ
2
I
m
λ
′
τ
0
ω
2
.This completes the proof.


## 5. Numerical Simulation

By selecting some related data from China environment protection database and doing the appropriate treatment, the following parameters can be obtained [19, 20]: a
=
2, b
=
1, c
=
1, p
=
1, θ
=
1, η
1
=
0.2, η
2
=
0.2, h
=
1, β
=
0.2, r
1
=
0.2, r
2
=
0.2, and m
=
0.2.Then, system (3) can be written as follows:(70)
x
˙
t
=
2
y
t
-
x
t
-
τ
-
0.2
x
t
-
0.2
u
t
x
t
,
y
˙
t
=
x
t
-
τ
-
0.3
y
t
-
0.2
y
2
t
-
E
t
y
t
-
0.2
u
t
y
t
,
u
˙
t
=
1
-
u
t
,
0
=
E
t
y
t
-
1
-
0.2
.
P
(
2.5649,2.821,1
,
0.1821
) is the unique positive equilibrium of system (70). ω
=
1.3544, τ
j
=
0.6043
+
5.4566
j, τ
0
=
0.6043, and τ
0 is the bifurcation parameter of system (70). When τ
=
0.2
<
0.6043, the dynamical responses of system (70) are shown by Figure 1, system (70) is stable at P, and the population and economic profits develop sustainably in this case; when τ
=
1
>
0.6043, the dynamical responses of system (70) are shown by Figure 2, system (70) is unstable, and the population and economic profits cannot develop sustainably in this case; Figure 3 shows that the dynamic behavior of the population changes with time delay and the Hopf bifurcation exists when τ
>
0.6043.Figure 1
Dynamical responses of system (70) when τ
=
0.2
<
0.6043.Figure 2
Dynamical responses of system (70) when τ
=
1
>
0.6043.Figure 3
The bifurcation diagram of system (70) when τ
∈
[
0,1.42
].


## 6. Conclusion

Based on the mathematical biology theory, the Hopf bifurcation theory of differential system, and the singular system theory, this paper considers a singular biological economic system with time delay in a polluted environment. The Hopf bifurcation occurs at the positive equilibrium with the change of time delay. We can proof that time delay has a great influence on the development of the population and economic development. In order to make the population development sustainable and ensure the maximization of economic benefits, the properties of Hopf bifurcation is necessary to be studied.



---

*Source: 5036305-2016-10-17.xml*