A Mathematical Model for the Dynamics of a Fish Algae Consumption Model with Impulsive Control Strategy

Abstract

A dynamic mathematical model of fish algae consumption with an impulsive control strategy is proposed and analyzed in detail. It is shown that the system has a globally asymptotically stable algae-eradication periodic solution which can be obtained using the Floquet theory of impulsive differential equations and small-amplitude perturbation techniques. The conditions for the permanence of the system can also be determined. Numerical results for impulsive perturbations show the rich dynamic behavior of the system. All these results may be useful in controlling eutrophication.

1. Introduction

One of the most common ecological and environmental problems of inland water bodies is eutrophication, which diminishes water quality by spurring excessive growth of algae and increasing amounts of suspended organic material [1–3]. Because of climate change and discharges of industrial wastewater into rivers, the self-purification capacity of the water body decreases. This provides favorable conditions for algae growth, leading to repeated outbreaks of abnormal proliferations of algae blooms in tributaries; recently, this phenomenon has been increasing [4]. Therefore, research on how to control the amount of algae in the environment is of great theoretical importance and practical significance.

There are many ways to control algae populations. The most widely used method is chemical control; this is not discussed in detail here because it has many negative impacts. More recently, another important method, called biological control, has been introduced. Biological control is the practice of using natural enemies (fish) to suppress algae, as has already been done for pest control [5–8]. These natural enemies can prevent algae populations from increasing to levels that cause eutrophication.

Mathematical models of ecological population dynamics have not only to account for growth and interactions, but also to improve the understanding of the functioning of food chains and webs and their dependence on some environmental conditions [9–12]. Impulsive differential equations have been widely used to study the mathematical properties of impulsive predator-prey models or three-species food-chain models [13–17]. Moreover, the theory of impulsive differential equations is being recognized, not only to be richer than the corresponding theory of differential equations without impulses, but also to represent a more natural framework for the mathematical modeling of real-world phenomena [18–20].

Recently, because of eutrophication, the nuisance algal (cyanobacteria and green algae) blooms frequently come forth in the Zeya reservoir which is located in Wenzhou. The algal blooms can cause clogging and blocking of the filtration system and result in millions of people without drinking water. In order to apply the biological principle to control algal blooms in the Zeya reservoir, we carry out good biological control agents for the algal blooms. Firstly, we bring up two species of filter-feeding fish, silver carp (Hypophthalmichthys molitrix) and bighead carp (Aristichthys nobilis), which are thought to be suited to control algal biomass directly in freshwater reservoir, especially cyanobacteria and green algae. Secondly, the decision maker must give cost-benefit and maintenance considerations; the fish should be harvested at an appropriate time. Finally, in order to prevent the rapid growth of the cyanobacteria and green algae, we must release a certain amount of filter-feeding fish at a fixed time. These processes are subject to short-term perturbations. Consequently, it is obvious to assume that these perturbations act instantaneously in the form of impulses. Hence, a fish algae consumption model can be described by the following differential equations:𝑑𝑥(𝑡)𝑑𝑡=𝑟¹𝑥(𝑡)𝐺⁰−𝑥(𝑡)𝐺₁−𝑥(𝑡)−𝑐¹𝑥²(𝑡)−𝑢¹𝑥(𝑡)𝑧(𝑡)𝑑𝑦(𝑡)𝑏₁+𝑥(𝑡),𝑑𝑡=𝑟²𝑦(𝑡)1−𝑦(𝑡)/𝑦^𝑚1−𝑦(𝑡)/𝑦_𝑛−𝑐²𝑦²(𝑡)−𝑢¹𝑎𝑦(𝑡)𝑧(𝑡)𝑥(𝑡)+𝑎𝑦(𝑡)+𝑏₂𝑑𝑧(𝑡),𝑡≠𝑛𝑇,𝑑𝑡=−𝑚𝑧(𝑡)+𝑢²𝑎𝑦(𝑡)𝑧(𝑡)𝑥(𝑡)+𝑎𝑦(𝑡)+𝑏₂+𝑢³𝑥(𝑡)𝑧(𝑡)ΔΔΔ𝑥𝑧𝑦(((𝑡𝑡𝑡)))===−00,,𝛿𝑧𝑡(𝑡=)+𝑛𝑇𝑝,,𝑏₁+𝑥(𝑡),(1.1) where 𝑥(𝑡), 𝑦(𝑡), 𝑧(𝑡) are the densities of the cyanobacteria population, green algae population, and fish population at time 𝑡, Δ𝑥(𝑡)=𝑥(𝑡⁺)−𝑥(𝑡);Δ𝑦(𝑡)=𝑦(𝑡⁺)−𝑦(𝑡), and Δ𝑧(𝑡)=𝑧(𝑡⁺)−𝑧(𝑡); 𝑟₁,𝑟₂ is a growth parameter which is related to the biological characteristics of populations and the rationalization of environmental resources; 𝑟₁𝐺₀  (0≤𝐺₀≤𝐺₁) is the carrying capacity of the cyanobacteria population 𝑥; 𝐺₁ is the limiting value of available resources; 𝑦_𝑚(0≤(𝑦_𝑚/𝑦_𝑛)≤1) is the maximum density of green algae population 𝑦 (i.e., the environmental carrying capacity); 𝑦_𝑛 is a nutritional parameter which is related to the resource conditions of the environment; 𝑢₁ is the cropping rate of fish; 𝑢₂ and 𝑢₃ are the average growth rate of fish that depends on the cyanobacteria and green algae for food; 𝑎 is the capture rate (as a proportion); 𝑚 is the average mortality rate for fish; 𝑏₁ and 𝑏₂ are the half-saturation constants; 𝑐₁and 𝑐₂ are the intraspecific competition rates of the cyanobacteria and green algae; 𝑇 is the period of the impulsive effect; 𝑁 is the set of all nonnegative integers; 𝛿(𝛿∈[0,1]) is the proportion of fish harvested at fixed times 𝑡=𝑛𝑇; and 𝑝>0 is the number of fish released at times 𝑡=𝑛𝑇.

The rest of this paper is arranged as follows: in the next section, some useful notations and basic properties of a system with impulsive effect are given, and then the main theorems and their proofs are developed using the Floquet theory. In Section 3, numerical results will be presented which exhibit the rich dynamic behaviors of system (1.1), and these results will be briefly discussed. Conclusions are given in the last section.

2. Preliminaries and Mathematical Analysis

Let 𝑅₊=[0,∞),𝑅⁺³={𝑋∈𝑅³∣𝑋>0}. Denote by 𝑓=(𝑓₁,𝑓₂,𝑓₃) the mapping defined by the right-hand sides of the first, second, and third equations of system (1.1). Let 𝑉∶𝑅₊×𝑅³⁺→𝑅₊; then 𝑉is said to belong to class 𝑉₀ if(1)𝑉 is continuous in (𝑛𝑇,(𝑛+1)𝑇]×𝑅³⁺, and for each 𝑋∈𝑅³⁺,𝑛∈𝑁, lim_{(𝑡,𝑦)→(𝑛𝑇}⁺_,𝑋)𝑉(𝑡,𝑦)=𝑉(𝑛𝑇⁺,𝑋) exists.(2)𝑉 is locally Lipschitzian in 𝑋.

Definition 2.1 (see [18]). Let 𝑉∈𝑉₀; then for (𝑡,𝑥)∈(𝑛𝑇,(𝑛+1)𝑇]×𝑅³⁺, the upper right derivative of 𝑉(𝑡) with respect to the impulsive differential system (1.1) is defined as 𝐷⁺𝑉(𝑡,𝑋)=lim_ℎ→0+sup1ℎ[𝑉(𝑡+ℎ,𝑋+ℎ𝑓(𝑡,𝑋))−𝑉(𝑡,𝑋)].(2.1)

Remark 2.2. (1) The solution of system (1.1) is a piecewise continuous function with 𝑋∶𝑅₊→𝑅³⁺; 𝑋(𝑡) is continuous on (𝑛𝑇,(𝑛+1)𝑇], and 𝑋(𝑛𝑇⁺)=lim_{(𝑡→𝑛𝑇}⁺₎𝑋(𝑛𝑇) exists.
(2) The smoothness properties of 𝑓 guarantee the global existence and uniqueness of solutions of system (1.1) (for details, see [18, 19]).
Because 𝑥^(𝑡)=𝑦′(𝑡)=𝑧′(𝑡)=0 whenever 𝑥(𝑡)=𝑦(𝑡)=𝑧(𝑡)=0   (𝑡≠𝑛𝑇) and 𝑧(𝑛𝑇⁺)=(1−𝛿)𝑧(𝑛𝑇)+𝑝, the following lemma can be stated.

Lemma 2.3. Assume 𝑋(𝑡) is a solution of system (1.1) such that(1)if𝑋(0⁺)≥0, then 𝑋(𝑡)≥0  for all 𝑡≥0,(2)if𝑋(0⁺)>0, then 𝑋(𝑡)>0  for all 𝑡>0.
Next, an important comparison theorem on impulsive differential equations [18] will be stated.

Lemma 2.4 (comparison theorem). Assume that 𝑉∈𝑉₀ and that 𝑉𝐷𝑡⁺,𝑋𝑉(𝑡𝑡,₊𝑋)≤𝑔(𝑡,𝑉(𝑡,𝑋)),𝑡≠𝑛𝑇,≤𝜑_𝑛(𝑉(𝑡,𝑋)),𝑡=𝑛𝑇,(2.2) where 𝑔∶𝑅₊×𝑅₊→𝑅 is continuous on (𝑛𝑇,(𝑛+1)𝑇]×𝑅³⁺ and for 𝑢∈𝑅₊, 𝑛∈𝑁, lim_{(𝑡,𝑦)→(𝑛𝑇}⁺_,𝑢)𝑔(𝑡,𝑣)=𝑔(𝑛𝑇⁺,𝑢) exists, with 𝜑_𝑛∶𝑅₊→𝑅₊ nondecreasing. Let 𝑟(𝑡) be the maximal solution of the scalar impulsive differential equation: 𝑑𝑢(𝑡)𝑢𝑑𝑡𝑡=𝑔(𝑡,𝑢(𝑡)),𝑡≠𝑛𝑇,₊𝑢0=𝜑_𝑛(𝑢(𝑡)),𝑡=𝑛𝑇,₊=𝑢₀,(2.3) existing on [0,∞). Then 𝑉(0⁺,𝑋₀)≤𝑢₀, which implies that 𝑉(𝑡,𝑋(𝑡))≤𝑟(𝑡), 𝑡≥0, where 𝑋(𝑡) is any solution of system (1.1).
If the cyanobacteria and green algae are eradicated, then system (1.1) will reduce to the following system: 𝑑𝑧(𝑡)𝑧𝑑𝑡𝑡=−𝑚𝑧(𝑡),𝑡≠𝑛𝑇,₊𝑧0=(1−𝛿)𝑧(𝑡)+𝑝,𝑡=𝑛𝑇,₊=𝑧₀.(2.4)

Therefore, the following result can be easily obtained.

Lemma 2.5. 𝑧^∗(𝑡)=𝑝exp(−𝑚(𝑡−𝑛𝑇))/(1−(1−𝛿)exp(−𝑚𝑇)) is a positive periodic solution of (2.4) with initial value 𝑧^∗(0⁺)=𝑝/(1−(1−𝛿)exp(−𝑚𝑇)),𝑡∈(𝑛𝑇,(𝑛+1)𝑇],𝑛∈𝑁.
𝑧(𝑡)=(𝑧(0⁺)−𝑝/(1−(1−𝛿)exp(−𝑚𝑇)))exp(−𝑚𝑡)+𝑧^∗(𝑡) is a general solution of (2.4) with initial value 𝑧₀=𝑧(0⁺)≥0, where 𝑡∈(𝑛𝑇,(𝑛+1)𝑇], 𝑛∈𝑁.
For a positive periodic solution 𝑧^∗(𝑡) of (2.4) and every solution 𝑧(𝑡) of (2.4) with 𝑧₀≥0, |𝑧(𝑡)−𝑧^∗(𝑡)|→0, 𝑡→∞.

From Lemma 2.5, the general solution 𝑧(𝑡) converges to the positive periodic solution 𝑧^∗(𝑡) of (2.4) when 𝑡→∞, and the complete expression for the cyanobacteria and green algae-eradication periodic solution (0,0,𝑧^∗(𝑡)) of system (1.1) becomes0,0,𝑧^∗(𝑡)=0,0,1𝑝−e(x1p−(−𝛿)𝑚e(x𝑡p−(−𝑛𝑇𝑚)𝑇)),(2.5) for 𝑡∈(𝑛𝑇,(𝑛+1)𝑇].

Now, the Floquet theorem for impulsive equations [20] will be introduced:𝑑𝑥(𝑡)𝑑𝑡=𝐴(𝑡)𝑥(𝑡),𝑡≠𝑡^𝑘𝑥𝑡₊,=𝑥(𝑡)+𝐵_𝑘𝑥(𝑡),𝑡=𝑡_𝑘,(2.6) with conditions as follows(1)𝐴(⋅)∈PC(𝑅,𝐶^𝑛×𝑛) and 𝐴(𝑡+𝑇)=𝐴(𝑡)(𝑡∈𝑅), where PC(𝑅,𝐶^𝑛×𝑛) is the set of all piecewise continuous matrix functions and 𝐶^𝑛×𝑛 is the set of 𝑛×𝑛 matrices.(2)𝐵_𝑘∈(𝐶^𝑛×𝑛) and det(𝐸+𝐵_𝑘)≠0, 𝑡_𝑘≤𝑇_𝑘+1.(3)There exists a 𝑞∈𝑁 such that 𝐵_𝑘+𝑞=𝐵_𝑘, 𝑡_𝑘+𝑞=𝑡_𝑘+𝑇.

Let 𝐻(𝑡) be a fundamental matrix of system (2.6), then there exists a unique nonsingular matrix 𝑀∈𝐶^𝑛×𝑛 such that 𝐻(𝑡+𝑇)=𝐻(𝑡)𝑀. 𝑀 is similar to the matrix of system (2.6) and has the same eigenvalues 𝑢₁,𝑢₂,…,𝑢_𝑛. These eigenvalues 𝑢₁,𝑢₂,…,𝑢_𝑛 of system (2.6) are called Floquet multipliers.

Lemma 2.6 (Floquet theorem [20]). Assume that (1)–(3) hold, then system (2.6) is(1)stable if and only if 𝑢_𝑖≤1(𝑖=1,…,𝑛), and to those 𝑢_𝑖=1, there correspond simple elementary divisors;(2)asymptotically stable if and only if all 𝑢_𝑖<1(𝑖=1,…,𝑛);(3)unstable if 𝑢_𝑖>1 for some 𝑖=1,…,𝑛.

After these preliminaries, the main theorems of this paper will now be presented, then it will be proven that each solution of system (1.1) is ultimately bounded.

Theorem 2.7. There exists a positive constant 𝑀 such that 𝑥(𝑡)≤𝑀, 𝑦(𝑡)≤𝑀, and 𝑧(𝑡)≤𝑀 for each solution of system (1.1) with all 𝑡 large enough.

Proof. Define 𝑉(𝑡,𝑥(𝑡)) such that 𝑉(𝑡,𝑥(𝑡))=𝑥(𝑡)+𝑢₂𝑦(𝑡)+𝑢₁𝑧(𝑡).(2.7)
Because (𝑑𝑥(𝑡)/𝑑𝑡)≤𝑟₁𝑥(𝑡)−𝑐₁𝑥²(𝑡), then 𝑥(𝑡)≤(𝑟₁/𝑐₁), and the upper right derivative of 𝑉(𝑡,𝑥(𝑡)) can be calculated along a solution of system (1.1) when 𝑡≠𝑛𝑇: 𝐷⁺(𝑉(𝑡))+𝐿𝑉(𝑡)≤−𝐿𝑐₂+𝑟₁𝑥(𝑡)−𝑐₁𝑥²(𝑡)+𝐿𝑢₂+𝑢₂𝑟₂𝑦(𝑡)𝑢₂𝑦²(𝑡)+𝐿𝑢₁−𝑢₁𝑚+𝑢¹𝑟₁𝑢₃𝑏₁𝑐₁𝑧(𝑡).(2.8)
Let0<𝐿<𝑚−(𝑟₁𝑢₃/𝑏₁𝑐₁), then (𝐿+𝑟₁)𝑥(𝑡)−𝑐₁𝑥²(𝑡), (𝐿𝑢₂+𝑢₂𝑟₂)𝑦(𝑡)−𝑐₂𝑢₂𝑦²(𝑡) are both bounded, so 𝐷⁺(𝑉(𝑡))+𝐿𝑉(𝑡) is bounded, and 𝐿₁ and 𝐿₂ can be chosen such that 𝑑𝐷⁺(𝑉(𝑡))𝑑𝑡≤−𝐿¹𝑉𝑡𝑉(𝑡)+𝐿₂,𝑡≠𝑛𝑇,₊≤𝑉(𝑡)+𝑢₁𝑝,𝑡=𝑛𝑇,(2.9) where 𝐿₁ and 𝐿₂ are positive constants. According to Lemma 2.4, 𝑉(𝑡)≤𝑉0⁺−𝐿²𝐿₁exp−𝐿¹𝑡+𝑢¹𝑝1−exp−𝐿₁𝑛𝑇1−exp−𝐿₁𝑇exp−𝐿¹(𝑡−𝑛𝑇)+𝐿²𝐿₁,(2.10) and hence 𝑉(𝑡)≤𝐿²𝐿₁+𝑢¹𝑝exp𝐿₁𝑇exp𝐿₁𝑇−1,(2.11) as 𝑡→∞, so 𝑉(𝑡,𝑥(𝑡)) is ultimately bounded. Therefore, each positive solution of system (1.1) is ultimately bounded, and the proof is completed.

Theorem 2.8. If 𝑟₁𝐺₀𝐺₁𝑇<𝑢¹𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₁(1−(1−𝛿)exp(−𝑚𝑇))₂,𝑇<𝑢¹𝑎𝑝(1−exp(−𝑚𝑇))𝑚𝑏₂(1−(1−𝛿)exp(−𝑚𝑇)),(2.12) then the periodic solution (0,0,𝑧^∗(𝑡)) is locally asymptotically stable.

Proof. The local stability of the periodic solution (0,0,𝑧^∗(𝑡)) may be determined by considering the behavior of small-amplitude perturbations of the solution.
Define 𝑥(𝑡)=𝑢(𝑡),𝑦(𝑡)=𝑣(𝑡),𝑧(𝑡)=𝑤(𝑡)+𝑧^∗(𝑡).
Then the linearization of system (1.1) becomes 𝑑𝑢(𝑡)𝑑𝑡=𝑟¹𝐺₀𝐺₁−𝑢¹𝑧^∗(𝑡)𝑑𝑣(𝑡)𝑏₁𝑢(𝑡),𝑑𝑡=𝑟²−𝑢¹𝑎𝑧^∗(𝑡)𝑑𝑤(𝑡)𝑏₂𝑣(𝑡),𝑡≠𝑛𝑇,𝑑𝑡=𝑢³𝑧^∗(𝑡)𝑏₁𝑢(𝑡)+𝑢²𝑎𝑧^∗(𝑡)𝑏₂ΔΔΔ𝑤𝑣𝑢(((𝑡𝑡𝑡)))===00−,,𝛿𝑣(𝑡),𝑡=𝑛𝑇,𝑣(𝑡)−𝑚𝑤(𝑡),(2.13) and, as a result, ⎛⎜⎜⎜⎜⎝𝑤𝑣𝑢(((𝑡𝑡𝑡)))⎞⎟⎟⎟⎟⎠=Φ(𝑡)⎛⎜⎜⎜⎜⎝𝑤𝑣𝑢(((000)))⎞⎟⎟⎟⎟⎠,0≤𝑡≤𝑇,(2.14) where Φ(𝑡) satisfies 𝑑Φ(𝑡)𝑑𝑡=⎛⎜⎜⎜⎜⎜⎜⎝𝑟¹𝐺₀𝐺₁−𝑢¹𝑧^∗(𝑡)𝑏₁0𝑟₂00−𝑢¹𝑎𝑧^∗(𝑡)𝑏₂𝑢₃0𝑧^∗(𝑡)𝑏₁𝑢₂𝑎𝑧^∗(𝑡)𝑏₂−𝑚⎞⎟⎟⎟⎟⎟⎟⎠Φ(𝑡),(2.15)Φ(0)=𝐼, the identity matrix, and ⎛⎜⎜⎜⎜⎝𝑣𝑢𝑛𝑛𝑇𝑇⁺₊𝑤𝑛𝑇₊⎞⎟⎟⎟⎟⎠=⎛⎜⎜⎜⎜⎝1000101−00𝛿⎞⎟⎟⎟⎟⎠⎛⎜⎜⎜⎜⎝𝑤𝑣𝑢(((𝑛𝑛𝑛𝑇𝑇𝑇)))⎞⎟⎟⎟⎟⎠.(2.16)
Therefore, the stability of the periodic solution (0,0,𝑧^∗(𝑡)) is determined by the eigenvalues 𝜃=⎛⎜⎜⎜⎜⎝1000101−00𝛿⎞⎟⎟⎟⎟⎠Φ(𝑇).(2.17)
Therefore, all eigenvalues of 𝜃 are given by 𝜆₁=exp^𝑇0𝑟¹𝐺₀𝐺₁−𝑢¹𝑧^∗(𝑡)𝜆𝑏₁𝑑𝑡₂,=exp^𝑇0𝑟²−𝑢¹𝑎𝑧^∗(𝑡)𝜆𝑏₂𝑑𝑡₃,=(1−𝛿)exp(−𝑚𝑇)<1.(2.18)
According to Lemma 2.6, (0,0,𝑧^∗(𝑡)) is locally asymptotically stable if 𝜆₁<1,𝜆₂<1, that is to say 𝑟₁𝐺₀𝐺₁𝑇<𝑢¹𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₁(1−(1−𝛿)exp(−𝑚𝑇))₂,𝑇<𝑢¹𝑎𝑝(1−exp(−𝑚𝑇))𝑚𝑏₂(1−(1−𝛿)exp(−𝑚𝑇)),(2.19) which completes the proof.

Theorem 2.9. If 𝑟₁𝐺₀𝐺₁𝑇<𝑢¹𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₁(1−(1−𝛿)exp(−𝑚𝑇))₂,𝑇<𝑢¹𝑎𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₂(1−(1−𝛿)exp(−𝑚𝑇))₁,−𝑢¹𝑟𝑏₁+𝑀1−(1𝑝e−xp𝛿)(e−x𝑚p𝑇(−)𝑚𝑇)<0,₂−𝑢¹𝑎𝑀+𝑎𝑀+𝑏₂1−(1𝑝e−xp𝛿)(e−x𝑚p𝑇(−)𝑚𝑇)<0,(2.20) then the periodic solution (0,0,𝑧^∗(𝑡)) is globally asymptotically stable.

Proof. From Theorem 2.8, the periodic solution (0,0,𝑧^∗(𝑡)) is locally asymptotically stable. Now it will be proved to be a global attractor. Let 𝑉(𝑡)=𝑥(𝑡)+𝑦(𝑡), then by Theorem 2.7, there exists a constant 𝑀>0 such that 𝑥(𝑡)≤𝑀, 𝑦(𝑡)≤𝑀, and 𝑧(𝑡)≤𝑀. Obviously, 𝑉^′||_(1.1)≤𝑥(𝑡)𝑟₁−𝑢¹𝑧(𝑡)𝑏₁+𝑀−𝑐¹𝑥²(𝑡)+𝑦(𝑡)𝑟₂−𝑢¹𝑎𝑧(𝑡)𝑀+𝑎𝑀+𝑏₂−𝑐²𝑦²(𝑡).(2.21)
However, for system (1.1), 𝑑𝑧(𝑡)𝑧𝑑𝑡𝑡≥−𝑚𝑧(𝑡),𝑡≠𝑛𝑇,₊=(1−𝛿)𝑧(𝑡)+𝑝,𝑡=𝑛𝑇.(2.22)
Therefore, there exists a 𝜏>0, and an 𝜀>0 small enough, such that 𝑧(𝑡)≥𝑧^∗1(𝑡)−𝜀 for all 𝑡>𝜏. This leads to 𝑧(𝑡)≥𝑧^∗1(𝑡)−𝜀=𝑝exp(−𝑚𝑇)1−(1−𝛿)exp(−𝑚𝑇)−𝜀.(2.23)
Let 𝛾=^Δ𝑝exp(−𝑚𝑇)/(1−(1−𝛿)exp(−𝑚𝑇))−𝜀; it is well known that 𝑟₁−(𝑢₁𝛾/(𝑏₁+𝑀))<0 and 𝑟₂−(𝑢₁𝑎𝛾/(𝑀+𝑎𝑀+𝑏₂))<0. Therefore, when 𝑡>𝜏, 𝑉^||_(1.1)≤𝑥(𝑡)𝑟₁−𝑢¹𝑧(𝑡)𝑏₁+𝑀−𝑐¹𝑥²(𝑡)+𝑦(𝑡)𝑟₂−𝑢¹𝑎𝑧(𝑡)𝑀+𝑎𝑀+𝑏₂−𝑐²𝑦²(𝑡)<0.(2.24)
Therefore, 𝑉(𝑡)→0, 𝑥(𝑡)→0, 𝑦(𝑡)→0 as 𝑡→∞, and it follows that the periodic solution (0,0,𝑧^∗(𝑡)) is a global attractor. This completes the proof.

Theorem 2.10. The system (1.1) is permanent, if 𝑟₁𝐺₀𝐺₁𝑇>𝑢¹𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₁(1−(1−𝛿)exp(−𝑚𝑇))₂,𝑇>𝑢¹𝑎𝑝(1−exp(−𝑚𝑇))𝑟𝑚𝑏₂(1−(1−𝛿)exp(−𝑚𝑇))₁,−𝑢¹𝑟𝑏₁+𝑀1−(1𝑝e−xp𝛿)(e−x𝑚p𝑇(−)𝑚𝑇)>0,₂−𝑢¹𝑎𝑀+𝑎𝑀+𝑏₂1−(1𝑝e−xp𝛿)(e−x𝑚p𝑇(−)𝑚𝑇)>0.(2.25)

Proof. Assume that 𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) is any solution of system (1.1) with 𝑋(0)>0. From Theorem 2.7, 𝑥(𝑡)≤𝑀, 𝑦(𝑡)≤𝑀, and 𝑧(𝑡)≤𝑀 with 𝑡≥0. From Lemma 2.5, 𝑧(𝑡)≥𝑧^∗(𝑡)−𝜀 for all 𝑡 large enough and some 𝑡 with 𝑧(𝑡)>𝑝exp(−𝑚𝑇)/(1−(1−𝛿)exp(−𝑚𝑇))−𝜀=𝜉₁. Therefore, it is necessary to find 𝜉₂ and 𝜉₃ such that 𝑥(𝑡)≥𝜉₂ and 𝑦(𝑡)≥𝜉₃ for 𝑡 large enough. Select an 𝜀₁ small enough so that 𝜑₁=exp^{(𝑛+1)𝑇𝑛𝑇}𝑟¹𝐺₀𝐺₁−𝑟¹𝑀𝐺₁−𝑀−𝑐¹𝑀−𝑢¹𝑏₁𝑣^∗(𝑡)+𝜀₁𝑑𝑡>1.(2.26)
Next, it is necessary to prove that 𝑥(𝑡)≥𝜉₂ for 𝑡 large enough, where 𝜉₂ is a positive constant.
It is easy to prove that there exists a 𝑡₁∈(0,∞) such that 𝑥(𝑡)≥𝜉₂. Otherwise, if 𝑥(𝑡)<𝜉₂ for all 𝑡>0, 𝑑𝑧(𝑡)𝑑𝑡≤𝑢²𝑎𝑀𝑏₂+𝑢³𝑀𝑏₁𝑧𝑡−𝑚𝑧(𝑡),𝑡≠𝑛𝑇,₊𝑧0=(1−𝛿)𝑧(𝑡)+𝑝,𝑡=𝑛𝑇,₊=𝑧₀,(2.27) and therefore 𝑧(𝑡)≤𝑣(𝑡) and 𝑣(𝑡)→𝑣^∗(𝑡), where 𝑣(𝑡) is the solution of 𝑑𝑣(𝑡)𝑑𝑡=𝑢²𝑎𝑀𝑏₂+𝑢³𝑀𝑏₁𝑣𝑡−𝑚𝑣(𝑡),𝑡≠𝑛𝑇,₊𝑣0=(1−𝛿)𝑣(𝑡)+𝑝,𝑡=𝑛𝑇,₊=𝑣₀,(2.28) with 𝑣^∗(𝑡)=𝑝exp−𝑚−𝑢²𝑎𝑀/𝑏₂−𝑢₃𝑀/𝑏₁(𝑡−𝑛𝑇)1−(1−𝛿)exp−𝑚−𝑢₂𝑎𝑀/𝑏₂−𝑢₃𝑀/𝑏₁𝑇,𝑡∈(𝑛𝑇,(𝑛+1)𝑇].(2.29)
Therefore, there exists a 𝑇₁>0 such that 𝑧(𝑡)≤𝑣(𝑡)≤𝑧^∗(𝑡)+𝜀₁, and 𝑑𝑥(𝑡)𝑑𝑡≥𝑥(𝑡)𝑟¹𝐺₀𝐺₁−𝑟¹𝑀𝐺₁−𝑀−𝑐¹𝑀−𝑢¹𝑏₁𝑣^∗(𝑡)+𝜀₁,(2.30) for 𝑡>𝑇₁. Let 𝑁₁∈𝑁 and 𝑁₁𝑇≥𝑇₂≥𝑇₁. Integrating (2.30) over (𝑛𝑇,(𝑛+1)𝑇]𝑛>𝑁₁ yields 𝑥((𝑛+1)𝑇)≥𝑥𝑛𝑇⁺exp^{(𝑛+1)𝑇𝑛𝑇}𝑟¹𝐺₀𝐺₁−𝑟¹𝑀𝐺₁−𝑀−𝑐¹𝑀−𝑢¹𝑏₁𝑣^∗exp_{(𝑛+1)𝑇𝑛𝑇}(𝑡)+𝜀₁𝑑𝑡=𝑥𝑛𝑇⁺,𝑟¹𝐺₀𝐺₁−𝑟¹𝑀𝐺₁−𝑀−𝑐¹𝑀−𝑢¹𝑏₁𝑣^∗(𝑡)+𝜀₁𝑑𝑡=𝑥𝑛𝑇⁺𝜑₁.(2.31)
Therefore, 𝑥((𝑁₁+𝑛)𝑇)≥𝑥(𝑁₁𝑇)𝜑^𝑛1→∞ as 𝑛→∞, but this leads to a contradiction because 𝑥(𝑡) is bounded. Therefore, there exists a 𝑡₁>0 such that 𝑥(𝑡₁)≥𝜉₂.
Second, if 𝑥(𝑡₁)≥𝜉₂ for all 𝑡≥𝑡₁, then the proof is complete. Hence, it is only necessary to consider those solutions which leave the region 𝑅={𝑥(𝑡)∶𝑥(𝑡)<𝜉₂} and reenter it. Let 𝑡^∗(𝑡)=inf_𝑡≥𝑡1{𝑥(𝑡)<𝜉₂}, then 𝑥(𝑡)≥𝜉₂, 𝑡∈(𝑡,𝑡^∗), and 𝑡^∗∈(𝑛₁𝑇,(𝑛₁+1)𝑇). However, obviously 𝑥(𝑡^∗)=𝜉₂ because 𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) is continuous. It is claimed here that there must exist a 𝑡₂∈((𝑛₁+1)𝑇,(𝑛₁+1)𝑇+𝑇₁) such that 𝑥(𝑡₂)≥𝜉₂. Otherwise, 𝑥(𝑡)<𝜉₂, for 𝑡∈((𝑛₁+1)𝑇,(𝑛₁+1)𝑇+𝑇₁),𝑇₁=𝑛₂𝑇+𝑛₃𝑇, 𝑛₂,𝑛₃∈𝑁, and then 𝑛₂−1𝑇>ln𝜀¹/(𝑀+𝑝)exp−𝜑₂𝑢₂𝑎𝑀/𝑏₂+𝑢₃𝑀/𝑏₁−𝑚,𝑛₂+1𝑇𝜑^𝑛3₁>1.(2.32)
Consider (2.28) with 𝑣(𝑡^∗+)=𝑧(𝑡^∗+); this leads to 𝑣(𝑡)=𝑣𝑛₁+1𝑇⁺−𝑝×exp−𝑢₂1−(1−𝛿)exp−𝑢₂𝑎𝑀/𝑏₂+𝑢₃𝑀/𝑏₁−𝑚𝑇𝑎𝑀𝑏₂+𝑢³𝑀𝑏₁−𝑚𝑡−𝑛₁+1𝑇+𝑣^∗(𝑡),(2.33) for 𝑡∈(𝑛𝑇,(𝑛+1)𝑇], 𝑛₁+1<𝑛<𝑛₁+𝑛₂+𝑛₃+1, and then ||𝑣(𝑡)−𝑣^∗(𝑡)||<(𝑀+𝑝)exp−𝑢²𝑎𝑀𝑏₂+𝑢³𝑀𝑏₁𝑧(𝑡)≤𝑣(𝑡)≤𝑣_∗−𝑚𝑡−𝑛₁+1𝑇<𝜀₁,(𝑡)+𝜀₁,(2.34) for 𝑡∈[(𝑛₁+𝑛₂+1)𝑇,(𝑛₁+1)𝑇+𝑇₁], which implies (2.30). Integrating (2.30) over 𝑡∈[(𝑛₁+𝑛₂+1)𝑇,(𝑛₁+1)𝑇+𝑇₁] yields 𝑥𝑛₁+𝑛₂+𝑛₃+1𝑇≥𝑥𝑛₁+𝑛₂+1𝑇𝜑^𝑛3₁.(2.35)
For every 𝑡∈(𝑡^∗,(𝑛₁+1)𝑇), there are only two possible cases: one, if 𝑥(𝑡)<𝜉₂ for 𝑡∈(𝑡^∗,(𝑛₁+1)𝑇), when 𝑡∈(𝑡^∗,(𝑛₁+𝑛₂+1)𝑇), then 𝑥(𝑡)<𝜉₂, and hence 𝑑𝑥(𝑡)𝑑𝑡≥𝑥(𝑡)𝑟¹𝐺₀𝐺₁−𝑟¹𝜉₂𝐺₁−𝑀−𝑐¹𝑀−𝑢¹𝑏₁𝑀=𝜑₂𝑥(𝑡).(2.36)
Integrating (2.36) over 𝑡∈(𝑡^∗,(𝑛₁+𝑛₂+1)𝑇),𝑥((𝑛₁+𝑛₂+1)𝑇)≥𝑥(𝑡^∗)exp(𝜑₂(𝑛₂+1)𝑇) is obvious.
Therefore,𝑥((𝑛₁+𝑛₂+𝑛₃+1)𝑇)≥𝜉₂exp(𝜑₂(𝑛₂+1)𝑇)𝜑^𝑛3₁>𝜉₂, but this leads to a contradiction. Hence, let 𝑡₃=inf_𝑡>𝑡∗{𝑥(𝑡)≥𝜉₂}, then 𝑥(𝑡₃)=𝜉₂, and (2.36) holds. Integrating (2.36) over [𝑡^∗,𝑡₃), 𝑥(𝑡)≥𝑥𝑡^∗exp𝜑₂𝑡−𝑡^∗≥𝜉₂exp𝜑₂𝑛₁+𝑛₃+1𝑇=𝜉₅.(2.37)
𝑥(𝑡₃)≥𝜉₂ is also true for 𝑡>𝑡₃. Hence, 𝑥(𝑡)≥𝜉₅, 𝑡>𝑡₃. Two, there exists a 𝑡₅∈(𝑡^∗,(𝑛₁+1)𝑇] such that 𝑥(𝑡₅)≥𝜉₂. Therefore, let 𝑡₄=inf_𝑡>𝑡∗{𝑥(𝑡)≥𝜉₂}, then 𝑥(𝑡)<𝜉₂ for 𝑡∈[𝑡^∗,𝑡₄) and 𝑥(𝑡₄)=𝜉₂. For 𝑡∈[𝑡^∗,𝑡₄), (2.36) holds, and integrating it over 𝑡∈[𝑡^∗,𝑡₄), 𝑥(𝑡)≥𝑥(𝑡^∗)exp(𝜑₂(𝑡−𝑡^∗))>𝜉₅. This process can be continued because 𝑥(𝑡₄)≥𝜉₂ and 𝑥(𝑡)≥𝜉₅, 𝑡>𝑡₄. This yields the two possible cases, 𝑥(𝑡)≥𝜉₅, 𝑡≥𝑡₁, and 𝑦(𝑡)≥𝜉₄, 𝑡≥𝑡₂. Let Ω={(𝑥,𝑦,𝑧)∶𝑥(𝑡)≥𝜉₂,𝑦(𝑡)≥𝜉₃,𝑧(𝑡)≥𝜉₁,𝑥(𝑡)+𝑦(𝑡)+𝑧(𝑡)≤3𝑀}. Then the set Ω is a global attractor. Moreover, every solution 𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (1.1) will eventually enter and remain in the set Ω. Therefore, system (1.1) is permanent, and the proof is completed.

3. Numerical Analysis

3.1. Bifurcation Analysis

In Section 2, it has been proved that the cyanobacteria and green algae-eradication periodic solution is locally asymptotically stable, and conditions have been established for the permanence of system (1.1). Now, to substantiate these theoretical results, the following parameters will be considered: 𝑢₁=0.5, 𝑢₂=0.6, 𝑢₃= 0.4, 𝑟₁=0.5, 𝑟₂=0.5, 𝐺₀=10, 𝐺₁=15, 𝑦_𝑚=15, 𝑦_𝑛=20, 𝑐₁=0.1, 𝑐₂=0.12, 𝑎=0.75, 𝑏₁=0.45, 𝑏₂=0.6, 𝛿=0.55, 𝑇=40 with initial values 𝑥₀=0.5, 𝑦₀=0.6, 𝑧₀=0.7.

From Theorem 2.8, it is known that the periodic solution (0,0,𝑧^∗(𝑡)) is locally asymptotically stable if 𝑝>𝑝_max=8.15; this periodic solution (0,0,𝑧^∗(𝑡)) is shown in Figure 1. It is clear that the variable fish 𝑧 oscillates in a stable cycle, but the cyanobacteria population 𝑥 and green algae population 𝑦 rapidly decrease to zero when 𝑝>𝑝_max. When 𝑝 is smaller than 𝑝_max, then the cyanobacteria and green algae-eradication periodic solution will become unstable, and it is possible that the cyanobacteria population, green algae population, and fish population can coexist.

(a)

(b)

(c)

Figure 1 

Dynamic behavior of system (1.1). When 𝑝>𝑝max=8.15, the cyanobacteria and green algae will become extinct. Time series evolving according to biological control system (1.1): (a) the cyanobacteria population 𝑥, (b) the green algae population 𝑦, and (c) the fish population 𝑧.

Figure 2 shows typical bifurcation diagrams for system (1.1) with respect to 𝑝 in the range 𝑝∈[0,10]. As 𝑝 increases, system (1.1) exhibits rich dynamic behavior, such as chaotic bands with periodic windows, crises, period-halving bifurcations, quasi-periodic oscillations, and narrow or wide periodic windows. When 𝑝<0.89, the cyanobacteria population 𝑥 enters a chaotic band. As 𝑝 increases beyond 0.89, the cyanobacteria population 𝑥 enters a chaotic band with periodic windows, crisis phenomena appear, and the chaotic attractor suddenly changes into a periodic attractor. Then the periodic attractor loses its stability, and system (1.1) again enters a chaotic band. Finally, the chaotic attractor changes into a periodic attractor again; details of these results are shown in Figure 3. When 1.15<𝑝<2.7, there is a cascade of period-halving bifurcation which leads system (1.1) to a 𝑇-periodic solution (Figure 4), where the period-halving bifurcation is the opposite of the bifurcation observed earlier. When 𝑝>2.75, the green algae population 𝑦 decreases to zero because of the increasing number of fish population 𝑧 and the competition between the cyanobacteria population 𝑥 and green algae population 𝑦. Then the cyanobacteria population 𝑥 and fish population 𝑧 can coexist in a stable cycle, but the numbers of the cyanobacteria population 𝑥 will decrease because of short supply of resources. When 𝑝>𝑝_max= 8.15 the cyanobacteria population 𝑥 will be eradicated because the number of fish population 𝑧 released is continually increasing, then the cyanobacteria and green algae-eradication periodic solution occurs. Once the cyanobacteria population𝑥 and green algae population 𝑦 have decreased to zero, fish population 𝑧 will be eradicated after a period of time because of lack of food. Therefore, it is apparent that the value of 𝑝 can effectively control the size of the cyanobacteria population and green algae population. All these numerical simulations are consistent with the theoretical proofs presented earlier.

(a)

(b)

(c)

Figure 2 

Bifurcation diagrams for system (1.1) showing the effect of 𝑝. (a) 𝑥 versus 𝑝, (b) 𝑦 versus 𝑝, and (c) 𝑧 versus 𝑝.

(a)

(b)

(c)

(d)

Figure 3 

Crisis. (a) Chaotic attractor when 𝑝=0.88, (b) phase diagram of periodic solution when𝑝=0.9, (c) chaotic attractor when 𝑝=1, and (d) phase diagram of periodic solution when𝑝=1.1.

(a)

(b)

(c)

(d)

Figure 4 

Period-halving bifurcation. (a) Chaotic attractor when 𝑝=1.3, (b) phase diagram of a 4T-periodic solution when 𝑝=1.55, (c) phase diagram of a 2T-periodic solution when𝑝=2, and (d) a period-T attractor when𝑝=2.5.

3.2. Strange Attractors and Power Spectra

Now the commonly used method called power spectra will be used to study the qualitative nature of strange attractors [21]. By calculating the largest Lyapunov exponent for strange attractor (a) (Figure 5(a)), this value is determined to be 0.18219. Obviously, the strange attractor is a chaotic attractor. Furthermore, the spectrum of chaotic attractor (a) (Figure 5(b)) is composed of dense wave bands and sparse wave peaks with low wave troughs. These results agree with the observation that the chaotic attractor comes into being because some cycles lose weak stability.

(a)

(b)

Figure 5 

Strange attractor and power spectrum. (a) Strange attractor when 𝑝=1.25, (b) power spectrum of attractor (a).

3.3. The Largest Lyapunov Exponent

Convincing evidence for deterministic chaos has come from several recent experiments [22–24]. From these results, it is clear that chaos plays a very important role in these studies, and therefore detecting and exploring chaos is very important [25–31]. Therefore, the largest Lyapunov exponent is considered to be the most useful diagnostic tool for chaotic systems. The largest Lyapunov exponent 𝜆 must be positive for a chaotic attractor, otherwise, if 𝜆 is negative, the system will enter a stable state or become a periodic attractor. The largest Lyapunov exponent for system (1.1) was calculated for various values of 𝑝, and Figure 6 shows the results for 𝑝 from 0 to 2.

Figure 6 

The largest Lyapunov exponents (LLE) for system (1.1) with 𝑝 varying between 0 and 2.

4. Conclusions

In this paper, the effects of impulsive perturbations on a fish algae model have been investigated. It has been proved by means of the Floquet theory of impulsive differential equations and small-amplitude perturbation techniques that system (1.1) has a stable algae-eradication periodic solution. Furthermore, the conditions for system (1.1) to be permanent have been determined using the comparison theorem. Typical bifurcation diagrams have been analyzed in detail, revealing that the system exhibits very rich dynamics. From Theorem 2.8, it is clear that the cyanobacteria and green algae-eradication periodic solution (0,0,𝑧^∗(𝑡)) is locally asymptotically stable if 𝑝>𝑝_max=8.15, but that the cyanobacteria population 𝑥 and green algae population 𝑦 rapidly decrease to zero. When the value of 𝑝>2.75, green algae population 𝑦 decreases to zero, but the cyanobacteria population 𝑥 and fish population 𝑧 can coexist in a cycle; when 2.25<𝑝<2.75, the cyanobacteria population, green algae population and fish population can coexist in a cycle. When 0<𝑝<2.25, system (1.1) undergoes complex dynamics. From the above analysis, it is apparent that the numbers of both algae species can be controlled effectively using an impulsive control strategy.

Acknowledgments

The authors would like to thank the editor and the anonymous referees for their valuable comments and suggestions on this paper. This work was supported by the National Natural Science Foundation of China (Grant no. 31170338 and Grant no. 30970305).