一类具有阿利效应的离散特化捕食系统的动力学研究

Dynamical Analysis of a Discrete Specialist Predator System with the Allee Effect

摘要: 基于一个具有阿利效应和Holling I型功能的连续特化捕食模型，首先采用半离散化方法建立了一类离散特化捕食系统，并分析了系统平衡点的性质及其动力学行为。其次，将捕食者种群的食物转化系数作为分岔参数，分析了该参数变化对正平衡点稳定性的影响。研究表明，系统在参数达到临界条件时会发生Neimark-Sacker分岔并出现准周期解。最后，通过数值模拟验证了上述理论结果的有效性，进一步揭示了离散化模型在捕食者–食饵系统研究中的动态特征。

Abstract: Based on a continuous specialist predator model characterized by the Allee effect and a Holling type I functional response, a discrete-time analogue is constructed via a semi-discretization method. The local stability and dynamical properties of the equilibria are systematically analyzed. Considering the predator’s food conversion factor as the bifurcation parameter, the effect of its variation on the stability of the positive equilibrium is analyzed in detail. The analysis demonstrates that a Neimark-Sacker bifurcation occurs when the bifurcation parameter reaches a critical value, leading to the emergence of quasiperiodic solutions. Finally, numerical simulations are further conducted to confirm the theoretical predictions and to highlight the dynamic features inherent in the discrete predator-prey model.

1. 引言

捕食者–食饵系统是数学生态学中经典且富有研究价值的模型之一，核心目标在于揭示捕食者与食饵之间相互作用对种群动态的影响。自Lotka和Volterra提出经典种群动力学模型以来，引入更多符合实际生态过程的机制[1] [2]以反映复杂的生物学现实，如密度依赖增长、功能反应以及种间竞争等[3]。

随着生态学对密度依赖效应理解的深化，Holling提出了功能反应理论[4]，May揭示了简单动力系统中复杂行为的可能性[5]。1931年W. C. Allee提出了一种重要的种群生态现象——Allee效应[6]，Allee指出小种群由于交配受限、防御能力不足等因素[7] [8]可能出现生长率下降的“低密度不利现象”。此后，大量研究表明，Allee效应会引起捕食系统产生灭绝门槛、不稳定性甚至周期振荡，显著影响生态系统的稳定性与可持续性。基于这一生态事实，Das和Mandal构建了一类具有Allee效应的特化捕食连续模型[11]，并提出了一种新的Allee函数形式：

$\alpha \left(X\right)=1-\exp \left(-\overline{a}X\right)$ (1)

该函数满足Allee效应的基本生态约束条件：当食饵密度较低时，增长率近似为零；而随着密度上升，繁殖能力逐渐恢复并趋于稳定。将该函数引入具有Holling I型功能反应的特化型捕食系统，且当猎物保持环境容纳量为$K$ 以及内禀增长率为$r$ 的logistic增长时，可得到以下连续模型：

$\begin{array}{l}\frac{\text{d}X}{\text{d}T}=\overline{r}X\left(1-{\text{e}}^{-\overline{a}X}\right)\left(1-\frac{X}{K}\right)-AXY\\ \frac{\text{d}Y}{\text{d}T}=\overline{m}AXY-DY\end{array}$ (2)

其中，$X$ 与$Y$ 分别表示食饵与捕食者种群，$a$ 为调节Allee效应效果的参数，$m$ 为食物转化系数，$A$ 为捕食者的捕食率，$D$ 为捕食者的内禀死亡率。

在上述模型中，当Allee效应较强，系统可能发生Hopf分岔，产生周期振荡行为；而当Allee效应减弱时，系统趋于稳定平衡。该结果表明Allee效应对生态系统的动态稳定性具有关键调节作用。由于离散系统能够呈现比连续模型更复杂的动力结构，Patra和Maitra [12]研究了双重Allee效应作用下的离散捕食–被捕食系统，离散化不仅可保留系统结构特征，还可能诱发更多复杂动力学现象[13]-[15]。基于上述研究背景，本文将具有Allee效应的特化捕食连续模型离散化，并系统研究了离散模型在正平衡点附近的局部动力学性质及Neimark-Sacker分岔发生的条件。通过平衡点分析、雅可比矩阵线性化、正则形归约及数值模拟等手段，揭示Allee效应与系统参数对离散捕食者–食饵系统复杂动力学的影响规律。

2. 模型建立及局部稳定性分析

2.1. 模型建立及不动点

通过对系统(2)进行无量纲化，变量可重新缩放为：

$x=\mu X,y=\nu Y,t=\beta T$ (3)

将(3)代入(2)中，可得：

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{\overline{r}}{\beta }x\left(1-{\text{e}}^{-\frac{\overline{a}x}{\mu }}\right)\left(1-\frac{x}{\mu K}\right)-\frac{A}{\beta \mu }xy\\ \frac{\text{d}y}{\text{d}t}=\overline{m}\frac{A}{\beta \mu }xy-\frac{D}{\beta }y\end{array}$(4)

对新模型(4)选择归一化参数，可得：

$\mu =\frac{1}{K},\beta =\overline{r},\nu =\frac{A}{\overline{r}}$ (5)

无量纲模型可推导为：

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(1-{\text{e}}^{-ax}\right)\left(1-x\right)-xy\\ \frac{\text{d}y}{\text{d}t}=mxy-dy\end{array}$ (6)

其中$a=\frac{\overline{a}}{\mu },m=\frac{\overline{m}AK}{\overline{r}},d=\frac{D}{\overline{r}}$ 。

采取了半离散化方法，将系统(6)转化为离散系统，设步长$h>0$ ，通过半离散化后，在区间$\left[n,t\right)$ 上，对于$t\to n+1$ ，令

$\begin{array}{l}{x}_{n+1}=x_n+h\left[x_n\left(1-{\text{e}}^{-a{x}_n}\right)\left(1-x_n\right)-x_n{y}_n\right]\\ y_{n+1}=y_n\cdot \exp \left[h\left(m{x}_n-d\right)\right]\end{array}$ (7)

通过观察可得，边界不动点为$E_0\left(0,0\right)$ 和$E_1\left(1,0\right)$ 。考虑到内部不动点，设$x_{n+1}=x_n=x^{*}$ ，$y_{n+1}=y_n=y^{*}$ ，则系统(7)的第二个方程可变为$y^{*}=y^{*}\exp \left(h\left(m{x}^{*}-d\right)\right)$ 。当$y^{*}\ne 0$ 时$\exp \left(h\left(m{x}^{*}-d\right)\right)=0$ ，得到：

$x^{*}=\frac{d}{m}\left(y^{*}\ne 0\right)$ (8)

记$g\left(x\right)=\left(1-{\text{e}}^{-ax}\right)\left(1-x\right)$ ，$y^{*}=g\left(\frac{d}{m}\right)=\left(1-{\text{e}}^{-\frac{ad}{m}}\right)\left(1-\frac{d}{m}\right)$ ，可以得到系统有两个边界不动点$E_0\left(0,0\right)$ ，$E_1\left(1,0\right)$ 和一个内部不动点$E^{*}=\left(\frac{d}{m},\left(1-{\text{e}}^{-\frac{ad}{m}}\right)\left(1-\frac{d}{m}\right)\right)$ 。

2.2. 局部稳定性分析

记$F_1=x+h\left(xg\left(x\right)-xy\right)$ ，$F_2=y\exp \left(h\left(mx-d\right)\right)$ ，求偏导得到一般形式的雅可比矩阵为：

$J\left(x,y\right)=\left(\begin{array}{cc}1+h\left(g+x{g}^{\prime }-y\right)& -hx\\ hmy{\text{e}}^{h\left(mx-d\right)}& {\text{e}}^{h\left(mx-d\right)}\end{array}\right)$ (9)

定理1边界不动点$E_0\left(0,0\right)$ 为非双曲不动点。

证明 在$E_0\left(0,0\right)$ 处有$g\left(0\right)=0$ ，于是$J\left(E_0\right)=\left(\begin{array}{cc}1& 0\\ 0& {\text{e}}^{-hd}\end{array}\right)$ ，特征值${\lambda }_1=1$ ，${\lambda }_2={\text{e}}^{-hd}\in \left(0,1\right)$ ，因此$E_0\left(0,0\right)$ 为非双曲不动点。

定理2

(i) 当$m<d$ 且$h<2/\left(1-{\text{e}}^{-a}\right)$ 时，$E_1$ 为吸引点；

(ii) 当$m>d$ 且$\left|{\lambda }_1\right|<1$ 时，$E_1$ 为鞍点；

(iii) 当$m=d$ 且${\lambda }_1=\pm 1$ 时，$E_1$ 为非双曲不动点。

证明 在$E_1\left(1,0\right)$ 处，

$J\left(E_1\right)=\left(\begin{array}{cc}1+h\left({\text{e}}^{-a}-1\right)& -h\\ 0& {\text{e}}^{h\left(m-d\right)}\end{array}\right)$ (10)

两特征值为对角元${\lambda }_1=1+h\left({\text{e}}^{-a}-1\right)$ ，${\lambda }_2={\text{e}}^{h\left(m-d\right)}$ 。

① 若$\left|{\lambda }_1\right|<1$ 且$\left|{\lambda }_2\right|<1$ ，则$E_1$ 是局部吸引点。

由$\left(1-{\text{e}}^{-a}\right)\in \left(0,1\right)$ 则${\lambda }_1<1$ 对任意$h>0$ 都成立，而通常只需要保证${\lambda }_1>-1$ 即$1-h\left(1-{\text{e}}^{-a}\right)>-1$ ，可以推得$h<2/\left(1-{\text{e}}^{-a}\right)$ ，通常取步长小的满足此不等式，使$\left|{\lambda }_1\right|<1$ ，而$\left|{\lambda }_2\right|<1$ 等价于${\text{e}}^{h\left(m-d\right)}<1$ ，即$m<d$ 。

② 若$m>d$ 即$\left|{\lambda }_2\right|>1$ ，而通常对于$h$ 有$\left|{\lambda }_1\right|<1$ ，则$E_1$ 是鞍点。

③ 若$m=d$ 或$h=\frac{2}{1-{\text{e}}^{-a}}$ 时，$E_1$ 为非双曲不动点。

定理3

(i) 若$F\left(-1\right)>0$ 且$C<1$ ，$E^{*}$ 为吸引点；

(ii) 若$F\left(-1\right)>0$ 且$C>1$ ，$E^{*}$ 为源点；

(iii) 若$F\left(-1\right)<0$ ，$E^{*}$ 为鞍点；

(iv) 若$-2<B<2$ 且$C=1$ ，$E^{*}$ 为非双曲点或可能产生N-S分岔。

证明 在$E^{*}=\left(x^{*},y^{*}\right)=\left(d/m,g\left(d/m\right)\right)$ 处，求平衡点可得：

$\{\begin{array}{l}0=x^{*}\left(1-{\text{e}}^{-a{x}^{*}}\right)\left(1-x^{*}\right)-x^{*}{y}^{*}\\ 0=m{x}^{*}{y}^{*}-dy\end{array}$ (11)

当$y^{*}=0$ 时，分别代入第一个方程$0=x^{*}\left(1-{\text{e}}^{-a{x}^{*}}\right)\left(1-x^{*}\right)$ 和第二个方程$0=y^{*}\left(m{x}^{*}-d\right)$ ，得到$E_0$ 和$E_1$ ，当$y^{*}>0$ 时必有$m{x}^{*}-d=0$ 可以得到$x^{*}=d/m$ ，代回第一个方程，因为$x^{*}>0$ ，方程可化简为$y^{*}=\left(1-{\text{e}}^{-a{x}^{*}}\right)\left(1-x^{*}\right)$ ，此时：

$J\left(E^{*}\right)=\left(\begin{array}{cc}1+h{x}^{*}{g}^{\prime }\left(x^{*}\right)& -h{x}^{*}\\ hm{y}^{*}& 1\end{array}\right)$ (12)

故$trJ=2+h{x}^{*}{g}^{\prime }\left(x^{*}\right)$ ，$\det J=1+h{x}^{*}{g}^{\prime }\left(x^{*}\right)+h^{2}m{x}^{*}{y}^{*}$ ，特征多项式${\lambda }^2-\left(trJ\right)\lambda +\det J=0$ ，根据Schur判据，$F\left(1\right)=1-trJ+\det J=1-\left(2+h{x}^{*}{g}^{\prime }\left(x^{*}\right)\right)+1+h{x}^{*}{g}^{\prime }\left(x^{*}\right)+h^{2}m{x}^{*}{y}^{*}>0$ 。

3. 分岔分析

将系统(7)重新梳理为：

$\{\begin{array}{l}{x}_{n+1}=F_1\left(x,y\right)=x+h\left(xg\left(x\right)-xy\right)\\ y_{n+1}=F_2\left(x,y\right)=y\exp \left(h\left(mx-d\right)\right)\end{array}$ (13)

其中，$g\left(x\right)=\left(1-{\text{e}}^{-ax}\right)\left(1-x\right)$ ，$E^{*}=\left(x^{*},y^{*}\right)=\left(d/m,g\left(d/m\right)\right)$ ，令$r=x^{*}=d/m,E={\text{e}}^{-ar},G=y^{*}=g\left(x^{*}\right)=\left(1-E\right)\left(1-r\right),G^{\prime }=g^{\prime }\left(r\right),G^{″}=g^{″}\left(r\right)$ ，故$trJ=2+hr{G}^{\prime },\det J=1+hrG+h^{2}mrG$ ，其特征多项式为$F\left(\lambda \right)={\lambda }^2-\left(trJ\right)\lambda +\det J$ 。

为进一步探讨系统在参数变化下的复杂动力学行为，本节选取食物转化系数$m$ 为分岔参数，研究系统在正平衡点$E^{*}$ 处的Neimark-Sacker分岔。

3.1. Neimark-Sacker分岔判据

系统要产生N-S分岔需要满足如下三个条件：

1、特征值条件

当系统的雅可比矩阵具有一对共轭复特征值为${\lambda }_{1,2}=r{\text{e}}^{\pm i\theta }$ ，且满足模长$\left|r\right|=1$ ，即${\left(trJ\right)}^2-4\det \left(J\right)<0$ 且$\det \left(J\right)=1$ ，其中${\left(trJ\right)}^2-4\det \left(J\right)<0$ ，化简可得$G^{\prime }+hmG=0$ ，$m$ 作为分岔参数，求临界$m_c$ ，满足$g^{\prime }\left(\frac{d}{m_c}\right)+h{m}_{c}g\left(\frac{d}{m_c}\right)=0$ 。

2、横截性条件

当分岔参数$m$ 穿过临界值时，特征值模随参数变化的导数不为0，即满足横截性条件：${\frac{\text{d}}{\text{d}m}\left(\det J\right)|}_{m=m_c}\ne 0$ ，化简可得${-h{x}^{*}\left(\frac{G^{\prime }+x^{*}{G}^{″}}{m}+h{x}^{*}{G}^{\prime }\right)|}_{m=m_c}\ne 0$ 。

3、非共振条件

保证特征值不满足低阶共振关系，避免产生次谐波分岔，当$m=m_c$ 时，${\lambda }^k\ne 1\left(k=1,2,3,4\right)$ ，此时，$trJ\left(m_c\right)=2+h\frac{d}{m_c}{g}^{\prime }\left(\frac{d}{m_c}\right)$ ，由此可得：

$\begin{array}{l}k=1:\lambda =1\iff trJ=2k=2:\lambda =-1\iff trJ=-2\\ k=3:\lambda ={\text{e}}^{\pm \frac{2\pi i}{3}}\iff trJ=-1k=4:\lambda =\pm i\iff trJ=0\end{array}$

3.2. 坐标平移与线性化

为便于研究正平衡点附近的局部动力学性质，现将平衡点平移到原点，令$u_n=x_n-x^{*},v_n=y_n-y^{*}$ ，代入系统可得：

$x^{*}+u_{n+1}=\left(x^{*}+u_n\right)+h\left(\left(x^{*}+u_n\right)g\left(x^{*}+u_n\right)-\left(x^{*}+u_n\right)\left(y^{*}+v_n\right)\right)$ (14)

利用平衡点条件$x^{*}g\left(x^{*}\right)-x^{*}{y}^{*}=0$ ，得

$\begin{array}{l}{u}_{n+1}=u_n+h\left[\left(x^{*}+u_n\right)g\left(x^{*}+u_n\right)-\left(x^{*}+u_n\right)g\left(y^{*}+v_n\right)\right]\\ v_{n+1}=v_n+h\left(y^{*}+v_n\right)+\frac{{\text{e}}^{hm{u}_n}-1}{h}\end{array}$ (15)

此时系统可写为

$\left(\begin{array}{c}{u}_{n+1}\\ v_{n+1}\end{array}\right)=J\left(\begin{array}{c}{u}_n\\ v_n\end{array}\right)+B\left(u_n,v_n\right)+C\left(u_n,v_n\right)+\cdots$ (16)

其中$J$ 是雅可比矩阵，$B\left(u_n,v_n\right)$ 包含了所有二阶非线性项，$C\left(u_n,v_n\right)$ 包含了所有三阶非线性项。

当分岔参数$m=m_c$ 时，系统的特征方程具有一对共轭复特征值${\lambda }_{1,2}=r{\text{e}}^{\pm i{\theta }_0}$ ，$\left|{\lambda }_{1,2}\right|=1$ ，令

$J=\left(\begin{array}{cc}{s}_{10}& s_{01}\\ t_{10}& t_{01}\end{array}\right)$

说明系统线性部分仅产生单位圆上的旋转而不引起轨道收敛或发散，取${\lambda }_0={\text{e}}^{i{\theta }_0}$ 表示系统的主特征值。

继续研究非线性影响，需将系统的线性部分转化为旋转形式。设矩阵的特征值为${\lambda }_{1,2}={\eta }_1\pm i{\eta }_2$ ，其中$\left|{\eta }_1+i{\eta }_2\right|=1$ ，考虑可逆矩阵：

$M=\left(\begin{array}{cc}{s}_{01}& 0\\ {\eta }_1-s_{10}& {\eta }_2\end{array}\right)$ ，令$\left(\begin{array}{c}{u}_n\\ v_n\end{array}\right)=M\left(\begin{array}{c}{X}_n\\ Y_n\end{array}\right)=\left(\begin{array}{c}{s}_{01}{X}_n\\ \left({\eta }_1-s_{10}\right)X_n+{\eta }_2{Y}_n\end{array}\right)$

其中$2{\eta }_1=s_{10}+t_{01}$ ，${\eta }_1^2+{\eta }_2^2=s_{10}{t}_{01}-s_{01}{t}_{10}$ ，可得${\eta }_1=1+hr{g}^{\prime }/2$ ，判别式$\Delta ={\left(trJ\right)}^2-4\det J=h^2\left(r^2{g^{\prime }}^2-4mrg\right)$ ，当$\Delta <0$ 时，${\eta }_2=\sqrt{-\Delta }/2=h\sqrt{4mrg-r^{2{g^{\prime }}^2}}/2$ ，代入具体数值，得到在正平衡点$E^{*}$ 处的临界条件为${\eta }_1=1-h^{2}dg/2$ ，${\eta }_2=h\sqrt{4dg-h^2{d}^2{g}^2}/2$ 。为进一步研究分岔方向和新生轨道的稳定性，对系统作线性变换：设$p$ 和$q$ 分别为Jacobian矩阵$J$ 在特征值${\lambda }_0$ 对应的右、左特征向量，即满足$Jp={\lambda }_{0}p,J^{\text{T}}q={\lambda }_{0}q$ 。为保证坐标变换的可逆性和分岔常数计算的一致性，需使二者满足归一化条件：$p^{\text{T}}q=1$ 。由$Jq=\lambda q$ 可以推得$\left(J-\lambda I\right)q=0$ ，取$q=\left(\begin{array}{c}-s_{01}\\ s_{10}-\lambda \end{array}\right)$ ，检验$\left(s_{10}-\lambda \right)\left(-s_{01}\right)+s_{01}\left(s_{10}-\lambda \right)=0$ 成立；由$p^{\prime }J=\lambda p^{\prime }$ 可以得$\left(J^{\prime }-\lambda I\right)p=0$ ，取$p=\left(\begin{array}{c}{t}_{10}\\ s_{10}-\lambda \end{array}\right)$ ，检验$p^{\prime }\left(J-\lambda I\right)=\left(\begin{array}{cc}{t}_{10}& s_{10}-\lambda \end{array}\right)\left(\begin{array}{cc}{s}_{10}-\lambda & s_{01}\\ t_{10}& t_{01}-\lambda \end{array}\right)$ ，第一行等于0成立，计算二者的复内积并归一化，即可确定唯一的一组左右特征向量。

计算复内积：${\Delta }_p=p^{\text{T}}q=p_1{q}_1+p_2{q}_2=\left(-t_{10}\right)\left(-s_{01}\right)+\left(t_{01}-\lambda \right)\left(s_{10}-\lambda \right)=t_{10}{s}_{01}+\left(t_{01}-\lambda \right)\left(s_{10}-\lambda \right)$

令$\widehat{q}=q$ ，$\widehat{p}=p/\Delta p$ ，使${\widehat{p}}^{\text{T}}\widehat{q}={\left(p/\Delta p\right)}^{\text{T}}q=1$，经过归一化处理系统在原点附近的线性部分被化为标准旋转形式：

$\left(\begin{array}{c}{u}_{n+1}\\ v_{n+1}\end{array}\right)=J_{{\theta }_0}\left(\begin{array}{c}{u}_n\\ v_n\end{array}\right)+B\left(u_n,v_n\right)+C\left(u_n,v_n\right)+\cdots$ (17)

其中$J_{{\theta }_0}$ 为旋转矩阵。此时线性部分仅产生纯旋转作用，系统的动力学性质由非线性项决定，为后续计算第一Lyapunov系数奠定基础。

3.3. 第一Lyapunov系数与分岔方向

为研究系统在分岔点附近的非线性行为，将(15)式展开至三阶项的泰勒级数，可得到如下模型：

$\begin{array}{l}{u}_{n+1}=s_{10}{u}_n+s_{01}{v}_n+s_{20}{u}_n^2+s_{11}{u}_n{v}_n+s_{02}{v}_n^2+s_{30}{u}_n^3+s_{21}{u}_n^2{v}_n+s_{12}{u}_n{v}_n^2+s_{03}{v}_n^3+O\left(\rho \right)\\ v_{n+1}=t_{10}{u}_n+t_{01}{v}_n+t_{20}{u}_n^2+t_{11}{u}_n{v}_n+t_{02}{v}_n^2+t_{30}{u}_n^3+t_{21}{u}_n^2{v}_n+t_{12}{u}_n{v}_n^2+t_{03}{v}_n^3+O\left(\rho \right)\end{array}$ (18)

$\begin{array}{l}{s}_{10}=1+hr{g}^{\prime },s_{01}=-hr,t_{10}=hmg,t_{01}=1,\\ s_{20}=\frac{h}{2}\left[2g^{\prime }\left(x^{*}\right)+x^{*}g\text{'}\text{'}\left(x^{*}\right)\right],s_{11}=-h,s_{02}=0\\ s_{30}=\frac{h}{6}\left(3g^{″}+x^{*}{g}^{‴}\right),s_{21}=\frac{h}{2}{g}^{″},s_{12}=s_{03}=0\\ t_{20}=\frac{1}{2}{h}^2{m}^2{y}^{*},t_{11}=hm,t_{02}=0\\ t_{30}=\frac{1}{6}{h}^3{m}^3{y}^{*},t_{21}=\frac{1}{2}{h}^2{m}^2,t_{12}=t_{03}=0\end{array}$

为精确计算Neimark-Sacker分岔的第一Lyapunov系数，需要对系统在平衡点处的非线性项进行多线性映射表达。由上文可知，$B\left(u_n,v_n\right)$ 包含所有二阶非线性项，其中$B\left(u,v\right)=\left(\begin{array}{c}{s}_{20}{u}^2+s_{11}uv+s_{02}{v}^2\\ t_{20}{u}^2+t_{11}uv+t_{02}{v}^2\end{array}\right)$ 。接下来与Hessian矩阵建立联系，对于系统中$u_{n+1}$ 方程的非线性部分，$F_1=s_{20}{u}^2+s_{11}uv+s_{02}{v}^2$ ，其Hessian矩阵为

$H_1=\left(\begin{array}{cc}\frac{{\partial }^{2}F}{\partial u^2}& \frac{{\partial }^{2}F}{\partial u\partial v}\\ \frac{{\partial }^{2}F}{\partial v\partial u}& \frac{{\partial }^{2}F}{\partial v^2}\end{array}\right)=\left(\begin{array}{cc}2s_{20}& s_{11}\\ s_{11}& 2s_{02}\end{array}\right)$ (19)

引入二阶多线性映射$D^{2}F\left(\xi ,\eta \right)$ ，

$D^{2}F\left(\xi ,\eta \right)=\left(\begin{array}{c}2s_{20}{\xi }_1{\eta }_1+s_{11}\left({\xi }_1{\eta }_2+{\xi }_2{\eta }_1\right)+2s_{02}{\xi }_2{\eta }_2\\ 2t_{20}{\xi }_1{\eta }_1+t_{11}\left({\xi }_1{\eta }_2+{\xi }_2{\eta }_1\right)+2t_{02}{\xi }_2{\eta }_2\end{array}\right)$ (20)

引入三阶多线性映射$D^{3}F\left(\xi ,\eta ,\varsigma \right)$ ，其中$F_1^{\left(3\right)}=s_{30}{u}^3+s_{21}{u}^{2}v+s_{12}u{v}^2+s_{02}{v}^3$

$\begin{array}{l}{D}^{3}F\left(\xi ,\eta ,\varsigma \right)=6s_{30}{\xi }_1{\eta }_1{\varsigma }_1+\left(2s_{21}\right)\left({\xi }_1{\eta }_1{\varsigma }_2+{\xi }_1{\eta }_2{\varsigma }_1+{\xi }_2{\eta }_1{\varsigma }_1\right)\\ +\left(2s_{12}\right)\left({\xi }_1{\eta }_2{\varsigma }_2+{\xi }_2{\eta }_1{\varsigma }_2+{\xi }_2{\eta }_2{\varsigma }_1\right)+6s_{03}{\xi }_2{\eta }_2{\varsigma }_2\end{array}$ (21)

对于任意二维映射

${\gamma }_{20}=\frac{1}{2}{p}^{\text{T}}{D}^{2}F\left(q,q\right),{\gamma }_{11}=\frac{1}{2}{p}^{\text{T}}{D}^{2}F\left(q,\widehat{q}\right),{\gamma }_{02}=\frac{1}{2}{p}^{\text{T}}{D}^{2}F\left(\widehat{q},\widehat{q}\right),{\gamma }_{21}=\frac{1}{2}{p}^{\text{T}}{D}^{3}F\left(q,q,\widehat{q}\right)$

为了确保N-S分岔发生，代入Lyapunov常数公式：

${\sigma }^{*}=-\mathrm{Re}\left[\frac{\left(1-2\lambda \right){\overline{\lambda }}^2}{1-\lambda }{\gamma }_{11}{\gamma }_{20}\right]-\frac{1}{2}{\left|{\gamma }_{11}\right|}^2-{\left|{\gamma }_{02}\right|}^2+\mathrm{Re}\left(\overline{\lambda }{\gamma }_{21}\right)$ (22)

经计算${\sigma }^{*}\approx 0.1739>0$ ，这表明系统在$m=m_c$ 处经历次临界Neimark-Sacker分岔，当$m<m_c$ 时正平衡点局部渐近稳定，当$m>m_c$ 时平衡点局部失稳，且其邻域出现一条不稳定小闭合曲线。

4. 数值模拟

为了验证理论分析的正确性，并进一步揭示系统在不同参数下的动力学特征，本节对模型进行数值模拟，绘制平衡点变化、分岔图及相图。

选取参数$a=2.0,d=0.3,h=0.1$ ，通过数值计算得到正平衡点$E^{*}=\left(\text{0}\text{.404554},\text{0}\text{.330321}\right)$ ，系统在$E^{*}$ 处发生Neimark-Sacker分岔时对应的捕食率参数$m_c=\text{0}\text{.74155679281}$ 。

4.1. 局部相图验证

为了验证上节中的局部分岔分析，本节进一步构造了临界值两侧的平衡点局部相图。中的(a)和(b)分别展示了$m=m_c-0.01$ 与$m=m_c+0.01$ 时，沿平衡点附近不同方向施加微小扰动后的轨道演化。

如中的(a)所示，小扰动轨道沿螺旋逐渐向内收缩，其半径随迭代逐步减小，最终靠近平衡点，说明此时平衡点为局部渐近稳定焦点。这与在临界值左侧特征值满足$\left|\lambda \right|<1$ 的结果完全一致。如中的(b)所示，小扰动轨道沿螺旋缓慢向外扩散，其半径随迭代逐步增大并最终离开平衡点邻域，说明平衡点在临界值右侧变为局部不稳定焦点，且其附近不存在稳定的小闭合不变曲线。

以上局部相图与第一Lyapunov系数${\sigma }^{*}>0$ 的结论一致，确认本系统在$m=m_c$ 处发生的是次临界Neimark-Sacker分岔：平衡点在分岔点处改变稳定性，产生的是不稳定的局部闭合曲线，系统随参数增大后最终落入的稳定振荡轨道属于远域全局吸引子，而非由正常形直接生成。

Figure 1. Local phase portraits near the equilibrium around the critical parameter ($m=m_c\pm 0.01$ )

图1. 临界参数附近平衡点的局部相图($m=m_c\pm 0.01$ )

为探索Neimark-Sacker分岔之后系统的全局动力学行为，本文在$m\in \left[0.70,1.20\right]$ 范围内计算了最大Lyapunov指数：

Figure 2. Variation of the maximum Lyapunov exponent over the extended range of m

图2. 参数m扩展区间内最大李雅普诺夫指数的变化

结果如显示，${\lambda }_{\max }\left(m\right)$ 在整个区间内均处于10⁻³量级，并围绕零点呈小幅正负波动，未出现显著的大正值区间，说明系统未进入典型的强混沌状态。Lyapunov指数并非严格为负或精确为零，而是在零附近振动，结合相图表明系统在分岔之后表现为准周期振荡及其轻微破裂形成的弱混沌结构。这种“准周期–弱混沌”型演化是离散系统的典型特征，与连续模型中超临界Hopf分岔后始终保持光滑极限环的情景不同。

综上，$m<m_c$ 时螺旋收紧，$m>m_c$ 时螺旋发散，表明平衡点局部稳定性变化与${\sigma }^{*}>0$ 的次临界Neimark-Sacker分岔一致，并在更大参数区间内主要表现为稳定闭合轨道、准周期结构以及弱混沌行为。相较于连续模型中典型的超临界Hopf分岔及其平滑的极限环振荡，离散模型能够在较大参数区间呈现更复杂的动力学结构，这是离散时间系统区别于连续系统的一个重要特征。

4.2. 全局分岔图、相图及连续模型的对比分析

下述三张图为离散系统在参数$m$ 变化下的相图，展示了Neimark-Sacker分岔的发生过程，其中，红点表示正平衡点$E^{*}$ 。如中的(a)所示，当$m<m_c=0.7416$ 时，正平衡点$E^{*}$ 为局部渐近稳定，系统轨迹沿螺旋路径逐渐收敛至平衡点。如图3中的(b)所示，当参数增大至临界值$m=m_c$ 时，平衡点的稳定性发生临界变化，其特征根模长约等于1，系统在平衡点附近出现一族近似闭合的准周期轨道，构成临界不变环，表明系统进入准周期振荡状态，灰度渐变的环形散点表示在临界点时出现的准周期不变环。如图3中的(c)所示，继续增大参数至$m>m_c=0.7416$ 时，平衡点由稳定转为不稳定，系统产生稳定的极限环并吸引附近轨道，其清晰的外部闭合曲线对应稳定的极限环。

Figure 3. Phase portraits of the system for $m<m_c,m=m_c$ and $m>m_c$

图3. 系统在$m<m_c,m=m_c,m>m_c$ 时的相图

这一连续变化过程清晰地揭示了系统从稳定焦点经临界环再到稳定极限环的演化特征，其局部性质对应次临界Neimark-Sacker分岔，即平衡点在临界值右侧失稳，邻域内不存在稳定小环，而稳定极限环位于较远区域。

下列两张图展示了系统在不同捕食效率参数m的分岔行为，其中(a)为x-m平面，(b)为y-m平面：

Figure 4. Bifurcation diagrams of prey and predator versus m

图4. 食饵和捕食者在不同参数m下的分岔图

如中的(a)所示，当$m<m_c\approx 0.7416$ 时，图中表现为平滑的单一分支，即系统收敛到稳定的正平衡点$E^{*}$ ；当$m=m_c$ 时，此时在分岔图中出现分支裂变，平衡点失去稳定性，产生闭合不动轨迹；当$m>m_c$ 时，图中出现的对称“喇叭形”区域系统轨道在平衡点附近出现准周期振荡，表明系统经历了Neimark-Sacker分岔。

中的(b)展示了相应的y-m图，可以观察到捕食者种群密度与食饵密度具有一致的变化趋势：在分岔点前保持稳定，在$m_c$ 之后出现波动。

综上所述，全局分岔图和相图显示，轨道在离开平衡点后在远处形成闭合环、准周期结构与弱混沌吸引子，这与次临界N-S分岔导致“局部失稳–全局吸引子”这一典型动力学图像相吻合，因此，局部分岔分析与全局数值模拟在分岔类型及后续演化路径上一致。

连续模型与离散模型采用了不同的无量纲化方式，因此，本节采用对应参数域的方式对两类模型的动力学结构进行对比，即连续模型的控制参数$\epsilon$ 与离散模型中的$m$ 在动力学作用上对应，均决定系统从稳定平衡向周期振荡的转变。

文献[11]表明，连续模型在$\epsilon ={\epsilon }_c$ 处发生超临界Hopf分岔，分岔后产生光滑稳定的极限环；随着$\epsilon$ 增大，极限环的幅度变化平稳，系统未出现混沌行为。对于离散模型，本文计算得到临界参数$m_c=0.7416$ ，系统在$m>m_c$ 发生次临界Neimark-Sacker分岔。分岔初期产生闭合轨道，与连续模型的极限环结构一致；当参数进一步增大时，离散模型的闭合轨道逐渐膨胀并出现厚度不均、局部破裂等准周期特征；在更大参数区间内，吸引子形成弱混沌结构，与最大李雅普诺夫指数在0附近的小幅正负波动相一致。二者分岔类型不同，但共同导致内部共存平衡失稳，并引发周期或准周期震荡，并且在离散系统中进一步出现了准周期和弱混沌等更复杂的全局动力学行为。

从生态学角度看，上述差异是合理的。连续模型对应连续繁殖、世代重叠的情形，因此动力学行为更为平滑；而离散模型反映季节性繁殖或世代不重叠的生态机制，离散时间结构本身更容易诱发准周期乃至弱混沌振荡。这解释了离散模型在更大参数区间中表现出比连续模型更复杂的动力学行为。

5. 结论

本文以具有阿利效应的连续特化捕食模型为基础，采用半离散化方法构建了一类具有Holling I型功能反应的离散特化捕食系统，首先系统分析了该系统的不动点性质、局部稳定性及分岔行为。通过雅可比矩阵分析与Schur判据，确定了系统在不同参数条件下的平衡特性：当食物转化系数较小时，正平衡点保持稳定；当食物转化系数增大至临界值时，系统在正平衡点处发生Neimark-Sacker分岔，产生准周期轨道或极限环，表现出由稳定态向振荡态的动力学转变。最后用MATLAB软件进行数值模拟，由局部分岔分析、全局相图、最大李雅普诺夫指数与参数扩展扫描的综合验证，本文的理论推导与数值模拟在分岔类型、平衡稳定性及分岔后动力学结构方面均表现出一致性。进一步验证了分岔条件的合理性与离散模型的动力学特征，揭示了离散动力系统在描述实际生态过程中的独特优势。该研究为进一步探讨阿利效应作用下生态系统的动态演化规律提供了新的理论参考。

基金项目

天津市教委科研计划项目(2021KJ009)。

NOTES

^*通讯作者。