一类食饵具有常数收获率和Holling第II类功能性反应的捕食者–食饵模型的分支分析

Bifurcation Analysis of a Predator-Prey Model with Constant Prey Harvesting Rate and Holling Type II Functional Response

摘要: 本文研究一类食饵种群具有常数收获率且有Holling第Ⅱ类功能性反应的捕食者–食饵模型的平衡点和分支问题。首先给出系统模型平衡点的存在条件，然后主要讨论了平衡点的类型及稳定性和正平衡点的Hopf分支，且得出了产生Hopf分支的条件；最后，对该模型做了数值仿真模拟实验得到了相关结论。

Abstract: In this paper, we study the equilibria and the bifurcation of a predator-prey model with constant prey harvesting rate and Holling type II functional response. Firstly, we analyze the existence condition of the equilibria of the system model. Then, we discuss the type and stability of the equilibria, Hopf bifurcation near positive equilibria and we obtain the conditions of Hopf bifurcation. Finally, the numerical simulation of the model is carried out and the relevant conclusions are obtained.

1. 引言

食饵捕食者之间的动态关系一直是种群动力学模型的一个重要研究课题，众多学者对食饵捕食者系统进行研究并得到了大量的结论 [1] - [9]。人们在研究生物种群时着重关注种群的演变规律和对如何采取措施对种群进行合理的开发利用和保护的问题进行研究。针对这些问题的研究，除了可以对种群的变化和发展进行探索分析和加以估计，学者们还可以建立具有食饵收获率的捕食者–食饵模型，然后对该模型进行分析，得出需要多少的食饵收获率才能使该捕食者食饵模型系统达到一个动态平衡又可以满足人们的需求。因此研究具有收获率的捕食者食饵模型对生态保护和可持续发展具有重大意义，而众多学者针对此类模型已进行了广泛的研究。例如，文 [5] 研究了一类有常数收获率的Lotka-Volterra模型，着重分析了系统模型平衡点的存在性和稳定性以及极限环的存在性。文 [6] 研究了食饵带有收获率参数的Holling II型捕食者–食饵模型的Bogdanov-Takens分支，证明了该模型的Bogdanov-Takens分支是退化的。文 [7] 用线性化方法，分析了具齐次Neumann 边界条件的Holling第二类功能性反应捕食者–食饵模型的稳定性和Hopf分支，但此模型没有考虑捕食者对食饵收获率的影响，受到对具有常数收获率的捕食者食饵模型 [2] [3] [5] [6] 的启发，我们在本文当中也引入捕食者对食饵的常数收获率，即考虑食饵具有常数收获率和Holling第二类功能性反应的捕食者–食饵模型。首先给出系统模型平衡点的存在条件，然后主要讨论了平衡点的类型及稳定性和正平衡点的Hopf分支，且得出了产生Hopf分支的条件；最后，对该模型做了数值仿真模拟实验得到了相关结论。

相应的常微分模型如下：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{k}\right)-\frac{emxy}{a+x}-h,\\ \frac{\text{d}y}{\text{d}t}=y\left(\frac{mx}{a+x}-d\right).\end{array}$ (1)

其中：x为食饵种群的密度，y为捕食者种群的密度； $r,a,e,m,k,d,h$ 都是常数且为正；r表示猎物的内增

长率；k表示环境承载力；e表示相对损失；d表示捕食者的死亡率；h为捕食者对食饵的常数收获率。 $$

表示捕食者的Holling II能反应。

首先对模型(1)作如下的变换

$u=x,v=\frac{em}{r}y,\text{d}\tau =r\text{d}t,$

则可将模型(1)变为

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}\tau }=u\left(1-\frac{u}{k}\right)-\frac{uv}{a+u}-h_0,\\ \frac{\text{d}v}{\text{d}\tau }=v\left(\frac{bu}{a+u}-c\right).\end{array}$ (2)

其中： $b=\frac{m}{r}$， $c=\frac{d}{r}$， $h_0=\frac{h}{r}$。

2. 平衡点类型及正平衡点的稳定性

考虑到模型(1)的生物意义，所以我们之后的讨论都在 $u\ge 0,v\ge 0$ 上进行。首先我们先考虑平衡点的存在性。令系统(2)右端为0，可得：

1) 若 $v=0$，则有

$u^2-ku+k{h}_0=0.$ (3)

记 $\Delta =k^2-4k{h}_0=k\left(k-4h_0\right)$。

a) 当 $k=4h_0$ 时，则系统(2)有唯一边界平衡点 $E_0=\left(\frac{k}{2},0\right)$。

b) 当 $k>4h_0$ 时，系统(2)有两个边界平衡点 $E_1=\left(\frac{k+\sqrt{\Delta }}{2},0\right)$ 和 $E_2=\left(\frac{k-\sqrt{\Delta }}{2},0\right)$。

2) 若 $v\ne 0$，则且当 $b>c$ 且 $h_0<-\frac{u\left[ac-\left(b-c\right)k\right]}{k\left(b-c\right)}$ 时，系统(2)存在唯一的正平衡点 $E_3\left(u^{*},v^{*}\right)$。其中， $u^{*}=\frac{ac}{b-c}$，。

接下来我们讨论平衡点的类型及其稳定性，系统(2)的Jacobi矩阵为：

$J=\left(\begin{array}{cc}1-\frac{2u}{k}-\frac{av}{{\left(a+u\right)}^2}& -\frac{u}{a+u}\\ \frac{abv}{{\left(a+u\right)}^2}& \frac{bu}{a+u}-c\end{array}\right).$ (4)

1) 判断边界平衡点 $E_0$ 的类型

将 $E_0=\left(\frac{k}{2},0\right)$ 代入Jacobi矩阵(4)中，可以得到系统(2)在 $E_0$ 处的Jacobi矩阵为

$J_{E_0}=\left(\begin{array}{cc}0& -\frac{k}{k+2a}\\ 0& \frac{bk}{k+2a}-c\end{array}\right).$

利用文 [10] [11]，通过计算，可得结论平衡点 $E_0$ 为鞍结点。

2) 判断边界平衡点 $E_1$ 和 $E_2$ 的类型

将 $E_1=\left(\frac{k+\sqrt{\Delta }}{2},0\right)$ 和代入Jacobi矩阵(4)中，可得到系统(2)在 $E_1$ 和 $E_2$ 处的Jacobi矩阵为

$J_{E_1}=\left(\begin{array}{cc}-\frac{\sqrt{\Delta }}{k}& -\frac{k+\sqrt{\Delta }}{2a+k+\sqrt{\Delta }}\\ 0& \frac{b\left(k+\sqrt{\Delta }\right)}{2a+k+\sqrt{\Delta }}-c\end{array}\right),J_{E_2}=\left(\begin{array}{cc}\frac{\sqrt{\Delta }}{k}& -\frac{k-\sqrt{\Delta }}{2a+k-\sqrt{\Delta }}\\ 0& \frac{b\left(k-\sqrt{\Delta }\right)}{2a+k-\sqrt{\Delta }}-c\end{array}\right)$

利用文 [10] [11]，通过计算，可得如下结论:

a) 若 $c>\frac{b\left(k+\sqrt{\Delta }\right)}{2a+k+\sqrt{\Delta }}$，那么 $E_1$ 为稳定结点；若 $c<\frac{b\left(k+\sqrt{\Delta }\right)}{2a+k+\sqrt{\Delta }}$，则 $E_1$ 为鞍点。

b) 若 $c<\frac{b\left(k-\sqrt{\Delta }\right)}{2a+k-\sqrt{\Delta }}$，那么 $E_2$ 为不稳定结点；若 $c>\frac{b\left(k-\sqrt{\Delta }\right)}{2a+k-\sqrt{\Delta }}$，则 $E_2$ 为鞍点。

3) 正平衡点 $E_3$ 的稳定性

将带入Jacobi矩阵(4)中，可以得到系统(2)在 $E_3\left(u^{*},v^{*}\right)$ 处的Jacobi矩阵为

$J_{E_3}=\left(\begin{array}{cc}{f}_1\left(v^{*}\right)& -\frac{c}{b}\\ \frac{v^{*}{\left(b-c\right)}^2}{ab}& 0\end{array}\right).$ (5)

其中 $f_1\left(v^{*}\right)=1-\frac{2ac}{k\left(b-c\right)}-\frac{v^{*}{\left(b-c\right)}^2}{a{b}^2}$。由于

$D=\det J_{E_3}=\frac{c{v}^{*}{\left(b-c\right)}^2}{a{b}^2},T=trace{J}_{E_3}=f_1\left(v^{*}\right)=1-\frac{2ac}{k\left(b-c\right)}-\frac{v^{*}{\left(b-c\right)}^2}{a{b}^2}$

若 $D>0$， $T<0$ ；我们可得到系统(2)在正平衡点 $E_3\left(u^{*},v^{*}\right)$ 处是局部渐近稳定的；

若 $D>0$， $T>0$ ；我们可得到系统(2)在正平衡点 $E_3\left(u^{*},v^{*}\right)$ 处是不稳定的。

令 $f_1\left(v^{*}\right)=0$，既有 $D>0$， $T=0$ ；由参考文献 [12]，则Jacobi矩阵 $J_{E_3}$ 的特征方程有一对纯虚的特征根 $\pm i{\omega }_0$ (其中i为虚数单位， ${\omega }_0=\frac{b-c}{b}\sqrt{\frac{v^{*}c}{a}}>0$ )。

将 $E_3\left(u^{*},v^{*}\right)$ 代入 $f_1\left(v^{*}\right)=0$ 中，解得 $h_0^{*}=\frac{a{c}^2\left[a\left(b+c\right)-k\left(b-c\right)\right]}{{\left(b-c\right)}^{3}k}$。

因此有如下结论：

定理1 当 $b>c$ 且 $h_0<-\frac{u\left[ac-\left(b-c\right)k\right]}{k\left(b-c\right)}$ 条件成立时，

a) 若 $h_0>h_0^{*}$，则系统(2)的正平衡点是局部渐进稳定的；

b) 若 $0<h_0<h_0^{*}$，则系统(2)的正平衡点是不稳定的；

c) 若 $h_0=h_0^{*}$，则系统(2)在正平衡点 $E_3\left(u^{*},v^{*}\right)$ 处出现Hopf分支。

3. Hopf分支

由上一节讨论我们可以知道，当 $b>c$ 且 $h_0<-\frac{u\left[ac-\left(b-c\right)k\right]}{k\left(b-c\right)}$ 时，系统(2)有唯一正平衡点 $E_3\left(u^{*},v^{*}\right)$。其中 $E_3$ 为初等焦点或者结点。下面我们讨论系统(2)在 $E_3\left(u^{*},v^{*}\right)$ 附近的Hopf分支。

系统(2)在 $E_3\left(u^{*},v^{*}\right)$ 处的雅可比矩阵为(5)，其特征方程为

${\lambda }^2-\lambda f_1\left(v^{*}\right)+\frac{v^{*}c{\left(b-c\right)}^2}{a{b}^2}=0.$ (6)

且当 $f_1\left(v^{*}\right)=0$ 时，Jacobi矩阵的特征方程(6)有一对纯虚的特征根 $\pm i{\omega }_0$。

我们取 $h_0$ 为分支参数，设(6)的特征根为 $\lambda =p\left(h_0\right)+i\omega \left(h_0\right)$，那么我们可以得到 $p\left(h_0^{*}\right)=0$， $\omega \left(h_0^{*}\right)={\omega }_0$。将 $u=\frac{ac}{b-c}$ 代入系统(2)的第一个方程可得：

$k{\left(b-c\right)}^2\left(cv+h_{0}b\right)+abc\left[ac-k\left(b-c\right)\right]=0.$

记 $f\left(v\right)=k{\left(b-c\right)}^2\left(cv+h_{0}b\right)+abc\left[ac-k\left(b-c\right)\right]$，则有 $f\left(v^{*}\right)=0$， $f^{\prime }\left(v^{*}\right)>0$。

由 $f\left(v^{*}\right)=0$ 解得 $v^{*}\left(h_0\right)=-\frac{\left\{ac\left[ac-k\left(b-c\right)\right]+{\left(b-c\right)}^2{h}_{0}k\right\}b}{{\left(b-c\right)}^{2}kc}$，且有 $v^{*}{}^{\prime }\left(h_0\right)=\frac{b}{c}>0$。

又 $p\left(h_0\right)=\frac{f_1\left(v^{*}\right)}{2}$。于是有

由 $p^{\prime }\left(h_0^{*}\right)<0$ 可知系统(2)满足产生Hopf分支的条件 [13]。因此系统(2)在 $h_0=h_0^{*}$ 附近可产生Hopf分支。接下来我们来讨论产生的Hopf分支的方向。

令 $w_1=u-u^{*}$， $w_2=v-v^{*}$。则系统(2)变成以下形式

$\{\begin{array}{l}\frac{\text{d}{w}_1}{\text{d}\tau }=-\frac{c}{b}{w}_2-\frac{a^2{b}^3-k{v}^{*}{\left(b-c\right)}^3}{a^2{b}^{3}k}{w}_1^2-\frac{{\left(b-c\right)}^2}{a{b}^2}{w}_1{w}_2-\frac{v^{*}{\left(b-c\right)}^4}{a^3{b}^4}{w}_1^3+\frac{{\left(b-c\right)}^3}{a^2{b}^3}{w}_1^2{w}_2+\cdots ,\\ \frac{\text{d}{w}_2}{\text{d}\tau }=\frac{v^{*}{\left(b-c\right)}^2}{ab}{w}_1-\frac{v^{*}{\left(b-c\right)}^3}{a^2{b}^2}{w}_1^2+\frac{{\left(b-c\right)}^2}{ab}{w}_1{w}_2+\frac{v^{*}{\left(b-c\right)}^4}{a^3{b}^3}{w}_1^3-\frac{{\left(b-c\right)}^3}{a^2{b}^2}{w}_1^2{w}_2+\cdots .\end{array}$ (7)

因 ${\omega }_0=\frac{b-c}{b}\sqrt{\frac{v^{*}c}{a}}$，再作变换：

$\xi =-\frac{v^{*}{\left(b-c\right)}^2}{ab}{w}_1,\eta =-{\omega }_0{w}_2,$

系统(7)变为

$\{\begin{array}{l}\frac{\text{d}\xi }{\text{d}\tau }=-{\omega }_0\eta +\frac{a^2{b}^3{\xi }^2-k{v}^{*}{\left(b-c\right)}^3{\xi }^2}{a{b}^2{v}^{*}k{\left(b-c\right)}^2}+\frac{{\left(b-c\right)}^2\xi \eta }{a{b}^2{\omega }_0}+\frac{{\xi }^3}{a{b}^2{v}^{*}}+\frac{\left(b-c\right){\xi }^2\eta }{a{b}^2{v}^{*}{\omega }_0}+\cdots ,\\ \frac{\text{d}\eta }{\text{d}\tau }={\omega }_0\xi +\frac{{\omega }_0{\xi }^2}{v^{*}\left(b-c\right)}-\frac{\xi \eta }{v^{*}}+\frac{{\omega }_0{\xi }^3}{v^{*}{}^2{\left(b-c\right)}^2}-\frac{{\xi }^2\eta }{v^{*}{}^2\left(b-c\right)}+\cdots .\end{array}$ (8)

引入复变量：，可以得到：

$\frac{\text{d}z}{\text{d}\tau }=i{\omega }_{0}z+{\displaystyle \underset{2\le k+l\le \text{3}}{\sum }{g}_{kl}}{z}^k{\bar{z}}^l,$

其中：

$g_{20}=\frac{1}{4}\left[\frac{ab}{v^{*}k{\left(b-c\right)}^2}-\frac{b-c}{a{b}^2}-\frac{1}{v^{*}}\right]+\frac{1}{4}i\left[\frac{{\omega }_0}{v^{*}\left(b-c\right)}-\frac{{\left(b-c\right)}^2}{a{b}^2{\omega }_0}\right],$

$g_{11}=\frac{1}{2}\left[\frac{ab}{v^{*}k{\left(b-c\right)}^2}-\frac{b-c}{a{b}^2}\right]+\frac{1}{2}i\frac{{\omega }_0}{v^{*}\left(b-c\right)},$

$g_{21}=\frac{1}{8v^{*}}\left[\frac{3}{a{b}^2{v}^{*}}-\frac{1}{v^{*}\left(b-c\right)}\right]+\frac{1}{8v^{*}}i\left[\frac{3w_0}{v^{*}{\left(b-c\right)}^2}-\frac{b-c}{a{b}^2{\omega }_0}\right].$

把以上计算结果代入到一阶Lyapunov系数计算公式 [13]

$l_1\left(0\right)=\frac{1}{2{\omega }_0^2}\mathrm{Re}\left(i{g}_{20}{g}_{11}+{\omega }_0{g}_{21}\right),$

我们可以得到

$l_1\left(0\right)=\frac{1}{8{\omega }_0}\left\{\frac{1}{2a{b}^2{\omega }_0^2}\left[\frac{ab{v}^{*}+5}{v^{*}{}^{2}k}-\frac{{\left(b-c\right)}^3}{a{b}^2}\right]-\frac{ab}{v^{*}{}^2{\left(b-c\right)}^{3}k}\right\}.$

因此，根据文献 [13] 我们可以得到以下结论：

定理2 由 $l_1\left(0\right)=\frac{1}{8{\omega }_0}\left\{\frac{1}{2a{b}^2{\omega }_0^2}\left[\frac{ab{v}^{*}+5}{v^{*}{}^{2}k}-\frac{{\left(b-c\right)}^3}{a{b}^2}\right]-\frac{ab}{v^{*}{}^2{\left(b-c\right)}^{3}k}\right\}$，

a) 当 $l_1\left(0\right)<0$ 时，模型(1)在 $$ 附近可以产生超临界的Hopf分支；

b) 当 $l_1\left(0\right)>0$ 时，模型(1)在 $E_3$ 附近可以产生次临界的Hopf分支。

4. 仿真算例

由上一节定理2可知，若取 $a=1.8$， $b=3.2$，， $k=3.6$，则由于计算较为复杂，通过maple软件可计算得 $h_0^{*}=0.0648$， $u^{*}=1.08$， $v^{*}=1.8432$， ${\omega }_0=0.692820323$，且我们可以算出 $l_1\left(0\right)=-0.00202732318<0$。那么系统(2)在 $h_0=h_0^{*}$ 的时候就会有超临界Hopf分支产生。当 $h_0$ 在大于的附近取值时，在平衡点 $E_3$ 附近则存在唯一的稳定极限环。如图1(a)所示。

(a) (b)

Figure 1. (a) When $a=1.8,b=3.2,c=1.2,k=3.6,h_0=0.066$, there is a stable limit cycle near $E_3$ ; (b) When $a=1.9,b=2.3,c=0.5,k=1.5,h_0=0.14$, there is a unstable limit cycle near $E_3$

图1. (a)当时，在 $E_3$ 附近存在一个稳定的极限环； (b)当 $a=1.9,b=2.3,c=0.5,k=1.5,h_0=0.14$ 时，在 $E_3$ 附近存在一个不稳定极限环

若取 $a=1.9$， $b=2.3$， $c=0.5$， $k=1.5$，则 $h_0^{*}=0.1422610882$， $u^{*}=0.5277777778$， $v^{*}=0.919158665$， ${\omega }_0=0.3849001794$，我们可以算出 $l_1\left(0\right)=0.5205964260>0$。那么系统(2)就会在 $h_0=h_0^{*}$ 处有次临界的Hopf分支产生。当 $h_0$ 在小于 $h_0^{*}$ 附近取值时，则在平衡点 $E_3$ 附近存在唯一的不稳定极限环。如图1(b)所示。

5. 结论

文 [7] 用线性化方法，分析了具齐次Neumann边界条件的Holling II型捕食者–食饵模型的稳定性和Hopf分支。研究具有收获率的捕食者食饵模型对生态保护和可持续发展具有重大意义，因此本文对食饵种群进一步增加常数收获率后进行分析新的模型的平衡点类型及其稳定性和正平衡点的Hopf分支，即讨论了食饵种群具有常数收获率和Holling II型功能反应的捕食者–食饵模型的平衡点类型及其稳定性和正平衡点的Hopf分支问题且通过运用Matlab软件和Matcont软件做出了相应的Hopf分支图。由上面的讨论我们可以知道，超临界的Hopf分支可以产生稳定的极限环，其中的参数和初值是通过适当的选取来获得，最终使得捕食者和食饵的种群数量呈现周期性的变化而不消失。次临界的Hopf分支可以产生不稳定的极限环，其中的参数和初值是通过适当的选取来获得，当初值处在不稳定的极限环里面时，轨线趋向于正平衡点 $E_3$，也就是说最终捕食者和食饵种群的数量趋于稳定而不消失

参考文献