毒素影响下捕食种群具有脉冲投放的捕食系统的动力学研究

Dynamical Research of a Predator-Prey System of Predator Population with Impulse Release Rate under Toxin Influence

摘要: 本文建立毒素影响下捕食者具有脉冲投放的捕食–食饵模型，并运用Bendixson-Dulac定理、Poincaré-Bendixson定理等方法分析了无脉冲系统的全局动力学；接着证明了不同条件下阶一周期解的存在性与唯一性；其次利用类Poincaré准则，得到了阶一周期解稳定的条件；最后，分别选取不同的脉冲系数、阈值以及毒素系数，讨论不同参数对鱼类的可持续发展的影响，最后数值仿真结果验证了结论的合理性。

Abstract: In this paper, we establish a predator-prey model with a constant feeding rate of predators under the influence of toxins. The existence of the equilibrium point and the stability of the positive equilibrium point of the system are studied, combined with the theory of pulsed differential equa-tions. The pulse control of the pulseless system is carried out according to the actual situation, the persistence of the predation system is discussed, and the uniqueness of the order 1 periodic solu-tion of the system is proved by the following function method and interval set theorem, and finally, the conclusion is verified by numerical simulation.

1. 引言

富营养化是指用氮和磷等营养物质过度污染海水，造成水生生态系统中物种分布的不平衡。它造成水生生态系统中物种分布的不平衡，影响系统内物质和能量的流动，破坏整个海洋生态系统。富营养化最常见的后果是形成藻类大量繁殖，大多数淡水藻类群包含形成有害藻华的物种，包括真核生物藻类，如绿藻、甲烷、隐藻和金藻，以及原核生物蓝藻。这些藻类有的产生臭味物质，有的则产生毒素，对鱼类健康构成最大风险。因此许多学者研究了在水环境中的捕食系统，得到了许多重要结论 [1] - [7] 。

本文我们考虑富营养化影响下鱼类(捕食者)受浮游植物(食饵)数量影响的捕食关系，建立了具有状态反馈脉冲收获的数学模型，通过分析其动力学性质研究这一生态系统的稳定性。

2. 建立模型

本文模拟了在富营养化下捕食者鱼类(捕食者)受浮游植物(食饵)影响而减少和Holling Ι功能反应的捕食系统模型，来反映鱼类的养殖与捕捞整个过程。模型如下

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(r-x\right)-\lambda xy\\ \frac{\text{d}y}{\text{d}t}=hxy-f{x}^{2}y-ky-b{y}^2\\ \Delta x=-\alpha x\\ \Delta y=-\beta y+u\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x\ne m\\ \begin{array}{l}\\ \end{array}\}& x=m\end{array}$ . (1)

其中x与y分别为浮游植物与鱼类在t时刻的数量，r表示浮游植物的生长率(出生率减去死亡率)，并且受到密度制约， $\lambda$ 为鱼类对浮游植物的消耗率，h为鱼类的捕食率，f表示为浮游植物过多影响下鱼类的死亡率，k表示鱼类的自然死亡率，b表示鱼类的种间竞争率。当浮游植物的数量达到m时，捕捞浮游植物，为了节省打捞成本和考虑生态平衡，同时捕捞鱼类，为了鱼类的可持续发展，还要投放一定数量的鱼苗进行补充来保证产量。 $0\le \alpha \le 1$ 为收获浮游植物的比率， $0\le \beta \le 1$ 为收获鱼类的比例，u为鱼苗的投放量。考虑种群实际情况，以上参数均为正常数且两种群具有正初始值情况，研究将在 $R_2^{+}=\left\{\left(x,y\right)|x\ge 0,y\ge 0\right\}$ 内进行。

3. 无脉冲系统分析

当不考虑脉冲时，系统(1)退化为

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(r-x\right)-\lambda xy\\ \frac{\text{d}y}{\text{d}t}=hxy-f{x}^{2}y-ky-b{y}^2\end{array}$ . (2)

3.1. 平衡点分析

根据

$\frac{\text{d}x}{\text{d}t}=\text{0},\frac{\text{d}y}{\text{d}t}=\text{0}$ . (3)

得垂直等倾线

$x=0,y_1=\frac{r-x}{\lambda }$ . (4)

与水平等倾线

$y=0,y_2=\frac{1}{b}\left(hx-f{x}^2-k\right)$ . (5)

容易计算出系统恒有平衡点 $\left(0,0\right),\left(r,0\right)$ 。

水平等倾线 $y_2=\frac{1}{b}\left(hx-f{x}^2-k\right)$ 与x轴相交情况分为三种：

当 $\Delta =\frac{1}{b}\left(h^2-4fk\right)>0$ ，等倾线 $y_2$ 与x轴交点有两个，并且在 $x_0=\frac{h}{2f}$ 时，水平等倾线 $y_2=\frac{1}{b}\left(hx-f{x}^2-k\right)$ 取得极大值；

当 $\Delta =\frac{1}{b}\left(h^2-4fk\right)=0$ ，等倾线 $y_2$ 与x轴交点只有一个；

当 $\Delta =\frac{1}{b}\left(h^2-4fk\right)<0$ ，等倾线 $y_2$ 与x轴没有交点。

由于第2、3种情况，系统不存在正平衡点，我们不考虑这两种情况。

垂直等倾线 $y_1=\frac{r-x}{\lambda }$ 与水平等倾线 $y_2=\frac{1}{b}\left(hx-f{x}^2-k\right)$ 相交情况分为三种：

1) 当 $H_1={\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)>0$ ，两线交于两点 $A_1\left(x_1^{*},y_1^{*}\right)$ 与 $A_2\left(x_2^{*},y_2^{*}\right)$ 。通过计算可得：

$x_1^{*}=\frac{\lambda h+b+\sqrt{{\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)}}{2\lambda f},y_1^{*}=\frac{r}{\lambda }-\frac{\lambda h+b+\sqrt{{\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)}}{2{\lambda }^{2}f}$ . (6)

$x_2^{*}=\frac{\lambda h+b-\sqrt{{\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)}}{2\lambda f},y_2^{*}=\frac{r}{\lambda }-\frac{\lambda h+b-\sqrt{{\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)}}{2{\lambda }^{2}f}$ . (7)

两点的分布情况将会分为以下三种情况

当 $y_1^{*}<y_2^{*}\le 0$ 时，系统不存在正平衡点。

当 $y_1^{*}\le 0,y_2^{*}>0$ 时，系统存在一个正平衡点 $A_2$ 。

当 $y_2^{*}>y_1^{*}>0$ 时，系统存在两个正平衡点 $A_1$ 与 $A_2$ 。

2) 当 $H_2={\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)=0$ ，两线相交于一点，系统有唯一的正平衡点 $A_3\left(x_3^{*},y_3^{*}\right)$ ，通过计算可得：

$x_1^{*}=\frac{\lambda h+b}{2\lambda f},y_1^{*}=\frac{r}{\lambda }-\frac{\lambda h+b}{2{\lambda }^{2}f}$ ， (8)

3) 当 $H_3={\left(\lambda h+b\right)}^2-4\lambda f\left(\lambda k+br\right)<0$ ，系统不存在正平衡点。

根据实际意义，下面重点讨论系统存在正平衡点的情况(如图1所示)。

3.2. 平衡点的稳定性

系统(2)的Jocabi矩阵

$J=\left[\begin{array}{cc}r-2x-ky& -\lambda x\\ hy-2fxy& hx-f{x}^2-k-2by\end{array}\right]$ ， (9)

对于平衡点 $\left(\text{0},0\right)$ 的特征根方程

$|\begin{array}{cc}r-\bar{\lambda }& 0\\ 0& -k-\bar{\lambda }\end{array}|=0$ ， (10)

可得 ${\bar{\lambda }}_1{\bar{\lambda }}_2<0$ ，则 $\left(\text{0},0\right)$ 是鞍点。

对于平衡点 $\left(r,0\right)$ 的特征根方程

$|\begin{array}{cc}-r-\bar{\lambda }& -\lambda r\\ 0& hr-f{r}^2-k-\bar{\lambda }\end{array}|=0$ ， (11)

因为

$hr-f{x}^2-k<0$ . (12)

${\bar{\lambda }}_1{\bar{\lambda }}_2=\det \left(J\right)>0$ ， ${\bar{\lambda }}_1+{\bar{\lambda }}_2=tr\left(J\right)<0$ ， (13)

$T^2-4D={\left[r-\left(hr-f{r}^2-k\right)\right]}^2+4r\left(hr-f{r}^2-k\right)={\left[r+\left(hr-f{r}^2-k\right)\right]}^2>0$ ， (14)

所以平衡点 $\left(r,0\right)$ 是稳定的结点。

对于正平衡点 $A_1\left(x_1^{*},y_1^{*}\right)$ 的特征根方程

$|\begin{array}{cc}-x_1^{*}-\bar{\lambda }& -\lambda x_1^{*}\\ h{y}_1^{*}-2f{x}_1^{*}{y}_1^{*}& -b{y}_1^{*}-\bar{\lambda }\end{array}|=0$ ， (15)

${\bar{\lambda }}_1{\bar{\lambda }}_2=\det \left(J\right)=b{x}_1^{*}{y}_1^{*}+\lambda x_1^{*}\left(h{y}_1^{*}-2f{x}_1^{*}{y}_1^{*}\right)=x_1^{*}{y}_1^{*}\left(b+\lambda h-2\lambda f{x}_1^{*}\right)<0$ . (16)

所以正平衡点 $A_1\left(x_1^{*},y_1^{*}\right)$ 是鞍点。

对于正平衡点 $A_2\left(x_2^{*},y_2^{*}\right)$ 的特征根方程

$|\begin{array}{cc}-x_2^{*}-\bar{\lambda }& -\lambda x_2^{*}\\ h{y}_2^{*}-2f{x}_2^{*}{y}_2^{*}& -b{y}_2^{*}-\bar{\lambda }\end{array}|=0$ ， (17)

${\bar{\lambda }}_1{\bar{\lambda }}_2=\det \left(J\right)=b{x}_2^{*}{y}_2^{*}+\lambda x_2^{*}\left(h{y}_2^{*}-2f{x}_2^{*}{y}_2^{*}\right)=x_2^{*}{y}_2^{*}\left(b+\lambda h-2\lambda f{x}_2^{*}\right)>0$ ， (18)

${\bar{\lambda }}_1+{\bar{\lambda }}_2=tr\left(J\right)=-x_2^{*}-b{y}_2^{*}<0$ . (19)

所以正平衡点 $A_2\left(x_2^{*},y_2^{*}\right)$ 是稳定的结点或交点。

系统(2)的轨线图如图1所示。

3.3. 极限环的不存在性

定理1：系统(2)在第一象限内不存在极限环。

Figure 1. The trajectory diagram of the system (2)

图1. 系统(2)的轨线图

证：令

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=P\left(x,y\right)\\ \frac{\text{d}y}{\text{d}t}=Q\left(x,y\right)\end{array}$ . (20)

取函数 $B\left(x,y\right)=\frac{1}{xy}$ ，则有

$\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}=-\frac{1}{y}-\frac{b}{x}<0$ ， (21)

根据Bendixson-Dulac定理可知，系统(2)在第一象限不存在极限环。

3.4. 无脉冲系统的一致有界性

定理2：系统(2)在第一象限内一致有界。

证：根据

$\frac{\text{d}x}{\text{d}t}=0,\frac{\text{d}y}{\text{d}t}=0$ ， (22)

得垂直等倾线

$x=0,y=\frac{r-x}{\lambda }$ ， (23)

与水平等倾线

$y=0,y_2=\frac{1}{b}\left(hx-f{x}^2-k\right)$ . (24)

根据环域定理，构造外境界线，令 $y_b>\max \{y_2,\frac{r}{\lambda }\}$ ， $y=y_B$ 与y轴相交于B。过点C作x轴垂线交x轴于点C，与 $y=y_B$ 相交于点D。

对于OB，OC，由于y轴与x轴是轨线，没有任何的轨线与它相交。

对于BD，当 $0<x<r$ 时，有 ${\frac{\text{d}y}{\text{d}t}|}_{y=y_B}=y_b\left(hx-f{x}^2-k-b{y}_B\right)<0$ ，故当轨线与直线BD相遇时，从它的上面穿向下面。

对于DC，由于 ${\frac{\text{d}x}{\text{d}t}|}_{x=r}=-\lambda ry<0$ ，故轨线与直线AB相遇时，均从它的右面穿向左面。

这样就围成了区域G的外境界线(如图2所示)。且G内除E₂外无其他奇点，奇点E₁是鞍点，所以从 $R_{+}^2$ 出发的解一致有界。

Figure 2. Boundedness of the system (2)

图2. 系统(2)的有界性

4. 阶1周期解的存在唯一性与稳定性

根据上述对该无脉冲系统的分析，当参数条件H₁成立时，系统存在两个正平衡点。从无脉冲轨线图可以看出，当浮游植物不断增加时，鱼类最终会走向灭亡，所以本文将讨论脉冲系统(1)周期解的存在性。

4.1. 阶1周期解的存在唯一性

稳定的正平衡点坐标为 $E\left(x^{*},y^{*}\right)$ ，当 $m=x^{*}$ 时，设垂直等倾线 $\frac{\text{d}x}{\text{d}t}=0$ 与相集N相交于点A，点A坐

标为 $\left(\left(1-\alpha \right)m,y_A\right)$ ，设轨迹 $L\left(A,t\right)$ 与脉冲集M交于点A₁，经过脉冲效应映射到相集N上的点 $A_1^{*}$ ，下面我们分情况讨论点 $A_1^{*}$ 的位置情况，并讨论此时阶一周期解的存在性。

定理3：当系统(2)脉冲系数处于以下几种情况时，系统(1)均存在一个唯一的阶1周期解。

证明如下：

情况1 当 $u=u_1$ 时

点A₁经过脉冲效应映射到相集N上的点 $A_1^{*}$ ，则坐标为 $A_1^{*}\left(\left(1-\alpha \right)m,\left(1-\beta \right)y_{A_1}+u\right)$ ，当 $u=u_1$ 时，使 $y_{A_1^{*}}=\left(1-\beta \right)y_A{}_{{}_1}+u_1=y_A$ 成立，可知点 $A_1^{*}$ 此时与点A重合，点A的后继函数为 $F\left(A\right)=y_A-y_{A_1^{*}}=0$ ，此时系统存在一个阶1周期解得证(如图3所示)。

情况2 当 $u>u_1$ 时

点A₁经过脉冲效应映射到相集N上的点 $A_1^{*}$ ，则坐标为 $A_1^{*}\left(\left(1-\alpha \right)m,\left(1-\beta \right)y_{A_1}+u\right)$ ，当 $u>u_1$ 时， $y_{A_1^{*}}=\left(1-\beta \right)y_{A_1}+u>y_A$ 成立，可知点 $A_1^{*}$ 必在点A上方。设过点 $A_1^{*}$ 的轨线与脉冲集M相交于点B，B点

Figure 3. $u=u_1^{*}$ , Schematic diagram of the existence of order-1 periodic solution

图3. $u=u_1^{*}$ 时，阶一周期解的存在性示意图

坐标为 $\left(m,y_B\right)$ ，点B经过脉冲映射与相集N相交于点 $B^{+}$ ，根据轨线方向和轨线不能相交的性质，我们可以得到点B在点A₁的下方，即 $y_B<y_{A_1}$ 。从而我们可以得到 $y_{B^{+}}<y_{A_1^{*}}$ ，所以点 $A_1^{*}$ 的后继函数为 $F\left(A_1^{*}\right)=y_{B^{+}}-y_{A_1^{*}}<0$ 。我们在点A上方选取另一个点C，使 $\epsilon >0$ 且足够小，得到点C的坐标为 $C\left(\left(1-\alpha \right)m,y_A+\epsilon \right)$ ，则过点C的轨线与脉冲集M相交与点C₁，点C₁坐标为 $\left(m,y_{C_1}\right)$ ，根据轨线方向和轨线不能相交的性质，我们可以得到点C₁在点A₁的下方，由于解对初值和时间的连续依赖性，点C₁非常接近点A₁，即 $y_{C_1}<y_{A_1}$ 。点C₁经过脉冲映射与相集N相交于点C₂，从而可以得到点C₂常接近点 $A_1^{*}$ ，并且在点 $A_1^{*}$ 的下方，即 $y_C<y_{C_2}<y_{A_1^{*}}$ 。所以点C的后继函数为 $F\left(C\right)=y_{C_2}-y_C>0$ 。则由引理1-1，引理1-2可知，在点C与点 $A_1^{*}$ 之间必存在一点，使这个点的后继函数为零，此时系统存在一个唯一阶1周期解得证(如图4所示)。

Figure 4. $u>u_1^{*}$ , Schematic diagram of the existence of order-1 periodic solution

图4. $u>u_1^{*}$ 时，阶一周期解的存在性示意图

接下来，我们考虑此时系统(1)的阶1周期解的唯一性。已知系统存在一个阶1周期解，它的初始点在相集上的点C和点 $A_1^{*}$ 之间。则我们在线段 $\overline{C{A}_1^{*}}$ 上随意选取两点E，G，坐标分别为 $E\left(\left(1-\alpha \right)m,x_E\right)$ ， $G\left(\left(1-\alpha \right)m,x_G\right)$ ，并使两点满足 $y_c<y_G<y_E<y_{A_1^{*}}$ 。那么从点E出发的轨线到达脉冲集M于点E₁，点G出发的轨线到达脉冲集M于点G₁，根据轨线方向和轨线不能相交的性质，我们可以得到点E₁必定在点G₁的下方，那么两点E，G在经过脉冲映射到达相集N，与相集N分别相交于点E₂，点G₂，可以得到点E₂必定在点G₂的下方，即 $y_{E_2}<y_{G_2}$ ，于是点E与点G的后继函数满足： $F\left(E\right)-F\left(G\right)=y_{E_2}-y_E-y_{G_2}+y_G<0$ ，这就说明了在线段 $\overline{C{A}_1^{*}}$ 上的后继函数是单调递减的，所以线段 $\overline{C{A}_1^{*}}$ 上仅存在一个唯一的阶1周期解，此时系统存在一个唯一的阶1周期解得证(如图5所示)。

Figure 5. Uniqueness of order-1 periodic solution

图5. 阶一周期解的唯一性

情况3 当 $u<u_1$ 时

点 $A_1$ 经过脉冲效应映射到相集N上的点 $A_1^{*}$ ，则坐标为 $A_1^{*}\left(\left(1-\alpha \right)m,\left(1-\beta \right)y_{A_1}+u\right)$ ，当 $u<u_1$ 时， $y_{A_1^{*}}=\left(1-\beta \right)y_{A_1}+u<y_A$ 成立，可知点 $A_1^{*}$ 必在点A下方。所以点 $A_1^{*}$ 的后继函数为 $F\left(A_1^{*}\right)=y_A-y_{A_1^{*}}<0$ 。我们在点A下方选取另一个点D，使 $\epsilon >0$ 且足够小，得到点D的坐标为 $D\left(\left(1-\alpha \right)m,\epsilon \right)$ ，则过点D的轨线与脉冲集M相交于点 $D_1$ ，点 $D_1$ 坐标为 $\left(m,y_{D_1}\right)$ 。点 $D_1$ 经过脉冲映射与相集N相交于点 $D_2$ ，点 $D_2$ 坐标为 $\left(\left(1-\alpha \right)m,y_{D_2}\right)$ ，由于 $\epsilon >0$ 且足够小，从而可以得到 $y_D<y_{D_2}$ 。所以点D的后继函数为 $F\left(D\right)=y_{D_2}-y_D>0$ 。则由引理1-1，引理1-2可知，在点D与点A之间必存在一点，使这个点的后继函数为零，此时系统存在一个阶1周期解得证(如图6所示)。

接下来，我们考虑此时系统(1)的阶1周期解的唯一性。由于正平衡点是渐近稳定的，且上述已经证明系统(1)的阶一周期解的初始点在点D与点A之间，那么仅考虑初始点落在直线段 $\overline{AD}$ 上，可以看出在线段 $\overline{AD}$ 上的点通过轨迹运动到线段 $\overline{A_1{D}_1}$ ，再通过脉冲作用被映射到线段 $\overline{A_1^{*}{D}_2}$ 。显然， $y_{D_1}<y_{A_1}$ ， $y_D<y_{D_2}<y_{A_1^{*}}<y_{A_1}$ ，根据系统(1)的向量场，重复上述过程，可以得到两个数列。根据这两列序列，因此有 $\left(x_{A_1^{*}},x_{D_2}\right)\subset \left(x_A,x_D\right)$ ， $d\left(A_1^{*},D_2\right)<d\left(A,D\right)$ 。根据区间套定理，，这意味着在点D与点A之间存在一个唯一的点，使这个点的后继函数为零。因此，系统(1)的阶一周期解具有唯一性。

定理4：当系统(2)的平衡点与脉冲集M的位置出于以下几种情况时，系统(1)不存在阶1周期解。

设此时系统存在两个正平衡点，坐标分别为 $E_1\left(x_1,y_1\right),E_2\left(x_2,y_2\right)$ ，且 $x_1<x_2$

$\left(1-\alpha \right)m<x_1<m<x_2$ ；

$x_1<\left(1-\alpha \right)m<m<x_2$ ；

$\left(1-\alpha \right)m<x_1<x_2<m$ ；

$x_1<\left(1-\alpha \right)m<x_2<m$ .

Figure 6. $u<u_1^{*}$ , Schematic diagram of the existence of order-1 periodic solution

图6. $u<u_1^{*}$ 时，阶一周期解的存在性示意图

证：类似于定理6的证明方法可知，根据系统(1)的向量场及轨线走势可知，当 $\left(1-\alpha \right)m<x_1<m<x_2$ 时，于是从相集上任意一点出发的轨线最多经过有限次脉冲效应趋于平衡点 $E_1\left(x_1,y_1\right)$ ，例如示意图3~7中的轨线的轨迹，此时系统(1)无周期解，当脉冲集与相集处于1)，2)，3)等不同位置时，用相同的思路可以证明剩下三种条件下，该系统均不存在阶一周期解(如图7(b)、图7(c)、图7(d)所示)。

(a) $\left(1-\alpha \right)m<x_1<m<x_2$ (b) $x_1<\left(1-\alpha \right)m<m<x_2$(c) $\left(1-\alpha \right)m<x_1<x_2<m$ (d) $x_1<\left(1-\alpha \right)m<x_2<m$

Figure 7. Schematic diagram of nonexistence of periodic solution of model (2); When there is two equilibrium point in model (1)

图7. 当系统(1)有两个平衡点时，模型(2)周期解的不存在性示意图

4.2. 阶1周期解的稳定性

定理5：当 $\left|\ln \frac{y_1}{y_2}\frac{r-\left(1-\alpha \right)m-\lambda \left(\left(1-\beta \right)y_1+u\right)}{r-\left(1-\alpha \right)m-\lambda y_0}\right|<1$ 时，系统(1)过C点的阶1周期解轨道渐近稳定。

证：设 $\Gamma \left(x=\varphi \left(t\right),y=\phi \left(t\right)\right)$ 是系统(1)的一个阶1周期解，有 $C_0=\left(\varphi \left(0\right),\phi \left(0\right)\right)$ ， $C_1=\left(\varphi \left(T\right),\phi \left(T\right)\right)$ ， $C_1^{+}=\left(\varphi \left(T^{+}\right),\phi \left(T^{+}\right)\right)$ ，则设 $C_1\left(m,y_1\right)$ 为脉冲点， $C^{+}\left(\left(1-\alpha \right)m,y_2\right)$ 为 $C_0$ 点的后继点。令

$P\left(x,y\right)=x\left(r-x\right)-\lambda xy$ ， $Q\left(x,y\right)=hxy-f{x}^{2}y-ky-b{y}^2$ ， (25)

$\alpha \left(x,y\right)=-\alpha x$ ， $\beta \left(x,y\right)=-\beta y+u$ ， $\Phi \left(x,y\right)=x-m$ ， (26)

计算得

$\frac{\partial P}{\partial x}=r-2x-\lambda y$ ， $\frac{\partial Q}{\partial y}=hx-f{x}^{2}y-k-2by$ ， $\frac{\partial \alpha }{\partial x}=-\alpha$ ， (27)

$\frac{\partial \beta }{\partial y}=-\beta$ ， $\frac{\partial \alpha }{\partial y}=\frac{\partial \beta }{\partial x}=\frac{\partial \Phi }{\partial y}=0$ ， $\frac{\partial \Phi }{\partial x}=1$ . (28)

$\begin{array}{c}{\Delta }_1=\frac{P_{+}\left(\frac{\partial \beta }{\partial y}\frac{\partial \Phi }{\partial x}-\frac{\partial \beta }{\partial x}\frac{\partial \Phi }{\partial y}+\frac{\partial \Phi }{\partial x}\right)+Q_{+}\left(\frac{\partial \alpha }{\partial x}\frac{\partial \Phi }{\partial y}-\frac{\partial \alpha }{\partial y}\frac{\partial \Phi }{\partial x}+\frac{\partial \Phi }{\partial y}\right)}{P\left(\frac{\partial \Phi }{\partial x}\right)+Q\left(\frac{\partial \Phi }{\partial y}\right)}\\ =\frac{\left(1-\beta \right)P\left(C^{+}\right)}{P\left(C_0\right)}\\ =\frac{\left(1-\beta \right)\left[r-\left(1-\alpha \right)m-\lambda \left(\left(1-\beta \right)y_1+u\right)\right]}{r-\left(1-\alpha \right)m-\lambda y_0}.\end{array}$ (29)

$\begin{array}{c}{\displaystyle {\int }_0^T\left(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}\right)\text{d}t}={\displaystyle {\int }_0^T\left(r-2x-ky+hx-f{x}^2-k-2by\right)\text{d}t}\\ ={\displaystyle {\int }_{\left(1-\alpha \right)m}^m\frac{1}{x}\text{d}x}+{\displaystyle {\int }_{y_2}^{y_1}\frac{1}{y}\text{d}y}+{\displaystyle {\int }_0^T\left(-\frac{u}{y}-x\right)\text{d}t}\\ =\ln \frac{1}{1-\alpha }+\ln \frac{y_1}{y_2}+{\displaystyle {\int }_0^T\left(-by-x\right)\text{d}t}.\end{array}$ (30)

显然 ${\displaystyle {\int }_0^T\left(-by-x\right)\text{d}t}<0$ ，进而我们可以得到：

$\begin{array}{c}{\mu }_2={\Delta }_1\exp \left\{{\displaystyle {\int }_0^T\left(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}\right)\text{d}t}\right\}\\ =\ln \frac{y_1}{y_2}\frac{r-\left(1-\alpha \right)m-\lambda \left(\left(1-\beta \right)y_1+u\right)}{r-\left(1-\alpha \right)m-\lambda y_0}\exp {\displaystyle {\int }_0^T\left(-by-x\right)\text{d}t},\end{array}$ (31)

当 $\left|\ln \frac{y_1}{y_2}\frac{r-\left(1-\alpha \right)m-\lambda \left(\left(1-\beta \right)y_1+u\right)}{r-\left(1-\alpha \right)m-\lambda y_0}\right|<1$ 时， $\left|{\mu }_2\right|<1$ 。由类Poincaré准则可知，系统(1)过C点的阶1周期解是轨道渐近稳定的。

5. 数值模拟和结论

基于前面对该脉冲控制系统的理论分析，通过数值模拟来验证污染环境下收获不同数量的食饵(浮游植物)对该系统的影响。研究模型

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(0.9-x\right)-0.26xy\\ \frac{\text{d}y}{\text{d}t}=0.5xy-0.1x^{2}y-0.4y-0.5y^2\\ \Delta x=-0.1x\\ \Delta y=-0.5y+0.1\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x<0.9\\ \begin{array}{l}\\ \end{array}\}& x=0.9\end{array}$ (32)

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(0.9-x\right)-0.26xy\\ \frac{\text{d}y}{\text{d}t}=0.5xy-0.1x^{2}y-0.4y-0.5y^2\\ \Delta x=-0.1x\\ \Delta y=-0.9y+0.1\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x<0.9\\ \begin{array}{l}\\ \end{array}\}& x=0.9\end{array}$ (33)

当我们保留其他参数相同时，使浮游植物的脉冲收获系数不同，模型(31)的捕食者与食饵的时间序列图如图8、图9所示。模型(32)的捕食者与食饵的时间序列图如图10、图11所示。

可以说明当富营养化时，食饵(浮游植物)占据大量的生存空间，影响捕食者(鱼类)的生存。当捕捞较少的浮游植物，浮游植物的种群数量一直维持在较高水平，使食饵(浮游植物)占据了更多的生存空间，使鱼类的持续生存受到威胁。但当我们加大浮游植物的收获量时，给了鱼类可发展的空间，减缓了鱼类数量的减少。我们改变浮游植物的收获系数时，也很难控制浮游植物持续迅速的生长，只能减缓这种情况，这就说明，食饵的脉冲收获系数发生改变会影响系统平衡处的状态值，并且收获系数适当大，有利于鱼类的可持续发展。

Figure 8. $\alpha =0.1$ , Time series diagram of x

图8. $\alpha =0.1$ , x的时间序列图

为了进一步研究经济阈值m对捕食者(鱼类)和食饵(浮游植物)可持续发展的影响，研究模型

Figure 9. $\alpha =0.1$ , Time series diagram of y

图9. $\alpha =0.1$ , y的时间序列图

Figure 10. $\alpha =0.9$ , Time series diagram of x

图10. $\alpha =0.9$ , x的时间序列图

Figure 11. $\alpha =0.9$ , Time series diagram of y

图11. $\alpha =0.9$ , y的时间序列图

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(0.9-x\right)-0.26xy\\ \frac{\text{d}y}{\text{d}t}=0.5xy-0.1x^{2}y-0.4y-0.5y^2\\ \Delta x=-0.1x\\ \Delta y=-0.5y+0.1\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x<1.5\\ \begin{array}{l}\\ \end{array}\}& x=1.5\end{array}$ . (34)

对比模型(32)，当我们保留其他参数相同时，使经济阈值m变大，模型(34)的解曲线如图12所示。当 $x_1<m$ 时，该系统(1)无周期解，也验证了定理4。当阈值m过大的时候，浮游植物的数量到达一定程度时，减少了鱼类的生存空间。大量藻类覆盖水面恶化了水中的通风、导致无法得到光照，光合作用受阻，使水体处于严重缺氧状态，致使鱼类大量死亡。

Figure 12. The solution curve of model (34)

图12. 模型(34)的解曲线

下面我们继续研究毒素系数f对捕食者(鱼类)和食饵(浮游植物)可持续发展的影响，研究模型

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(0.9-x\right)-0.26xy\\ \frac{\text{d}y}{\text{d}t}=0.5xy-0.3x^{2}y-0.4y-0.5y^2\\ \Delta x=-0.1x\\ \Delta y=-0.5y+0.1\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x<0.9\\ \begin{array}{l}\\ \end{array}\}& x=0.9\end{array}$ . (35)

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(0.9-x\right)-0.26xy\\ \frac{\text{d}y}{\text{d}t}=0.5xy-0.6x^{2}y-0.4y-0.5y^2\\ \Delta x=-0.1x\\ \Delta y=-0.5y+0.1\end{array}\begin{array}{cc}\begin{array}{c}\\ \\ \end{array}\}& x<0.9\\ \begin{array}{l}\\ \end{array}\}& x=0.9\end{array}$ . (36)

当我们保留其他参数相同时，使毒素系数f不同，模型(35)的捕食者与食饵的时间序列图如图13、图14所示。模型(36)的捕食者与食饵的时间序列图如图15、图16所示。

Figure 13. $f=0.3$ , Time series diagram of x

图13. $f=0.3$ ，x的时间序列图

Figure14. $f=0.3$ , Time series diagram of y

图14. $f=0.3$ ，y的时间序列图

Figure 15. $f=0.6$ , Time series diagram of x

图15. $f=0.6$ ，x的时间序列图

Figure 16. $f=0.6$ , Time series diagram of y

图16. $f=0.6$ ，y的时间序列图

当毒素系数增大时，显然看到，捕食者(鱼类)和食饵(浮游植物)的收获周期明显缩短，可以说明，毒素系数发生改变会影响系统平衡处的状态值，并且毒素系数越大，其影响程度越大。所以当富营养化情况越严重时，养殖鱼需要耗费大量的人力物力，不利渔业的发展，这就需要我们根据毒素对系统的影响来制定更加有效的脉冲收获策略。

NOTES

^*通讯作者。