DOI:10.12677/AAM.2022.118636,PDF,HTML,XML,被引量下载: 305浏览: 538科研立项经费支持
作者:张 琰^*,张 蒙：北京建筑大学理学院，北京
关键词:毒素；常数投放；阶一周期解；类Poincaré准则；数值模拟；Toxin；Constant Release；Order-1 Periodic Solution；The Analogue of Poincaré’s Criterion；Numerical Simulation

文章引用：张琰, 张蒙. 毒素影响下捕食种群具有常数投放率的捕食–食饵模型的动力学研究[J]. 应用数学进展, 2022, 11(8): 6044-6054. https://doi.org/10.12677/AAM.2022.118636

1. 引言

随着现代化工业的迅速发展，在促进国民经济发展的同时，造成的环境污染问题也是不容忽视的，这严重影响了空气、土壤以及水源。环境污染会对生态系统造成直接的影响和破坏。同时环境污染问题已是人类面临的重大课题之一。

近些年来，许多学者发现脉冲微分方程的理论和方法在生命科学领域中有很好的应用前景。例如陈 [1]、Liu [2]、黄立状 [3] 等考虑了在害虫防治时，大量使用杀虫剂会使环境恶化，因此采取脉冲投放天敌的方法抑制害虫生长。毒素对生物种群的影响已有大量文献研究 [4] - [10]。SKB [6] 等研究了毒素环境下的阶段结构渔业模型，利用Pontryagins极大值原理得到了最优收获策略。蓝桂杰 [7] 等研究了在污染环境下毒素脉冲输入和心理效应对随机捕食–食饵系统的影响，得到系统持续生存与灭绝的充分条件。张茜 [8] 等研究了污染物排放导致捕食者减少的捕食食饵Gompertz模型的动力学性质。

本文考虑水环境污染导致鱼类的减少，对捕食者(鱼类)进行常量投放，可以防止浮游植物(食饵)占据大量生存空间，同时利于鱼类的可持续发展。因此建立了在毒素影响下捕食者具有常数投放和状态反馈脉冲控制模型。

2. 建立模型

本文模拟了在毒素影响下捕食者具有常数投放和HollingΙ功能反应的捕食系统模型，来反映鱼类的养殖与捕捞整个过程。模型如下

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(r-x\right)-kxy\\ \frac{\text{d}y}{\text{d}t}=hxy-fy+u\\ \Delta x=-\alpha x\\ \Delta y=-\beta y\end{array}\begin{array}{cc}\begin{array}{l}\\ \end{array}\}& y<m\\ & \\ \begin{array}{l}\\ \end{array}\}& y=m\end{array}$.(1)

其中x与y分别为浮游植物与鱼类在t时刻的数量，r表示浮游植物的生长率(出生率减去死亡率)，k为鱼类对浮游植物的消耗率，h为鱼类的捕食率，f表示为毒素影响下鱼类的死亡率，u表示鱼类的投放率。当鱼类的数量达到m时，收获鱼类，为了避免浮游植物大量繁殖，也要捕捞一定量的浮游植物。 $0\le \alpha \le 1$为收获浮游植物的比率， $0\le \beta \le 1$为收获鱼类的比例。考虑种群实际情况，以上参数均为正常数且两种群具有正初始值情况研究将在 $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)-kxy\\ \frac{\text{d}y}{\text{d}t}=hxy-fy+u\end{array}$.(2)

3.1. 平衡点分析

求方程组

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\text{0}\\ \frac{\text{d}y}{\text{d}t}=\text{0}\end{array}$.(3)

当 $y=0$时，x无解。则在横坐标轴上不存在奇点。

当 $x=0$时， $y=\frac{u}{f}$，故存在奇点 $E_{\text{1}}\left(\text{0,}\frac{u}{f}\right)$。

当 $x>0,y>0$时，有 $r-x-ky=0$，将其代入式子 $hxy-fy+u=0$中，得

$-hk{y}^2+\left(hr-f\right)y+u=0$.(4)

由韦达定理可得

$y_1{y}^{*}=\frac{u}{-hk}<0$，(5)

所以必有一正根 $y^{*}$，由求根公式可得

$y^{*}=\frac{f-hr-\sqrt{{\left(hr-f\right)}^2+4ukh}}{-2hk}$，(6)

$x^{*}=r-k{y}^{*}=r+\frac{f-hr-\sqrt{{\left(hr-f\right)}^2+4ukh}}{2h}$.(7)

由于当 $x^{*}>0$时，有

$r+\frac{f-hr-\sqrt{{\left(hr-f\right)}^2+4ukh}}{2h}>0$，(8)

$2hr+f-hr>\sqrt{{\left(hr-f\right)}^2+4ukh}$，(9)

${\left(2hr+f-hr\right)}^{\text{2}}>{\left(hr-f\right)}^2+4ukh$，(10)

$4h^2{r}^2+4hr\left(f-hr\right)>4ukh$，(11)

可以得出

$fr-uk>0$.(12)

则在第一象限存在一个正平衡点 $E_2\left(x^{*},y^{*}\right)$。根据实际情况，下面重点讨论系统存在正平衡点的情况。

3.2. 平衡点的稳定性

系统(2)的Jocabi矩阵为

$J=\left[\begin{array}{cc}r-2x-ky& -kx\\ hy& hx-f\end{array}\right]$，(13)

对于 $E_{\text{1}}\left(\text{0},\frac{u}{f}\right)$的特征方程

$\left|\begin{array}{cc}r-\frac{ku}{f}-\bar{\lambda }& 0\\ \frac{hu}{f}& -f-\bar{\lambda }\end{array}\right|=0$.(14)

因为 $fr-ku>0$，所以 ${\bar{\lambda }}_1{\bar{\lambda }}_2<0$，则 $E_{\text{1}}\left(\text{0},\frac{u}{f}\right)$是鞍点。

通过计算，得到 $E_{\text{2}}\left(x^{*},y^{*}\right)$的特征方程为

$\left|\begin{array}{cc}-x^{*}-\bar{\lambda }& -k{x}^{*}\\ \frac{h\left(r-x^{*}\right)}{k}& h{x}^{*}-f-\bar{\lambda }\end{array}\right|=0$，(15)

${\bar{\lambda }}_1{\bar{\lambda }}_2=\det \left(J\right)=\left(f-h{x}^{*}\right)x^{*}+h\left(r-x^{*}\right)x^{*}=\left(\frac{u}{y^{*}}+hk{y}^{*}\right)x^{*}>0$，(16)

${\bar{\lambda }}_1+{\bar{\lambda }}_2=tr\left(J\right)=-x^{*}+h{x}^{*}-f$，(17)

因为 $hxy-fy+u=0$，则 $hx-f=\frac{-u}{y}$，有 $-x^{*}+h{x}^{*}-f<0$，即

${\bar{\lambda }}_1+{\bar{\lambda }}_2<0$.(18)

所以 $E_{\text{2}}\left(x^{*},y^{*}\right)$是稳定的结点或者焦点。

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

3.3. 全局稳定性

定理1：当 $\lambda =\frac{k}{h}$时，系统(2)在第一象限是全局渐近稳定。

证：建立一个Lyapunov函数

$V\left(x,y\right)=\left[\left(x-x^{*}\right)-x^{*}\ln \left(\frac{x}{x^{*}}\right)\right]+\lambda \left[\left(y-y^{*}\right)-y^{*}\ln \left(\frac{y}{y^{*}}\right)\right]$，(19)

$\frac{\text{d}V}{\text{d}t}=\frac{x-x^{*}}{x}\frac{\text{d}x}{\text{d}t}+\frac{y-y^{*}}{y}\frac{\text{d}y}{\text{d}t}=\left(x-x^{*}\right)\left(r-x-ky\right)+\lambda \left(y-y^{*}\right)\left(hx-f+\frac{u}{y}\right)$，(20)

考虑到平衡点方程

$\{\begin{array}{l}r-x^{*}-k{y}^{*}=0\\ h{x}^{*}{y}^{*}-f{y}^{*}+u=0\end{array}$，(21)

$\frac{\text{d}V}{\text{d}t}=\left(x-x^{*}\right)\left[\left(x^{*}-x\right)-k\left(y-y^{*}\right)\right]+\lambda \left(y-y^{*}\right)\left[h\left(x-x^{*}\right)+u\left(\frac{1}{y}-\frac{1}{y^{*}}\right)\right]$，(22)

当 $\lambda =\frac{k}{h}$时，显然 $\frac{\text{d}V}{\text{d}t}<0$，所以系统(2)在第一象限是全局渐近稳定的。

3.4. 极限环的存在性

$\{\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}$.(23)

令 $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{u}{x{y}^2}<0$，(24)

则无脉冲系统在第一象限不存在极限环。

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

根据

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

得垂直等倾线

$x=0,y=\frac{r-x}{k}$，(26)

与水平等倾线

$y=\frac{u}{f-hx}$.(27)

根据环域定理，构造外境界线，过点C作x轴平行线交水平等倾线于点A，过点A作AB垂直于x轴交x轴于点B。

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

对于OB，当 $y=0$时，有 ${\frac{\text{d}y}{\text{d}t}|}_{y=0}=u>0$，故当轨线与x轴相遇时，均从它的下面穿向上面。

对于AC，由于 $\text{0}<x<x_A$时，有 ${\frac{\text{d}y}{\text{d}t}|}_{y=y_A}=hxy-fy+u<0$，故轨线与直线AC相遇时，从上面穿向下面。

可知 $B\left(\frac{fr-uk}{hr},\frac{r}{k}\right)$，当 $\frac{fr-uk}{hr}>r$时，可知B点在R点右侧。

对于AB，由于 ${\frac{\text{d}x}{\text{d}t}|}_{x=\frac{fr-uk}{hr}}=x\left(r-x\right)-kxy<0$，故当轨线与直线AB相遇时，均从它的右面穿向左面。

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

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

本节利用文献 [11] 中的半连续动力系统的理论和后继函数法分析系统(1)阶一周期解的存在唯一性与稳定性。

定义1：系统为半连续动力系统，则

脉冲集为 $M=\left\{\left(x,y\right)\in R_2^{+}|y=m,x\ge 0\right\}$，

脉冲映射为 $\phi \left(x,y\right)\in M\to \left(\left(\left(1-\alpha \right)x,\left(1-\beta \right)y\right)\right)\in R_2^{+}$，

相集为 $N=\phi \left(M\right)=\left\{\left(x,y\right)\in R_2^{+}|y=\left(1-\beta \right)m,x\ge 0\right\}$。

定义2：假设系统有一条从点 $D\left(x_D,y_D\right)\in N$出发的轨线，到达 $M$后与其相交于点 $D\left(x_{D_1},h\right)\in M$，经过脉冲作用后，跳跃至点 $D_2\left(x_{D_2},y_{D_2}\right)\in N$，称 $D_2$为 $D$的后继点，且记 $D$的后继函数为 $f\left(D\right)=x_{D_2}-x_D$。

定义3：对于系统来说，如果存在一点 $D\in N$和 $t_0>0$，使得 $f\left(D,t_0\right)=D_1\in M$和 $\phi \left(f\left(D,t_0\right)\right)=D\in N$成立，那么 $f\left(D,t_0\right)$称为系统的阶一周期解。如图3所示。

引理1：对于系统来说，如果存在两点 $E\in N$和 $F\in N$，使得 $f\left(E\right)f\left(F\right)<0$，那么在点 $E$和点 $F$之间一定存在一点 $D\in N$，使得 $f\left(D\right)=0$。

4.1. 阶1周期解的存在性

当 $m=y^{*}$时， $\Delta y=y^{+}-y^{*}=-\beta y^{*},y=\left(1-\beta \right)y^{*},0<\beta <1$。相集N与垂直等倾线相交于点A，过A微分方程的轨线与脉冲集与一点 $A_1$，微分方程的轨线段记为 $\pi \left(A,A_1\right)$， $A_1$经对应参数脉冲映射 $\phi$到相集上一点 $A^{\prime }$，映射线段记为 $\phi \left(A_1,A^{\prime }\right)$， $\pi \left(A,A_1\right)+\phi \left(A_1,A^{\prime }\right)$称为脉冲参数为 $-\alpha$的阶一轨线段， $A^{\prime }$也称为A脉冲参数为 $-\alpha$的阶一后继点。设点A坐标为 $\alpha$， $A^{\prime }$点坐标为 $a^{\prime }$，则点A的后继函数 $F_{-\alpha }^1\left(A\right)=a^{\prime }-a<0$。

在相集N存在一点B，使B点坐标为 $\left(\tau ,\left(1-\beta \right)y^{*}\right)$， $\tau$为充分小的正数。过B微分方程的轨线与脉冲集与一点 $B_1$，微分方程的轨线段记为 $\pi \left(B,B_1\right)$， $B_1$经对应参数脉冲映射 $\phi$到相集上一点 $B^{\prime }$，映射线段记为 $\phi \left(B_1,B^{\prime }\right)$， $\pi \left(B,B_1\right)+\phi \left(B_1,B^{\prime }\right)$称为脉冲参数为 $-\alpha$的阶一轨线段， $B^{\prime }$也称为B脉冲参数为 $-\alpha$的阶一后继点。设点B坐标为b， $B^{\prime }$点坐标为 $b^{\prime }$，则点B的后继函数 $F_{-\alpha }^1\left(B\right)=b^{\prime }-b>0$(如图4所示)。

则点A和点B之间必存在一点D，设D的后继点为 $D^{\prime }$，设D点坐标为d， $D^{\prime }$点坐标为 $d^{\prime }$，则点D的后继函数为 ${F^{\prime }}_{-\alpha }\left(D\right)=d^{\prime }-d=0$。这就证明了系统(1)存在阶1周期解 $D{D}^{\prime }$(如图3所示)。

4.2. 阶1周期解的唯一性

已经知道A₂是A₀的后继点，且A₂在A₀的右侧，那么从A₂出发的轨线将会到达脉冲集，并形成交点A₃，在脉冲作用后，轨线到达相集上的点A₄，A₄是A₂的后继点，且A₄位于A₂的右侧。接着从A₄出发的轨线也将与脉冲集相交于点A₅，在脉冲作用后，轨线又达到相集上的点A₆，显然A₆是A₄的后继点，且A₆位于A₄的右侧。以此类推，可以得到两个序列

$x_{A_1},x_{A_3},x_{A_5},x_{A_7},\cdots ,x_{A_{2n-1}},\cdots \in \overline{A_1{B}_1}$，(28)

$x_{A_0},x_{A_2},x_{A_4},x_{A_6},\cdots ,x_{A_{2n}},\cdots \in \overline{A_0{B}_0}$.(29)

同理，还可以得到轨线从B₀出发的两个序列

$x_{B_1},x_{B_3},x_{B_5},x_{B_7},\cdots ,x_{B_{2n-1}},\cdots \in \overline{A_1{B}_1}$，(30)

$x_{B_0},x_{B_2},x_{B_4},x_{B_6},\cdots ,x_{B_{2n}},\cdots \in \overline{A_0{B}_0}$.(31)

根据这两列序列，可以看出线段 $\overline{A_1{B}_1}$通过脉冲作用被映射到线段 $\overline{A_0{B}_0}$，因此有 $\left(x_{A_1},x_{B_1}\right)\subset \left(x_{A_0},x_{B_0}\right)$， $d\left(A_1,B_1\right)<d\left(A_0,B_0\right)$。根据区间套定理，可以得到系统(1)在阶一周期解存在时是唯一的，且有 $\underset{n\to \infty }{\lim }{x}_{A_{2n}}=x_C$， $\underset{n\to \infty }{\lim }{x}_{B_{2n}}=x_C$成立(如图5所示)。

4.3. 阶1周期解的稳定性

定理2：当 $\alpha h{x}_{1}m>\frac{\beta u}{1-\beta }$时，系统(1)过 $A$点的阶1周期解轨道渐近稳定。

证：不妨设 $A\left(x_1,m\right)$点为脉冲点， $A^{+}\left(\left(1-\alpha \right)x_1,\left(1-\beta \right)m\right)$为A点的后继点。令

$P\left(x,y\right)=x\left(r-x\right)-kxy$， $Q\left(x,y\right)=hxy-fy+u$，(32)

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

计算得

$\frac{\partial P}{\partial x}=r-2x-ky$， $\frac{\partial Q}{\partial y}=hx-f$， $\frac{\partial \alpha }{\partial x}=-\alpha$， $\frac{\partial \beta }{\partial y}=-\beta$， $\frac{\partial \alpha }{\partial y}=\frac{\partial \beta }{\partial x}=\frac{\partial \Phi }{\partial x}=0$， $\frac{\partial \Phi }{\partial y}=1$，(34)

$\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-\alpha \right)Q\left(A^{+}\right)}{Q\left(A\right)}=\frac{\left(1-\alpha \right)\left[h\left(1-\alpha \right)x_1\left(1-\beta \right)m-f\left(1-\beta \right)m+u\right]}{h{x}_{1}m-fm+u}\end{array}$.(35)

$\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\right)\text{d}t}\\ ={\displaystyle {\int }_{\left(1-\alpha \right)x_1}^{x_1}\frac{1}{x}\text{d}x}+{\displaystyle {\int }_{\left(1-\beta \right)m_1}^{m_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{1}{1-\beta }+{\displaystyle {\int }_0^T\left(-\frac{u}{y}-x\right)\text{d}t}\end{array}$，(36)

${\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\}=\frac{h\left(1-\alpha \right)x_{1}m-fm+\frac{u}{1-\beta }}{h{x}_{1}m-fm+u}\exp {\displaystyle {\int }_0^T\left(-\frac{u}{y}-x\right)\text{d}t}$.(37)

当 $\alpha h{x}_{1}m>\frac{\beta u}{1-\beta }$时， $\left|{\mu }_2\right|<1$。由类Poincaré准则可知，系统(1)过A点的阶一周期解是轨道渐近稳定的。

5. 数值模拟和结论

基于前面对系统(1)的理论分析，通过数值模拟来验证污染环境下投放不同数量的捕食者(鱼类)对该系统的影响。当系统存在正平衡点时，取参数值为 $r=0.9,k=0.25,h=0.3,f=0.5,u=0.1$，此时，系统(1)的捕食者与食饵的时间序列图如图6、图7所示。当我们保留其他参数不变，使 $u=0.5$，系统(1)的捕食者与食饵的时间序列图如图8、图9所示。

可以说明当污染环境下，当我们投放少量的捕食者(鱼类)，食饵(浮游植物)占据了更多的生存空间，使鱼类的持续生存受到威胁，使收获量减小。但当我们加大投放量时，浮游动物适量减小，而鱼类的数量和收获量大幅度增加。我们通过人工对投放鱼类进行合理控制，这样既增加了资源，也能够维持渔业的正常生存。

由此可知，人工投放资源是维护生态平衡的重要手段。同时我们更应该为保护环境付诸行动，合理利用自然资源，防止环境的污染和破坏，这样才能使自然环境、人文环境、经济环境共同平衡可持续发展，扩大有用资源的再生产，保障社会的稳定可持续发展。

基金项目

本研究由北京建筑大学研究生创新项目资助，项目号：PG2022141。

NOTES

^*通讯作者。

参考文献

[1]	陈兰荪. 害虫治理与半连续动力系统几何理论[J]. 北华大学学报(自然科学版), 2011, 12(1): 1-9.
[2]	Liu, B., Tian, Y. and Kang, B. (2012) Dynamics on a Holling II Predator-Prey Model with State-Dependent Impulsive Control. International Journal of Biomathematics, 5, Article ID: 1260006. https://doi.org/10.1142/S1793524512600066
[3]	黄立壮, 刘琼, 陈武大仁, 等. 一类具有状态脉冲反馈控制的害虫治理模型[J]. 广西大学学报(自然科学版), 2020, 45(4): 911-919.
[4]	Das, T., Mukherjee, R.N. and Chaudhuri, K.S. (2009) Harvesting of a Prey-Predator Fishery in the Presence of Toxicity. Applied Mathematical Model-ling, 33, 2282-2292. https://doi.org/10.1016/j.apm.2008.06.008
[5]	Truscott, J.E. and Brindley, J. (1994) Ocean Plankton Populations as Excitable Media. Bulletin of Mathematical Biology, 56, 981-998.
[6]	Bhatia, S.K., Chauhan, S. and Agarwal, A. (2017) A Stage-Structured Prey-Predator Fishery Model in the Presence of Toxicity with Taxation as a Control Parameter of Harvesting Effort. Journal of Nonlinear Analysis and Application, 2017, 83-104.
[7]	蓝桂杰, 陈哲文, 魏春金, 等. 污染环境下毒素脉冲输入和心理效应对随机捕食-食饵系统的影响[J]. 系统科学与数学, 2019, 39(12): 2070-2092.
[8]	张茜, 刘兵, 石丽云. 在脉冲污染环境中具有阶段结构的捕食食饵Gompertz模型的动力学性质[J]. 生物数学学报, 2010(2): 9.
[9]	何伟莲. 具有反馈控制和有毒物质的竞争系统的持久性和全局吸引性[J]. 商丘师范学院学报, 2007, 23(12): 31-39.
[10]	Wang, F., Meng, X.-Z. and Cheng, H.-D. (2010) Dy-namical Analysis of a Chemostat-Type Competition Model with Pulsed Input in a Polluted Environment. Mathematica Applicata, 23, 204-212.
[11]	陈兰荪, 程惠东. 害虫综合防治建模驱动“半连续动力系统理论”兴起[J]. 数学建模及其应用, 2021, 10(1): 1-16. https://doi.org/10.19943/j.2095-3070.jmmia.2021.01.01