1. 模型的组成

近年来，关于捕食–食饵模型的动力学分析作为一个重要课题，吸引了许多学者的广泛研究。1932年，Allee通过金鱼种群增长的实验研究中发现，金鱼密度较大时有利于金鱼种群的增长，提出当种群的密度减少到一定量时，将会保持在一个很低的水平并趋于灭绝，即Allee效应[1] [2]。Allee效应主要分为两类，分别为弱Allee效应和强Allee效应[3]-[7]。人们普遍认为Allee效应可能会增加低密度种群的灭绝风险。因此，Allee效应的种群生态学调查对保护生物学很重要。如果用$x\left(t\right)$ 表示食饵种群的数量；$y\left(t\right)$ 表示捕食者种群的数量，食饵具有Allee效应的增长模型表示为：

$\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{k}\right)\left(x-v\right)$

r表示食饵的自然增长率，k表示环境容纳量，v是表示Allee效应的阈值总体水平，且$0<v<k$ 。

此外学者们提出多种功能反应，其中Holling类型I，II，III被广泛讨论[8]-[10]。

Holling I：

$f\left(x\right)=\{\begin{array}{l}\frac{b}{a}x,0<x<a\\ b,x>a\end{array}$

Holling II：

$f\left(x\right)=\frac{\alpha x}{\gamma +x}$

Holling III：

$f\left(x\right)=\frac{\alpha x^2}{\gamma +x^2}$

在现实生活中，食饵种群大部分是有防御能力的，随着食饵种群数量的不断增加，其防御能力也会增强，对捕食者就会起到相对抑制的作用。因此，根据这类现象，Andrews提出了Holling IV型功能性反应函数。

$f\left(x,y\right)=\frac{mxy}{\alpha +\beta x+x^2}$

最近几年，Cui等[11]考虑了具有强Allee效应和的Holling II型功能反应函数的扩散捕食者–食饵系统。但是对于研究Allee效应和Holling IV功能反应的扩散捕食者–食饵模型比较少见。因此本文研究具有Allee效应的Holling IV功能反应的捕食者–食饵模型，具体模型为虑相互作用的物种食饵和捕食者的系统，$x\left(t\right)$ 为食饵种群的数量；$y\left(t\right)$ 为捕食者种群的数量。

2. 解的性质

首先我们证明解的正性。事实上，设

${\phi }_1\left(x,y\right)=r\left(1-\frac{x}{k}\right)\left(x-v\right)-\frac{my}{\alpha +\beta x+x^2}$ ，${\phi }_2\left(x,y\right)=\frac{emx}{\alpha +\beta x+x^2}-d$

由模型(1.1)可得：

$\frac{\text{d}x}{\text{d}t}=x{\phi }_1\left(x,y\right)$ ，$\frac{\text{d}y}{\text{d}t}=y{\phi }_2\left(x,y\right)$

因此：

$x\left(t\right)=x\left(0\right){\text{e}}^{{\displaystyle \underset{0}{\overset{t}{\int }}{\phi }_1\left(x\left(s\right),y\left(s\right)\right)\text{d}s}}>0$ ，$y\left(t\right)=y\left(0\right){\text{e}}^{{\displaystyle \underset{0}{\overset{t}{\int }}{\phi }_2\left(x\left(s\right),y\left(s\right)\right)\text{d}s}}>0$ .

这意味着解是正的。

引理2.2 若系统初值$\left(x_0,y_0\right)\in R_{+}^2$ ，且满足初值条件的解都是有界的。

证：构造一个关于x和y的函数$w\left(t\right)$ ：

$w\left(t\right)=ex\left(t\right)+y\left(t\right)$

所以

$\begin{array}{c}\frac{\text{d}w\left(t\right)}{\text{d}t}=er\left(1-\frac{x}{k}\right)\left(x-v\right)-dy\\ =ervx\left(1-\frac{x}{k}\right)\left(\frac{x}{v}-1\right)-dy\\ =-ervx-dy+\frac{er}{k}{x}^2\left(-x+k+v\right)\end{array}$

令

$\wp =\frac{er}{k}{x}^2\left(-x+k+v\right)$ ，

经过计算可得$\wp$ 在$x=\left(2k+2v\right)/3$ 处取得最大值M，则可得到：

$\frac{\text{d}w\left(t\right)}{\text{d}t}\le M-ervx-dy\le M-\delta w$

其中：

$\delta =\min \left\{rv,d\right\}$ ，$M=\frac{4}{27}\frac{ev}{k}{\left(k+v\right)}^3$ .

根据微分比较定理可得：

$0\le w\left(x\left(t\right),y\left(t\right)\right)\le \frac{\delta \left(1-{\text{e}}^{-Mt}\right)}{M}+w\left(x\left(0\right),y\left(0\right)\right){\text{e}}^{-Mt}$

当$t\to +\infty$ ，有$0\le w\left(t\right)\le \delta /M$ 。

所以系统(1.1)的解是有界的。

3. 模型的动力学分析

本节先分析平衡点的存在性，再分析稳定性：

1) 显然系统存在平衡点$E_0=\left(0,0\right)$ ；

2) 当$x\ne 0,y=0$ 时，令

$rx\left(1-\frac{x}{k}\right)\left(x-v\right)-\frac{mxy}{\alpha +\beta x+x^2}=0$

可得：

$rx\left(1-\frac{x}{k}\right)\left(x-v\right)=0$

可得$x=k$ 或$x=v$ ，此时系统存在两个平衡点$E_1=\left(k,0\right)$ ，$E_2=\left(v,0\right)$ 。

3) 当$x\ne 0,y\ne 0$ 时，令

$\frac{emxy}{\alpha +\beta x+x^2}-dy=0$

可得：

$x_i^{*}=\frac{em-\beta d\pm \sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d}$

再由

$\frac{emxy}{\alpha +\beta x+x^2}-dy=0$

可得：

$\frac{mxy}{\alpha +\beta x+x^2}=\frac{dy}{e}$

将其代入

$rx\left(1-\frac{x}{k}\right)\left(x-v\right)-\frac{mxy}{\alpha +\beta x+x^2}=0$

可得：

$rx\left(1-\frac{x}{k}\right)\left(x-v\right)=\frac{dy}{e}$

$y_i^{*}=\frac{e}{d}r{x}_i^{*}\left(1-\frac{x_i^{*}}{k}\right)\left(x_i^{*}-v\right),i=1,2$ ，

系统(1.1)的Jacobian矩阵为：

$J=\left(\begin{array}{cc}r\left(1-\frac{x}{k}\right)\left(x-v\right)-\frac{rx\left(x-v\right)}{k}+rx\left(1-\frac{x}{k}\right)+\frac{m{\left(\alpha -x\right)}^{2}y}{{\left(\alpha +\beta x+x^2\right)}^2}& -\frac{mx}{\alpha +\beta x+x^2}\\ \frac{em{\left(\alpha -x\right)}^{2}y}{{\left(\alpha +\beta x+x^2\right)}^2}& \frac{emx}{\alpha +\beta x+x^2}-d\end{array}\right)$

于是通过计算各平衡点对应的Jacobian矩阵来分析局部稳定性。

3.1. $E_0=\left(0,0\right)$ 的稳定性

设J₀是系统(1.1)在平衡点E₀处的Jacobian矩阵，所以：

$J_0=\left(\begin{array}{cc}-rv& 0\\ 0& -d\end{array}\right)$

所以两特征根分别为${\lambda }_1=-rv,{\lambda }_2=-d$ ，系统在E₀是稳定的，否则系统(1.1)是不稳定的。

3.2. $E_1=\left(k,0\right)$ 的稳定性

设J₁是系统(1.1)在平衡点E₁处的Jacobian矩阵，所以

$J_1=\left(\begin{array}{cc}-r\left(k-v\right)& \frac{mk}{\alpha +\beta x+x^2}\\ 0& \frac{emk}{\alpha +\beta k+k^2}-d\end{array}\right)$

所以两特征根分别为：

${\lambda }_1=-r\left(k-v\right),{\lambda }_2=\frac{emk}{\alpha +\beta k+k^2}-d$ ，

如果

$d>\frac{emk}{\alpha +\beta k+k^2}$

系统在E₁是稳定的，否则系统(1.1)是不稳定的。

3.3. $E_2=\left(v,0\right)$ 的稳定性

设J₂是系统(1.1)在平衡点E₂处的Jacobian矩阵，所以

$J_2=\left(\begin{array}{cc}rv\left(1-\frac{v}{k}\right)& \frac{mv}{\alpha +\beta x+x^2}\\ 0& \frac{mv}{\alpha +\beta v+v^2}-d\end{array}\right)$

所以两特征根分别为：

${\lambda }_1=rv\left(1-\frac{v}{k}\right),{\lambda }_2=\frac{emv}{\alpha +\beta v+v^2}-d$

由于${\lambda }_1$ 为正，所以根据Routh-Hurwitz的稳定性条件可得，E₂是不稳定的。

3.4. $E_3=\left(x^{\ast },y^{\ast }\right)$ 的稳定性

设J₃是系统(1.1)在平衡点E₃处的Jacobian矩阵，所以

$J_3=\left(\begin{array}{cc}r\left(1-\frac{x^{*}}{k}\right)\left(x^{*}-v\right)-\frac{r{x}^{*}(x^{*}-v)}{k}+r{x}^{*}\left(1-\frac{x^{*}}{k}\right)+\frac{m{\left(\alpha -x^{*}\right)}^2{y}^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}& -\frac{m{x}^{*}}{\alpha +\beta x^{*}+x^{*}{}^2}\\ \frac{em{\left(\alpha -x^{*}\right)}^2{y}^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}& \frac{em{x}^{*}}{\alpha +\beta x^{*}+x^{*}{}^2}-d\end{array}\right)$

当参数H₁满足条件：

$r\left(1-\frac{x^{*}}{k}\right)\left(x^{*}-v\right)-\frac{r{x}^{*}\left(x^{*}-v\right)}{k}+r{x}^{*}\left(1-\frac{x^{*}}{k}\right)+\frac{m{\left(\alpha -x^{*}\right)}^2{y}^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}<0$

时，正平衡点$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 是局部渐近稳定的。

设$A=r\left(1-\frac{x_1^{*}}{k}\right)\left(x_1^{*}-v\right)-\frac{r{x}_1^{*}\left(x_1^{*}-v\right)}{k}+r{x}_1^{*}\left(1-\frac{x_1^{*}}{k}\right)+\frac{m{\left(\alpha -x_1^{*}\right)}^2{y}_1^{*}}{{\left(\alpha +\beta x_1^{*}+x_1^{*}{}^2\right)}^2}$

$B=-\frac{m{x}_1^{*}}{\alpha +\beta x_1^{*}+x_1^{*}{}^2}$ ，$C=\frac{em{\left(\alpha -x_1^{*}\right)}^2{y}_1^{*}}{{\left(\alpha +\beta x_1^{*}+x_1^{*}{}^2\right)}^2}$

$D=0$ 所以特征方程为${\lambda }^2-A\lambda +BC=0$ 。

所以根据根与系数的关系，可得方程的两个特征根${\lambda }_1,{\lambda }_2$ 满足：

${\lambda }_1+{\lambda }_2=A,{\lambda }_1{\lambda }_2=BC=\frac{e{m}^2{\left(\alpha -x_1^{*}\right)}^2{x}_1^{*}{y}_1^{*}}{{\left(\alpha +\beta x_1^{*}+x_1^{*}{}^2\right)}^3}$ (1.2)

令

$x_1^{*}=\frac{em-\beta d-\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d}$

$x_2^{*}=\frac{em-\beta d+\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d}$

所以

$\begin{array}{c}\alpha -x_1^{*2}=\alpha -{\left[\frac{em-\beta d-\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d}\right]}^2\\ =\alpha -\frac{{\left(em-\beta d\right)}^2+\left(e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha \right)-2\left(em-\beta d\right)\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d}\\ =\frac{8\alpha d^2-2{\left(em-\beta d\right)}^2+2\left(em-\beta d\right)\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{4d^2}\\ =\frac{4\alpha d^2-{\left(em-\beta d\right)}^2+\left(em-\beta d\right)\sqrt{{\left(em-\beta d\right)}^2-4d^2\alpha }}{2d^2}>0\end{array}$

因为

$\left(em-\beta d\right)>\sqrt{{\left(em-\beta d\right)}^2-4d^2\alpha }$

所以

$\begin{array}{l}\frac{4\alpha d^2-{\left(em-\beta d\right)}^2+\left(em-\beta d\right)\sqrt{e^2{m}^2-2e\beta dm+{\beta }^2{d}^2-4d^2\alpha }}{2d^2}\\ >\frac{4\alpha d^2-{\left(em-\beta d\right)}^2+{\left(em-\beta d\right)}^2-4\alpha d^2}{2d^2}=0\end{array}$

所以${\lambda }_1{\lambda }_2>0$ ，因此若${\lambda }_1+{\lambda }_2<0$ ，两个特征根具有负实根，此时$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 是局部渐近稳定的

同理对于平衡点$E_{32}=\left(x_2^{*},y_2^{*}\right)$ 可得${\lambda }_1{\lambda }_2<0$ 所以不稳定。

3.5. Hopf分支

由3.4可知，系统(1.1)正平衡点$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 的局部稳定性需要满足一定的条件。因此，本小节主要分析系统(1.1)在平衡点$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 发生在Hopf分支的条件和Hopf分支的方向。

记${\eta }_1={\lambda }_1+{\lambda }_2$ ，${\eta }_2={\lambda }_1{\lambda }_2$ 其中${\lambda }_1$ 和${\lambda }_2$ 是(1.2)中定义的，且${\eta }_2>0$ 选取捕食者自然死亡率d作为分支参数。如果存在$d=d^{*}>0$ ，使得H₂：${\eta }_1\left(d^{*}\right)=0$ ，${\eta }_{12}\left(d^{*}\right)<0$ ，$2{\eta }_{12}{\eta }_2+{\eta }_1{\eta }_{21}\ne 0$ 成立，其中${\eta }_{12}$ ${\eta }_{21}$ 表示对d的求导。

假设条件H₂成立，当$d=d^{*}$ ，系统会在正平衡点$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 处发生Hopf分支。

证：因为正平衡点E₃₁处的Jacobian矩阵对应的特征方程如下：

${\lambda }^2-{\eta }_1\lambda +{\eta }_2=0$ (1.3)

当$d=d^{*}$ 时，即${\eta }_1=0$ 特征方程可以为：

${\lambda }^2+{\eta }_2=0$

此时方程会出现一对纯虚根为：

${\lambda }_1=i\sqrt{{\eta }_2}$ ，${\lambda }_2=-i\sqrt{{\eta }_2}$

现在验证横截条件，对特征方程(1.3)关于d进行求导，有：

$2\lambda \frac{\text{d}\lambda }{\text{d}d}-\lambda {\eta }_{12}-\frac{\text{d}\lambda }{\text{d}d}{\eta }_1+{\eta }_{21}=0$

整理可得：

$\frac{\text{d}\lambda }{\text{d}d}=\frac{\lambda {\eta }_{12}-{\eta }_{21}}{2\lambda -{\eta }_1}$

则有：

$\frac{\text{d}\lambda }{\text{d}d}=\frac{2\lambda {\eta }_{12}{\eta }_2-{\eta }_1{\eta }_{21}}{4{\eta }_2+{\eta }_1{}^2}+i\frac{2{\eta }_{21}\sqrt{{\eta }_2}-{\eta }_1\sqrt{{\eta }_2}{\eta }_{21}}{4{\eta }_2+{\eta }_1{}^2}$ 当$\lambda =i\sqrt{{\eta }_2}$

因此此时系统(1.1)会在正平衡点$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 处发生Hopf分支。

下面通过计算第一系数Lyapunov系数来分析Hopf分支的方向。

首先做线性替换所以令$u=x-x^{*}$ ，$v=y-y^{*}$ 对系统(1.1)泰勒展开可得：

$\frac{\text{d}u}{\text{d}t}=a_{10}u+a_{01}v+a_{11}uv+a_{20}{u}^2+a_{02}{v}^2+a_{30}{u}^3+a_{21}{v}^3+a_{12}{u}^{2}v+a_{03}u{v}^2+P\left(u,v\right)$

$\frac{\text{d}v}{\text{d}t}=b_{10}u+b_{01}v+b_{11}uv+b_{20}{u}^2+b_{02}{v}^2+b_{30}{u}^3+b_{21}{v}^3+b_{12}{u}^{2}v+b_{03}u{v}^2+Q\left(u,v\right)$

$a_{10}=r\left(1-\frac{x^{*}}{k}\right)\left(x^{*}-v\right)-\frac{r{x}^{*}\left(x^{*}-v\right)}{k}+r{x}^{*}\left(1-\frac{x^{*}}{k}\right)+\frac{m{\left(\alpha -x^{*}\right)}^2{y}^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}$

$a_{01}=-\frac{m{x}^{*}}{\alpha +\beta x^{*}+x^{*}{}^2}$ ，$a_{11}=\frac{m\left(\alpha -x^{*}{}^2\right)}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}$

$a_{20}=-\frac{r\left(x^{*}-v\right)}{k}+r\left(1-\frac{x^{*}}{k}\right)-\frac{r{x}^{*}}{k}-\frac{m\left(\alpha -x^{*}\right)y^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}-\frac{m\left(\beta +2x^{*}\right)\left(\alpha -x^{*2}\right)y^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}$ ，$a_{02}=0$

$\begin{array}{c}{a}_{30}=-\frac{r}{k}+\frac{m{y}^{*}}{3{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}+\frac{2m\left(\beta +2x^{*}\right)\left(\alpha -x^{*2}\right)y^{*}}{3{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^3}-\frac{2m\left(\alpha -x^{*2}\right)y^{*}}{3{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^3}\\ +\frac{2m\left(\beta +2x^{*}\right)x^{*}{y}^{*}}{3{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^3}+\frac{m{\left(\beta +2x^{*}\right)}^2\left(\alpha -x^{*2}\right)y^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^4}\end{array}$

$a_{03}=0$ ，$a_{21}=-\frac{m\left(\alpha -x^{*}\right)}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}+\frac{m\left(\beta +2x^{*}\right)\left(\alpha -x^{*2}\right)}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^3}$ ，$a_{12}=0$

$b_{10}=\frac{em{\left(\alpha -x^{*}\right)}^2{y}^{*}}{{\left(\alpha +\beta x^{*}+x^{*}{}^2\right)}^2}=e{y}^{*}{a}_{11}$ ，$b_{01}=0$ ，$b_{11}=e{a}_{11}$ ，$b_{20}=e{y}^{*}{a}_{21}$ ，$b_{02}=0$ ，

$b_{30}=e\left(\frac{r}{k}+a_{30}\right)$ ，$b_{03}=0$ ，$b_{21}=e{a}_{21}$ ，$b_{12}=0$

$P\left(u,v\right)={\displaystyle \sum _{i+j=4}^{+\infty }{a}_{ij}{u}^i{v}^j}$ ，$Q\left(u,v\right)={\displaystyle \sum _{i+j=4}^{+\infty }{b}_{ij}{u}^i{v}^j}$

则由第一Lyapunov系数$\sigma$ 计算可得：

$\begin{array}{c}\sigma =-\frac{3\pi }{2a_{01}{\Delta }^{\frac{3}{2}}}\{[a_{10}{b}_{10}\left(a_{11}{b}_{02}+a_{02}{b}_{11}+a_{11}^2\right)+a_{10}{a}_{01}\left(a_{20}{b}_{11}+a_{11}{b}_{02}+b_{11}^2\right)+b_{10}^2\left(a_{11}{a}_{02}+2a_{02}{b}_{02}\right)\\ -2a_{10}{b}_{10}\left(b_{02}^2-a_{20}{a}_{02}\right)-2a_{10}{a}_{01}\left(a_{20}^2-b_{20}{b}_{02}\right)-a_{01}^2\left(2a_{20}{b}_{20}+b_{11}{b}_{20}\right)+\left(a_{10}{b}_{01}-2a_{10}^2\right)\left(b_{11}{b}_{02}-a_{11}{a}_{20}\right)]\\ -\left(a_{10}^2+a_{01}{b}_{10}\right)\left[3\left(b_{10}{b}_{03}-a_{01}{a}_{30}\right)+2a_{10}\left(a_{21}+b_{12}\right)+\left(b_{10}{a}_{21}-a_{01}{b}_{21}\right)\right]\}\end{array}$

其中$\Delta =a_{10}{b}_{01}-a_{01}{b}_{10}$ 。

若$\sigma >0$ 时：$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 产生次临界Hopf分支；

当$\sigma <0$ 时：$E_{31}=\left(x_1^{*},y_1^{*}\right)$ 产生超临界Hopf分支。

4. 结语

本文研究了具有Allee效应和具有Holling IV功能反应的捕食者–食饵模型。对其有界性、解的性质、平衡点的存在性以及其局部稳定性和分支进行了分析。通过构造函数得出了解的有界性。最后对模型线性化处理，通过计算Lyapunov系数来分析Hopf分支的方向。

参考文献