食饵具有扩散和恐惧效应的捕食–食饵模型的稳定性分析

Stability Analysis of Predator-Prey Model with Dispersal and Fear Effect

摘要: 本文提出了一个食饵具有恐惧效应和扩散的捕食–食饵模型，通过定性分析，得到系统详尽的稳定性性态。我们还研究了恐惧效应对系统动力学行为的影响，同时还发现扩散对捕食者和食饵的持久性也有很大的影响。数值模拟进一步验证理论分析的正确性。

Abstract: In this paper, we propose a predator-prey model with fear effect and diffusion. By making full use of qualitative analysis, we obtain the detailed dynamic behavior of the system. We also study the influence of fear effect on the system and find that diffusion has a large effect on the persistence of predator and prey. Numerical simulation further demonstrates the feasibility of our theoretical conclusion.

1. 引言

在现实的生态圈中，由于自然环境以及人类在陆地上的开发，已经使许多生态系统的栖息地支离破碎，连续的种群会破裂成斑块种群。这种碎片化严重威胁着许多濒危物种的持久和生存。种群在斑块间的扩散现象对防止种群灭绝以及种群的持续发展起到了至关重要的作用。20世纪50年代至今，学者们关于扩散在斑块种群作用的研究取得很大进展。1951，Skellem [1] 提出斑块种群扩散纯理论研究。1974年，Levin [2] 建立单种群在n斑块环境下的扩散模型来研究扩散系数对种群的平衡点的稳定性和持续生存的影响。后来，Holt [3] 等学者通过建立单物种扩散模型来研究扩散对物种密度的影响。 [4] [5] [6] [7] [8] 讨论了两斑块坏境下的捕食–食饵扩散模型，给出了使食饵、捕食者共存的最佳斑块分布。2018年Liu [9] 提出了一个食饵庇护所对斑块环境下Leslie-Gower捕食食饵模型

$\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}\tau }=\alpha x_2-\beta x_1+\epsilon \left(r_1-b_1{x}_1-a_{1}y\right)x_1,\\ \frac{\text{d}{x}_2}{\text{d}\tau }=\beta x_1-\alpha x_2+\epsilon \left(-d{x}_2\right),\\ \frac{\text{d}y}{\text{d}\tau }=\epsilon \left(r_2-\frac{a_{2}y}{x_1}\right)y,\end{array}$

研究表明，系统的唯一正平衡点总是全局渐近稳定的，虽然扩散系数对捕食者和食饵的密度有影响，但并不会改变系统的稳定性。常数食饵庇护所对系统的持续生存没有影响，但是增加了食饵平衡密度；在一定条件下，捕食者种群密度会随着躲进庇护所的食饵的数量增加而降低。

在现实的生物种群中，食饵对捕食者的恐惧表现在很多方面，比如栖息地利用、觅食行为、繁殖和心理变化。在 [10] 中，Zantte等人首次提出了以下具有恐惧效应的捕食–食饵模型

$\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\frac{r_{0}u}{1+kv}-du-a{u}^2-puv,\\ \frac{\text{d}v}{\text{d}t}=cpuv-mv,\end{array}$

这里的参数k反映了食饵反捕食者行为的恐惧水平。此外，该作者还考虑了HoIIing-II型功能反应函数的情况，发现恐惧效应可以导致从超临界Hopf分叉到亚临界Hopf分叉的周期振荡。这与以前没有恐惧效应的经典捕食–食饵模型有很大的不同。在 [10] 之后，对食饵考虑恐惧效应的研究引起了广泛关注，详见文献 [11] - [17]。

2. 模型的建立

据笔者所知，至今为止还没有学者研究过具有恐惧效应的捕食食饵斑块模型。实际上，一旦食饵遭受了恐惧效应，它通常会导致斑块之间的迁移。因此，研究具有恐惧效应和扩散的捕食者–食饵系统是非常有意义的。受前人研究成果的启发，本文将考虑如下模型：

$\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=\alpha x_2-\beta x_1+x_1\left(\frac{r_1}{1+k_{1}y}-d_1-b_1{x}_1\right)-a_1{x}_{1}y,\\ \frac{\text{d}{x}_2}{\text{d}t}=\beta x_1-\alpha x_2-d_2{x}_2,\\ \frac{\text{d}y}{\text{d}t}=\left(a_2{x}_1-d_3-b_{2}y\right)y,\end{array}$ (1)

其中 $\alpha ,\beta ,r_1,k_1,b_1,b_2,a_1,a_2,h,k,e$ 均为正常数。

做下述无量纲变换

其中 $l=r_1,g=\frac{b_1}{r_1},\varphi =\frac{\alpha b_1}{r_1^2},\phi =\frac{a_1{r}_1}{b_1}$。仍用 $t,x_1,x_2,y$ 表示 $\tau ,{\bar{x}}_1,{\bar{x}}_2,\bar{y}$。于是得到如下系统

$\begin{array}{l}\frac{\text{d}{x}_1}{\text{d}t}=x_2-a{x}_1+x_1\left(\frac{\text{1}}{1+ky}-d-x_1\right)-x_{1}y,\\ \frac{\text{d}{x}_2}{\text{d}t}=b\left(c{x}_1-x_2\right),\\ \frac{\text{d}y}{\text{d}t}=f\left(x_1-e-hy\right)y,\end{array}$ (2)

这里

可以证明，系统(2)的解在初值条件下是非负有界的。本文将对系统(2)平衡点的存在性和稳定性，以及扩散和恐惧效应对系统的影响进行研究。

3. 平衡点的局部稳定性

本节将讨论系统(2)平衡点的存在性和局部渐近稳定性。首先令(2)式右端为零，可得系统(2)的平衡点满足下列方程组

$\{\begin{array}{l}-a{x}_1+x_2+x_1\left(\frac{1}{1+ky}-d_1-x_1\right)-x_{1}y=0,\\ b\left(c{x}_1-x_2\right)=0,\\ f\left(x-e-hy\right)y=0,\end{array}$

通过简单的计算可知，不论参数取何值，系统(2)总是存在灭绝平衡点 $E_0\left(0,0,0\right)$。当 $a+d<c+1$ 时，存在边界平衡点 $E_1\left(x_{11},x_{22},0\right)$，其中 $x_{11}=c+1-a-d,x_{22}=c{x}_{11}$。对于可能的正平衡点，我们关注下列方程的正根

$k\left(h+1\right)y^2+\left(h+1+k\left(a+d+e-c\right)\right)y+a+d+e-c-1=0.$

定义上述方程的判别式

$\begin{array}{c}\Delta ={\left(h+1+k\left(a+d+e-c\right)\right)}^2-4k\left(h+1\right)\left(a+d+e-c-1\right)\\ ={\left(h+1-k\left(a+d+e-c\right)\right)}^2+4k\left(h+1\right)>0.\end{array}$

注意到 $a+d-c>0$，因此可以分析出当 $c+1>a+d+e$ 时，系统(2)存在唯一的正平衡点 $E^{*}\left(x_1^{*},c{x}_1^{*},y^{*}\right)$，其中 $y^{*}=\frac{-\left(h+1+k\left(a+d+e-c\right)\right)+\sqrt{\Delta }}{2k\left(h+1\right)},x_1^{*}=e+h{y}^{*}$。总结上述分析可以得到如下结果。

定理 1

(1) 系统(2)始终存在一个灭绝平衡点 $E_0\left(0,0,0\right)$。

(a) 若 $c+1<a+d$，那么 $E_0\left(0,0,0\right)$ 是局部渐近稳定的；

(b) 若 $c+1=a+d$， $E_0\left(0,0,0\right)$ 是鞍结点；

(c) 若 $c+1>a+d$， $E_0\left(0,0,0\right)$ 是不稳定的。

(2) 当 $c+1>a+d$，那么系统(2)有唯一个边界平衡点 $E_1\left(x_{11},x_{22},0\right)$。

(a) 若 $a+d<c+1<a+d+e$，它是局部渐近稳定的；

(b) 若 $c+1=a+d+e$， $E_1\left(x_{11},x_{22},0\right)$ 是鞍结点；

(c) 若 $c+1>a+d+e$， $E_1\left(x_{11},x_{22},0\right)$ 是不稳定的。

(3) 当 $c+1>a+d+e$，那么系统(2)有唯一一个正平衡点 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$，且是局部渐近稳定的。

证明：系统(2)在任意平衡点 $E\left(x_1,x_2,y\right)$ 的雅可比矩阵为

$J_E=\left(\begin{array}{ccc}-a-d-2x-y+\frac{1}{1+ky}& 1& -x\left(1+k\right)\\ bc& -b& 0\\ fy& 0& f\left(-2hy-e+x\right)\end{array}\right).$

(1) 系统在平衡点 $E_0\left(0,0,0\right)$ 处的雅可比矩阵为

$J_{E_0}=\left(\begin{array}{ccc}-a-d+1& 1& 0\\ bc& -b& 0\\ 0& 0& -fe\end{array}\right).$

$J_{E_0}$ 的特征方程是

$\left(\lambda +fe\right)\left({\lambda }^2+\left(a+b+d-1\right)\lambda +\left(a+d-1-c\right)b\right)=0.$

显然 $J_{E_0}$ 有一个特征值是 ${\lambda }_1=-fe$，另外两个特征值 ${\lambda }_2,{\lambda }_3$ 计算如下

$\begin{array}{l}{\lambda }_2=\frac{-\left(a+b+d-1\right)+\sqrt{{\left(a+b+d-1\right)}^2-4\left(a+d-1-c\right)}}{2},\\ {\lambda }_3=\frac{-\left(a+b+d-1\right)-\sqrt{{\left(a+b+d-1\right)}^2-4\left(a+d-1-c\right)}}{2}.\end{array}$

通过简单的计算可得到如下结果：

当 $c+1<a+d$ 时，有 ${\lambda }_2<0,{\lambda }_3<0$，因此 $E_0\left(0,0,0\right)$ 是局部渐近稳定的；当 $c+1>a+d$ 时，有 ${\lambda }_2>0,{\lambda }_3<0$，因此 $E_0\left(0,0,0\right)$ 是不稳定的；当 $c+1=a+d$，此时 ${\lambda }_2=0,{\lambda }_3<0$。为了确定 $E_0\left(0,0,0\right)$ 的类型，将系统(2)在 $E_0$ 处线性化，并将线性部分对角化。做变换

$\left(\begin{array}{c}{x}_1\\ x_2\\ y\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{c}& 0& -\frac{1}{b}\\ 1& 0& 1\\ 0& 1& 0\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\\ Y\end{array}\right),$

得到如下系统

$\{\begin{array}{l}{\stackrel{\dot{}}{X}}_1=-\frac{b{X}_1^2}{\left(b+c\right)c}-\frac{c{Y}^2}{\left(b+c\right)c}+\frac{2X_{1}Y}{b+c}-\frac{b\left(1+k\right)X_1{X}_2}{b+c}+\frac{c\left(1+k\right)X_{2}Y}{b+c}+O\left({\left|X_1,X_2,Y\right|}^3\right),\\ {\stackrel{\dot{}}{X}}_2=-fe{X}_2-fh{X}_2^2-\frac{f}{b}{X}_{2}Y+\frac{f}{c}{X}_1{X}_2+O\left({\left|X_1,X_2,Y\right|}^3\right),\\ \stackrel{\dot{}}{Y}=-\left(b+c\right)Y+\frac{b{X}_1^2}{\left(b+c\right)c}+\frac{c{Y}^2}{\left(b+c\right)c}-\frac{2X_{1}Y}{b+c}-\frac{c\left(1+k\right)X_{2}Y}{b+c}+\frac{b\left(1+k\right)X_1{X}_2}{b+c}+O\left({\left|X_1,X_2,Y\right|}^3\right),\end{array}$

由于 $-\frac{b}{\left(b+c\right)c}\ne 0$，根据 [18] 可以得到 $E_0\left(0,0,0\right)$ 是鞍结点。

(2) 系统(2)在平衡点 $E_1\left(x_{11},x_{22},0\right)$ 处的雅可比矩阵为

$J_{E_0}=\left(\begin{array}{ccc}-x_{11}-c& 1& -x_{11}\left(1+k\right)\\ bc& -b& 0\\ 0& 0& f\left(x_{11}-e\right)\end{array}\right).$

经过简单计算，可得相应的特征值 ${\lambda }_1,{\lambda }_2,{\lambda }_3$ 如下

$\begin{array}{l}{\lambda }_1=\frac{-\left(b+c+x_{11}\right)+\sqrt{{\left(b+c+x_{11}\right)}^2-4b{x}_{11}}}{2}<0,\\ {\lambda }_2=\frac{-\left(b+c+x_{11}\right)-\sqrt{{\left(b+c+x_{11}\right)}^2-4b{x}_{11}}}{2}<0,\\ {\lambda }_3=f\left(x_{11}-e\right),\end{array}$

当 $c+1-e<a+d<c+1$，有 ${\lambda }_3<0$，故 $E_1\left(x_{11},x_{22},0\right)$ 是局部渐近稳定的；当 $c+1-e>a+d$ 时， ${\lambda }_3>0$，因此 $E_1\left(x_{11},x_{22},0\right)$ 是不稳定的；当 $c+1-e=a+d$ 时，有 ${\lambda }_3=0$。类似于(1)中的方法，先做变换 $\left(x_1,x_2,y\right)=\left(X_1,X_2,Y\right)+\left(e,ce,0\right)$，将 $E_1$ 平移到原点并作泰勒展开得到如下系统

$\{\begin{array}{l}{\stackrel{\dot{}}{X}}_1=-\left(c+e\right)X_1\text{+}{X}_2-e\left(k+1\right)Y-X_1^2-\left(k+1\right)X_{1}Y+e{k}^2{Y}^2+O\left({\left|X_1,X_2,Y\right|}^3\right),\\ {\stackrel{\dot{}}{X}}_2=bc{X}_1-b{X}_2+O\left({\left|X_1,X_2,Y\right|}^3\right),\\ \stackrel{\dot{}}{Y}=f{X}_{1}Y-fh{Y}^2+O\left({\left|X_1,X_2,Y\right|}^3\right),\end{array}$ (3)

再次执行下列变换

$\left(\begin{array}{c}{X}_1\\ X_2\\ Y\end{array}\right)=\left(\begin{array}{ccc}1& 1& -\left(k+1\right)\\ {\lambda }_1+c+e& {\lambda }_2+c+e& -c\left(k+1\right)\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{U}_1\\ U_2\\ V\end{array}\right),$

系统(3)变为

$\{\begin{array}{l}{\stackrel{\dot{}}{U}}_1={\lambda }_1{U}_1+O\left({\left|U_1,U_2,V\right|}^2\right),\\ {\stackrel{\dot{}}{U}}_2={\lambda }_2{U}_2+O\left({\left|U_1,U_2,V\right|}^2\right),\\ \stackrel{\dot{}}{V}=-f\left(h+k+1\right)V^2+f{U}_{1}V+f{U}_{2}V+O\left({\left|U_1,U_2,V\right|}^3\right),\end{array}$

因为 $V^2$ 的系数是 $-f\left(h+k+1\right)\ne 0$，由 [18] 可以分析出 $E_1\left(x_1^{*},c{x}_1^{*},0\right)$ 是鞍结点。

(3) 系统(2)在 $E_2\left(x_1^{*},c{x}_1^{*},y\right)$ 处的雅可比矩阵是

$J_{E_2}=\left(\begin{array}{ccc}-c-x_1^{*}& 1& -x_1^{*}\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)\\ bc& -b& 0\\ f{y}^{*}& 0& -fh{y}^{*}\end{array}\right),$

那么相应的 $J_{E_2}$ 的特征方程是

${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0,$

其中

$\begin{array}{l}{a}_1=fh{y}^{*}+b+c+x_1^{*}>0,\\ a_2=fh{y}^{*}\left(b+c+x_1^{*}\right)+b{x}_1^{*}+f{x}_1^{*}y\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)>0,\\ a_3=bfh{x}_1^{*}{y}^{*}+bf{x}_1^{*}{y}^{*}\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)>0,\end{array}$

那么

$\begin{array}{c}{a}_1{a}_2-a_3=\left(fh{y}^{*}+b+c+x_1^{*}\right)\left(fh{y}^{*}\left(b+c+x_1^{*}\right)+b{x}_1^{*}+f{x}_1^{*}y\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)\right)\\ -bfh{x}_1^{*}{y}^{*}-bf{x}_1^{*}{y}^{*}\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)\\ =\left(b+c+x_1^{*}\right)b{x}_1^{*}+fh{y}^{*}\left(b+c+x_1^{*}\right)\left(b+c+x_1^{*}+fh{y}^{*}\right)\\ +f{x}_1^{*}y\left(1+\frac{k}{{\left(1+k{y}^{*}\right)}^2}\right)\left(fh{y}^{*}+c+x_1^{*}\right)\\ >0.\end{array}$

根据Rourth-Hurwitz判据 [19]，可得 $J_{E_2}$ 的所有特征值的实部均为非负的，因此 $E_2\left(x_2^{*},c{x}_2^{*},y\right)$ 是局部渐近稳定的。定理1证毕。

基于以上的定性分析，将系统(2)平衡点的存在性和稳定性的结果总结在如下的表1中。

Table 1. The existence and stability of equilibrium of system (2)

表1. 系统((2)的平衡点的存在性及稳定性

4. 全局稳定性

定理2 若 $a+d>c+1$，则 $E_0\left(0,0,0\right)$ 是全局渐近稳定的。

证明：构造李雅普诺夫函数

$V=x_1+\frac{1}{b}{x}_2+\frac{1}{f+1}y,$

显然，在第一象限V是恒正的， $V\left(0,0,0\right)=0$。沿着系统(2)的解轨线计算关于t的全导数，得到

$\begin{array}{c}{D}^{+}V\left(t\right)=-x_1^2-\frac{f}{f+1}\left(e+hy\right)y+\left(\frac{1}{1+ky}-a-d+c-\frac{1}{f+1}y\right)x_1\\ \le -x_1^2-\frac{f}{f+1}\left(e+hy\right)y+\left(1+c-a-d-\frac{1}{f+1}y\right)x_1\\ \le 0.\end{array}$

当 $1+c-a-d<0$ 时，对任意的 $x_1,x_2,y$，都有 $D^{+}V\left(t\right)<0$，此外，当且仅当在边界平衡点 $E_0\left(0,0,0\right)$ 处有 $D^{+}V\left(t\right)=0$，这满足李雅普诺夫渐近稳定性定理 [20]，故 $E_0$ 是全局渐近稳定的。

定理3 如果 $a+d<c+1<a+d+e$，并且满足 $k\le \max \left\{2,2h\right\}$，则 $E_1\left(x_{11},x_{22},0\right)$ 是全局渐近稳定的。

证明：构造李雅谱诺夫函数

$V=\frac{1}{x_{22}}\left(x_1-x_1^{*}-x_1^{*}\ln \frac{x_1}{x_1^{*}}\right)+\frac{1}{bc{x}_{11}}\left(x_2-x_2^{*}-x_2^{*}\ln \frac{x_2}{x_2^{*}}\right)+\frac{1}{f{x}_{22}}y,$

则V沿着系统(2)关于t的导数为

$\begin{array}{c}\stackrel{\dot{}}{V}=\frac{1}{x_{22}}\cdot \frac{x_1-x_{11}}{x_1}\cdot {\stackrel{\dot{}}{x}}_1+\frac{1}{bc{x}_{11}}\cdot \frac{x_2-x_{22}}{x_2}\cdot {\stackrel{\dot{}}{x}}_2+\frac{1}{f{x}_{22}}\cdot \stackrel{\dot{}}{y}\\ =\frac{1}{x_{22}}\left(x_1-x_{11}\right)\left(\frac{x_2}{x_1}-\frac{x_{22}}{x_{11}}-\left(x_1-x_{11}\right)+\frac{1}{1+ky}-1-y\right)+\frac{1}{x_{11}}\left(x_2-x_{22}\right)\left(\frac{x_1}{x_2}-\frac{x_{11}}{x_{22}}\right)+\frac{1}{x_{22}}\left(x_1-e-hy\right)y\\ =\left(2-\frac{x_{11}{x}_2}{x_1{x}_{22}}-\frac{x_1{x}_{22}}{x_{11}{x}_2}\right)+\frac{1}{x_{22}}\left(x_{11}-e\right)y-\frac{1}{x_{22}}{\left(x_1-x_{11}\right)}^2-\frac{h}{x_{22}}{y}^2+\frac{k}{x_{22}\left(1+ky\right)}\left(x_1-x_{11}\right)\left(-y\right)\\ \le \frac{k}{2x_{22}\left(1+ky\right)}\left\{{\left(x_1-x_{11}\right)}^2+y^2\right\}-\frac{1}{x_{22}}{\left(x_1-x_{11}\right)}^2-\frac{h}{x_{22}}{y}^2+\left(2-\frac{x_{11}{x}_2}{x_1{x}_{22}}-\frac{x_1{x}_{22}}{x_{11}{x}_2}\right)+\frac{1}{x_{22}}\left(x_{11}-e\right)y\\ \le \frac{k}{2x_{22}}\left\{{\left(x_1-x_{11}\right)}^2+y^2\right\}-\frac{1}{x_{22}}{\left(x_1-x_{11}\right)}^2-\frac{h}{x_{22}}{y}^2+\left(2-\frac{x_{11}{x}_2}{x_1{x}_{22}}-\frac{x_1{x}_{22}}{x_{11}{x}_2}\right)+\frac{1}{x_{22}}\left(x_{11}-e\right)y\\ =\frac{1}{x_{22}}\left\{\left(\frac{k}{2}-1\right){\left(x_1-x_{11}\right)}^2+\left(\frac{k}{2}-h\right)y^2\right\}+\left(2-\frac{x_{11}{x}_2}{x_1{x}_{22}}-\frac{x_1{x}_{22}}{x_{11}{x}_2}\right)+\frac{1}{x_{22}}\left(x_{11}-e\right)y.\end{array}$

当 $c+1-e<a+d$，即 $x_{11}<e$，又当 $a>0$ 时， $2-a-\frac{1}{a}\le 0$ 恒成立。故对任意的 $x_1,x_2,y$ 都有 $\stackrel{\dot{}}{V}\le \text{0}$，当且仅当在 $E_1\left(x_{11}^{*},c{x}_{11}^{*},0\right)$ 处等号取到。利用李雅普诺夫渐近稳定性定理 [20]，故 $E_1$ 是全局渐近稳定的。

定理4 如果 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 存在，并且满足条件 $k\le \max \left\{2,2h\right\}$，则 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 是全局渐近稳定的。

证明：构造李雅普诺夫函数

$V=k_1\left(x_1-x_1^{*}-x_1^{*}\ln \frac{x_1}{x_1^{*}}\right)+k_2\left(x_2-x_2^{*}-x_2^{*}\ln \frac{x_2}{x_2^{*}}\right)+k_3\left(y-y^{*}-y^{*}\ln \frac{y}{y^{*}}\right)$

沿着系统(2)的解，计算V的导数，可得

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\frac{\partial V}{\partial x_1}\cdot \frac{\text{d}{x}_1}{\text{d}t}+\frac{\partial V}{\partial x_2}\cdot \frac{\text{d}{x}_2}{\text{d}t}+\frac{\partial V}{\partial y}\cdot \frac{\text{d}y}{\text{d}t}\\ =k_1\left(x_1-x_1^{*}\right)\left(\frac{x_2}{x_1}-a+\frac{1}{1+ky}-d-x_1-y\right)+k_2\left(x_2-x_2^{*}\right)\cdot \left(bc\cdot \frac{x_1}{x_2}-b\right)\\ +k_3\left(y-y^{*}\right)\cdot f\cdot \left(x_1-e-hy\right)\\ =k_1\left(x_1-x_1^{*}\right)\left(\frac{x_2}{x_1}-\frac{x_2^{*}}{x_1^{*}}\right)+k_2\left(x_2-x_2^{*}\right)\cdot b\cdot c\cdot \left(\frac{x_1}{x_2}-\frac{x_1^{*}}{x_2^{*}}\right)\\ +k_1\left(x_1-x_1^{*}\right)\left(\frac{1}{1+ky}-\frac{1}{1+k{y}^{*}}+x_1^{*}-x_1+y^{*}-y\right)+k_{3}f\left(y-y^{*}\right)\cdot \left(x_1-x_1^{*}-h\left(y-y^{*}\right)\right)\end{array}$

这里取 $k_1=\frac{1}{x_2^{*}},k_2=\frac{1}{bc{x}_1^{*}},k_3=\frac{1}{f{x}_2^{*}},k\le \max \left\{2,2h\right\}$。又当 $a>0$ 时， $2-a-\frac{1}{a}\le 0$ 恒成立，得到

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\left(x_1-x_1^{*}\right)\left(\frac{x_2}{x_2^{*}{x}_1}-\frac{1}{x_1^{*}}\right)+\left(x_2-x_2^{*}\right)\left(\frac{x_1}{x_1^{*}{x}_2}-\frac{1}{x_2^{*}}\right)-\frac{{\left(x_1-x_1^{*}\right)}^2}{x_2^{*}}+\frac{k}{x_2^{*}}\frac{\left(x_1-x_1^{*}\right)\left(y^{*}-y\right)}{\left(1+ky\right)\left(1+k{y}^{*}\right)}-\frac{h}{x_2^{*}}{\left(y-y^{*}\right)}^2\\ =\left(2-\frac{x_1^{*}{x}_2}{x_2^{*}{x}_1}-\frac{x_2^{*}{x}_1}{x_1^{*}{x}_2}\right)+\frac{1}{x_2^{*}}\left(-{\left(x_1-x_1^{*}\right)}^2+\frac{k\left(x_1-x_1^{*}\right)\left(y^{*}-y\right)}{\left(1+ky\right)\left(1+k{y}^{*}\right)}-h{\left(y-y^{*}\right)}^2\right)\\ \le \frac{1}{x_2^{*}}\left\{\left(\frac{k}{2}-1\right){\left(x_1-x_1^{*}\right)}^2+\left(\frac{k}{2}-h\right){\left(y-y^{*}\right)}^2\right\},\end{array}$

显然，对任何 $x_1>0,x_2>0,y>0$，除了正平衡点 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 处 $\frac{\text{d}V}{\text{d}t}=0$，其它地方均是 $\frac{\text{d}V}{\text{d}t}<0$，利用李雅普诺夫渐近稳定定理，此时系统(2)在 $E_2\left(x_2^{*},c{x}_2^{*},y^{*}\right)$ 处是全局渐近稳定的。

5. 扩散和恐惧效应对系统(1)的影响

在这部分，我们通过分析稳定的正平衡点 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 得到恐惧效应和扩散对系统(2)中食饵与捕食者种群密度(TPA)的影响。利用无量纲变换，还原到系统(1)中，分析出扩散和恐惧效应对种群密度的影响。我们用 $T_1$ 表示在 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 处具有扩散的食饵种群密度， ${\tilde{T}}_1,\tilde{y}$ 分别表示系统(1)对应的正平衡点存在时食饵与捕食者的种群密度，定义 $T_1=x_1^{*}+c{x}_1^{*}$。

5.1. 恐惧效应对食饵与捕食者密度的影响

通过一些简单的计算，我们得到

$\begin{array}{l}\frac{\text{d}\tilde{y}}{\text{d}{k}_1}=\frac{b_1}{r_1{a}_1}\cdot \frac{\text{d}{y}^{*}}{\text{d}{k}_1}=\frac{\text{d}{y}^{*}}{\text{d}k}=\frac{\sqrt{\Delta }+\left(a+d+e-c-2\right)k-h-1}{2\sqrt{\Delta }{k}^2},\\ \frac{\text{d}{\tilde{T}}_1}{\text{d}{k}_1}=\frac{r_1}{b_1}\cdot \frac{\text{d}{T}_1}{\text{d}{k}_1}=\frac{r_1^2{a}_1}{b_1^2}\cdot \frac{\text{d}{T}_1}{\text{d}k}=\frac{r_1^2{a}_1}{b_1^2}\cdot \frac{\text{d}\left(x_1^{*}+c{x}_1^{*}\right)}{\text{d}k}\\ =\frac{r_1^2{a}_1}{b_1^2}\cdot \frac{h\left(c+1\right)\left(\sqrt{\Delta }+\left(a+d+e-c-2\right)k-h-1\right)}{2\sqrt{\Delta }{k}^2},\end{array}$

注意到当 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 存在时，可推出 $h+1+k\left(a+d+e-c\right)<\sqrt{\Delta }<h+1+k$，通过放缩法可得到 $\frac{\text{d}{y}^{*}}{\text{d}{k}_1}<\frac{a+d+e-c-1}{2\sqrt{\Delta }k}<0,\frac{\text{d}{T}_1}{\text{d}{k}_1}<0$。因此，我们可以得出食饵种群密度与捕食者种群密度均随着恐惧效应的增大而减小。

5.2. 扩散对食饵与捕食者密度的影响

接下来，我们分析扩散系数 $\alpha ,\beta$ 对系统(1)的影响。

$\begin{array}{l}\frac{\text{d}\tilde{y}}{\text{d}\beta }=\frac{b_1}{a_1{r}_1}\cdot \frac{\text{d}{y}^{*}}{\text{d}\beta }=\frac{b_1}{a_1{r}_1}\cdot \frac{\text{d}{y}^{*}}{\text{d}a}\cdot \frac{\text{d}a}{\text{d}\beta }=\frac{b_1}{a_1{r}_1^2}\cdot \frac{\left(a+d+e-c\right)k-h-1-\sqrt{\Delta }}{2\sqrt{\Delta }\left(h+1\right)},\\ \frac{\text{d}{\tilde{T}}_1}{\text{d}\beta }=\frac{r_1}{b_1}\cdot \frac{\text{d}{T}_1}{\text{d}\beta }=\frac{r_1}{b_1}\cdot \frac{\text{d}{T}_1}{\text{d}a}\cdot \frac{\text{d}a}{\text{d}\beta }=\frac{1}{b_1}\cdot \frac{h\left(c+1\right)\left(\left(a+d+e-c\right)k-h-1-\sqrt{\Delta }\right)}{2\sqrt{\Delta }\left(h+1\right)},\\ \frac{\text{d}\tilde{y}}{\text{d}\alpha }=\frac{b_1}{a_1{r}_1}\cdot \frac{\text{d}{y}^{*}}{\text{d}\alpha }=\frac{b_1}{a_1{r}_1}\cdot \frac{1}{2\left(h+1\right)}\frac{\beta d_2}{r_1{\left(\alpha +d_2\right)}^2}\cdot \frac{\sqrt{\Delta }+h+1-\left(a+d+e-c\right)k}{\sqrt{\Delta }},\\ \frac{\text{d}{\tilde{T}}_1}{\text{d}\alpha }=\frac{r_1}{b_1}\cdot \frac{\text{d}{T}_1}{\text{d}\alpha }=\frac{r_1}{b_1}\cdot \left(1+c\right)\cdot h\cdot \frac{\text{d}{y}^{*}}{\text{d}\alpha }=\frac{\left(1+c\right)h}{2a_1{r}_1\left(h+1\right)}\frac{\beta d_2}{{\left(\alpha +d_2\right)}^2}\cdot \frac{\sqrt{\Delta }+h+1-\left(a+d+e-c\right)k}{\sqrt{\Delta }},\end{array}$

当 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 存在时，可推出条件 $h+1+k\left(a+d+e-c\right)<\sqrt{\Delta }<h+1+k$ 成立，通过放缩法可以得到 $\frac{\text{d}\tilde{y}}{\text{d}\beta }<-\frac{b_1}{a_1{r}_1^2}\cdot \frac{1}{\sqrt{\Delta }}<0,\frac{\text{d}{\tilde{T}}_1}{\text{d}\beta }<0,\frac{\text{d}\tilde{y}}{\text{d}\alpha }>\frac{b_1}{a_1{r}_1}\cdot \frac{\beta d_2}{r_1{\left(\alpha +d_2\right)}^2}\cdot \frac{1}{\sqrt{\Delta }}>0,\frac{\text{d}{\tilde{T}}_1}{\text{d}\alpha }>0$。

因此，食饵(捕食者)种群密度随着 $\alpha$ 的增大而增大，随着 $\beta$ 的增大而减小。

6. 数值模拟

为了验证以上理论分析结果的正确性，运用MATLAB软件进行了数值模拟。通过数值模拟可直观的显示出系统的全局动力学行为以及种群密度随着时间变化而变化的规律。

图1分别展示了模型(2)中灭绝平衡点 $E_0\left(0,0,0\right)$ 的全局渐近稳定性，食饵与捕食者均灭绝。其中参数选取如下 $a=0.85,b=1,c=0.1,d=0.28,e=0.6,f=1,h=1,k=0.2$。

图2分别展示了模型(2)中边界平衡点 $E_1\left(x_{11},c{x}_{11},0\right)$ 的全局渐近稳定性，食饵持续生存，捕食者灭绝。其中参数选取如下 $a=0.65,b=1,c=0.5,d=0.29,e=0.6,f=1,h=1,k=0.2$。

图3分别展示了模型(2)中正平衡点 $E_2\left(x_1^{*},c{x}_1^{*},y^{*}\right)$ 的全局渐近稳定性，我们发现食饵与捕食者均能够持续生存。参数选取如下 $a=0.2,b=1,c=0.5,d=0.1,e=0.1,f=1,h=0.8,k=0.2$。

Figure 1. Global asymptotic stability and the time series of E₀

图1. E₀的全局相图和时间序列图

Figure 2. Global asymptotic stability and the time series of $E_1$ with $k\le \max \left\{2,2h\right\}$

图2. 当 $k\le \max \left\{2,2h\right\}$ 时 $E_1$ 的全局相图和时间序列图

Figure 3. Global asymptotic stability and the time series of $E_2$ with $k\le \max \left\{2,2h\right\}$

图3. 当 $k\le \max \left\{2,2h\right\}$ 时 $E_2$ 的全局相图和时间序列图

7. 结论

本文提出了一个食饵具有恐惧效应和扩散的捕食–食饵模型。主要是利用微分方程的定性和稳定性理论详细地分析了平衡点的存在性和稳定性。此外通过构造合适的Lyapunov函数，证明了平衡点 $E_0,E_1,E_2$ 在一定条件下是全局渐近稳定的。最后，通过数值模拟我们还分析出当恐惧效应增大时，食饵和捕食者的种群密度均减小。也得到了食饵和捕食者的种群密度随着扩散系数 $\alpha$ 的增大而增大，随着 $\beta$ 的增大而减小。以上分析对生态系统可持续发展问题的研究及发展具有一定的理论价值。

致谢

在此特别感谢福建省自然科学基金项目(项目代码：2021J01614，2021J01613)的资助。