ABSTRACT

In this paper, we study the bifurcation of a predator-prey model with constant prey harvesting rate. By using the qualitative and bifurcation theory of ordinary differential equations, we analyze the existence conditions of equilibria, and discuss the Hopf and Bogdanov-Takens bifurcations near positive equilibria, respectively. The conditions and conclusions of corresponding bifurcations are obtained.

Keywords:Predator-Prey Model, Harvesting Rate, Hopf Bifurcation, Bogdanov-Takens Bifurcation

一类具有常数收获率的捕食者–食饵模型的 分支分析

邹桂华，刘宣亮

华南理工大学数学学院，广东 广州

收稿日期：2019年6月19日；录用日期：2019年6月29日；发布日期：2019年7月12日

摘 要

本文研究一类食饵种群具有常数收获率的捕食者–食饵模型的分支问题。利用常微分方程的定性与分支理论，分析了模型平衡点的存在条件，主要讨论了正平衡点的Hopf分支和Bogdanov-Takens分支，并得出了相应分支产生的条件及其结论。

关键词 :捕食者–食饵模型，收获率，Hopf分支，Bogdanov-Takens分支

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

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

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

1. 引言

近几十年来，人们对种群动力学模型进行了广泛的研究，特别是关于捕食者–食饵模型的研究，得到了大量的结论 [1] - [8] 。例如，对具有常数收获率的捕食者–食饵模型，文 [2] [3] [4] 进行了研究，讨论了模型的各种分支现象。文 [5] [6] [7] 考虑了疾病在捕食者种群或食饵种群中传播的捕食者–食饵模型，研究了模型平衡点的稳定性和分支行为。虽然对捕食者–食饵模型有很多的研究工作，但经典的捕食者–食饵模型仍存在一些不足的地方。例如文 [8] 考虑了由于食饵种群对捕食者种群的恐惧会造成食饵种群出生率的下降，于是对食饵种群的出生率引入了恐惧效应因子 $f\left(k,y\right)$，其中y为捕食者种群密度，k为与恐惧水平有关的常数，即将原来食饵种群的出生率 $r_0$ 修改为与捕食者种群密度有关的出生率： $f\left(k,y\right)r_0$，然后研究了捕食者–食饵模型的动力学性质。

受文 [8] 的启发，对食饵种群具有常数收获率的Leslie-Gower类捕食者–食饵模型 [2] [3] ，我们在食饵种群的出生率中也引入文 [8] 中考虑过的恐惧效应因子： $f\left(k,y\right)=1/\left(1+ky\right)$，即考虑如下的捕食者–食饵模型：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{r_{0}x}{1+ky}-dx-c{x}^2-axy-h,\\ \frac{\text{d}y}{\text{d}t}=r_{2}y\left(1-\frac{y}{bx}\right),\end{array}$ (1)

其中：x为食饵种群的密度，y为捕食者种群的密度， $\frac{r_0}{1+ky}$ 为食饵种群的出生率， $r_2$ 为捕食者种群的内

禀增长率，ax为每一个捕食者捕食食饵的数量，bx为捕食者种群中与食饵有关的容纳量，c为食饵种群内部竞争的死亡率，d为食饵种群的自然死亡率，h为捕食者对食饵的常数收获率。所有参数都为正参数。

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

$x=\frac{u}{bk},y=\frac{v}{k},\text{d}t=\frac{1}{r_2}\text{d}\tau ,$

则可将模型(1)变为

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}\tau }=\frac{ru}{1+v}-\alpha u^2-\beta uv-\theta u-h_0,\\ \frac{\text{d}v}{\text{d}\tau }=\frac{v}{u}\left(u-v\right).\end{array}$ (2)

其中： $r=\frac{r_0}{r_2},\alpha =\frac{c}{r_{2}bk},\beta =\frac{a}{r_{2}k},h_0=\frac{bkh}{r_2},\theta =\frac{d}{r_2}$。

2. 平衡点及其类型

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

1) 若 $v=0$，则有

$\alpha u^2+\left(\theta -r\right)u+h_0=0.$ (3)

记 $\Delta ={\left(\theta -r\right)}^2-4\alpha h_0$。

i) 当 $h_0=\frac{{\left(\theta -r\right)}^2}{4\alpha }$ 且时，则系统(2)有唯一边界平衡点 $E_0=\left(\frac{r-\theta }{2\alpha },0\right)$。

ii) 当 $h_0<\frac{{\left(\theta -r\right)}^2}{4\alpha }$ 且 $r>\theta$ 时，则系统(2)有两个边界平衡点

$E_{\text{1}}=\left(\frac{r-\theta -\sqrt{\Delta }}{2\alpha },0\right)$ 和 $E_{\text{2}}=\left(\frac{r-\theta +\sqrt{\Delta }}{2\alpha },0\right).$

2) 若 $v\ne 0$，则 $u=v$，且

$\left(\alpha +\beta \right)u^3+\left(\alpha +\beta +\theta \right)u^2+\left(\theta +h_0-r\right)u+h_0=0.$ (4)

记 $L\left(u\right)=\left(\alpha +\beta \right)u^3+\left(\alpha +\beta +\theta \right)u^2+\left(\theta +h_0-r\right)u+h_0$，则有

$L^{\prime }\left(u\right)=3\left(\alpha +\beta \right)u^2+2\left(\alpha +\beta +\theta \right)u+\left(\theta +h_0-r\right),$$L\left(0\right)=h_0>0$

当 $r>\theta +h_0$ 时，记

$u^{*}=\frac{-\left(\alpha +\beta +\theta \right)+\sqrt{{\left(\alpha +\beta +\theta \right)}^2-3\left(\alpha +\beta \right)\left(\theta +h_0-r\right)}}{3\left(\alpha +\beta \right)},$

则 $u^{*}$ 为 $L^{\prime }\left(u\right)=0$ 的唯一正根。如下结论成立：

i) 当 $L\left(u^{*}\right)=0$ 时，系统(2)有唯一正平衡点 $E_{\text{3}}\left(u^{*},u^{*}\right)$。

ii) 当 $L\left(u^{*}\right)<0$ 时，系统(2)有两个正平衡点 $E_4\left(u_4,u_4\right)$ 和 $E_5\left(u_5,u_5\right)$，其中 $u_4<u^{*}<u_5,L\left(u_4\right)=L\left(u_5\right)=0,L^{\prime }\left(u_4\right)<0,L^{\prime }\left(u_5\right)>0$。

通过讨论系统(2)在平衡点处的Jacobi矩阵的特征值，利用文 [9] [10] ，可得如下结论： $E_{\text{0}}$ 为一个鞍结点， $E_{\text{1}}$ 为不稳定的非鞍初等奇点， $E_{\text{2}}$ 为鞍点， $E_3$ 为高次奇点， $E_4$ 为鞍点， $E_5$ 为非鞍初等奇点。

3. Hopf分支

由上一节可知，当 $r>\theta +h_0$ 且 $L\left(u^{*}\right)<0$ 时，系统(2)有两个正平衡点 $E_4\left(u_4,u_4\right)$ 和 $E_5\left(u_5,u_5\right)$，其中 $E_{\text{4}}$ 是鞍点， $E_5$ 为初等焦点或结点。下面讨论系统(2)在 $E_5\left(u_5,u_5\right)$ 附近的Hopf分支。

系统(2)在 $E_5\left(u_5,u_5\right)$ 处的雅可比矩阵为：

$A=\left(\begin{array}{cc}{f}_1\left(u_5\right)& f_2\left(u_5\right)\\ 1& -1\end{array}\right),$

其中 $f_1\left(u_5\right)=\frac{r}{1+u_5}-\left(2\alpha +\beta \right)u_5-\theta$， $f_2\left(u_5\right)=-\left(\beta +\frac{r}{{\left(1+u_5\right)}^2}\right)u_5$。则矩阵A的特征方程为

${\lambda }^2+\left(1-f_1\left(u_5\right)\right)\lambda -\left(f_1\left(u_5\right)+f_2\left(u_5\right)\right)=0$ . (5)

当 $f_1\left(u_5\right)=1$ 时，矩阵A的特征方程(5)有一对纯虚根 $\pm i{\omega }_0$ (i为虚数单位， ${\omega }_{\text{0}}=\sqrt{-1-f_2\left(u_5\right)}>0$ )。由 $E_{\text{5}}$ 是系统(2)的平衡点以及 $f_1\left(u_5\right)=1$ 可得

$\{\begin{array}{l}\left(\alpha +\beta \right)u_5^3+\left(\alpha +\beta +\theta \right)u_5^2+\left(\theta +h_0-r\right)u_5+h_0=0,\\ \left(2\alpha +\beta \right)u_5^2+\left(2\alpha +\beta +\theta +1\right)u_5+\theta -r+1=0.\end{array}$ (6)

由辗转相除法，可得

$u_5=-\frac{\left(4\alpha h_0-r+1\right)\left(\alpha +\beta \right)+\alpha \theta \left(r-\theta \right)+\beta \theta +{\beta }^2{h}_0}{2\left(2h_0-r\right){\alpha }^2+\left(4h_0-r\right)\alpha \beta -\alpha {\theta }^2+{\beta }^2{h}_0+\beta \theta +\alpha +\beta }.$

将 $u_5$ 代入 $\left(2\alpha +\beta \right)u_5^2+\left(2\alpha +\beta +\theta +1\right)u_5+\left(\theta -r+1\right)=0$，得到A的特征根有一对纯虚根的条件： $\theta <\frac{\beta }{\alpha }+1-2\alpha -\beta$ 且

$\begin{array}{c}{h}_0=\frac{1}{2{\left(2\alpha +\beta \right)}^2}[\left(4{\alpha }^2-4\alpha -1\right)\left(\alpha +\beta \right)-\beta \left(\beta +\theta \right)+\alpha \left({\beta }^2+{\theta }^2\right)+2\alpha \left(2\alpha +\beta \right)r\\ -\left(\alpha \left(2\alpha +\beta +\theta -1\right)-\beta \right)\sqrt{{\left(2\alpha +\beta -\theta -1\right)}^2+4r\left(2\alpha +\beta \right)}]\\ \equiv h_0^{*}\end{array}$

取 $h_0$ 为分支参数。设特征方程(5)的特征根为： $\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$。由于 $L\left(u_5\right)=0,L^{\prime }\left(u_5\right)>0$，则由 $L\left(u_5\right)=0$ 确定的隐函数 $u_5=u_5\left(h_0\right)$ 的导数： ${u^{\prime }}_5\left(h_0\right)=-{L^{\prime }}_{h_0}\left(u_5\right)/{L^{\prime }}_{u_5}\left(u_5\right)=-\left(u_5+1\right)/L^{\prime }\left(u_5\right)<0$，又 $p\left(h_0\right)=\frac{1}{2}\left[f_1\left(u_5\right)-1\right]$，于是有

$p^{\prime }\left(h_0\right)=\frac{1}{2}\frac{\partial f_1}{\partial u_5}{u^{\prime }}_5\left(h_0\right)=\frac{1}{2}\left[-\frac{r}{{\left(1+u_5\right)}^2}-\left(2\alpha +\beta \right)\right]{u^{\prime }}_5\left(h_0\right)>0$ .

由 $$ 知系统(2)产生Hopf分支的横截条件 [11] 满足，因此系统(2)在 $h_0=h_0^{*}$ 附近可产生Hopf分支。而且，当 $h_0<h_0^{*}$ 时，平衡点 $E_5$ 是稳定的，当 $h_0>h_0^{*}$ 时，平衡点 $E_5$ 是不稳定的。下面讨论Hopf分支的方向。

令 $w_{\text{1}}=u-u_5$， $w_2=v-u_5$，则系统(2)变成以下形式

$\{\begin{array}{l}\frac{\text{d}{w}_1}{\text{d}\tau }=w_1-\left(\frac{r}{{\left(u_5+1\right)}^2}+\beta \right)u_5{w}_2-\alpha w_1^2-\left(\frac{r}{{\left(u_5+1\right)}^2}+\beta \right)w_1{w}_2+\frac{r{u}_5}{{\left(u_5+1\right)}^3}{w}_2^2+\frac{r{w}_1{w}_2^2}{{\left(u_5+1\right)}^3}-\frac{r{u}_5}{{\left(u_5+1\right)}^4}{w}_2^3+\cdots ,\\ \frac{\text{d}{w}_2}{\text{d}\tau }=w_1-w_2-\frac{1}{u_5}{w}_1^2+\frac{2}{u_5}{w}_1{w}_2-\frac{1}{u_5}{w}_2^2+\frac{1}{u_5^2}{w}_1^3-\frac{2}{u_5^2}{w}_1^2{w}_2+\frac{1}{u_5^2}{w}_1{w}_2^2+\cdots .\end{array}$ (7)

再作变换：

$\xi =w_2-w_1,$ $\eta =-{\omega }_0{w}_2,$

系统(7)变为

$\{\begin{array}{l}\frac{\text{d}\xi }{\text{d}\tau }=-{\omega }_0\eta +\left(\alpha -\frac{1}{u_5}\right){\xi }^2+\left(\frac{2\alpha }{{\omega }_0}+\frac{1+{\omega }_0^2}{u_5{\omega }_0}\right)\xi \eta +\frac{1}{{\omega }_0^2}\left(\alpha +\frac{{\omega }_0^2+1}{u_5}-\frac{r{u}_5}{{\left(u_5+1\right)}^3}\right){\eta }^2\\ -\frac{1}{u_5^2}{\xi }^3-\frac{1}{u_5^2{\omega }_0}{\xi }^2\eta +\frac{r}{{\left(u_5+1\right)}^3{\omega }_0^2}\xi {\eta }^2+\frac{r{\eta }^3}{{\left(u_5+1\right)}^4{\omega }_0^3}+\cdots ,\\ \frac{\text{d}\eta }{\text{d}\tau }={\omega }_0\xi +\frac{{\omega }_0}{u_5}{\xi }^2+\frac{{\omega }_0}{u_5^2}{\xi }^3+\frac{1}{u_5^2}{\xi }^2\eta +\cdots .\end{array}$

引入复变量： $z=\xi +i\eta$，可以得到：

$\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[\alpha -\frac{2}{u_5}-\frac{\alpha }{{\omega }_0^2}-\frac{1}{u_5{\omega }_0^2}+\frac{r{u}_5}{{\omega }_0^2{\left(u_5+1\right)}^3}\right]+\frac{1}{4}i\left(\frac{{\omega }_0}{u_5}-\frac{2\alpha }{{\omega }_0}-\frac{1+{\omega }_0^2}{u_5{\omega }_0}\right),$

$g_{11}=\frac{1}{2}\left[\alpha +\frac{\alpha }{{\omega }_0^2}+\frac{1}{u_5{\omega }_0^2}-\frac{r{u}_5}{{\left(u_5+1\right)}^3{\omega }_0^2}\right]+\frac{{\omega }_0}{2u_5}i,$

$g_{21}=\frac{1}{8}\left[-\frac{2}{u_5^2}+\frac{r}{{\left(u_5+1\right)}^3{\omega }_0^2}\right]+\frac{1}{8}i\left[\frac{3{\omega }_0^2+1}{u_5^2{\omega }_0}-\frac{3r}{{\left(u_5+1\right)}^4{\omega }_0^3}\right].$

将上述结果代入一阶Lyapunov系数的计算公式 [11]

$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}{16{\omega }_0^2}\left\{\left[\alpha +\frac{\alpha }{{\omega }_0^2}+\frac{1}{u_5{\omega }_0^2}-\frac{r{u}_5}{{\left(u_5+1\right)}^3{\omega }_0^2}\right]\frac{2\alpha u_5+1+{\omega }_0^2}{u_5{\omega }_0}-2\left(\alpha -\frac{1}{u_5}\right)\frac{{\omega }_0}{u_5}-\frac{2{\omega }_0}{u_5^2}+\frac{r}{{\left(u_5+1\right)}^3{\omega }_0}\right\}.$

因此，由文 [11] 可得以下结论：

i) 当 $l_1\left(0\right)<0$ 时，模型(1)在 $E_{\text{5}}$ 附近有超临界的Hopf分支产生；

ii) 当 $l_1\left(0\right)>0$ 时，模型(1)在 $E_{\text{5}}$ 附近有次临界的Hopf分支产生。

例如，取 $\alpha =\text{0}\text{.6},\beta =\text{5},r=\text{12},\theta =0.\text{2}$，则 $h_0^{*}=1.286933942$， $u_5=0.851699139$， ${\omega }_0=2.497849632$，可算出 $l_1\left(0\right)=-0.00309846244<0$。于是系统(2)在 $h_0=h_0^{*}$ 时产生超临界的Hopf分支。当 $h_0$ 在大于 $h_0^{*}$ 附近取值时，在平衡点 $E_5$ 附近存在唯一稳定极限环。如图1(a)所示。若取 $\alpha =0.\text{6},\beta =\text{5},r=\text{3},\theta =0.\text{2}$，则 $h_0^{*}=0.2330270625$， $u_5=0.\text{2072544242}$， ${\omega }_0=0.6803518274$，可算出 $l_1\left(0\right)=\text{18}\text{.74767753}>0$，于是系统(2)在 $$ 处产生次临界的Hopf分支。当 $h_0$ 在小于 $h_0^{*}$ 附近取值时，在平衡点 $E_5$ 附近存在唯一不稳定极限环。如图1(b)所示。

4. Bogdanov-Takens分支

当 $f_1\left(u_5\right)=1$ 且 $f_2\left(u_5\right)=-1$ 时，系统(2)在正平衡点 $E_{\text{5}}$ 处的特征方程(5)有两个零特征值，此时 $u_5$ 满足方程组(6)及

$\beta u_5^3+\left(2\beta -1\right)u_5^2+\left(\beta +r-2\right)u_5-1=0.$

利用辗转相除法，最终可得到特征方程(5)有两个零特征值的参数条件：

$\{\begin{array}{l}r=\frac{1}{8{\left(\alpha +\beta \right)}^2}\{-2\alpha \beta \left(2\alpha +2\beta -13\right)+2\theta \left(2{\alpha }^2+6\alpha \beta -\beta \right)-{\beta }^3-\theta \left(-6{\beta }^2+\beta \theta +2\alpha \right)\\ +16{\alpha }^2+10{\beta }^2-\left[{\left(2\alpha +\beta \right)}^2+\beta \theta +2\left(\alpha +\beta \right)\right]\sqrt{{\left(\beta +\theta \right)}^2+8\left(\alpha +\beta \right)}\}\triangleq r_0,\\ h_0=\frac{\alpha \left(\theta -r_0\right)\left(\beta +\theta \right)+\left(2\alpha +2r_0-3\right)\left(\alpha +\beta \right)-2\beta \theta }{{\left(2\alpha +\beta \right)}^2+\beta \theta +2\left(\alpha +\beta \right)}\triangleq h_{00}.\end{array}$

图1. (a) 当 $\alpha =0.6,\beta =5,r=12,\theta =0.2,h_0=1.29$ 时，在 $E_5$ 附近有一个稳定极限环；(b) 当 $\alpha =0.6,\beta =5,r=3,\theta =0.2,h_0=0.23$ 时，在 $E_5$ 附近有一个不稳定极限环

我们选取 $h_0,r$ 作为分支参数，讨论系统(2)在 $\left(h_0,r\right)=\left(h_{00},r_0\right)$ 附近的Bogdanov-Takens分支问题。

记 $X={\left(u,v\right)}^{\text{T}}$， $\mu ={\left(h_0,r\right)}^{\text{T}}$， $X_0={\left(u_5,u_5\right)}^{\text{T}}$， ${\mu }_0={\left(h_{00},r_0\right)}^{\text{T}}$。于是系统(2)可以写成

$\frac{\text{d}X}{\text{d}\tau }=F\left(X,\mu \right)=\left(\begin{array}{c}\frac{ru}{1+v}-\alpha u^2-\beta uv-\theta u-h_0\\ \frac{v}{u}\left(u-v\right)\end{array}\right),$

其中 $X\in R^2,\mu \in R^2$。 $F\left(X,\mu \right)$ 在 $\left(X_0,{\mu }_0\right)$ 处的泰勒展开式为

$\begin{array}{c}F\left(X,\mu \right)=DF\left(X_0,{\mu }_0\right)\left(X-X_0\right)+F_{\mu }\left(X_0,{\mu }_0\right)\left(\mu -{\mu }_0\right)+\frac{1}{2}{D}^{2}F\left(X_0,{\mu }_0\right)\left(X-X_0,X-X_0\right)\\ +F_{\mu X}\left(X_0,{\mu }_0\right)\left(X-X_0,\mu -{\mu }_0\right)+\cdots \end{array}$

其中

$DF\left(X_0,{\mu }_0\right)=\left(\begin{array}{cc}{f}_1\left(u_5\right)& f_2\left(u_5\right)\\ 1& -1\end{array}\right)=\left(\begin{array}{cc}1& -1\\ 1& -1\end{array}\right).$

设有二阶可逆矩阵 $P=\left(p_1,p_2\right)$，使得 $P^{-1}DF\left(X_0,{\mu }_0\right)P=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$。可取

$p_1=\left(\begin{array}{c}1\\ 1\end{array}\right),$ $p_2=\left(\begin{array}{c}2\\ 1\end{array}\right)$

再设 $P^{-1}={\left(q_1,q_2\right)}^{\text{T}}$，可得

$q_1=\left(\begin{array}{c}-1\\ 2\end{array}\right),$ $q_2=\left(\begin{array}{c}1\\ -1\end{array}\right).$

利用文献 [12] 中的记号，通过计算，不难得到：

$\tilde{a}\equiv \frac{1}{2}{p}_1^{\text{T}}\left(q_2\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_1=-\left(\alpha +\beta \right)-\frac{r}{{\left(1+u_5\right)}^3},$

$\tilde{b}\equiv p_1^{\text{T}}\left(q_1\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_1+p_1^{\text{T}}\left(q_2\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_2=-2\alpha -\beta -\frac{r}{{\left(1+u_5\right)}^2},$

$S_1^{\text{T}}\equiv q_2^{\text{T}}{F}_{\mu }\left(X_0,{\mu }_0\right)=\left(-1,\frac{u_5}{1+u_5}\right),$

$\begin{array}{c}{S}_2\equiv [\frac{2\tilde{a}}{\tilde{b}}\left(p_1^{\text{T}}\left(q_1\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_2+p_2^{\text{T}}\left(q_2\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_2\right)\\ \begin{array}{c}\\ \end{array}-p_{\text{1}}^{\text{T}}\left(q_2\cdot D^{2}F\left(X_0,{\mu }_0\right)\right)p_2]\times F_{\mu }^{\text{T}}\left(X_0,{\mu }_0\right)q_1\\ -\frac{2\tilde{a}}{\tilde{b}}{\displaystyle \underset{i=1}{\overset{2}{\sum }}\left(q_i\cdot F_{\mu x}\left(X_0,{\mu }_0\right)\right)p_i}+\left(q_2\cdot F_{\mu x}\left(X_0,{\mu }_0\right)\right)p_1\\ ={\left(s_{21},s_{22}\right)}^{\text{T}}\end{array}$

其中

$s_{21}=\frac{4\alpha u_5^4+\left(4\alpha r+12\alpha +5\right)u_5^3+\left(12\alpha -4r+15\right)u_5^2+\left(4\alpha +15\right)u_5+5}{u_5{\left(1+u_5\right)}^3\left(2\alpha u_5+1\right)}.$

$s_{22}=-\frac{6\alpha u_5^4+\left(4\alpha r+16\alpha +7\right)u_5^3+\left(14\alpha -6r+20\right)u_5^2+\left(4\alpha +19\right)u_5+6}{{\left(1+u_5\right)}^4\left(2\alpha u_5+1\right)}.$

记 ${\lambda }_1=h_0-h_{00}$， ${\lambda }_2=r-r_0$，则

${\beta }_1\equiv S_1^{\text{T}}\left(\mu -{\mu }_0\right)=-{\lambda }_1+\frac{{\lambda }_2{u}_5}{1+u_5},$ ${\beta }_{\text{2}}\equiv S_{\text{2}}^{\text{T}}\left(\mu -{\mu }_0\right)=s_{21}{\lambda }_1+s_{22}{\lambda }_2.$

由文献 [12] 可知，当 $0<\Vert \left({\lambda }_1,{\lambda }_2\right)\Vert \ll 1$ 时，系统(2)在正平衡点 $E_{\text{5}}$ 附近的中心流形上的轨线行为局部拓扑等价于以下系统

$\{\begin{array}{l}\frac{\text{d}{z}_1}{\text{d}\tau }=z_2,\\ \frac{\text{d}{z}_2}{\text{d}\tau }={\beta }_1+{\beta }_2{z}_1+\tilde{a}{z}_1^2+\tilde{b}{z}_1{z}_2.\end{array}$ (8)

由于 $\tilde{a}<0,\tilde{b}<0$，再作变量变换

$\tau =-\frac{\tilde{b}}{\tilde{a}}{t}_1,$ $z_1=\frac{\tilde{a}}{{\tilde{b}}^2}{\eta }_1,$ $z_2=-\frac{{\tilde{a}}^2}{{\tilde{b}}^3}{\eta }_{\text{2}},$

于是系统(8)变为

$\{\begin{array}{l}\frac{\text{d}{\eta }_1}{\text{d}{t}_1}={\eta }_2,\\ \frac{\text{d}{\eta }_2}{\text{d}{t}_1}={\gamma }_1+{\gamma }_2{\eta }_1+{\eta }_1^2-{\eta }_1{\eta }_2.\end{array}$ (9)

其中： ${\gamma }_1=\frac{{\tilde{b}}^4}{{\tilde{a}}^3}{\beta }_1$， ${\gamma }_2=\frac{{\tilde{b}}^2}{{\tilde{a}}^2}{\beta }_2$。对系统(9)，利用文 [11] 中关于Bogdanov-Takens分支的讨论，可知，当 $0<\Vert \left({\gamma }_1,{\gamma }_2\right)\Vert \ll 1$ 时，系统(9)在原点附近有以下分支现象：

i) 有鞍结点分支曲线：；

ii) 有Hopf分支曲线： $H=\left\{\left({\gamma }_1,{\gamma }_2\right):{\gamma }_1=0,{\gamma }_2<0\right\}$ ；

iii) 有同宿分支曲线： $HL=\left\{\left({\gamma }_1,{\gamma }_2\right):25{\gamma }_1+6{\gamma }_2^2=o\left({\gamma }_2^2\right),{\gamma }_2<0\right\}$。

例如，取参数 $\alpha =\text{0}\text{.6},\beta =5,\theta =0.2$ 时，可算出 $r_0=2.414799684$， $h_{00}=0.1590754631$，当 $\left(r,h_0\right)=\left(r_0,h_{00}\right)$ 时，系统(2)在正平衡点 $E_5\left(0.146243216,0.146243216\right)$ 处有两个零特征值。当 $0<\Vert \left(r-r_0,h_0-h_{00}\right)\Vert \ll 1$ 时，在 $E_5$ 附近有Bogdanov-Takens分支产生。出现的鞍结点分支曲线SN，Hopf分支曲线H与同宿分支曲线HL将整个参数平面 $O{\gamma }_1{\gamma }_2$ 分为四个区域，分别为图2(a)中区域①~④。下面在 $\left(r_0,h_0{}_0\right)$ 附近适当选取几组参数值 $\left(r,h_0\right)$，使得 $\left({\gamma }_1,{\gamma }_2\right)$ 分别位于每一个区域上，并作出系统(2)对应的相图。

取 $h_0=0.1690754631$。当 $r=2.464799684$ 时， $\left({\gamma }_1,{\gamma }_2\right)=\left(0.0404349372,0.0870823038\right)\in$ ①，此时系统(2)的相图如图2(b)所示，系统(2)无正平衡点；当 $r=2.492799684$ 时， $\left({\gamma }_1,{\gamma }_2\right)=\left(0.0005403447,-0.0761137454\right)\in$ ②，此时系统(2)的相图如图2(c)所示；当 $r=2.493799684$ 时， $\left({\gamma }_1,{\gamma }_2\right)=\left(-0.0008844622,-0.0819421757\right)\in$ ③，此时系统(2)的相图如图2(d)所示，系统(2)有一个不稳定极限环(红色闭曲线)出现；当 $r=2.514799684$ 时， $\left({\gamma }_1,{\gamma }_2\right)=\left(-0.0308054066,-0.2043392125\right)\in$ ④，此时系统(2)的相图如图2(e)所示；当 $r=2.494410328$ 时， $\left({\gamma }_1,{\gamma }_2\right)=\left(-0.0017545115,-0.0855012698\right)\in HL$，此时系统(2)的相图如图2(f)所示，系统(2)有一条同宿轨(红色曲线)出现。

5. 结论

文 [8] 中给出了具有恐惧效应的捕食者–食饵模型，讨论了模型的平衡点的稳定性，极限环的存在性和余维1的Hopf分支。本文对食饵种群进一步增加常数收获率后进行分析，发现模型除了有余维1的Hopf分支外，还会产生余维2的Bogdanov-Takens分支。本文对这两种分支现象分别进行了研究，并作出了相应的相图和分支图，作图时使用了Matlab软件和Matcont软件。由本文的结论可看出，通过Hopf分支可产生极限环，其中超临界Hopf分支可产生稳定极限环，表明参数和初值的适当选取，可使捕食者和食饵种群数量最终呈周期性变化而不消失。次临界Hopf分支可产生不稳定极限环，表明当初值位于不稳

图2. 当α = 0.6, β = 5, θ = 0.2, h₀ = 0.1690754631时，分别取r = 2.464799684, 2.492799684, 2.493799684, 2.514799684, 2.494410328，系统(2)对应的相图分别为 (b), (c), (d), (e), (f), 其中在图(d)中有一个不稳定极限环(红色闭曲线)出现，在图(f)中有一条同宿轨(红色曲线)出现

定极限环内时，轨线趋向正平衡点E₅，即捕食者和食饵种群数量最终趋于稳定而不消失。由Bogdanov-Takens分支可产生三种分支曲线SN，H和HL的结论，表明除了由Hopf分支产生极限环外，还可由鞍结点分支产生两个正平衡点，以及由同宿分支出现同宿轨和分支出不稳定极限环，当初值位于同宿轨或极限环内时，轨线趋于正平衡点E₅，即参数和初值的适当选取，可使捕食者和食饵种群数量最终趋于稳定而不消失。

基金项目

国家自然科学基金(11572127)。

文章引用

邹桂华,刘宣亮. 一类具有常数收获率的捕食者–食饵模型的分支分析
Bifurcation Analysis of a Predator-Prey Model with Constant Prey Harvesting[J]. 动力系统与控制, 2019, 08(03): 181-190. https://doi.org/10.12677/DSC.2019.83020

参考文献