一类具有相互干扰和Holling-III功能反应模型的定性分析

Qualitative Analysis on a Class of Holling-III Function Response Models with Mutual Interference

摘要: 通过分析模型在平衡点处的线性近似系统，发现该模型在某些条件下具有独特的正平衡点。当给定模型具有正平衡点时，运用Bendixson-Dulac定理获得系统在第一象限中的闭合轨道。借助Liénard变换，利用张芷芬的唯一性定理得到该系统极限环的存在性和唯一性。最后，利用Mathematica软件对模型进行了四组数值模拟，验证了所得结论的合理性。

Abstract: In this paper, by analyzing the linear approximation system of the model at the equilibrium point, we first find that the model has a unique positive equilibrium point under some conditions. When the given model has positive equilibrium point, the closed orbit of the system in the first quadrant is obtained via Bendixson-Dulac theorem. By using the Liénard transformation and the uniqueness theorem of Zhang Zhifen, the existence and uniqueness of the limit cycle of the system is obtained. Finally, four numerical simulations of the model are given to illustrate the validity of our results by using Mathematica software.

1. 前言

1965年C. S. Holling在实验的基础上，对于不同类型的物种，提出了三种形式的功能性反应函数，具有这三类功能反应函数的食饵–捕食系统已经有大量的研究 [1] [2] [3] [4] [5] 。然而人们又发现，当营养资源达到一个较高水平时，对某些个体的增长率呈抑制作用，基于这样的生物背景，文献 [6] 提出了Monod-Haldane型功能反应函数。后来人们将Monod-Haldane型功能反应函数称之为Holling-IV函数 [7] 。

在一般情况下，具有Holling类型的功能反应的捕食–被捕食系统可表述为下面的微分方程组模型 [8] ：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=xg\left(x\right)-y\varphi \left(x\right),\\ \frac{\text{d}y}{\text{d}t}=-yq\left(y\right)+ky\varphi \left(x\right),\end{array}$

其中 $g\left(x\right)$ 为被捕食者种群的相对增长率(当没有捕食者时)， $\varphi \left(x\right)$ 为捕食者功能反应函数， $q\left(y\right)$ 为捕食者种群的死亡率，k为常数。由于食饵种群的生存环境可能会因为各种原因而遭到不同程度的破坏，所以食饵种群的增长率遵循的函数是非线性函数 $g\left(x\right)$ ；又因为捕食者种群在捕食过程中相互之间可能会存在干扰，设定 $m\left(0<m\le 1\right)$ 为系统中的干扰系数，同时捕食者种群的功能性反应函数 $\varphi \left(x\right)$ 也是非线性的。基于上述文献的讨论提出如下具相互干扰和Holling-III功能反应模型：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)-\frac{k{x}^2{y}^m}{1+x^2},\\ \frac{\text{d}y}{\text{d}t}=y\left(-d+\frac{ek{x}^2}{1+x^2}\right),\end{array}$ (1)

其中：x，y分别表示在t这个时刻食饵和捕食者种群密度，a，b，c分别表示食饵种群的内禀增长率、密度制约系数和移除的比例系数，函数 $a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)$ 表示该食饵种群内部增长率；e，d分别表示捕

食者种群内的变换系数(捕获食饵后转化为生育率的比例系数)和死亡率，而 $\frac{k{x}^2}{1+x^2}$ 表示Holling-III功能性反应函数。显然系数 $a,b,c,d,e,k,\alpha ,m$ 均是正的常数，且 $0<m\le 1$ ， $\alpha \ge 2$ ， $h\left(x\right)$ 需要满足条件：

$h\left(0\right)=0,h^{\left(n\right)}\left(x\right)\ge 0,n=1,2,3,4,x>0.$

同时，根据模型的生态意义，需在第一象限内讨论种群模型，记

$R_{\text{+}}=\left\{\left(x,y\right)|x>0,y>0\right\},{\bar{R}}_{\text{+}}=\left\{\left(x,y\right)|x\ge 0,y\ge 0\right\}.$

2. 预备知识

定义1 [9] 设常系数齐次线性系统

$\{\begin{array}{l}\stackrel{\dot{}}{x}=a_{11}{x}_1+a_{12}{x}_2,\\ \stackrel{\dot{}}{y}=a_{21}{x}_1+a_{22}{x}_2,\end{array}$

的向量形式： $\stackrel{\dot{}}{x}=Ax$ ，其中：

$x=\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right),A=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right)$ ，

用T表示矩阵A的迹，D表示矩阵A的行列式，设 $\Delta =T^2-4D$ 。对于系统的平衡点 $O\left(0,0\right)$ ，此时 $D\ne 0$ ，那么称 $O\left(0,0\right)$ 为初等奇点；当 $D\equiv 0$ 时，称 $O\left(0,0\right)$ 为高次奇点。对于初等奇点，有如下的分类：

1) 当 $D>0$ 时，平衡点为鞍点；

2) 当 $D>0$ ， $\Delta <0$ ， $T<0\left(T>0\right)$ 时，平衡点为稳定(不稳定)结点；

3) 当 $D>0$ ， $\Delta >0$ ， $T<0\left(T>0\right)$ 时，平衡点为稳定(不稳定)焦点；

4) 当 $D>0$ ， $T=0$ 时，平衡点为中心。

定义2 [9] 考察非线性振动方程

$\frac{{\text{d}}^{2}x}{\text{d}{t}^2}+f\left(x\right)\frac{\text{d}x}{\text{d}t}+g\left(x\right)=0$ ，

可通过变换 $x=x,y=-{\displaystyle {\int }_0^{t}g\left(x\right)\text{d}t}$ ，化为一阶方程组：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=y-F\left(x\right),\\ \frac{\text{d}y}{\text{d}t}=-g\left(x\right).\end{array}$ (2)

该方程组称为Liénard方程组，而该方程组对应的方程称作Liénard方程，变换 $x=x$ ， $y=-{\displaystyle {\int }_0^{t}g\left(x\right)\text{d}t}$ 称为Liénard变换。

引理1 [8] 考虑Liénard方程组，若曲线 $F\left(x_1\right)=F\left(x_2\right)$ 与曲线 $G\left(x_1\right)=G\left(x_2\right)$ 在区域： $D=\left\{\left(x_1,x_2\right),x_{02}<x_2<0,0<x_1<x_{01},x_{02}<x_{01}\right\}$ 内无交点，这里 $G\left(x\right)={\displaystyle {\int }_0^{x}g\left(t\right)\text{d}t}$ ，则Liénard方程在带域 $x_{02}<x<x_{01}$ 中没有极限环。

引理2 (Poincare-Bendixson环域定理 [9] )设Γ是由两条单闭曲线L₁与L₂所构成的环域，并且在Γ内不存在系统的奇点；当t增加时，从L₁与L₂出发的轨线都进入(或者离开) Γ，则Γ内存在闭轨线L，且位置为： $L_1\subset L\subset L_2$ 。

3. 模型的定性分析

作时间变换 $\left(1+x^2\right)\text{d}\tau =\text{d}t$ ，仍然记 $\tau$ 为t，化为

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)-k{x}^2{y}^m,\\ \frac{\text{d}y}{\text{d}t}=y\left[\left(ek-d\right)x^2-d\right].\end{array}$ (3)

再令 $x=\bar{x}$ ， $k{y}^m=\bar{y}$ ， $\left(ek-d\right)=\omega$ ，仍然记 $\bar{x},\bar{y}$ 为 $x,y$ ，则(3)化为

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)-x^{2}y,\\ \frac{\text{d}y}{\text{d}t}=my\left(\omega x^2-d\right).\end{array}$ (4)

其中 $a,b,c,d$ 均为正数， $\omega$ 符号不定。

3.1. 模型平衡点的稳定性分析

令系统(4)右边等于零，得到两个平衡点： $O\left(0,0\right)$ ， $P_0\left(x_0,y_0\right)$ ，其中： $x_0=\sqrt{d/\omega }$ ，

$y_0=\frac{1}{x_0}\left(1+x_0^2\right)\left(a-b{x}_0^{\alpha }-c{x}_0^{2\alpha }-h\left(x_0\right)\right)$ 。当且仅当 $\omega >0$ ， $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ 时， $P_0\left(x_0,y_0\right)$ 是 $R_{+}$ 内唯一的正平衡点。

因为系统在平衡点 $O\left(0,0\right)$ 处的线性近似系统为 $\frac{\text{d}x}{\text{d}t}=ax$ ， $\frac{\text{d}y}{\text{d}t}=-\text{d}y$ ，该线性近似系统对应的矩阵为：

$A_1=\left(\begin{array}{cc}a& 0\\ 0& -d\end{array}\right)$ ，

易知 $\det \left(A_1\right)\ne 0$ 且 $\det \left(A_1\right)<0$ ，由此即知 $O\left(0,0\right)$ 是鞍点。令

$Q\left(x\right)=b{x}^{\alpha }+c{x}^{2\alpha }+h\left(x\right)-a,$

由 $Q\left(0\right)=-a<0$ 及 $Q^{\prime }\left(x\right)=\alpha b{x}^{\alpha }+2c\alpha x^{2\alpha }+h^{\prime }\left(x\right)>0$ ，当 $x>0$ 时， $Q\left(x\right)$ 有唯一的正零点 $x_1$ 。当 $\omega >0$ ， $x_1>x_0$ 时，系统在平衡点 $E\left(x_1,0\right)$ 处的近似系统对应的矩阵 $A_2$ 为：

$A_2=\left(\begin{array}{cc}{x}_1\left(1+x_1^2\right)\left(-\alpha b{x}_1^{\alpha -1}-2c\alpha x_1^{2\alpha -1}-h^{\prime }\left(x_1\right)\right)& -x_1^2\\ 0& m\left(\omega x_1^2-d\right)\end{array}\right)$ 。

显然 $x_1\left(1+x_1^2\right)\left(-\alpha b{x}_1^{\alpha -1}-2c\alpha x_1^{2\alpha -1}-h^{\prime }\left(x_1\right)\right)<0$ ，因为 $x_1>x_0$ ， $\omega x_0^2-d=0$ ，所以有 $m\left(\omega x_1^2-d\right)>0$ ，得到 $\det \left(A_2\right)<0$ ，即 $E\left(x_1,0\right)$ 是鞍点。当 $\omega >0$ ， $0<x_1<x_0$ 时， $m\left(\omega x_1^2-d\right)<0$ ，得到 $\det \left(A_2\right)>0$ ， $tr\left(A_2\right)<0$ ，即 $E\left(x_1,0\right)$ 是稳定结点，并在 ${\bar{R}}_{+}$ 内是渐近稳定的。同理可知当 $\omega \le 0$ 时， $E\left(x_1,0\right)$ 在 ${\bar{R}}_{+}$ 内是渐近稳定结点。

令

$\begin{array}{c}P=a\left(x_0^2-1\right)-b\left[\alpha -1+\left(\alpha +1\right)x_0^2\right]x_0^{\alpha }-c\left[2\alpha -1+\left(2\alpha +1\right)x_0^2\right]x_0^{2\alpha }\\ -\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)-\left(x_0^2-1\right)h\left(x_0\right).\end{array}$

由定义1可以得到，若 $\omega >0$ ， $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ， $P<0\left(>0\right)$ ，则 $P_0$ 是系统(1)的稳定(不稳定)焦点或结点；若 $\omega >0$ ， $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ， $P=0$ ，则 $P_0$ 是系统(1)的中心–焦点型奇点。

3.2. 闭轨的不存在性

定理1 当 $\omega <0$ 时，系统(1)不存在全部位于 ${\bar{R}}_{+}$ 内的闭轨线。

证明 作Dulac函数 $B\left(x,y\right)=x^r{y}^{\theta }$ ，这里 $r,\theta$ 是待定常数。根据系统(1)可得

$BP=x^{r+1}{y}^{\theta }\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)-x^{r+2}{y}^{\theta +1}$ ， $BQ=m{x}^r{y}^{\theta +1}\left(\omega x^2-d\right)$ 。

对 $\forall \left(x,y\right)\in {\bar{R}}_{+}$ ，总有

$\begin{array}{c}\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}=B\{\left[\left(r+1\right)a-\left(\theta +1\right)md\right]+\left[\left(r+3\right)a+\left(\theta +1\right)m\omega \right]x^2-\left[\left(r+1\right)b+\alpha b\right]x^{\alpha }\\ -\left[\left(r+3\right)b+\alpha b\right]x^{\alpha +2}-\left[\left(r+1\right)c+2\alpha c\right]x^{2\alpha }-\left[\left(r+3\right)c+2\alpha c\right]x^{2\alpha +2}\\ -\left(r+1\right)h\left(x\right)-\left(r+3\right)x^{2}h\left(x\right)-x\left(1+x^2\right)h^{\prime }\left(x\right)-\left(r+2\right)xy\}.\end{array}$

令 $r+1=0$ ， $a\left(r+3\right)+\left(\theta +1\right)m\omega =0$ 。解得 $r=-1$ ， $\theta =-1-\frac{2a}{m\omega }$ 。由 $\omega <0$ 知

$\begin{array}{l}\left(r+1\right)b+\alpha b=\alpha b>0,\left(r+3\right)b+\alpha b=2b+\alpha b>0,\left(r+1\right)c+2\alpha c=2\alpha c>0,\\ \left(r+3\right)c+2\alpha c=2c+2\alpha c>0,\left(r+1\right)a-\left(\theta +1\right)md=md\frac{2a}{m\omega }=\frac{2ad}{w}<0,\end{array}$

从而有 $\frac{\partial \left(BP\right)}{\partial x}+\frac{\partial \left(BQ\right)}{\partial y}<0$ 。据Bendixson-Dulac判别法，定理1成立。

3.3. 极限环的存在性和唯一性分析

为证明系统(1)极限环的不存在性和存在唯一性，将系统(4)化为Liénard系统，令 $u=\ln y$ ， $-x^2\text{d}t=\text{d}\tau$ ， $x=v$ ，仍然记 $v,u,\tau$ 为 $x,y,t$ ，则系统(4)化为

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}={\text{e}}^y-l\left(x\right)\equiv \phi \left(y\right)-F\left(x\right),\\ \frac{\text{d}y}{\text{d}t}=-\frac{m\left(\omega x^2-d\right)}{x^2}\equiv -g\left(x\right),\end{array}$ (5)

其中：

$l\left(x\right)=\frac{\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)}{x}$ 。

奇点从 $\left(x_0,y_0\right)$ 变为 $\left(x_0,ln{y}_0\right)$ 。

定理2 如果 $P_0$ 是稳定的焦点或是结点或是中心–焦点型焦点，即

$\omega >0,a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right),P\le 0,$

且函数 $h\left(x\right)$ 满足条件

$h\left(0\right)=0,h^{\left(n\right)}\left(x\right)\ge 0,n=1,2,3,4,x>0,$

则 $P_0\left(x_0,y_0\right)$ 在 $R_{+}$ 内是全局渐近稳定。

证明 仅须证明系统(5)在 $R_{+}$ 内无极限环，为此记

$G\left(x\right)={\displaystyle \int g\left(x\right)\text{d}x}={\displaystyle \int \frac{m\left(\omega x^2-d\right)}{x^2}\text{d}x}=m\omega x+\frac{md}{x}+l.$

由引理1仅须证明两曲线

$\{\begin{array}{l}F\left(u\right)=F\left(v\right),\\ G\left(u\right)=G\left(v\right),0<u<x_0<v,\end{array}$ (6)

无交点。假设这两曲线有交点，那么满足(6)式，有：

$\begin{array}{c}0=a-a{x}_0^2+b{x}_0^2\left(\frac{v^{\alpha -1}-u^{\alpha -1}}{v-u}\right)+c{x}_0^2\left(\frac{v^{2\alpha -1}-u^{2\alpha -1}}{v-u}\right)+b{x}_0^2\left(\frac{v^{\alpha +1}-u^{\alpha +1}}{v-u}\right)\\ +c{x}_0^2\left(\frac{v^{2\alpha +1}-u^{2\alpha +1}}{v-u}\right)+x_0^2\left(\frac{\frac{h\left(v\right)}{v}-\frac{h\left(u\right)}{u}}{v-u}\right)+x_0^2\left(\frac{vh\left(v\right)-uh\left(u\right)}{v-u}\right).\end{array}$

由祁锋不等式 [10] ：

$\frac{v^{\alpha \text{+}1}-u^{\alpha \text{+}1}}{v-u}>\left(\alpha +1\right){\left(uv\right)}^{\frac{\alpha }{2}},$

可得：

$\begin{array}{c}0>a\left(1-x_0^2\right)+b{x}_0^2\left(\left(\alpha -1\right)x_0^{\alpha -2}+\left(\alpha +1\right)x_0^{\alpha }\right)+cuv\left(\left(\text{2}\alpha -1\right)x_0^{2\alpha -2}+\left(\text{2}\alpha +1\right)x_0^{2\alpha }\right)\\ +x_0^2\left(\frac{\frac{h\left(v\right)}{v}-\frac{h\left(u\right)}{u}}{v-u}\right)+x_0^2\left(\frac{vh\left(v\right)-uh\left(u\right)}{v-u}\right).\end{array}$ (7)

下面证明：

$\frac{vh\left(v\right)-uh\left(u\right)}{v-u}\ge {\left(xh\left(x\right)\right)}^{\prime }{}_{x=x_0},$ (8)

$\frac{\frac{h\left(v\right)}{v}-\frac{h\left(u\right)}{u}}{v-u}\ge {\left(\frac{h\left(x\right)}{x}\right)}^{\prime }{}_{x=x_0}$ (9)

先证明式(8)成立。记 $H\left(x\right)=xh\left(x\right)$ 。由泰勒公式有：

$\frac{H\left(v\right)-H\left(u\right)}{v-u}=H^{\prime }\left(x_0\right)+\frac{1}{2\left(v-u\right)}\left[H^{″}\left({\xi }_2\right){\left(v-x_0\right)}^2-H^{″}\left({\xi }_1\right){\left(u-x_0\right)}^2\right].$

由 $v-x_0=\frac{x_0}{u}\left(x_0-u\right)>\left(x_0-u\right)>0$ ，可得：

$\begin{array}{c}{H}^{″}\left({\xi }_2\right){\left(v-x_0\right)}^2-H^{″}\left({\xi }_1\right){\left(u-x_0\right)}^2>{\left(u-x_0\right)}^2\left(H^{″}\left({\xi }_2\right)-H^{″}\left({\xi }_1\right)\right)\\ ={\left(u-x_0\right)}^2{H}^{\left(3\right)}\left({\xi }_0\right)\left({\xi }_2-{\xi }_1\right),{\xi }_1<{\xi }_0<{\xi }_2.\end{array}$

因为

$H^{\left(3\right)}\left({\xi }_0\right)=3h^{″}\left({\xi }_0\right)+x{h}^{\left(3\right)}\left({\xi }_0\right)\ge 0,$

从而

$\frac{H\left(v\right)-H\left(u\right)}{v-u}>H^{\prime }\left(x_0\right)+\frac{1}{2\left(v-u\right)}{\left(u-x_0\right)}^2{H}^{\left(3\right)}\left({\xi }_0\right)\left({\xi }_2-{\xi }_1\right)>H^{\prime }\left(x_0\right),$

所以在定理条件下式(8)得证。

为证式(9)成立，由式(8)的证明知仅须证

${\left(\frac{h\left(x\right)}{x}\right)}^{\left(n\right)}\ge 0,x>0,n=1,2,3.$

因为

${\left(\frac{h\left(x\right)}{x}\right)}^{\prime }=\frac{x{h}^{\prime }\left(x\right)-h\left(x\right)}{x^2}\equiv \frac{E_1\left(x\right)}{x^2},$

${\left(\frac{h\left(x\right)}{x}\right)}^{\prime \text{}\prime }=\frac{x^2{h}^{″}\left(x\right)-2x{h}^{\prime }\left(x\right)+2h\left(x\right)}{x^3}\equiv \frac{E_2\left(x\right)}{x^3},$

${\left(\frac{h\left(x\right)}{x}\right)}^{\left(3\right)}=\frac{x^3{h}^{\left(3\right)}\left(x\right)-3x^2{h}^{″}\left(x\right)+6x{h}^{\prime }\left(x\right)-6h\left(x\right)}{x^4}\equiv \frac{E_3\left(x\right)}{x^4}.$

又因为 $E_1\left(0\right)=0$ ， $E_1{}^{\prime }\left(x\right)=x{h}^{″}\left(x\right)\ge 0$ ，所以 $E_1\left(x\right)\ge 0$ 。同理可得

$E_2\left(0\right)=0,E_2{}^{\prime }\left(x\right)=x^2{h}^{\left(3\right)}\left(x\right)\ge 0,$

所以 $E_2\left(x\right)\ge 0$ 。同理有 $E_3\left(0\right)=0$ ， $E_3{}^{\prime }\left(x\right)=x^3{h}^{\left(4\right)}\left(x\right)\ge 0$ ，所以 $E_3\left(x\right)\ge 0$ 。从而式(9)得证。

将式(8)和式(9)代入式(7)得：

$\begin{array}{c}0>a\left(1-x_0^2\right)+b{x}_0^2\left(\left(\alpha -1\right)x_0^{\alpha -2}+\left(\alpha +1\right)x_0^{\alpha }\right)+cuv\left(\left(2\alpha -1\right)x_0^{2\alpha -2}+\left(2\alpha +1\right)x_0^{2\alpha }\right)\\ +x_0^2{\left(xh\left(x\right)\right)}^{\prime }{}_{x=x_0}+x_0^2{\left(\frac{h\left(x\right)}{x}\right)}^{\prime }{}_{x=x_0}\\ =-P.\end{array}$

但由已知条件有 $P\le \text{0}$ ，产生矛盾。从而式(6)代表的两曲线无交点，定理证毕。

定理3 如果P₀是不稳定的焦点或者结点，即 $\omega >0$ ， $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ，当 $P>0$ 时，且 $h\left(x\right)$ 满足条件：

$h\left(0\right)=0,h^{\left(n\right)}\left(x\right)\ge 0,n=1,2,3,4,x>0.$

则系统(4)在 $R_{+}$ 内存在唯一极限环。

证明 首先证明系统极限环的存在性。作直线 $L_1:x=x_1$ ，则当 $y>0$ 时，

${\frac{\text{d}x}{\text{d}t}|}_{L_1}=-x_1^{2}y<0,$

所以L₁是无切直线，系统(4)轨线自右而左穿过L₁。再作直线

$L_2:y+m\omega x-n=0,0<x<x_1,$

当n充分大时，则

${\frac{\text{d}{L}_2}{\text{d}t}|}_{\left(3\right)}={\left[\frac{\text{d}y}{\text{d}t}+m\omega \frac{\text{d}x}{\text{d}t}\right]|}_{y=n-m\omega x}=m\omega x\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right)-md\left(n-m\omega x\right)<0,$

即L₂是无切直线，系统(4)轨线自上而下穿过L₂。由于 $x=0$ ， $y=0$ 均是系统(3)的轨线，而P₀是不稳定的，据引理2可知P₀外围存在极限环。

下面证明系统极限环的唯一性，因为

$F\left(x\right)=\frac{1}{x}\left(1+x^2\right)\left(a-b{x}^{\alpha }-c{x}^{2\alpha }-h\left(x\right)\right),$

所以：

$\begin{array}{l}f\left(x\right)=F^{\prime }\left(x\right)=\frac{1}{x^2}[a{x}^2-a-b{x}^{\alpha +2}-\alpha b{x}^{\alpha +2}+b{x}^{\alpha }-\alpha b{x}^{\alpha }+c{x}^{2\alpha }-2\alpha c{x}^{2\alpha }-c{x}^{2\alpha +2}\\ -2\alpha c{x}^{2\alpha }-c{x}^{2\alpha +2}-2\alpha c{x}^{2\alpha +2}-x^3{h}^{\prime }\left(x\right)-x{h}^{\prime }\left(x\right)-x^{2}h\left(x\right)+h\left(x\right).\end{array}$

从而：

$\begin{array}{c}{\left(\frac{f}{g}\right)}^{\prime }m\frac{{\left(\omega x^2-d\right)}^2}{x}=\left(\omega x^2-d\right)[2a-{\alpha }^{2}b{x}^{\alpha }-2b{x}^{\alpha }-3\alpha b{x}^{\alpha }-{\alpha }^{2}b{x}^{\alpha -2}+\alpha b{x}^{\alpha -2}\\ -4{\alpha }^{2}c{x}^{2\alpha -2}-4{\alpha }^{2}c{x}^{2\alpha -2}+2\alpha c{x}^{2\alpha -2}-2c{x}^{2\alpha }-6\alpha c{x}^{2\alpha }+2\omega [a-a{x}^2\\ +2\omega [a-a{x}^2-b{x}^{\alpha }+b{x}^{\alpha +2}+\alpha b{x}^{\alpha }+\alpha b{x}^{\alpha +2}-c{x}^{2\alpha }+c{x}^{2\alpha +2}+2\alpha c{x}^{2\alpha }\\ +2\alpha c{x}^{2\alpha +2}+\left(\omega x^2-d\right)[-x^2{h}^{″}\left(x\right)-h^{″}\left(x\right)-4x{h}^{\prime }\left(x\right)-2h\left(x\right)]\\ +2\omega \left[x{h}^{\prime }\left(x\right)+x^3{h}^{\prime }\left(x\right)-h\left(x\right)+x^{2}h\left(x\right)\right]\\ \equiv {\Delta }_1+{\Delta }_2.\end{array}$

其中：

${\Delta }_2=\left(\omega x^2-d\right)[-x^2{h}^{″}\left(x\right)-h^{″}\left(x\right)-4x{h}^{\prime }\left(x\right)-2h\left(x\right)]+2\omega \left[x{h}^{\prime }\left(x\right)+x^3{h}^{\prime }\left(x\right)-h\left(x\right)+x^{2}h\left(x\right)\right],$

$\begin{array}{c}{\Delta }_1=\left(\omega x^2-d\right)[2a-{\alpha }^{2}b{x}^{\alpha }-2b{x}^{\alpha }-3\alpha b{x}^{\alpha }-{\alpha }^{2}b{x}^{\alpha -2}+\alpha b{x}^{\alpha -2}-4{\alpha }^{2}c{x}^{2\alpha -2}\\ -4{\alpha }^{2}c{x}^{2\alpha -2}+2\alpha c{x}^{2\alpha -2}-2c{x}^{2\alpha }-6\alpha c{x}^{2\alpha }]+2\omega [a-a{x}^2-b{x}^{\alpha }+b{x}^{\alpha +2}\\ +\alpha b{x}^{\alpha }+\alpha b{x}^{\alpha +2}-c{x}^{2\alpha }+c{x}^{2\alpha +2}+2\alpha c{x}^{2\alpha }+2\alpha c{x}^{2\alpha +2}\\ =2\omega \left(a-a\frac{d}{\omega }\right)+{\alpha }^{2}bd{x}^{\alpha }+2bd{x}^{\alpha }+3\alpha bd{x}^{\alpha }+{\alpha }^{2}bd{x}^{\alpha -2}-\alpha bd{x}^{\alpha -2}-{\alpha }^{2}b\omega x^{\alpha }\\ -{\alpha }^{2}b\omega x^{\alpha +2}-2b\omega x^{\alpha }+3\alpha b\omega x^{\alpha }-\alpha b\omega x^{\alpha +2}+4{\alpha }^{2}cd{x}^{2\alpha -2}+4{\alpha }^{2}cd{x}^{2\alpha }-2\alpha cd{x}^{2\alpha -2}\\ +2cd{x}^{2\alpha }+6\alpha cd{x}^{2\alpha }-4{\alpha }^{2}c\omega x^{2\alpha }-4{\alpha }^{2}c\omega x^{2\alpha +2}-2c\omega x^{2\alpha }+6\alpha c\omega x^{2\alpha }-2\alpha c\omega x^{2\alpha +2},\end{array}$

由 $P>0$ 及 $x_0^2=\frac{d}{\omega }$ ，有

$\begin{array}{c}{\Delta }_1<2\omega \{-b\left[\alpha -1+\left(\alpha +1\right)x_0^2\right]x_0^{\alpha }-c\left[2\alpha -1+\left(2\alpha +1\right)x_0^2\right]x_0^{2\alpha }-\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)\\ -\left(x_0^2-1\right)h\left(x_0\right)\}+{\alpha }^{2}bd{x}^{\alpha }+2bd{x}^{\alpha }+3\alpha bd{x}^{\alpha }+{\alpha }^{2}bd{x}^{\alpha -2}-\alpha bd{x}^{\alpha -2}-{\alpha }^{2}b\omega x^{\alpha }\\ -{\alpha }^{2}b\omega x^{\alpha +2}-2b\omega x^{\alpha }+3\alpha b\omega x^{\alpha }-\alpha b\omega x^{\alpha +2}+4{\alpha }^{2}cd{x}^{2\alpha -2}+4{\alpha }^{2}cd{x}^{2\alpha }-2\alpha cd{x}^{2\alpha -2}\\ +2cd{x}^{2\alpha }+6\alpha cd{x}^{2\alpha }-4{\alpha }^{2}c\omega x^{2\alpha }-4{\alpha }^{2}c\omega x^{2\alpha +2}-2c\omega x^{2\alpha }+6\alpha c\omega x^{2\alpha }-2\alpha c\omega x^{2\alpha +2}\\ =S_1\left(x\right)+S_2\left(x\right)+2\omega \left[-\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)-\left(x_0^2-1\right)h\left(x_0\right)\right],\end{array}$

其中：

$\begin{array}{c}{S}_2\left(x\right)=2\omega \left\{-b\left[\alpha -1+\left(\alpha +1\right)x_0^2\right]x_0^{\alpha }\right\}+c[4{\alpha }^{2}d{x}^{2\alpha -2}+4{\alpha }^{2}d{x}^{2\alpha }-2\alpha d{x}^{2\alpha -2}\\ -\alpha d{x}^{\alpha -2}-{\alpha }^2\omega x^{\alpha }-{\alpha }^2\omega x^{\alpha +2}-2\omega x^{\alpha }+3\alpha \omega x^{\alpha }-\alpha \omega x^{\alpha +2}],\end{array}$

$S_1\left(x\right)=2\omega s_1\left(x\right),$

其中：

$\begin{array}{c}{s}_1\left(x\right)=2\alpha x_0^{2\alpha }-x_0^{2\alpha }+2\alpha x_0^{2\alpha +2}+x_0^{2\alpha +2}-2{\alpha }^2{x}_0^2{x}^{2\alpha -2}-2{\alpha }^2{x}_0^2{x}^{2\alpha }+\alpha x_0^2{x}^{2\alpha -2}\\ -x_0^2{x}^{2\alpha }-3\alpha x_0^2{x}^{2\alpha }+2{\alpha }^2{x}^{2\alpha }+2{\alpha }^2{x}^{2\alpha +2}+\alpha x^{2\alpha +2}+x^{2\alpha }-3\alpha x^{2\alpha }.\end{array}$

经计算可得 $s_1\left(x_0\right)=0$ ， ${s^{\prime }}_1\left(x_0\right)=0$ ，且 $\alpha \ge 2$ 知当 $x>x_0$ 时， ${s^{\prime }}_1\left(x\right)\ge 0$ ；当 $0<x<x_0$ 时， ${s^{\prime }}_1\left(x\right)\le 0$ ，因此

$s_1\left(x\right)\ge s_1\left(x_0\right)=0,S_1\left(x\right)\le S_1\left(x_0\right)=0,$

同理可得 $S_2\left(x\right)\le S_2\left(x_0\right)=0$ 。即

${\Delta }_1\le 2\omega \left[-\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)-\left(x_0^2-1\right)h\left(x_0\right)\right].$

因此

$\begin{array}{c}{\left(\frac{f}{g}\right)}^{\prime }m\frac{{\left(\omega x^2-d\right)}^2}{x}<2\omega \left[-\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)-\left(x_0^2-1\right)h\left(x_0\right)\right]\\ +\left(\omega x^2-d\right)[-x^2{h}^{″}\left(x\right)-h^{″}\left(x\right)-4x{h}^{\prime }\left(x\right)-2h\left(x\right)]\\ +2\omega \left[x{h}^{\prime }\left(x\right)+x^3{h}^{\prime }\left(x\right)-h\left(x\right)+x^{2}h\left(x\right)\right].\end{array}$

令

$\begin{array}{c}W\left(x\right)=2\omega \left[-\left(x_0^2+1\right)x_0{h}^{\prime }\left(x_0\right)-\left(x_0^2-1\right)h\left(x_0\right)\right]+\left(\omega x^2-d\right)[-x^2{h}^{″}\left(x\right)-h^{″}\left(x\right)-4x{h}^{\prime }\left(x\right)-2h\left(x\right)]\\ +2\omega \left[x{h}^{\prime }\left(x\right)+x^3{h}^{\prime }\left(x\right)-h\left(x\right)+x^{2}h\left(x\right)\right].\end{array}$

显然有 $W\left(x_0\right)=0$ ，及

$W^{\prime }\left(x\right)\text{=}\omega \left(x^2-x_0^2\right)\left[-x^2{h}^{\left(3\right)}\left(x\right)-h^{\left(3\right)}\left(x\right)-6h^{\prime }\left(x\right)-6x{h}^{″}\left(x\right)\right].$

由于 $h^{\left(n\right)}\left(x\right)\ge 0$ ， $n=1,2,3$ ， $x>0$ 。所以当 $x>x_0$ 时， $W^{\prime }\left(x\right)\le 0$ ；当 $0<x<x_0$ 时， $W^{\prime }\left(x\right)\ge 0$ ，从而当 $0<x<+\infty$ 时，有

${\left(\frac{f}{g}\right)}^{\prime }m\frac{{\left(\omega x^2-d\right)}^2}{x}<W\left(x\right)<W\left(x_0\right)=0,$

所以 ${\left(\frac{f}{g}\right)}^{\prime }<0$ ，由引理2知唯一性得证，证毕。

通过对系统的稳定点和极限环的定性分析可以推出一些特殊形式的食饵捕食者系统如下，同样可以得出类似的结论。

推论1 系统

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(a-b{x}^{\alpha }-c{x}^{2\alpha }\right)-\frac{k{x}^2{y}^m}{1+x^2}\\ \frac{\text{d}y}{\text{d}t}=y\left(-d+\frac{ek{x}^2}{1+x^2}\right)\end{array}$

其中 $a,b,c,d,e,k,\alpha ,m$ 所有的系数均为正数，且 $0<m\le 1$ ， $\alpha \ge 2$ 时，系统在第一象限存在极限环的充要条件是该系统在第一象限存在唯一的不稳定奇点。显然当极限环存在时则必唯一。

推论2 系统

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(c_0-c_1{x}^{\alpha }-\cdots -c_n{x}^{n\alpha }\right)-\frac{k{x}^2{y}^m}{1+x^2}\\ \frac{\text{d}y}{\text{d}t}=y\left(-d+\frac{ek{x}^2}{1+x^2}\right)\end{array}$

其中 $a,b,c,d,e,k,\alpha ,m$ 所有的系数均为正数，且 $0<m\le 1$ ， $\alpha \ge 2$ ， $n\ge 3$ 时，系统在第一象限存在极限环的充要条件是该系统在第一象限存在唯一的不稳定奇点。显然当极限环存在时则必唯一。

4. 数值仿真

下面分四种情况进行数值仿真，以验证所得结果的正确性。

第1种情况：取 $a=6$ ， $b=0.5$ ， $c=0.5$ ， $d=2$ ， $e=0.5$ ， $m=1$ ， $k=2$ ， $\alpha =4$ ， $h\left(x\right)=0$ ，于是有 $\omega =-1<0$ ，此时系统(3)从第一象限出发的初值不同的轨线均无限趋近于边界平衡点 $E\left(\sqrt[4]{3},0\right)$ (如图1所示)。

Figure 1. Boundary equilibrium point E asymptotic stable phase diagram

图1. 边界平衡点E点渐近稳定相图

第2种情况：取 $a=10$ ， $b=\frac{1}{4}$ ， $c=\frac{1}{16}$ ， $d=2$ ， $e=1$ ， $m=1$ ， $k=3$ ， $\alpha =4$ ， $h\left(x\right)=0$ ，于是有 $\omega =1>0$ ，

$x_0=\sqrt{2}$ ，且 $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ， $P\le 0$ ，此时系统(3)从第一象限出发的初值不同的轨线均无限趋近于正平衡点 $P_0\left(\sqrt{2},12\sqrt{2}\right)$ (如图2所示)。

第3种情况：取 $a=16$ ， $b=0.25$ ， $c=0.25$ ， $d=1$ ， $e=0.5$ ， $m=1$ ， $k=3$ ， $\alpha =2$ ， $h\left(x\right)=0$ ，于是有 $\omega =ek-d=0.5>0$ ， $x_0=\sqrt{2}$ ，且 $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ， $P>0$ ，即使取不同的初值，系统(3)始终在

正平衡点 $P_0\left(\sqrt{2},\frac{87\sqrt{2}}{4}\right)$ 的外围存在唯一稳定的极限环(如图3所示)。

第4种情况：取 $a=30$ ， $b=\frac{1}{4}$ ， $c=\frac{1}{8}$ ， $d=1$ ， $e=0.5$ ， $m=1$ ， $k=3$ ， $\alpha =2$ ， $h\left(x\right)=\frac{1}{8}{x}^{3\alpha }$ ，于是

有 $\omega =0.5>0$ ， $x_0=\sqrt{2}$ ，且 $a>b{x}_0^{\alpha }+c{x}_0^{2\alpha }+h\left(x_0\right)$ ， $P>0$ ，即使取不同的初值，系统(3)在正平衡点 $P_0\left(\sqrt{2},42\sqrt{2}\right)$ 的外围存在唯一稳定的极限环(如图4所示)。

Figure 2. Positive equilibrium point P₀ asymptotic stable phase diagram

图2. 正平衡点P₀点渐近稳定相图

Figure 3. Positive equilibrium point P₀ unique phase diagram of the outer limit ring of the point

图3. 正平衡点P₀点外围极限环唯一相图

Figure 4. Unique phase diagram of the limit ring outside the positive equilibrium point P₀

图4. 正平衡点P₀点外围极限环唯一相图

5. 结语

基于生态环境的复杂性以及捕食者种群在捕食过程中相互之间存在一定的干扰，该模型考虑的因素相对来说要更加全面，得到的结果也更加的会符合实际捕食–食饵种群的数量变化情况。通过对生态模型(1)在平衡点处线性近似系统的分析，首先得出该模型在一定的条件下存在唯一的正平衡点，当模型(1)的正平衡点存在时，利用Bendixson-Dulac定理分析得出系统存在全部位于第一象限的闭轨线；然后在Liénard变换的基础上，再利用张芷芬唯一性定理得到该系统极限环的存在性及唯一性；最后利用Mathematica软件对模型的不同情况进行了四组具有代表性的数值仿真，数值模拟结果表明所得定性结论的有效性。同时，数值模拟所得图形表明，理论分析较好地反映了该类具相互干扰和Holling-III功能反应模型在食饵种群内部增长率、捕食者种群内的变换系数和死亡率满足一定条件时，食饵和捕食者的比例最后会趋于一个常量。

基金项目

湖南省大学生创新训练项目(编号：202211535022)；湖南省教育厅科研课题资助项目(编号：19C0556)。