北京建筑大学理学院，北京

收稿日期：2022年7月24日；录用日期：2022年8月17日；发布日期：2022年8月29日

摘要

基于水环境污染日益严重，本文建立了一类毒素影响下鱼类具有常数投放的捕食–食饵模型，研究了系统平衡点的存在性与稳定性，利用后继函数法证明了系统阶一周期解存在性，其次利用类Poincaré准则，得到了阶一周期解稳定的条件，最后通过数值模拟验证结论的重要性。

关键词

毒素，常数投放，阶一周期解，类Poincaré准则，数值模拟

Dynamical Research of a Predator-Prey Model of Predator Population with Constant Release Rate under Toxin Influence

Yan Zhang^*, Meng Zhang

School of Science, Beijing University of Civil Engineering and Architecture, Beijing

Received: Jul. 24^th, 2022; accepted: Aug. 17^th, 2022; published: Aug. 29^th, 2022

ABSTRACT

In this paper, a predator-prey model of fish with constant release under the influence of toxins is established based on the increasing pollution of water environment. The existence and stability of the equilibrium point of the system are studied, and the existence of the order-1 periodic solution is proved by using the method of successor function. Secondly, the condition of the stability of the order-1 periodic solution is obtained by using the analogue of Poincaré’s criterion. Finally, the importance of the conclusion is verified by numerical simulation.

Keywords:Toxin, Constant Release, Order-1 Periodic Solution, The Analogue of Poincaré’s Criterion, Numerical Simulation

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

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所示。

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

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₂出发的解一致有界。

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

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所示。

图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所示)。

图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所示)。

图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所示。

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

图6. $u=0.1$，x的时间序列图

图7. $u=0.1$，y的时间序列图

图8. $u=0.5$，x的时间序列图

图9. $u=0.5$，y的时间序列图

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

基金项目

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

文章引用

张 琰,张 蒙. 毒素影响下捕食种群具有常数投放率的捕食–食饵模型的动力学研究
Dynamical Research of a Predator-Prey Model of Predator Population with Constant Release Rate under Toxin Influence[J]. 应用数学进展, 2022, 11(08): 6044-6054. https://doi.org/10.12677/AAM.2022.118636

参考文献