长春理工大学理学院，吉林 长春

收稿日期：2020年12月14日；录用日期：2021年1月3日；发布日期：2021年1月15日

摘要

本文考虑了脉冲控制下的一个带有毒素效应功能反应的植物–食植动物模型。理论分析得到了系统关于植物灭绝的平凡周期解的存在性以及它的局部和全局渐近稳定性条件。数值模拟验证了系统存在非常复杂的动力学行为，包括双稳态以及非平凡周期解分支等。结果表明毒素效应不会改变平凡周期解的局部动力学，但它改变了系统的全局稳定性条件。虽然脉冲控制是可行的，但是毒素效应确实增加了控制的难度。

关键词

植物–食植动物模型，毒素功能反应，周期解，脉冲控制

Impulsive Control of One Plant-Herbivore Model with Toxin-Determined Functional Response

Tianyu Liu, Shuai Wang^*

School of Science, Changchun University of Science and Technology, Changchun Jilin

Received: Dec. 14^th, 2020; accepted: Jan. 3^rd, 2021; published: Jan. 15^th, 2021

ABSTRACT

In this paper, one impulsive plant-herbivore model with toxin-determined functional response is considered. The trivial periodic solution corresponding to plant extinction, and its local and global asymptotic stability conditions are established. Some dynamical behavior, such as bistability and bifurcation of non-trivial periodic solutions are verified by numerical simulation. The results show the toxin effect does not change the local dynamics of the trivial periodic solutions in the plant-herbivore model with impulsive control, but it changes their global stability conditions. Although impulsive control is feasible in this model, the toxin effect does increase the difficulty of impulsive control.

Keywords:Plant-Herbivore Model, Toxin-Determined Functional Response, Periodic Solution, Impulsive Control

Copyright © 2021 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. 引言

传统的Holling-II功能反应通常用来描述食植动物对植物的消耗，但是考虑到某些植物能够释放一些化学毒素抵制食植动物，传统的捕食模型不能精确描述植物与食植动物之间的相互作用。Feng和DeAngelis等人 [1] 构建了一个由毒素效应决定的功能反应，其可以描述植物–食植动物在毒素效应下的相互作用。这种由毒性决定的功能反应模型为植物的这种化学毒素防御理论提供了非常有力的理论基础，因此受到了广泛的关注。截至目前所有关于毒素效应的工作都是在连续系统中完成的 [2] [3] [4]，但是脉冲微分方程这种半连续的系统更适合模拟实际生产中的脉冲控制现象，因此在本文中我们主要研究在植物–食植动物Beddington-DeAngelis模型中实施脉冲控制的可行性。

Holling II型功能反应可以对一般情形下的食草动物和植物的摄食率进行了恰当的描述，但现实中也可能存在这样的情况：由于食植动物个体之间的相互干扰，食植动物种群密度的增加意味着食植动物的摄食率降低 [5] [6] [7]，这就是Beddington-DeAngelis功能反应被提出的主要原因。实际上，在 [8] [9] 中，Beddington-DeAngelis功能反应已经被用来描述植物的毒素效应。本文中研究的带有毒素效应功能反应植物–食植动物脉冲控制的Beddington-DeAngelis (B-D)模型如下所示

$\{\begin{array}{l}\begin{array}{l}\stackrel{\dot{}}{x}=rx\left(1-\frac{x}{K}\right)-D\left(x,y\right)y,\\ \stackrel{\dot{}}{y}=\gamma D\left(x,y\right)y-dy,\end{array}\}t\ne nT,\\ y\left(n{T}^{+}\right)=y\left(nT\right)+\mu Tt=nT.\end{array}$ (1.1)

这里

$D\left(x,y\right)=B\left(x,y\right)\left(1-\frac{B\left(x,y\right)}{4G}\right),B\left(x,y\right)=\frac{cx}{1+hcx+\beta uy}.$

在方程(1.1)中 $x,y$ 分别代表植物和食植动物的种群数量。连续部分描述的是在毒素效应下植物与食植动物的数量变化规律，离散部分描述的是食植动物的脉冲释放，这里T代表释放周期，μ代表单位时间的释放率。在本文中我们主要通过研究系统的脉冲控制下系统对应于植物灭绝的周期解的存在性以及稳定性来验证脉冲控制的可行性。

2. 动力学分析

定理2.1系统(1.1)存在一个对应于植物灭绝的平凡周期解

$\left(\bar{x}\left(t\right),\bar{y}\left(t\right)\right)=\left(0,y^{\ast }{\text{e}}^{-d\left(tmodT\right)}\right)$ (2.1)

这里

$y^{\ast }=\frac{\mu T}{1-{\text{e}}^{-dT}}.$ (2.2)

该周期解是局部渐近稳定的当且仅当

$r<c/\beta$

并且

$\mu >{\mu }_1=\frac{1}{\beta T}\frac{\left(1-{\text{e}}^{-\left(r\beta /c\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{{\text{e}}^{-\left(r\beta /c\right)dT}-{\text{e}}^{-dT}}.$ (2.3)

证明：考虑x(0) = 0，那么系统(1.1)可以简化为

$\{\begin{array}{l}\stackrel{\dot{}}{x}=0,\\ \stackrel{\dot{}}{y}=-dy,\\ y\left(n{T}^{+}\right)=y\left(nT\right)+\mu T.\end{array}$ (2.4)

则是 $\bar{x}\left(t\right)=0$ 平凡的。由上面方程的后两个式子可知

$y\left(\left(n+1\right)T^{+}\right)=y\left(n{T}^{+}\right){\text{e}}^{-dT}+\mu T,$ (2.5)

它是指数稳定的离散系统，并且随着 $n\to \infty$ 其收敛到

$y^{\ast }=\frac{\mu T}{1-{\text{e}}^{-dT}}.$ (2.6)

为了讨论局部渐近稳定性，我们定义一个小的摄动 $\tilde{x}=x-\bar{x}=x,\tilde{y}=y-\bar{y}$。那么原系统等价于下面系统

$\{\begin{array}{l}\stackrel{\dot{}}{\tilde{x}}=r\tilde{x}\left(1-\tilde{x}/K\right)-cD\left(\tilde{x},\bar{y}+\bar{y}\right)\left(\bar{y}+\tilde{y}\right),\\ \stackrel{\dot{}}{\tilde{y}}=\gamma cD\left(\tilde{x},\bar{y}+\tilde{y}\right)\left(\bar{y}+\tilde{y}\right)-d\tilde{y}.\end{array}$ (2.7)

这里脉冲部分消失是因为

$\tilde{y}\left(n{T}^{+}\right)=y\left(nT\right)+\mu T+\bar{y}\left(nT\right)-\mu T=\tilde{y}\left(nT\right).$

接下来我们研究周期系统(2.2)。我们通过研究该周期系统零解的稳定性来研究原系统平凡周期解的稳定性 [10] [11]。考虑如下矩阵方程

$\left(\begin{array}{c}\tilde{x}\left(nT\right)\\ \tilde{y}\left(nT\right)\end{array}\right)=\Phi \left(t\right)\left(\begin{array}{c}\tilde{x}\left(0\right)\\ \tilde{y}\left(0\right)\end{array}\right),$ (2.8)

和

$\frac{\text{d}\Phi \left(t\right)}{\text{d}t}=A\left(t\right)\Phi \left(t\right).$ (2.9)

这里

$A\left(t\right)=\left(\begin{array}{c}r-\frac{2\tilde{x}}{K}-c\frac{\partial D\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{x}}\left(\bar{y}+\tilde{y}\right)-c\left(D\left(\tilde{x},\bar{y}+\tilde{y}\right)-c\left(\bar{y}+\tilde{y}\right)\frac{\partial D\left(\tilde{x},\bar{y}\right)}{\partial \tilde{y}}\right)\\ \frac{\partial D\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{x}}\left(\bar{y}+\tilde{y}\right)D\left(\tilde{x},\bar{y}+\tilde{y}\right)+\frac{\partial D\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{y}}\left(\bar{y}+\tilde{y}\right)-d\end{array}\right).$

由于 ${D\left(\tilde{x},\bar{y}+\tilde{y}\right)|}_{\tilde{x}=0}=0,{\frac{\partial B\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{y}}|}_{\tilde{x}=0,\tilde{y}=0}=0$，并且

$\frac{\partial D\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{x}}=\left(1-\frac{B\left(\tilde{x},\bar{y}+\tilde{y}\right)}{2G}\right)\frac{\partial B\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{x}},$

$\frac{\partial D\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{y}}=\left(1-\frac{B\left(\tilde{x},\bar{y}+\tilde{y}\right)}{2G}\right)\frac{\partial B\left(\tilde{x},\bar{y}+\tilde{y}\right)}{\partial \tilde{y}},$

那么

$A\left(t\right)=\left(\begin{array}{c}r-\frac{c\bar{y}}{1+\beta y}0\\ \frac{1}{1+\beta y}-d\end{array}\right),$

又

$\left(\begin{array}{c}\tilde{x}\left(n{T}^{+}\right)\\ \tilde{y}\left(n{T}^{+}\right)\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}\tilde{x}\left(nT\right)\\ \tilde{y}\left(nT\right)\end{array}\right),$ (2.10)

那么基本解矩阵满足

$\Gamma =\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\Phi \left(T\right)=\Phi \left(T\right).$ (2.11)

我们有 $\Phi \left(T\right)=\Phi \left(0\right)\exp \left(A\left(t\right)\right)$，$\Phi \left(0\right)$ 是单位矩阵。设 ${\lambda }_1,{\lambda }_2$ 是 $\exp \left(A\left(t\right)\right)$ 的特征根，那么有

${\lambda }_1={\text{e}}^{{\displaystyle {\int }_0^T\left(r-\frac{c\bar{y}}{1+\beta \bar{y}}\right)\text{d}t}},{\lambda }_2={\text{e}}^{-{\displaystyle {\int }_0^{T}d\text{d}t}}.$

由Floquet定理，我们可以得到平凡周期解 $\left(\bar{x}\left(t\right),\bar{y}\left(t\right)\right)=\left(0,y^{\ast }{\text{e}}^{-d\left(tmodT\right)}\right)$ 的局部稳定性条件是

${\lambda }_1={\text{e}}^{{\displaystyle {\int }_0^T\left(r-\frac{c\bar{y}}{1+\beta \bar{y}}\right)\text{d}t}}<1,{\lambda }_2={\text{e}}^{{\displaystyle {\int }_0^T\left(-d\right)\text{d}t}}<1.$

第二个条件是平凡的，对于第一个条件有

${\displaystyle {\int }_0^{T}c{y}^{\ast }{\text{e}}^{-dt}/\left(1+\beta y^{\ast }{\text{e}}^{-dt}\right)\text{d}t}>rT$

$\iff \left({\text{e}}^{-\left(r\beta /c\right)dT}-{\text{e}}^{-dT}\right)\beta y^{\ast }>\left(1-{\text{e}}^{-\left(r\beta /c\right)dT}\right)$ (2.12)

不等式右端恒大于0，因此我们只需

${\text{e}}^{-\left(r\beta /c\right)dT}-{\text{e}}^{-dT}>0$

也就是

$r<c/\beta$

将 代入公式可得

$\mu >{\mu }_1=\frac{\left(1-{\text{e}}^{-\left(r\beta /c\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{\beta T\left({\text{e}}^{-\left(r\beta /c\right)dT}-{\text{e}}^{-dT}\right)}.$

因此结论得证。

平凡周期解的全局稳定性条件

由于毒素作用，系统的动力学变得复杂。这可能会增加脉冲控制的难度。在这一部分中，我们研究了通过研究系统平凡周期解的全局稳定性条件来讨论毒素对脉冲控制的影响。

我们首先讨论系统(1.1)解的有界性。为了保证 $D\left(x,y\right)$ 总是正数，我们对G有如下要求：

$G>\frac{\max \left\{B\left(x,y\right)\right\}}{4}=\underset{x\to \infty }{\lim }\frac{B\left(x,y\right)}{4}=\frac{1}{4h}.$

那么有

定理2.2 如果 $G>\frac{1}{4h}$，系统(1.1)的解总是一致有界的。

证明：如果 $G>\frac{1}{4h}$，那么x有上界K。定义计算 $V\left(t\right)=\gamma x\left(t\right)+y\left(t\right)$，$V\left(t\right)$ 在 $t\in \left(nT,\left(n+1\right)T\right)$ 上的右导数有

$D^{+}V\left(t\right)=\gamma rx\left(1-x/K\right)-dy.$

那么存在常数 $0<s<d$ 满足：

$\begin{array}{c}{D}^{+}V\left(t\right)+sV\left(t\right)\le \gamma rx\left(1-\frac{x}{K}\right)+s\gamma x\\ =\gamma x\left(r+s\right)-\gamma r{x}^2/K\\ \le L,L=\frac{\gamma {\left(r+s\right)}^{2}K}{4r}\end{array}$

或者

$D^{+}V\left(t\right)\le -sV\left(t\right)+L.$

由于 $t=nT$

$V\left(n{T}^{+}\right)=\gamma x\left(n{T}^{+}\right)+y\left(n{T}^{+}\right)=V\left(nT\right)+\mu T.$

那么有

$V\left(\left(n+1\right)T^{+}\right)\le V\left(nT\right){\text{e}}^{-sT}+L/s+\mu T$

或者

$V^{\ast }\left(t\right)\le \frac{L/s+\mu T}{1-{\text{e}}^{-sT}}=B$ 当 $n\to \infty$。

所以V(t)是一致有界的。

定理2.3平凡周期解

$\left(\bar{x}\left(t\right),\bar{y}\left(t\right)\right)=\left(0,y^{\ast }{\text{e}}^{-d\left(tmodT\right)}\right)$ (2.13)

是全局吸引的若下面条件成立

(H1) $G>\frac{1}{4h}$，

(H2) $r<\frac{cM}{\beta }$，

(H3) $\mu >{\mu }_2=\frac{\left(1+K\right)\left(1-{\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{\beta T\left({\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}-{\text{e}}^{-dT}\right)}$。

这里 $M=1-\frac{K}{4G\left(1+\alpha K\right)}$。

证明：为了讨论全局吸引性，我们接着考虑系统(2.2) [11]。设系统在时刻 $t=t_0$ 的初值为 $\left({\tilde{x}}_0,{\tilde{y}}_0\right)$ 则关

于 $\stackrel{\dot{}}{x}$ 我们将方程(2.7)的第一式除以 $\frac{c\tilde{x}}{1+\alpha \tilde{x}}$ 并积分得

${\displaystyle {\int }_{t_0}^t\frac{\stackrel{\dot{}}{\tilde{x}}}{\tilde{x}}\text{d}\tau }={\displaystyle {\int }_{t_0}^t\left[r/c\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\left(1-\alpha \tilde{x}\right)}{c\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }$

为了得到 ${\lim }_{t\to |+\infty }\tilde{x}\left(t\right)=0$，只需满足当 $t\to +\infty$ 时， ${\displaystyle {\int }_{t_0}^t\frac{\stackrel{\dot{}}{\tilde{x}}}{c\tilde{x}}\text{d}\tau }\to -\infty$。

我们将时间区间划分如下：

$\begin{array}{c}{\displaystyle {\int }_{t_0}^t\frac{\stackrel{\dot{}}{\tilde{x}}}{\tilde{x}}\text{d}\tau }={\displaystyle {\int }_{t_0}^{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }\\ +{\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }\\ +{\displaystyle {\int }_{\left(\lfloor \frac{t}{T}\rfloor \right)T}^t\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }\end{array}$

易知第一、三项是有界的，现只需考虑第二项满足

${\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }\to -\infty$

随着 $t\to +\infty$。

接着我们做一些必要的放缩，由于

$\stackrel{\dot{}}{\tilde{y}}>-d\tilde{y}.$

接着， $\tilde{y}\left(t\right)>-\bar{y}\left(t\right)>-y^{\ast }$，那么

$\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\ge y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right).$

从而

$\begin{array}{l}{\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-\frac{D\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{\tilde{x}}\left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)\right]\text{d}\tau }\\ ={\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-c\left(1-\frac{B\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}{4G}\frac{\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)}{1+hc\tilde{x}+\beta \left(\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}\right)\right]\text{d}\tau }\\ \le {\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-cM\frac{\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)}{1+hc\tilde{x}+\beta \left(y\left(\tau \right)+\tilde{y}\left(\tau \right)\right)}\right]\text{d}\tau }\\ \le {\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-cM\frac{y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)}{1+hc\tilde{x}+\beta \left(y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)\right)}\right]\text{d}\tau }\\ \le {\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-cM\frac{y^{\ast }\left({\text{e}}^{-d\left(tmod\text{}T\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)}{1+hcK+\beta \left(y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)\right)}\right]\text{d}\tau },\end{array}$

这里

$\begin{array}{c}M=\min \left\{1-B\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)/4G\right\}\\ =1-\max \left\{B\left(\tilde{x}\left(\tau \right),\bar{y}\left(\tau \right)+\tilde{y}\left(\tau \right)\right)/4G\right\}\\ =1-\frac{K}{4G\left(1+\alpha K\right)}.\end{array}$ (2.14)

我们重新考虑

$\begin{array}{l}{\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-cM\frac{y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)}{1+hc\tilde{x}+\beta \left(y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)\right)}\right]\text{d}\tau }\\ ={\displaystyle {\sum }_{n=\left(\lfloor \frac{t_0}{T}\rfloor +1\right)}^{\lfloor \frac{t}{T}\rfloor -1}I}\left(n\right),\end{array}$ (2.15)

这里

$\begin{array}{c}I\left(n\right)={\displaystyle {\int }_{\left(\lfloor \frac{t_0}{T}\rfloor +1\right)T}^{\left(\lfloor \frac{t}{T}\rfloor \right)T}\left[r\left(1-\tilde{x}\left(\tau \right)/K\right)-cM\frac{y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)}{1+hc\tilde{x}+\beta \left(y^{\ast }\left({\text{e}}^{-d\left(tmodT\right)}-{\text{e}}^{-d\left(t-t_0\right)}\right)\right)}\right]\text{d}\tau }\\ =rT+\frac{cM}{\beta d}\cdot \ln \frac{1+hcK+\beta y^{\ast }{\text{e}}^{-dT}\left(1-{\text{e}}^{-d\left(nT-t_0\right)}\right)}{1+hcK+\beta y^{\ast }\left(1-{\text{e}}^{-d\left(nT-t_0\right)}\right)}.\end{array}$ (2.16)

若要满足 ${\lim }_{t\to |+\infty }\tilde{x}\left(t\right)=0$，只需

$\underset{n\to +\infty }{\lim }I\left(n\right)=rT+\frac{cM}{\beta d}\cdot \ln \frac{1+hcK+\beta y^{\ast }{\text{e}}^{-dT}}{1+hcK+\beta y^{\ast }}<0.$

将 $y^{\ast }$ 带入解出 $\mu$ 得

$\mu >{\mu }_2=\frac{\left(1+K\right)\left(1-{\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{\beta T\left({\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}-{\text{e}}^{-dT}\right)},$

这里 $r<\frac{cM}{\beta }$。

下面考虑 $\tilde{y}$。由于 $\tilde{x}$ 收敛到0，那么存在一个时刻 $t_m$ 使得

$\forall t>t_m,\gamma cD\left(\tilde{x},\bar{y}+\tilde{y}\right)<\frac{d}{2},$

也就是

$\forall t>t_m,\stackrel{\dot{}}{\tilde{y}}\le \gamma cD\left(\tilde{x},\bar{y}+\tilde{y}\right)\bar{y}-\frac{d}{2}\tilde{y}.$

那么由于 $\gamma cD\left(\tilde{x},\bar{y}+\tilde{y}\right)$ 趋于0随着 $t\to +\infty$，同样 $\tilde{y}$ 也趋于0随着 $t\to +\infty$。故证明了对于系统(2.7) (0,0)是全局吸引的，因此对于系统(1.1)平凡周期解也是全局吸引的。根据周期解的局部渐近稳定性和全局吸引性，我们可以得到系统 (1.1)平凡周期解的全局渐近稳定性。

注2.1 我们接下来与无毒素效应的植物–食植动物模型做对比。因为 $M=1-\frac{K}{4G\left(1+\alpha K\right)}$，所以有

$\underset{G\to +\infty }{\lim }M=1,$

并且

$\underset{G\to +\infty }{\lim }{\mu }_2=\frac{\left(1+K\right)\left(1-{\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{\beta T\left({\text{e}}^{-\left(r\beta /\left(cM\right)\right)dT}-{\text{e}}^{-dT}\right)}=\frac{\left(1+K\right)\left(1-{\text{e}}^{-\left(r\beta /c\right)dT}\right)\left(1-{\text{e}}^{-dT}\right)}{\beta T\left({\text{e}}^{-\left(r\beta /c\right)dT}-{\text{e}}^{-dT}\right)}\stackrel{\dot{}}{=}{\mu }_3$ (2.17)

易知 $\mu >{\mu }_3$ 正是无毒素效应的植物–食植动物模型在相同脉冲控制下平凡周期解的全局渐近稳定性条件。并且，由于 ${\mu }_2$ 是关于M的单调递减函数，所以 ${\mu }_2>{\mu }_3$ 如果 $G\ne +\infty$ 因此我们有结论，和无毒素效应的植物–食植动物脉冲模型相比，毒素效应下脉冲模型需要更大的释放率才能得到平凡周期解的全局稳定性，这增加了脉冲控制的难度。

3. 数值模拟

这部分我们主要展示模型(1.1)的动力学的复杂性。首先图1展示了系统存在着非平凡周期解分支。考虑如下参数取值： $r=2$，$K=8$，$c=1$，$h=0.45$，$\beta =0.0\text{2}$，$\gamma =0.\text{9}$，$d=0.25$，$T=1$。此时，系统平凡周期解得局部渐近稳定性条件得临界值是 ${\mu }_1=1.22$。在图1中，我们在关于参数 $\mu$ 的分支图中证

图1. 本图展示的是在不同的G的取值下系统关于参数μ的非平凡周期解分支。图中X轴是参数μ，Y轴是植物的数量(密度)。此时系统平凡周期解的局部渐近稳定性条件的临界值为μ = 1.22。(a) 当G = 0.6时，非平凡周期解的存在区间为 $\mu \in \left[1.2,1.22\right]$ ；(b) 当 $G=1$ 时，非平凡周期解的存在区间为 $\mu \in \left[1.17,1.22\right]$ ；(c) 当 $G=2$ 时，非平凡周期解的存在区间为 $\mu \in \left[1.14,1.22\right]$ ；(d) 当 $G=+\infty$ 时，非平凡周期解的存在区间为 $\mu \in \left[1.12,1.22\right]$。随着G取值的增加，非平凡周期解的存在区间也随着增大

实了上述结果，并且发现在平凡周期解失去稳定性以后系统会产生非平凡的周期解。在不同的G的取值下，图1给出了系统的关于参数 $\mu$ 的非平凡周期解分支，其中图1(d)对应的是没有毒素反应下的脉冲B-D模型的分支图，其对应的G的取值为 $G=+\infty$。具体来看，在图1(a)中，非平凡周期解的存在区间为[1.2, 1.22]；在图1(b)中，非平凡周期解的存在区间为[1.17, 1.22]；在图1(c)中，非平凡周期解的存在区间为[1.14, 1.22]；在图1(d)中，非平凡周期解的存在区间为[1.11, 1.22]。通过对比，我们发现随着G取值的增加，分支图中非平凡周期解的存在区间也随之变大；相反，G值的减小也造成了非平凡周期解的存在区间变小。

图2. 本图展示的是当 $\mu =1.25$ 时系统的双稳现象。(a) (b) 当初值为(0.001, 0.1)时，植物和食植动物的种群数量变化轨迹如图所示，此时对应于系统的平凡周期解；(c) (d) 当初值为(1, 0.1)时，植物和食植动物的种群数量变化轨迹如图所示，此时对应于系统的非平凡周期解。此时系统有两个稳定的周期解，其中一个是平凡的，另一个是非平凡的

图3. 本图展示的是当 $\mu =1.2$ 时系统的双稳现象。(a) (b) 当初值为(0.001, 0.1)时，植物和食植动物的种群数量变化轨迹如图所示，此时对应于系统的一个非平凡周期解；(c) (d) 当初值为(4, 0.1)时，植物和食植动物的种群数量变化轨迹如图所示，此时对应于系统的另一个非平凡周期解。此时系统有两个稳定的非平凡周期解

在图2和图3中，我们展示了系统的两个稳定的T-周期解的共存，即双稳现象。在图2中，我们选择参数取值 $G=0.48$，$\mu =1.25$。此时，由于参数 $\mu >{\mu }_1=1.22$。故系统有一个稳定的平凡T-周期解 $\left(0,y^{\ast }\right)$，如图2(a)，图2(b)所示。然而如果我们选择一个新的初值 $x=4,y=1$。由于这个初值远离平凡周期解，我们得到了一个稳定的非平凡的T周期解，如图2(c)，图2(d)所示。在图3中，我们选择参数取值为 $G=0.48$，$\mu =1.2<{\mu }_1=1.22$，此时由于非平凡周期解分支的存在性，原先的平凡周期解丧失稳定性转变成一个非平凡的周期解，此时我们可以得到两个不同的稳定的非平凡周期解。

4. 结论

本文研究了一类具有毒素效应功能反应的脉冲植物–食植动物B-D模型，得到了脉冲控制下系统的平凡周期解的局部和全局稳定性条件。结果表明在带有毒素环境下脉冲控制在B-D模型中是可行的，然而毒素效应改变了平凡周期解的全局稳定性条件，增加了脉冲控制的难度。

文章引用

刘天宇,王 帅. 带有毒素效应功能反应的植物–食植动物模型的脉冲控制研究
Impulsive Control of One Plant-Herbivore Model with Toxin-Determined Functional Response[J]. 应用数学进展, 2021, 10(01): 37-47. https://doi.org/10.12677/AAM.2021.101004

参考文献