一类具退化承载力捕食者–食饵模型的Hopf分支

Hopf Bifurcation of a Predator-Prey Model with Degenerate Carrying Capacity

摘要: 本文研究带Holling-II型功能反应的具有退化承载能力的捕食者–食饵模型的稳定性与Hopf分支。首先讨论平衡点的局部渐近稳定性，然后以退化系数ρ为分支参数，给出Hopf分支存在的条件。最后利用规范型理论和中心流形定理分析Hopf分支的方向及分支周期解的稳定性。

Abstract: In this paper, we study the stability and the Hopf bifurcation of apredator-prey model with degenerate carrying capacity with Holling-II type of functional response. Firstly, the local asymptotic stability of the equilibrium point is discussed. Then, the existence condition of the Hopf bifurcation is given by taking the degradation coefficient ρ as the branching parameter. Finally, the direction of the Hopf bifurcation and the stability of its periodic solution are analyzed by means of canonical theory and central manifold theorem.

1. 引言

高原鼠兔是高山草甸生态系统中常见的一种小型狐形哺乳动物。高山草甸是中国和西藏重要的生态系统组成部分。高山草甸提供空气和水净化、生物多样性等重要的生态服务维护，是生态系统中非常重要的组成部分。近几十年里，人们已经观察到忽略高山草甸生态系统对环境造成了一些严重的负面影响 [1] [2] [3]。鼠兔在高山草甸经常被视为害虫，因为它们被观察到与牲畜争夺食物 [4] 并在土壤中制造洞，增加土壤侵蚀 [5] [6] [7]。各种研究表明，高原鼠兔密度高，导致草甸退化 [4] [8]。鼠兔成为植被生长区域减少的原因(通过挖坑和扔土堆)，即减少对猎物的承载能力 [9]。刘等人 [10] 在他们的模型中考虑了这一想法，并分析了它的存在性和稳定性以及平衡点的分岔和极限环的存在性。高原鼠兔与高山草甸植被之间存在捕食–被捕食系统的关系，考虑承载能力下降的事实以及植被高度的影响。在文献 [11] 中作者建立了非线性退化承载力的捕食者–食饵模型，模型如下：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=gx\left(1-\frac{x}{\frac{k}{1+\rho y}}\right)-\frac{P}{\beta +x}xy,\\ \frac{\text{d}y}{\text{d}t}=-dy+\frac{qxy}{\beta +x}-\eta xy,\end{array}$ (1)

其中x和y分别表示食饵和捕食者种群的密度。参数g代表内禀增长率。 $\frac{P}{\beta +x}$ 是Holling圆盘函数 [12]，代表捕食者的捕食率。参数k代表没有捕食者时的环境容纳量。参数 $\rho$ 代表由于高原鼠兔挖洞造成承载力退化系数。d为捕食者种群的自然死亡率，q为转化率， $\eta x$ 为捕食者的衰减率，与食饵高度有关。

本文主要研究模型(1)平衡点的稳定性，其次以参数 $\rho$ 为分支参数，讨论Hopf分支的存在性。接下来通过中心流行定理，规范型理论和文献 [13] 的方法分析Hopf分支的方向以及分支周期解的稳定性，最后用数值模拟验证所得到的结论。

2. 平衡点的存在性和稳定性

2.1. 平衡点的存在性

对于模型(1)，

(i) 总存在平凡平衡点 $E_0=\left(0,0\right)$ 和半平凡平衡点 $E_1=\left(k,0\right)$。

(ii) 假设(H₁) $\eta \beta +d-q<0,{\left(\eta \beta +d-q\right)}^2-4\eta \beta d\ge 0$ 成立。

当且仅当(H₁₁) $x^{*}<k<x_1^{*}$，系统(1)存在一个正常数平衡点 $E^{*}=\left(x^{*},y^{*}\right)$。

当且仅当(H₁₂) $k>x_1^{*}$，系统(1)存在两个正常数平衡点 $E^{*}=\left(x^{*},y^{*}\right)$， $E_1^{*}=\left(x_1^{*},y_1^{*}\right)$ 其中

$x^{*}=\frac{-\left(\eta \beta +d-q\right)-\sqrt{{\left(\eta \beta +d-q\right)}^2-4\eta \beta d}}{2\eta },x_1^{*}=\frac{-\left(\eta \beta +d-q\right)+\sqrt{{\left(\eta \beta +d-q\right)}^2-4\eta \beta d}}{2\eta }$

$y^{*}=\frac{g\left(k-x^{*}\right)\left(\beta +x^{*}\right)}{g\rho \left(\beta +x^{*}\right)x^{*}+Pk},y_1^{*}=\frac{g\left(k-x_1^{*}\right)\left(\beta +x_1^{*}\right)}{g\rho \left(\beta +x_1^{*}\right)x_1^{*}+Pk}$

2.2. 平衡点的局部稳定性

系统(1)在 $\left(x,y\right)$ 处的Jacobi矩阵如下

$J=\left(\begin{array}{cc}g-\frac{2gx}{k}\left(1+\rho y\right)-\frac{P\beta y}{{\left(\beta +x\right)}^2}& -\frac{g\rho x^2}{k}-\frac{Px}{\beta +x}\\ \frac{q\beta y}{{\left(\beta +x\right)}^2}-\eta y& -d+\frac{qx}{\beta +x}-\eta x\end{array}\right)$ (2)

下面通过计算系统(1)在每个平衡点处的Jacobi矩阵的特征值，来确定这些平衡点的稳定性。

定理1 (i) 平凡平衡点 $E_0=\left(0,0\right)$ 是无条件不稳定的。

(ii) 若 $\eta k+\frac{d\beta }{k}>0$ 成立，则半平凡平衡点 $E_1=\left(k,0\right)$ 不稳定；否则 $E_1$ 是局部渐近稳定的。

证明(i) 系统(1)在平衡点 $E_0$ 处的Jacobi矩阵为

$J_{E_0}=\left(\begin{array}{cc}g& 0\\ 0& -d\end{array}\right)$ (3)

矩阵(3)的特征值为 ${\lambda }_1=g,{\lambda }_2=-d$ 因此，平衡点 $E_0$ 是鞍点，是不稳定的。

(ii) 系统(1)在平衡点 $E_1$ 处的Jacobi矩阵为

$J_{E_1}=\left(\begin{array}{cc}-g& -g\rho k-\frac{Pk}{\beta +k}\\ 0& -d+\frac{qk}{\beta +k}-\eta k\end{array}\right)$ (4)

若 $-d+\frac{qk}{\beta +k}-\eta k>0$ 即 $\eta k+\frac{d\beta }{k}>0$ 成立，则 $E_1$ 是鞍点，不稳定。

若 $-d+\frac{qk}{\beta +k}-\eta k<0$ 即 $\eta k+\frac{d\beta }{k}<0$ 成立，则 $E_1$ 渐近稳定。

定理2假设(H₁) (H₁₁)成立。若

(H₂) $q>\frac{\eta {\left(\beta +x^{*}\right)}^2}{\beta },\rho >\frac{P\left(k-\beta -2x^{*}\right)}{g{\left(\beta +x^{*}\right)}^2}$ 且 $\left(k-\beta -2x^{*}>0\right)$，则 $E^{*}$ 是局部渐近稳定的。

当(H₂)中， $\rho <\frac{P\left(k-\beta -2x^{*}\right)}{g{\left(\beta +x^{*}\right)}^2}$ 且 $\left(k-\beta -2x^{*}>0\right)$ 成立，则 $E^{*}$ 是不稳定的。

当(H₁) (H₁₂) (H₂)成立，经计算 $E_1^{*}$ 是鞍点，不稳定。此时 $E^{*}$ 稳定性同上。

证明系统(1)在正常数平衡点 $E^{*}$ 处的Jacobi矩阵为

$J_{E^{*}}=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right)$ (5)

其中

$a_{11}=g-\frac{2g{x}^{*}}{k}\left(1+\rho y^{*}\right)-\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^2},a_{12}=-\frac{g\rho x^{*2}}{k}-\frac{P{x}^{*}}{\beta +x^{*}}$

$a_{21}=\frac{q\beta y^{*}}{{\left(\beta +x^{*}\right)}^2}-\eta y^{*},a_{22}=0$

矩阵(5)的特征方程为

${\lambda }^2-Tr\left[J\left(E^{*}\right)\right]\lambda +Det\left[J\left(E^{*}\right)\right]=0$ (6)

其中

$\begin{array}{l}Tr\left[J\left(E^{*}\right)\right]=g-\frac{2g{x}^{*}}{k}\left(1+\rho y^{*}\right)-\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^2}\\ Det\left[J\left(E^{*}\right)\right]=\left(\frac{g\rho x^{*}{}^2}{k}+\frac{P{x}^{*}}{\beta +x^{*}}\right)\left(\frac{q\beta y^{*}}{{\left(\beta +x^{*}\right)}^2}-\eta y^{*}\right)\end{array}$

则当

$Tr\left[J\left(E^{*}\right)\right]<0,Det\left[J\left(E^{*}\right)\right]>0$

即条件(H₂)成立时，特征方程(6)的根均具有负实部，进而系统(1)的正常数平衡点 $E^{*}$ 是局部渐近稳定的。反之，特征方程(6)的根不全有负实部且实部不为零，因此正常数平衡点 $E^{*}$ 是不稳定的。

3. Hopf分支的存在性

本节选参数 $\rho$ 来研究系统(1)在正常数平衡点 $E^{*}$ 处的Hopf分支的存在性。

引理1假设条件(H₁) (H₁₁)成立，则存在 ${\rho }^{*}>0$，使得当 $\rho ={\rho }^{*}$ 时

$g-\frac{2g{x}^{*}}{k}\left(1+{\rho }^{*}{y}^{*}\right)-\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^2}=0$ (7)

证明 令 $g-\frac{2g{x}^{*}}{k}\left(1+\rho y^{*}\right)-\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^2}=0$。整理得

${\rho }^{*}=\frac{P\left(k-\beta -2x^{*}\right)}{g{\left(\beta +x^{*}\right)}^2}$ (8)

且 $\left(k-\beta -2x^{*}>0\right)$，因此，存在 ${\rho }^{*}>0$，使得(7)成立，且 ${\rho }^{*}>0$ 是唯一的。

引理2若条件 $\rho ={\rho }^{*}$，则模型(1)在 $E^{*}$ 处产生Hopf分支。

证明计算可知：

${\lambda }_{1,2}=\pm i\sqrt{D\text{e}t\left[J\left(E^{*}\right)\right]}$ (9)

下面验证横截性条件。即如果 ${\lambda }_{1,2}=\alpha \left(\rho \right)\pm i\omega \left(\rho \right)$ 是特征方程(6)的根，则

$\alpha \left(\rho \right)={\frac{1}{2}Tr\left[J\left(E^{*}\right)\right]|}_{\rho ={\rho }^{*}}=0$

${\alpha }^{\prime }\left(\rho \right)={\frac{\text{d}}{\text{d}\rho }Tr\left[J\left(E^{*}\right)\right]|}_{\rho ={\rho }^{*}}=-\frac{g^{2}k{x}^{*}\left(k-x^{*}\right)\left(\beta +2x^{*}\right)}{P{\left(\frac{\left(k-2x^{*}-\beta \right)x^{*}}{\beta +x^{*}}+k\right)}^2}<0$

表明满足横截性条件。因此，由文献 [14] 的Poincaré-Andronov-Hopf分支定理可得，当 $\rho$ 经过 ${\rho }^{*}$ 时，系统(1)在 $E^{*}$ 处产生Hopf分支。

4. Hopf分支的方向与稳定性

这一节我们主要研究当 $\rho ={\rho }^{*}$ 时，模型(1)在 $E^{*}$ 处发生Hopf分支的方向及由分支产生的周期解的稳定性。作变换 $\bar{x}=x-x^{*},\bar{y}=y-y^{*}$，则平衡点 $E^{*}=\left(x^{*},y^{*}\right)$ 平移到原点，为了方便，变换后仍用x和y表示 $\bar{x}$ 和 $\bar{y}$，则模型变为

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=g\left(x+x^{*}\right)\left(1-\frac{x+x^{*}}{\frac{k}{1+\rho (y+y^{*})}}\right)-\frac{P}{\beta +\left(x+x^{*}\right)}\left(x+x^{*}\right)\left(y+y^{*}\right)\\ \frac{\text{d}y}{\text{d}t}=-d\left(y+y^{*}\right)+\frac{q\left(x+x^{*}\right)\left(y+y^{*}\right)}{\beta +\left(x+x^{*}\right)}-\eta \left(x+x^{*}\right)\left(y+y^{*}\right)\end{array}$ (10)

对于模型(10)用Taylor展式可化简为

$\left(\begin{array}{c}\frac{\text{d}x}{\text{d}t}\\ \frac{\text{d}y}{\text{d}t}\end{array}\right)=J\left(E^{*}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}f\left(x,y,\rho \right)\\ g\left(x,y,\rho \right)\end{array}\right)$ (11)

其中

$\begin{array}{l}f\left(x,y,\rho \right)=a_1{x}^2+a_{2}xy+a_3{y}^2+a_4{x}^3+a_5{x}^{2}y+\cdots \\ g\left(x,y,\rho \right)=b_1{x}^2+b_{2}xy+b_3{y}^2+b_4{x}^3+b_5{x}^{2}y+\cdots \end{array}$

并且

$\begin{array}{l}{a}_1=-\frac{g\left(1+\rho y^{*}\right)}{k}+\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^3},a_2=-\frac{2\rho g{x}^{*}}{k}-\frac{P\beta }{{\left(\beta +x^{*}\right)}^2},a_3=b_3=0,\\ a_4=-\frac{P\beta y^{*}}{{\left(\beta +x^{*}\right)}^4},a_5=-\frac{\rho g}{k}+\frac{P\beta }{{\left(\beta +x^{*}\right)}^3},b_1=-\frac{q\beta y^{*}}{{\left(\beta +x^{*}\right)}^3},\\ b_2=\frac{q\beta }{{\left(\beta +x^{*}\right)}^2}-\eta ,b_4=\frac{q\beta y^{*}}{{\left(\beta +x^{*}\right)}^4},b_5=-\frac{q\beta }{{\left(\beta +x^{*}\right)}^3}\end{array}$

定义矩阵

$T=\left(\begin{array}{cc}N& 1\\ M& 0\end{array}\right)$

其中， $N=-\frac{\alpha \left(\rho \right)-a_{22}}{\omega \left(\rho \right)}$， $M=-\frac{a_{21}}{\omega \left(\rho \right)}$，可得

$T^{-1}J\left({\rho }^{*}\right)T=\phi \left(\rho \right)=(\begin{array}{cc}\alpha \left(\rho \right)& -\omega \left(\rho \right)\\ \omega \left(\rho \right)& \alpha (\; \rho \; ))\end{array}$

则

$T^{-1}=\left(\begin{array}{cc}0& \frac{1}{M}\\ 1& -\frac{N}{M}\end{array}\right)$

当 $\rho ={\rho }^{*}$ 时，有

$N_0={N|}_{\rho ={\rho }^{*}}=\frac{a_{22}}{\omega \left({\rho }^{*}\right)}=0,M_0={M|}_{\rho ={\rho }^{*}}=-\frac{a_{21}}{\omega \left({\rho }^{*}\right)},\omega \left({\rho }^{*}\right)={\omega \left(\rho \right)|}_{\rho ={\rho }^{*}}$

通过变换 ${\left(x,y\right)}^{\text{T}}=T{\left(u,v\right)}^{\text{T}}$，模型可重新写为

$\left(\begin{array}{c}\frac{\text{d}u}{\text{d}t}\\ \frac{\text{d}v}{\text{d}t}\end{array}\right)=\phi \left(\rho \right)\left(\begin{array}{c}u\\ v\end{array}\right)+\left(\begin{array}{c}{F}^1\left(u,v,\rho \right)\\ F^2\left(u,v,\rho \right)\end{array}\right)$

其中

$\begin{array}{c}{F}^1\left(u,v,\rho \right)=\frac{1}{M}g\left(Nu+v,Mu,\rho \right)\\ =\left(\frac{N^2}{M}{b}_1+N{b}_2\right)u^2+\left(\frac{2N}{M}{b}_1+b_2\right)uv+\frac{b_1}{M}{v}^2+\left(\frac{N^3}{M}{b}_4+N^2{b}_5\right)u^3\\ +\left(\frac{3N^2}{M}{b}_4+2N{b}_5\right)u^{2}v+\left(\frac{3N}{M}{b}_4+b_5\right)u{v}^2+\frac{b_4}{M}{v}^3\end{array}$

$\begin{array}{c}{F}^2\left(u,v,\rho \right)=f\left(Nu+v,Mu,\rho \right)-\frac{N}{M}g\left(Nu+v,Mu,\rho \right)\\ =\left(N{a}_1+M{a}_2-\frac{N^2}{M}{b}_1-N{b}_2\right)N{u}^2+\left(2N{a}_1+M{a}_2-\frac{2N^2}{M}{b}_1-N{b}_2\right)uv\\ +\left(a_1-\frac{N}{M}{b}_1\right)v^2+\left(N{a}_4+M{a}_5-\frac{N^2}{M}{b}_4-N{b}_5\right)N^2{u}^3\\ +\left(3N{a}_4+2M{a}_5-\frac{3N^2}{M}{b}_4-2N{b}_5\right)N{u}^{2}v\\ +\left(3N{a}_4+M{a}_5-\frac{3N^2}{M}{b}_4-N{b}_5\right)u{v}^2+\left(a_4-\frac{N}{M}{b}_4\right)v^3\end{array}$

再用极坐标 $u=r\cos \theta ,v=r\sin \theta$ 系统等价于

$\{\begin{array}{l}\stackrel{\dot{}}{r}=\alpha \left(\rho \right)r+a\left(\rho \right)r^3+\cdots \\ \stackrel{\dot{}}{\theta }=\omega \left(\rho \right)+c\left(\rho \right)r^2+\cdots \end{array}$

对于模型在 $\rho ={\rho }^{*}$ 处作Taylor展式有

$\{\begin{array}{l}\stackrel{\dot{}}{r}={\alpha }^{\prime }\left({\rho }^{*}\right)\left(\rho -{\rho }^{*}\right)r+a\left({\rho }^{*}\right)r^3+o\left({\left(\rho -{\rho }^{*}\right)}^{2}r,\left(\rho -{\rho }^{*}\right)r^3,r^5\right)+\cdots \\ \stackrel{\dot{}}{\theta }=\omega \left({\rho }^{*}\right)+{\omega }^{\prime }\left({\rho }^{*}\right)\left(\rho -{\rho }^{*}\right)+c\left({\rho }^{*}\right)r^2+o\left({\left(\rho -{\rho }^{*}\right)}^2,\left(\rho -{\rho }^{*}\right)r^2,r^4\right)+\cdots \end{array}$

为了得到Hopf分支方向和分支周期解的稳定性，需要计算Liapunov系数 $h\left({\rho }^{*}\right)$ 的符号：

$h\left({\rho }^{*}\right):=\frac{1}{16}\left(F_{uuu}^1+F_{uvv}^1+F_{uuv}^2+F_{vvv}^2\right)+\frac{1}{16\omega \left({\rho }^{*}\right)}\left[F_{uv}^1\left(F_{uu}^1+F_{vv}^1\right)-F_{uv}^2\left(F_{uu}^2+F_{vv}^2\right)-F_{uu}^1{F}_{uu}^2+F_{vv}^1{F}_{vv}^2\right]$

所有的偏导数都取值于分支点 $\left(u,v,\rho \right)=\left(u,v,{\rho }^{*}\right)$，并且

$F_{uuu}^1\left(0,0,{\rho }^{*}\right)=6\left(\frac{N_0^3}{M_0}{b}_4+N_0^2{b}_5\right)$

$F_{uvv}^1\left(0,0,{\rho }^{*}\right)=2\left(\frac{3N_0}{M_0}{b}_4+b_5\right)$

$F_{uuv}^2\left(0,0,{\rho }^{*}\right)=2\left(3N_0^2{a}_4+2N_0{M}_0{a}_4-\frac{3N_0^3}{M_0}{b}_4-2N_0^2{b}_5\right)$

$F_{vvv}^2\left(0,0,{\rho }^{*}\right)=6\left(a_4-\frac{N_0}{M_0}{b}_4\right)$

$F_{uu}^1\left(0,0,{\rho }^{*}\right)=2\left(\frac{N_0^2}{M_0}{b}_1+N_0{b}_2\right)$

$F_{uv}^1\left(0,0,{\rho }^{*}\right)=\frac{2N_0}{M_0}{b}_1+b_2$

$F_{vv}^1\left(0,0,{\rho }^{*}\right)=\frac{2b_1}{M_0}$

$F_{uu}^2\left(0,0,{\rho }^{*}\right)=2N_0\left(N_0{a}_1+M_0{a}_2-\frac{N_0^2}{M_0}{b}_1-N_0{b}_2\right)$

$F_{uv}^2\left(0,0,{\rho }^{*}\right)=2N_0{a}_1+M_0{a}_2-\frac{2N_0^2}{M_0}{b}_1-N_0{b}_2$

$F_{vv}^2\left(0,0,{\rho }^{*}\right)=2\left(a_1-\frac{N_0}{M_0}{b}_1\right)$

故 $h\left({\rho }^{*}\right)=\frac{1}{8}\left(b_5-3b_4\right)+\frac{1}{8\omega \left({\rho }^{*}\right)}\left(\frac{b_1{b}_2}{M_0}-a_1{a}_2{M}_0+\frac{2a_1{b}_2}{M_0}\right)$。

因此，可以得到

${\mu }_2=-\frac{h\left({\rho }^{*}\right)}{{\alpha }^{\prime }(\; \rho \; *\; )}$

由Poincaré-Andronov-Hopf分支定理， ${{\alpha }^{\prime }\left(\rho \right)|}_{\rho ={\rho }^{*}}$ 和上面关于 $h\left({\rho }^{*}\right)$ 的计算可得下列结论。

定理3假设(H₁) (H₂)成立，则当 $\rho ={\rho }^{*}$ 时，系统(1)在正平衡点 $E^{*}=\left(x^{*},y^{*}\right)$ 处产生Hopf分支。若 $h\left({\rho }^{*}\right)<0$ 则Hopf分支方向是亚临界的，分支周期解是渐近稳定的；若 $h\left({\rho }^{*}\right)>0$ 则Hopf分支方向是超临界的，分支周期解是不稳定的。

5. 数值模拟

本节通过数值模拟来验证前面所得到的结果。在系统(1)中选取参数 $g=0.1$， $k=1.6$， $P=0.1$， $\beta =0.1$， $q=0.65$， $d=0.45$， $\eta =0.1$，满足系统有正平衡点 $E^{*}$ 经计算， ${\rho }^{*}\approx 6.64$。当 $\rho =7>{\rho }^{*}$ 时，正平衡点 $E^{*}=\left(0.27732,0.21397\right)$ 是渐近稳定的，“如图1”所示。当 $\rho$ 减小到临界值 ${\rho }^{*}$，正平衡点 $E^{*}=\left(0.27732,0.22402\right)$ 失稳并且经历Hopf分支，进一步计算可得 $h\left({\rho }^{*}\right)\approx -6.3<0$，由定理3，Hopf分支是亚临界的，分支周期解是渐近稳定的。当 $\rho =6<{\rho }^{*}$ 时，正平衡点 $E^{*}=\left(0.27732,0.22402\right)$ 是不稳定的，且系统(1)存在一个稳定的极限环，“如图2”所示。

Figure 1. When $\rho =7>{\rho }^{*}$, the positive equilibrium $E^{*}$ of system (1) is asymptotically stable

图1. 当 $\rho =7>{\rho }^{*}$ 时，系统(1)的正平衡点 $E^{*}$ 是渐近稳定的

Figure 2. When $\rho =6<{\rho }^{*}$, the positive equilibrium $E^{*}$ of system (1) is unstable

图2. 当 $\rho =6<{\rho }^{*}$ 时，系统(1)的正平衡点 $E^{*}$ 是不稳定的