近年来, 捕食者-食饵生态经济系统受到数学工作者的广泛关注和深入研究.文献[1-4]研究了建立在微分方程基础上具有线性人为收获的捕食者-食饵系统, 得到了复杂的动力学行为, 如平衡点的稳定性^[1-4], Hopf分支^[1-4], 周期解^[3]等.文献[5]研究了一类具有非线性收获率的捕食者-食饵生态经济系统的稳定性和Hopf分支.它们都是用微分代数方程来描述的.然而, 在捕食者-食饵生态经济系统中, 人们常以等间隔时间周期来统计食饵种群数量, 捕食者种群数量以及收获努力量.对于离散型变量, 差分方程是研究它们之间变化规律的有效方法.本文的研究基于下面具有非线性收获率的捕食者-食饵生态经济系统

$ \begin{equation} \left\{ \begin{aligned} &\dot{x}=x\left(a-kx-y\right)-e\frac{x}{1+mx}, \\ &\dot{y}=y\left(-s+x\right), \\ &0=e\left(p\frac{x}{1+mx}-c\right)-v, \end{aligned} \right. \end{equation} $

(1.1)

这是受限形式的微分代数方程, 其中$x=x\left( t \right)$, $y=y\left(t \right)$分别表示食饵种群的密度, 捕食者种群密度, $e=e\left( t \right)$表示渔民的收获努力量, 在渔业中通常用出动的渔船数或渔网数来度量, 相当于系统中的第二个捕食者. $a$, $s$, $k$为正常数, 分别表示食饵种群的内禀增长率, 捕食者种群的死亡率, 食饵种群的死亡率.第三个方程中, $p$表示单位捕获物的价格, $c$表示单位收获努力的收获成本, $v$表示收获经济利润, $m$为非负参数, 参见文献[5, 6, 7].通过Euler方法将(1.1) 离散化, 得到一个差分代数方程, 再讨论系统不动点的稳定性和Neimark-Sacker分支, 最后通过Matlab数值仿真, 并将结果和系统(1.1) 中的平衡点的稳定性, Hopf分支值以及分支方向进行对比分析.

对(1.1) 式采用Euler方法进行离散化, 得

$ \begin{equation*} \left\{ \begin{aligned} &{{x}_{n+1}}={{x}_{n}}+\delta {{x}_{n}}\left( a-k{{x}_{n}}-{{y}_{n}}-\frac{{{e}_{n}}}{1+m{{x}_{n}}} \right), \\ &{{y}_{n+1}}={{y}_{n}}+\delta {{y}_{n}}(-s+{{x}_{n}}), \\ &v={{e}_{n}}\left( p\frac{{{x}_{n}}}{1+m{{x}_{n}}}-c \right), \end{aligned} \right. \end{equation*} $

其中$\delta$为步长.这是一个差分代数方程, 映射形式为

$ \begin{equation} \left\{ \begin{aligned} &x\mapsto x+\delta x\left( a-kx-y-\frac{e}{1+mx} \right), \\ &y\mapsto y+\delta y(-s+x), \\ &v=e\left( p\frac{x}{1+mx}-c \right). \end{aligned}\right. \end{equation} $

(2.1)

为了方便, 令

$ \begin{eqnarray*} && f(x, y, e)=\left( \begin{matrix} {{f}_{1}}(x, y, e) \\ {{f}_{2}}(x, y, e) \\ \end{matrix} \right), \\ && g(x, y, e)=e\left( p\frac{x}{1+mx}-c \right)-v, \\ && X={{(x, y, e)}^{T}}, \end{eqnarray*} $

其中

$ \begin{eqnarray*}&&{{f}_{1}}(x, y, e)=x+\delta x\left( a-kx-y-\frac{e}{1+mx} \right), \\ &&{{f}_{2}}(x, y, e)=y+\delta y(-s+x).\end{eqnarray*} $

在区域$\text{R}_{+}^{3}=\left\{ (x, y, e)|x>0, y>0, e>0 \right\}$内讨论系统(2.1).类似文献[8, 9], 可以得到当$v>0$时, 系统(2.1) 以$\text{R}_{+}^{3}$中点为初值的任意解具有正不变性和有界性.

系统(2.1) 中, 不动点满足如下方程组

$ \left\{ \begin{aligned} &a-kx-y-\frac{e}{1+mx}=0, \\ &-s+x=0, \\ &v=e\left( p\frac{x}{1+mx}-c \right), \\ \end{aligned}\right. $

容易求得不动点为

$ {{X}_{0}}={{({{x}_{0}}, {{y}_{0}}, {{e}_{0}})}^{T}}={{\left( s, a-ks-\frac{{{e}_{0}}}{1+ms}, \frac{(1+ms)v}{(p-cm)s-c} \right)}^{T}}, $

由实际生态经济意义, 假定

$ a-ks-\frac{{{e}_{0}}}{1+ms}>0, (p-cm)s-c>0.$

(H1)

系统(2.1) 在不动点$X_{0}$处的广义Jacobi矩阵为

$ J({{X}_{0}})=\left( \begin{matrix} 1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)&-{{x}_{0}}\delta &-\frac{{{x}_{0}}\delta }{1+m{{x}_{0}}} \\ {{y}_{0}}\delta &1&0 \\ \frac{{{e}_{0}}p}{{{\left( 1+m{{x}_{0}} \right)}^{2}}}&0&\frac{p{{x}_{0}}}{1+m{{x}_{0}}}-c \\ \end{matrix} \right), $

$J({{X}_{0}})$的相应广义特征方程为

$ \det \left( \begin{matrix} 1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)-\lambda &-{{x}_{0}}\delta &-\frac{{{x}_{0}}\delta }{1+m{{x}_{0}}} \\ {{y}_{0}}\delta &1-\lambda &0 \\ \frac{{{e}_{0}}p}{{{\left( 1+m{{x}_{0}} \right)}^{2}}}&0&\frac{p{{x}_{0}}}{1+m{{x}_{0}}}-c \\ \end{matrix} \right)=0. $

即

$ {{\lambda }^{2}}+P\lambda +Q=0, $

其中

$ \begin{eqnarray*}&&P=-2-{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)-\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}, \\ &&Q=1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}.\end{eqnarray*} $

令

$ F(\lambda )={{\lambda }^{2}}+P\lambda +Q, $

则

$ \begin{eqnarray*}&&F(1)={{x}_{0}}{{y}_{0}}{{\delta }^{2}}, \\ &&F(-1)=4+2{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{2{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}.\end{eqnarray*} $

定义1  假设${{\lambda }_{1}}$和${{\lambda}_{2}}$是系统(2.1) 在正不动点${{X}_{0}}$处广义Jacobi矩阵的广义特征方程的两个根, 则

1) 若$\left| {{\lambda }_{1}} \right|<1$, $\left| {{\lambda}_{2}} \right|<1$, 则${{X}_{0}}$为汇, 且局部渐近稳定;

2) 若$\left| {{\lambda }_{1}} \right|>1$, $\left| {{\lambda}_{2}} \right|>1$, 则${{X}_{0}}$为源, 且不稳定;

3) 若$\left| {{\lambda }_{1}} \right|<1$, $\left| {{\lambda}_{2}} \right|>1$(或$\left| {{\lambda }_{1}} \right|>1$, $\left|{{\lambda }_{2}} \right|<1$), 则${{X}_{0}}$为鞍点;

4) 若$\left| {{\lambda }_{1}} \right|=1$或$\left| {{\lambda}_{2}} \right|=1$, 则${{X}_{0}}$为非双曲不动点.

由定义1看出, 系统不动点的稳定性和Jacobi矩阵的特征根有关.接下来需要以下引理, 参见文献[10].

引理1  令$F(\lambda )={{\lambda }^{2}}+P\lambda +Q$, 若$F(1)>0$, 则

1) $\left| {{\lambda }_{1}} \right|<1$, $\left| {{\lambda }_{2}} \right|<1$, 当且仅当$F(-1)>0, Q<1$;

2) $\left| {{\lambda }_{1}} \right|>1$, $\left| {{\lambda }_{2}} \right|>1$, 当且仅当$F(-1)>0, Q>1$;

3) $\left| {{\lambda }_{1}} \right|<1$, $\left| {{\lambda }_{2}} \right|>1$(或$\left| {{\lambda }_{1}} \right|>1$, $\left| {{\lambda }_{2}} \right|<1$), 当且仅当$F(-1)<0$;

4) ${{\lambda }_{1}}=-1$或$\left| {{\lambda }_{2}} \right|\ne 1$, 当且仅当$F(-1)=0$和$P\ne 0, 2$;

5) ${{\lambda }_{1}}$和${{\lambda }_{2}}$为复数, 且$\left| {{\lambda }_{1}} \right|=\left| {{\lambda }_{2}} \right|=1$, 当且仅当${{P}^{2}}-4Q<0$, $Q=1$.

由定义1和引理1容易得到如下关于不动点的稳定性定理.

定理1  当系统(2.1) 满足(H1) 时, 有唯一的正不动点${{X}_{0}}$, ${{X}_{0}}$为

1) 汇, 当且仅当$2{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{2{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}+4>0$, 且${{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}<0$;

2) 源, 当且仅当$2{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{2{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}+4>0$, 且${{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}>0$;

3) 鞍点, 当且仅当$2{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{2{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}+4<0$;

4) 非双曲不动点, 当且仅当

$ {{P}^{2}}-4Q<0, Q=1\Leftrightarrow v={{v}_{0}}. $

(H2)

事实上, 由$Q=1$, 即

${{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}=0, $

又${{e}_{0}}=\frac{(1+ms)v}{(p-cm)s-c}$, ${{y}_{0}}=a-ks-\frac{{{e}_{0}}}{1+ms}$.取满足(H1), (H2) 的参数值$a$, $c$, $k$, $m$, $s$, $p$, $\delta $, 从而可以找到相应的$v$值, 记为${{v}_{0}}$.

由以上不动点稳定性分析中的定理1, 可以得到下面关于分支存在条件的定理.

定理2  系统(2.1) 当参数$v$在${{v}_{0}}$的小邻域内变化时, 不动点${{X}_{0}}$处可能会出现Neimark-Sacker分支.

取满足(H1), (H2) 的参数, 系统(2.1) 可表为

$\left\{ \begin{aligned} &x\mapsto x+\delta x\left( a-kx-y-\frac{e}{1+mx} \right), \\ &y\mapsto y+\delta y(-s+x), \\ &{{v}_{0}}=e\left( p\frac{x}{1+mx}-c \right). \\ \end{aligned}\right.$

(4.1)

选取${{v}^{*}}$作为变化参数, 系统(4.1) 在微小扰动下变为

$\left\{ \begin{aligned} &x\mapsto x+\delta x\left( a-kx-y-\frac{e}{1+mx} \right), \\ &y\mapsto y+\delta y(-s+x), \\ &{{v}_{0}}+{{v}^{*}}=e\left( p\frac{x}{1+mx}-c \right), \\ \end{aligned}\right.$

(4.2)

其中$\left| {{v}^{*}} \right|\ll 1$为一个很小的扰动参数.

定义系统(2.1) 的局部参数化$\psi $:

$ \begin{eqnarray*}&&X=\psi (Z)={{X}_{0}}+{{U}_{0}}Z+{{V}_{0}}H(Z), \\ &&g(v, \psi (Z))=0, \end{eqnarray*} $

其中

$ {{U}_{0}}=\left( \begin{matrix} 1&0 \\ 0&1 \\ 0&0 \\ \end{matrix} \right), {{V}_{0}}=\left( \begin{matrix} 0 \\ 0 \\ 1 \\ \end{matrix} \right), Z={{({{z}_{1}}, {{z}_{2}})}^{T}}. $

$H:{{\text{R}}^{\text{2}}}\to \text{R}$是一个光滑映射, 参见文献[11, 12].基于上述函数$\psi$的定义, 有

$ D\psi (Z)=\left( \begin{matrix} 1&0 \\ 0&1 \\ -\frac{ep}{(1+mx)((p-cm)x-c)}&0 \\ \end{matrix} \right)=\left( \begin{matrix} {{D}_{{{Z}_{1}}}}\psi (Z),&{{D}_{{{Z}_{2}}}}\psi (Z) \\ \end{matrix} \right), $

系统(2.1) 可降阶变为形如

$ Z\mapsto f(\psi (Z)), Z\in A\subset {{\rm R}^{2}}$

(4.3)

的系统, 其中$A={{\psi }^{-1}}(B({{X}_{0}}))$, $B({{X}_{0}})$为包含${{X}_{0}}$的一个开域.

系统(4.3) 在正不动点${{Z}_{0}}=0$处的Jacobi矩阵为

$ \begin{aligned} {{D}_{Z}}f(\psi (0))&=Df({{X}_{0}})D\psi (0) \\ &=\left( \begin{matrix} 1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right) &-{{x}_{0}}\delta &-\frac{{{x}_{0}}\delta }{1+m{{x}_{0}}} \\ {{y}_{0}}\delta &1&0 \\ \end{matrix} \right)\left( \begin{matrix} Dg(X) \\ U_{0}^{T} \\ \end{matrix} \right)_{x={{x}_{0}}}^{-1}\left( \begin{matrix} 0&0 \\ 1&0 \\ 0&1 \\ \end{matrix} \right) \\ &=\left( \begin{matrix} 1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{(1+m{{x}_{0}})((p-cm){{x}_{0}}-c)}&-{{x}_{0}}\delta \\ {{y}_{0}}\delta &1 \\ \end{matrix} \right), \end{aligned}$

其中隐含的$v$为$v={{v}_{0}}+{{v}^{*}}$.

由(4.2) 式得

$\left\{ \begin{aligned} {{z}_{1}}\mapsto &{{a}_{1}}{{z}_{1}}+{{a}_{2}}{{z}_{2}}+{{a}_{11}}z_{1}^{2}+{{a}_{12}}{{z}_{1}}{{z}_{2}}+{{a}_{22}}z_{2}^{2}+{{a}_{111}}z_{1}^{3}+{{a}_{112}}z_{1}^{2}{{z}_{2}}+{{a}_{122}}{{z}_{1}}z_{2}^{2} \\ &+{{a}_{222}}z_{2}^{3}+O({{(\left| {{z}_{1}} \right|+\left| {{z}_{2}} \right|)}^{4}}), \\ {{z}_{2}}\mapsto &{{b}_{1}}{{z}_{1}}+{{b}_{2}}{{z}_{2}}+{{b}_{11}}z_{1}^{2}+{{b}_{12}}{{z}_{1}}{{z}_{2}}+{{b}_{22}}z_{2}^{2}+{{b}_{111}}z_{1}^{3}+{{b}_{112}}z_{1}^{2}{{z}_{2}}+{{b}_{122}}{{z}_{1}}z_{2}^{2} \\ &+{{b}_{222}}z_{2}^{3}+O({{(\left| {{z}_{1}} \right|+\left| {{z}_{2}} \right|)}^{4}}). \end{aligned}\right.$

(4.4)

为了求出(4.4) 式中的系数, 需要计算

$ \begin{aligned} {{f}_{1{{Z}_{1}}}}(X)&=D{{f}_{1}}(X){{D}_{{{Z}_{1}}}}\psi (Z) \\ &=1+\delta a-2\delta kx-\delta y-\frac{\delta e}{1+mx}+\frac{\delta exm}{{{\left( 1+mx \right)}^{2}}}+\frac{\delta xep}{{{\left( 1+mx \right)}^{2}}\left( \left( p-cm \right)x-c \right)}, \end{aligned}\\ {{f}_{1{{Z}_{2}}}}(X)=D{{f}_{1}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=-\delta x, {{f}_{2{{Z}_{1}}}}(X)=D{{f}_{2}}(X){{D}_{{{Z}_{1}}}}\psi (Z)=\delta y, \\ {{f}_{2{{Z}_{2}}}}(X)=D{{f}_{2}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=1+\delta x-\delta s, \\ {{f}_{1{{Z}_{1}}{{Z}_{1}}}}(X)=D{{f}_{1{{Z}_{1}}}}(X){{D}_{{{Z}_{1}}}}\psi (Z) \\ \begin{aligned} =\frac{2\left(-{{x}^{3}}kc{{m}^{3}}+{{x}^{3}}kp{{m}^{2}}-3{{x}^{2}}kc{{m}^{2}}+2{{x}^{2}}kpm-3xkcm+xkp+ecm-kc-ep\right)\delta}{{{\left(1+mx\right)}^{2}}\left(-xp+xcm+c\right)}, \end{aligned}\\ {{f}_{1{{Z}_{1}}{{Z}_{2}}}}(X)=D{{f}_{1{{Z}_{1}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=-\delta, {{f}_{1{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{1{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, \\ {{f}_{2{{Z}_{1}}{{Z}_{1}}}}(X)=D{{f}_{2{{Z}_{1}}}}(X){{D}_{{{Z}_{1}}}}\psi (Z)=0, {{f}_{2{{Z}_{1}}{{Z}_{2}}}}(X)=D{{f}_{2{{Z}_{1}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=\delta, \\ {{f}_{2{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{2{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, \\ \begin{aligned} {{f}_{1{{Z}_{1}}{{Z}_{1}}{{Z}_{1}}}}(X)&=D{{f}_{1{{Z}_{1}}{{Z}_{1}}}}(X){{D}_{{{Z}_{1}}}}\psi (Z) =-\frac{2\delta e\left( -p+cm \right)\left( 3c{{m}^{2}}x+3cm-3pmx-2p \right)}{{{\left( 1+mx \right)}^{3}}{{\left( -xp+xcm+c \right)}^{2}}}, \end{aligned} \\ {{f}_{1{{Z}_{1}}{{Z}_{1}}{{Z}_{2}}}}(X)=D{{f}_{1{{Z}_{1}}{{Z}_{1}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, {{f}_{1{{Z}_{1}}{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{1{{Z}_{1}}{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, \\ {{f}_{1{{Z}_{2}}{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{1{{Z}_{2}}{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, {{f}_{2{{Z}_{1}}{{Z}_{1}}{{Z}_{1}}}}(X)=D{{f}_{2{{Z}_{1}}{{Z}_{1}}}}(X){{D}_{{{Z}_{1}}}}\psi (Z)=0, \\ {{f}_{2{{Z}_{1}}{{Z}_{1}}{{Z}_{2}}}}(X)=D{{f}_{2{{Z}_{1}}{{Z}_{1}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, {{f}_{2{{Z}_{1}}{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{2{{Z}_{1}}{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0, \\ {{f}_{2{{Z}_{2}}{{Z}_{2}}{{Z}_{2}}}}(X)=D{{f}_{2{{Z}_{2}}{{Z}_{2}}}}(X){{D}_{{{Z}_{2}}}}\psi (Z)=0.$

将${{X}_{0}}$分别代入可得

$ {{f}_{1{{Z}_{1}}}}({{X}_{0}})=1-\delta k{{x}_{0}}+\frac{\delta e{{x}_{0}}m}{{{\left( 1+m{{x}_{0}} \right)}^{2}}}+\frac{\delta {{x}_{0}}{{e}_{0}}p}{{{\left( 1+m{{x}_{0}} \right)}^{2}}\left( \left( p-cm \right){{x}_{0}}-c \right)}, {{f}_{1{{Z}_{2}}}}({{X}_{0}})=-\delta {{x}_{0}}, {{f}_{2{{Z}_{1}}}}({{X}_{0}})=\delta {{y}_{0}}, \\ {{f}_{2{{Z}_{2}}}}({{X}_{0}})=1, \\ {{f}_{1{{Z}_{1}}{{Z}_{1}}}}({{X}_{0}})=\frac{2\left( -{{x}_{0}}^{3}kc{{m}^{3}}+{{x}_{0}}^{3}kp{{m}^{2}}-3{{x}_{0}}^{2}kc{{m}^{2}}+2{{x}_{0}}^{2}kpm-3{{x}_{0}}kcm+{{x}_{0}}kp+{{e}_{0}}cm-kc-{{e}_{0}}p \right)\delta }{{{\left( 1+m{{x}_{0}} \right)}^{2}}\left( -{{x}_{0}}p+{{x}_{0}}cm+c \right)}, \\ {{f}_{1{{Z}_{1}}{{Z}_{2}}}}({{X}_{0}})=-\delta, {{f}_{1{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0, {{f}_{2{{Z}_{1}}{{Z}_{1}}}}({{X}_{0}})=0, {{f}_{2{{Z}_{1}}{{Z}_{2}}}}({{X}_{0}})=\delta, {{f}_{2{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0, \\ {{f}_{1{{Z}_{1}}{{Z}_{1}}{{Z}_{1}}}}({{X}_{0}})=-\frac{2\delta {{e}_{0}}\left( -p+cm \right)\left( 3c{{m}^{2}}{{x}_{0}}+3cm-3pm{{x}_{0}}-2p \right)}{{{\left( 1+m{{x}_{0}} \right)}^{3}}{{\left( -{{x}_{0}}p+{{x}_{0}}cm+c \right)}^{2}}}, {{f}_{1{{Z}_{1}}{{Z}_{1}}{{Z}_{2}}}}({{X}_{0}})=0, {{f}_{1{{Z}_{1}}{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0, \\ {{f}_{1{{Z}_{2}}{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0, {{f}_{2{{Z}_{1}}{{Z}_{1}}{{Z}_{1}}}}({{X}_{0}})=0, \\ {{f}_{2{{Z}_{1}}{{Z}_{1}}{{Z}_{2}}}}({{X}_{0}})=0, {{f}_{2{{Z}_{1}}{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0, {{f}_{2{{Z}_{2}}{{Z}_{2}}{{Z}_{2}}}}({{X}_{0}})=0. $

所以

$ {{a}_{1}}=1-\delta k{{x}_{0}}+\frac{\delta {{e}_{0}}{{x}_{0}}m}{{{\left( 1+m{{x}_{0}} \right)}^{2}}}+\frac{\delta {{x}_{0}}{{e}_{0}}p}{{{\left( 1+m{{x}_{0}} \right)}^{2}}\left( \left( p-cm \right){{x}_{0}}-c \right)}, {{a}_{2}}=-\delta {{x}_{0}}, \\ {{a}_{11}}=\frac{2\left( -{{x}_{0}}^{3}kc{{m}^{3}}+{{x}_{0}}^{3}kp{{m}^{2}}-3{{x}_{0}}^{2}kc{{m}^{2}}+2{{x}_{0}}^{2}kpm-3{{x}_{0}}kcm+{{x}_{0}}kp+{{e}_{0}}cm-kc-{{e}_{0}}p \right)\delta }{{{\left( 1+m{{x}_{0}} \right)}^{2}}\left( -{{x}_{0}}p+{{x}_{0}}cm+c \right)}, \\ {{a}_{12}}=-\delta, {{a}_{22}}=0, {{a}_{111}}=-\frac{2\delta {{e}_{0}}\left( -p+cm \right)\left( 3c{{m}^{2}}{{x}_{0}}+3cm-3pm{{x}_{0}}-2p \right)}{{{\left( 1+m{{x}_{0}} \right)}^{3}}{{\left( -{{x}_{0}}p+{{x}_{0}}cm+c \right)}^{2}}}, \\ {{a}_{112}}=0, {{a}_{122}}=0, {{a}_{222}}=0, {{b}_{1}}=\delta {{y}_{0}}, {{b}_{2}}=1, {{b}_{11}}=0, {{b}_{12}}=\delta, \\ {{b}_{22}}=0, {{b}_{111}}=0, {{b}_{112}}=0, {{b}_{122}}=0, {{b}_{222}}=0. $

系统(4.4) 的线性化部分在$(0, 0)$处的特征方程为${{\lambda}^{2}}+P({{v}^{*}})\lambda +Q({{v}^{*}})=0, $其中

$ \begin{eqnarray*}&&P({{v}^{*}})=-2-{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)-\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}, \\ &&Q({{v}^{*}})=1+{{x}_{0}}\delta \left( -k+\frac{{{e}_{0}}m}{{{(1+m{{x}_{0}})}^{2}}} \right)+\frac{{{x}_{0}}\delta {{e}_{0}}p}{{{(1+m{{x}_{0}})}^{2}}((p-cm){{x}_{0}}-c)}+{{x}_{0}}{{y}_{0}}{{\delta }^{2}}, \\ &&{{y}_{0}}=a-ks-\frac{{{v}_{0}}+{{v}^{*}}}{(p-cm)s-c}, \\ &&{{e}_{0}}=\frac{1+ms}{(p-cm)s-c}({{v}_{0}}+{{v}^{*}}).\end{eqnarray*} $

因此

$ {{\lambda }_{1}}=-\frac{P({{v}^{*}})}{2}+\frac{i}{2}\sqrt{4Q({{v}^{*}})-{{P}^{2}}({{v}^{*}})}=\alpha +i\omega, \\ {{\lambda}_{2}}=-\frac{P({{v}^{*}})}{2}-\frac{i}{2}\sqrt{4Q({{v}^{*}})-{{P}^{2}}({{v}^{*}})}=\alpha -i\omega, $

其中$\alpha =-\frac{P({{v}^{*}})}{2}, \omega=\frac{\sqrt{4Q({{v}^{*}})-{{P}^{2}}({{v}^{*}})}}{2}, $且${{\left| {{\lambda }_{1, 2}}\right|}_{{{v}^{*}}=0}}=\sqrt{Q(0)}=1$, ${{\left. \frac{d\left|{{\lambda }_{1, 2}} \right|}{d{{v}^{*}}}\right|}_{{{v}^{*}}=0}}\ne 0$.当选择参数满足(H1) 和(H2), 且$P(0)\ne 0$时, ${{({{\lambda }_{1, 2}})}^{n}}\ne 1(n=1, 2, 3, 4)$.此时, 系统(4.4) 在其不动点$(0, 0)$处的特征值不在单位圆和坐标轴的交点上, 即系统产生的是Neimark-Sacker分支, 而不是折叠分支或者Flip分支, 参见文[13].

下面求(4.4) 式在${{v}^{*}}=0$处的标准型

$T=\left( \begin{matrix} {{a}_{2}}&0 \\ \alpha -{{a}_{1}}&-\omega \\ \end{matrix} \right)=\left( \begin{matrix} {{a}_{2}}&0 \\ \frac{{{b}_{2}}-{{a}_{1}}}{2}&-\omega \\ \end{matrix} \right)$

是可逆的.通过变换$\left( \begin{matrix} {{z}_{1}} \\ {{z}_{2}} \\ \end{matrix} \right)=T\left( \begin{matrix} u \\ v \\ \end{matrix} \right), $系统(4.1) 变为如下形式

$\left\{ \begin{aligned} &u\mapsto \alpha u-\omega v+{{{\bar{f}}}_{1}}(u, v), \\ &v\mapsto \omega u+av+{{{\bar{f}}}_{2}}(u, v), \\ \end{aligned}\right.$

(4.5)

在(4.5) 式中

$ \begin{eqnarray*}&&{{\bar{f}}_{1}}(u, v)={{\tilde{a}}_{11}}{{u}^{2}}+{{\tilde{a}}_{12}}uv+{{\tilde{a}}_{22}}{{v}^{2}}+{{\tilde{a}}_{111}}{{u}^{3}}+{{\tilde{a}}_{112}}{{u}^{2}}v+{{\tilde{a}}_{122}}u{{v}^{2}}+{{\tilde{a}}_{222}}{{v}^{3}}+O({{(\left| u \right|+\left| v \right|)}^{4}}), \\ &&{{\bar{f}}_{2}}(u, v)={{\tilde{b}}_{11}}{{u}^{2}}+{{\tilde{b}}_{12}}uv+{{\tilde{b}}_{22}}{{v}^{2}}+{{\tilde{b}}_{111}}{{u}^{3}}+{{\tilde{b}}_{112}}{{u}^{2}}v+{{\tilde{b}}_{122}}u{{v}^{2}}+{{\tilde{b}}_{222}}{{v}^{3}}+O({{(\left| u \right|+\left| v \right|)}^{4}}).\end{eqnarray*} $

令

$ \begin{eqnarray*} &&\xi =\frac{\alpha -{{a}_{1}}}{{{a}_{2}}}, {{\tilde{a}}_{11}}={{a}_{2}}({{a}_{11}}+\xi {{a}_{12}}+{{\xi }^{2}}{{a}_{22}}), {{\tilde{a}}_{12}}=-({{a}_{12}}+2\xi {{a}_{22}})\omega, \\ &&{{\tilde{a}}_{22}}=\frac{{{a}_{22}}{{\xi }^{2}}}{{{a}_{2}}}, {{\tilde{a}}_{111}}=a_{2}^{2}{{a}_{111}}+\xi a_{2}^{2}{{a}_{112}}+{{\xi }^{2}}a_{2}^{2}{{a}_{122}}+\xi {{a}_{222}}, \\ &&{{\tilde{a}}_{112}}=-\omega \left( \frac{{{a}_{112}}}{{{a}_{2}}}+2\xi {{a}_{2}}{{a}_{122}}+3{{\xi }^{2}}{{a}_{2}}{{a}_{222}} \right), {{\tilde{a}}_{122}}={{\omega }^{2}}({{a}_{112}}+3\xi {{a}_{222}}), \\ &&{{\tilde{a}}_{222}}=-\frac{{{a}_{222}}{{\xi }^{3}}}{{{a}_{2}}}, {{\tilde{b}}_{11}}=\frac{a_{2}^{2}}{\omega }\left( \xi {{a}_{11}}-{{b}_{11}}+\xi (\xi {{a}_{12}}-{{b}_{12}})+{{\xi }^{2}}(\xi {{a}_{22}}-{{b}_{22}}) \right), \\ &&{{\tilde{b}}_{12}}=-{{a}_{2}}\left( \xi {{a}_{12}}-{{b}_{12}}+2\xi (\xi {{a}_{22}}-{{b}_{22}}) \right), {{\tilde{b}}_{22}}=\omega (\xi {{a}_{22}}-{{b}_{22}}), \\ &&{{\tilde{b}}_{111}}=\frac{{{a}_{2}}}{\omega }\left( a_{2}^{2}(\xi {{a}_{111}}-{{b}_{111}})+\xi a_{2}^{2}(\xi {{a}_{112}}-{{b}_{112}})+{{\xi }^{2}}a_{2}^{2}(\xi {{a}_{122}}-{{b}_{122}})+\xi (\xi {{a}_{222}}-{{b}_{222}}) \right), \\ &&{{\tilde{b}}_{112}}=-\left( \xi {{a}_{112}}-{{b}_{112}}+2\xi a_{2}^{2}(\xi {{a}_{122}}-{{b}_{122}})+3{{\xi }^{2}}a_{2}^{2}(\xi {{a}_{222}}-{{b}_{222}}) \right), \\ &&{{\tilde{b}}_{122}}={{a}_{2}}\omega \left( \xi {{a}_{122}}-{{b}_{122}}+3\xi (\xi {{a}_{222}}-{{b}_{222}}) \right), {{\tilde{b}}_{222}}=-{{\omega }^{2}}(\xi {{a}_{222}}-{{b}_{222}}). \end{eqnarray*} $

由${{\bar{f}}_{1}}(u, v), {{\bar{f}}_{2}}(u, v)$可计算其在$(0, 0)$处的二阶, 三阶导数分别为

$ \begin{eqnarray*}&&{{\bar{f}}_{1uu}}=2{{\tilde{a}}_{11}}, {{\bar{f}}_{1uv}}=2{{\tilde{a}}_{12}}, {{\bar{f}}_{1vv}}=2{{\tilde{a}}_{22}}, {{\bar{f}}_{1uuu}}=6{{\tilde{a}}_{111}}, \\ && {{\bar{f}}_{1uuv}}=6{{\tilde{a}}_{112}}, {{\bar{f}}_{1uvv}}=6{{\tilde{a}}_{122}}, {{\bar{f}}_{1vvv}}=6{{\tilde{a}}_{222}}, \end{eqnarray*} $

$ \begin{eqnarray*} &&{{\bar{f}}_{2uu}}=2{{\tilde{b}}_{11}}, {{\bar{f}}_{2uv}}=2{{\tilde{b}}_{12}}, {{\bar{f}}_{2vv}}=2{{\tilde{b}}_{22}}, \\ &&{{\bar{f}}_{2uuu}}=6{{\tilde{b}}_{111}}, {{\bar{f}}_{2uuv}}=6{{\tilde{b}}_{112}}, {{\bar{f}}_{2uvv}}=6{{\tilde{b}}_{122}}, {{\bar{f}}_{2vvv}}=6{{\tilde{b}}_{222}}.\end{eqnarray*} $

为使系统(4.5) 存在Neimark-Sacker分支, 要求当${{v}^{*}}=0$时下面的判定值不为零

$\beta =\left( -\text{Re}\left( \frac{(1-2{{\lambda }_{1}})\lambda _{2}^{2}}{1-{{\lambda }_{1}}}{{\xi }_{1}}{{\xi }_{2}} \right)-\frac{1}{2}|{{\xi }_{2}}{{|}^{2}}-|{{\xi }_{3}}{{|}^{2}}+\text{Re}({{\lambda }_{2}}{{\xi }_{4}}) \right), $

其中

$ \begin{eqnarray*}&&{{\xi }_{1}}=\frac{1}{8}({{\bar{f}}_{1uu}}-{{\bar{f}}_{1vv}}+2{{\bar{f}}_{2uv}}+i({{\bar{f}}_{2uu}}-{{\bar{f}}_{2vv}}-2{{\bar{f}}_{1uv}})), \\ &&{{\xi }_{2}}=\frac{1}{4}({{\bar{f}}_{1uu}}+{{\bar{f}}_{1vv}}+i({{\bar{f}}_{2uu}}+{{\bar{f}}_{2vv}})), \\ &&{{\xi }_{3}}=\frac{1}{8}({{\bar{f}}_{1uu}}-{{\bar{f}}_{1vv}}-2{{\bar{f}}_{2uv}}+i({{\bar{f}}_{2uu}}-{{\bar{f}}_{2vv}}+2{{\bar{f}}_{1uv}})), \\ &&{{\xi }_{4}}=\frac{1}{16}({{\bar{f}}_{1uuu}}+{{\bar{f}}_{1uvv}}+{{\bar{f}}_{2uuv}}+{{\bar{f}}_{2vvv}}+i({{\bar{f}}_{2uuu}}+{{\bar{f}}_{2uvv}}-{{\bar{f}}_{1uuv}}-{{\bar{f}}_{1vvv}})).\end{eqnarray*} $

由以上分析和Neimark-Sacker分支理论, 参见文[13], 有以下定理

定理3  假设条件(H1) 和(H2) 成立, 且满足$P(0)\ne 0$.若$\beta \ne 0$, 则当参数$v$在${{v}_{0}}$的小邻域内变化时, 系统(2.1) 在正不动点${{X}_{0}}$处存在Neimark-Sacker分支.若$\beta <0~({\rm resp.}~ \beta >0)$, 则当$v>{{v}_{0}}~({\rm resp.}~ v<{{v}_{0}})$时, 系统(2.1) 在正不动点${{X}_{0}}$处产生一个渐近稳定(不稳定)的周期轨道.

利用MATLAB来展示一个具体的例子.

设步长$\delta \text{=0}\text{.0}1$, 并且以文献[5]中相同的初值出发.容易验证它们满足(H1) 和(H2), 那么系统(2.1) 变成

$ \begin{equation*} \left\{ \begin{aligned} &x\mapsto x+0.01x\left( 4-x-y-\frac{e}{1+0.1x} \right), \\ & y\mapsto y+0.01y(-2+x), \\ & v=e\left( \frac{x}{1+0.1x}-1 \right). \\ \end{aligned}\right. \end{equation*} $

(5.1)

计算得出系统(5.1) 存在一个正不动点${{X}_{0}}=(2, 1. 1210767, 1. 05470795)$, 根据定理1可算出${{v}_{0}}=0. 70313863$, 并求得广义Jacobi矩阵的特征值为

${{\lambda }_{1}}\text{=}0.9999+0.01497i, {{\lambda }_{2}}=0.9999-0.01497i, $

且

$ \left| {{\lambda }_{1, 2}} \right|=1, \beta \text{=}-715. 5163<0. $

根据定理3, 当$v>{{v}_{0}}$时, 不动点${{X}_{0}}=(2, 1.1210767, 1.05470795)$是不稳定的, 且有一个渐近稳定的周期轨道; 当$v<{{v}_{0}}$时, 不动点${{X}_{0}}=(2, 1.1210767, 1.05470795)$是渐近稳定的, 且无周期轨道.

图 1展示了当$v=0. 67<{{v}_{0}}$时系统(5.1) 的不动点${{X}_{0}}$是渐近稳定的; 图 3展示了当$v=0. 72074>{{v}_{0}}$时系统(5.1) 的不动点${{X}_{0}}$是不稳定的.说明当经济利润$v$由小变大通过${{v}_{0}}$时, 生态系统(5.1) 将失去稳定性.以上事实说明, 渔民在作业时不能一味地追求经济利润最大化, 否则可能导致生态系统的不平衡, 甚至会引发灾害.渔民们要有意识地保持经济利润在某个合理的区间内. 图 2显示出当$v$在${{v}_{0}}$的小邻域内, 系统在不动点处产生了一个渐近稳定的周期轨道.

与文献[5]中的结果进行比较发现, 当连续系统(1.1) 经过Euler方法离散化得到系统(2.1) 后, 平衡点和分支值都发生了变化.步长$\delta $越小, 则变化越小; 反之步长$\delta $越大, 则变化越大.这里取较小的值$\delta \text{=}0.01$, 计算得平衡点由$(2, 1.111111, 1.0666665)$微小地变化到$(2, 1.1210767, 1.05470795)$, 且分支值由0.711111微小地变化到0.70313863.综合两次实验发现, 系统(1.1) 经过Euler离散化后在其Hopf分支值的小邻域内产生Neimark-Sacker分支, 且系统平衡点的稳定性, 分支方向和产生的极限环稳定性均保持一致, 即系统的Hopf分支性质得以保持.因此Neimark-Sacker分支可以看做映射的Hopf分支.但在实验过程中也发现, 连续系统经过Euler方法离散化后, 必须经过多次迭代(实验中为10000次)才能看出系统的渐近行为, 而原连续系统在较小的时间间隔(实验中为[0,250])就可看出系统的渐近行为, 这是应该注意的问题.