0 引言

时空斑图是在空间或时间上具有某种规律性的非均匀宏观结构, 在化学^[1]、物理学^[2]、神经网络^[3]、生态学^[4]等领域有着重要的应用.特别地, 在生态学中, 为了掌握物种空间分布的多样性, 国内外学者对物种的空间斑图进行了大量而深入的研究^[5-10].文献[5]研究了带有自扩散项的系统中Allee效应对斑图形成的影响.文献[6]讨论了带有交叉扩散项的水生生物模型的斑图的形成.文献[7]研究了带有负交叉扩散项的捕食者-食饵模型中斑图的生成与选择, 并得到了不同类型的空间斑图. Sambath等人研究了以下带有Holling-Ⅱ型功能反应和双曲死亡率的捕食者-食饵模型^[11]

$ \left\{ \begin{array}{l} \dfrac{{{\rm{d}}}U}{{{\rm{d}}}\tau}=rU\Big(1-\dfrac{U}{K}\Big)- \dfrac{a_{1}UV}{n+U}, \\ \dfrac{{{\rm{d}}}V}{{{\rm{d}}}\tau}=\dfrac{a_{2}UV}{n+U}- \dfrac{\delta V^{2}}{\epsilon+\eta V}.\\ \end{array} \right. $

(1)

其中$U(\tau)$和$V(\tau)$分别表示在$\tau$时刻食饵和捕食者的密度. $r$, $K$, $a_{1}$, $n$, $a_{2}$, $\delta$均是正常数, 分别表示食饵的内禀增长率, 食饵的环境容纳量, 捕食者对食饵的捕食率, 捕食者对食饵的最大摄入比率, 捕食者的增长率以及死亡率. $\epsilon$, $\eta$也是正常数, 分别表示在浮游生物死亡的背景下由于水和自我保护导致的光衰减的系数.文献[11]主要讨论了常微分系统的局部稳定性和自扩散效应下Hopf分支的方向以及产生周期解的稳定性.然而, 自然界中捕食者-食饵系统中的生物资源最有可能被收获以获得经济利益, 人类需要开发生物资源和收获一些生物物种, 如渔业, 林业和野生动植物^[12].但收获对种群分布将产生重要的影响.因此, 下面我们考虑的是具有非线性食饵收获效应的捕食者-食饵的反应扩散系统

$ \left\{\!\! \begin{array}{l} \dfrac{\partial U}{\partial\tau}=D_{11}\Delta U+D_{12}\Delta V+rU\Big(1-\dfrac{U}{K}\Big) -\dfrac{a_{1}UV}{n+U}-\dfrac{qEU}{m_{1}E+m_{2}U}, \\ \dfrac{\partial V}{\partial\tau}=D_{21}\Delta U+D_{22}\Delta V+ \dfrac{a_{2}UV}{n+U}-\dfrac{\delta V^{2}}{\epsilon+\eta V}. \end{array} \right. $

(2)

其中$D_{11}, D_{22}$是自扩散系数且是正常数, $D_{12}, D_{21}$是交叉扩散系数, $q$和$E$均是正常数, 分别表示收获能力, 收获努力量, $m_{1}$, $m_{2}$表示合适的正常数, $\Delta\equiv\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}$为拉普拉斯算子.

为了简化模型, 作如下的无量纲化:

$ u=\dfrac{U}{K}, v=\dfrac{a_{1}}{K}V, t=r\tau, h=\dfrac{qE}{rm_{2}K}, \rho=\dfrac{m_{1}E}{m_{2}K}, \beta=\dfrac{n}{K}, s=\dfrac{1}{r}, \delta=a_{2}\eta, \epsilon=\eta, \alpha=\dfrac{a_{2}}{r}, $

$ \gamma=\dfrac{K}{a_{1}}, d_{11}=\dfrac{D_{11}}{r}, d_{12}=\dfrac{D_{12}}{ra_{1}}, d_{21}=\dfrac{a_{1}D_{21}}{r}, d_{22}=\dfrac{D_{22}}{r}. $

这样, 系统(2)就转化为

$ \left\{\!\! \begin{array}{l} \dfrac{\partial u(X, t)}{\partial t}=d_{11}\Delta u+d_{12}\Delta v+u(1-u)-\dfrac{suv}{\beta+u}-\dfrac{hu}{\rho+u}, X\in\Omega, t>0, \\[2mm] \dfrac{\partial v(X, t)}{\partial t}=d_{21}\Delta u+d_{22}\Delta v+ \dfrac{\alpha uv}{\beta+u}-\dfrac{\alpha \gamma v^{2}}{1+\gamma v}, X\in\Omega, t>0, \\[2mm] \dfrac{\partial u}{\partial \nu}=\dfrac{\partial v}{\partial \nu}=0, X\in\partial \Omega, t>0, \\[2mm] u(X, 0)=u_{0}(X)\geq0, v(X, 0)=v_{0}(X)\geq0, X\in\Omega.\\ \end{array} \right. $

(3)

其中$(X, t)\in\Omega\times {\textbf{R}}^{+}$, $\Omega\subset {{\bf R}}^{2}$为边界光滑的有界区域, $\nu$为$\partial\Omega$上单位外法向量.

对于上述系统(3), 文献[13]研究了在自扩散效应下时滞对Hopf分支产生的周期解的稳定性及其方向的影响, 文献[14]讨论了当$h=0$, $d_{12}=d_{21}=0$时, 时滞效应下系统周期解的存在性.然而, 在生态学中, 物种的空间分布规律是重要的研究课题, 特别地, 人类的收获效应对种群的空间斑图的影响尚未被研究.因此本文将讨论系统(3)中非线性食饵收获效应对种群斑图生成的影响.

1 稳定性分析

考虑系统(3)对应的常微分系统

$ \left\{\!\! \begin{array}{l} \dfrac{{{\rm{d}}}u}{{{\rm{d}}}t}=u(1-u)-\dfrac{suv}{\beta+u}- \dfrac{hu}{\rho+u}, \\ \dfrac{{{\rm{d}}}v}{{{\rm{d}}}t}=\dfrac{\alpha uv}{\beta+u}- \dfrac{\alpha \gamma v^{2}}{1+\gamma v}. \end{array} \right. $

(4)

考虑到生态学意义, 本文只研究系统(3)和系统(4)的正平衡点.通过求解, 易知当$h < \rho < \beta$时, 系统(3)和系统(4)存在唯一的正平衡点$E_{\ast}=(u_{\ast}, v_{\ast})$, 其中$v_{\ast}=\frac{u_{\ast}}{\beta\gamma}$且$u_{\ast}$满足如下的三次方程:

$ A_{0}u^{3}+A_{1}u^{2}+A_{2}u+A_{3}=0, $

(5)

其中

$ A_{0}=\beta\gamma, A_{1}=\beta^{2}\gamma+\beta\gamma \rho+s-\beta\gamma, A_{2}=\beta^{2}\gamma\rho+\beta\gamma h+\rho s-\beta^{2}\gamma-\beta\gamma\rho, A_{3}=\beta^{2}\gamma(h-\rho). $

图 1中, 通过食饵零增长的等倾线和捕食者零增长的等倾线的交点来标记正平衡点$E_{\ast}$, 其中参数值$\alpha=1.5$, $\beta=0.8$, $\gamma=0.8$, $\rho=0.4$, $h=0.2$, $s=1$.食饵的等倾线包括$u=0$和曲线$1-u-\frac{sv}{\beta+u}-\frac{h}{\rho+u}=0$, 捕食者的等倾线包括$v=0$和曲线$v=\frac{u}{\beta\gamma}$.

将系统(4)用如下形式来表示:

$ \dot{Y}=F(Y, \Lambda)=(P(u, v), Q(u, v))^{\rm T}, $

(6)

其中$Y=(u, v)^{\rm T}$, $\Lambda$表示所有参数的集合.系统(4)在平衡点$(u_{\ast}, v_{\ast})$处的Jacobi矩阵为

$ J=\left( \begin{array}{ll} P_{u}(u_{\ast}, v_{\ast}) & P_{v}(u_{\ast}, v_{\ast})\\ Q_{u}(u_{\ast}, v_{\ast}) & Q_{v}(u_{\ast}, v_{\ast}) \end{array} \right) =\left( \begin{array}{ll} a_{11} & a_{12}\\ a_{21} & a_{22} \end{array} \right), $

(7)

其中

$ a_{11}=\dfrac{\partial P}{\partial u}\Big|_{E_{\ast}}=-u_{\ast}+ \dfrac{su_{\ast}^{2}}{\beta\gamma(\beta+u_{\ast})^{2}}+\dfrac{hu_{\ast}} {(\rho+u_{\ast})^{2}}, \quad a_{12}=\dfrac{\partial P}{\partial v}\Big|_{E_{\ast}}= -\dfrac{su_{\ast}}{\beta+u_{\ast}}, \\ a_{21}=\dfrac{\partial Q}{\partial u}\Big|_{E_{\ast}}= \dfrac{\alpha u_{\ast}}{\gamma(\beta+u_{\ast})^{2}}, \quad a_{22}= \dfrac{\partial Q}{\partial v}\Big|_{E_{\ast}}= -\dfrac{\alpha\beta u_{\ast}}{(\beta+u_{\ast})^{2}}. $

为了便于下面的讨论, 我们作如下假设.

$H1: \frac{su_{\ast}}{\beta\gamma(\beta+u_{\ast})^{2}}+\frac{h}{(\rho+u_{\ast})^{2}}-\frac{\alpha\beta}{(\beta+u_{\ast})^{2}} < 1;$

$H2: (\frac{su_{\ast}}{\beta\gamma(\beta+u_{\ast})^{2}}+\frac{h}{ (\rho+u_{\ast})^{2}}-1)\beta < \frac{s}{(\beta+u_{\ast})\gamma}; $

$H3: (1-\frac{su_{\ast}}{\beta\gamma(\beta+u_{\ast})^{2}} -\frac{h}{(\rho+u_{\ast})^{2}})d_{22}+\frac{\alpha\beta}{ (\beta+u_{\ast})^{2}}d_{11}>0;$

$H4: (\frac{su_{\ast}}{\beta\gamma(\beta+u_{\ast})^{2}}+\frac{h}{(\rho+u_{\ast})^{2}}-1)d_{22}-\frac{\alpha\beta}{(\beta+u_{\ast})^{2}}d_{11}-\frac{\alpha}{(\beta+u_{\ast})^{2}\gamma}d_{12}+\frac{s}{\beta+u_{\ast}}d_{21}>0.$

定理1.1  对于系统(4), 有如下陈述成立.

(ⅰ)若$H1$和$H2$成立, 则系统的平衡点$(u_{\ast}, v_{\ast})$是局部渐近稳定的;

(ⅱ)若$h=h_{H}$成立, 则系统在平衡点$(u_{\ast}, v_{\ast})$处发生了Hopf分支, 其中

$ h_{H}=\dfrac{\alpha\beta^{2}\gamma(\rho+u_{\ast})^{2}+ \beta\gamma(\beta+u_{\ast})^{2}(\rho+u_{\ast})^{2}-su_{\ast}(\rho+ u_{\ast})^{2}}{\beta\gamma(\beta+u_{\ast})^{2}}. $

证明  (ⅰ)矩阵$J$的特征方程为

$ \lambda^{2}-{\rm tr}_{0}\lambda+{\rm det}J=0, $

(8)

其中,

$ {\rm tr}_{0}=-u_{\ast}+\dfrac{su_{\ast}^{2}}{\beta\gamma(\beta+u_{\ast})^{2}}+ \dfrac{hu_{\ast}}{(\rho+u_{\ast})^{2}}-\dfrac{\alpha\beta u_{\ast}}{(\beta+u_{\ast})^{2}}, $

$ {\rm det}J=\dfrac{\alpha\beta u_{\ast}^{2}}{(\beta+u_{\ast})^{2}}\Big(1-\dfrac{su_{\ast}}{\beta \gamma(\beta+u_{\ast})^{2}}-\dfrac{h}{(\rho+u_{\ast})^{2}}\Big)+ \dfrac{\alpha su_{\ast}^{2}}{(\beta+u_{\ast})^{3}\gamma}, $

tr$_{0}$和det$J$分别表示$J$的迹和行列式.根据$H1$和$H2$, 易知tr$_{0} < 0$和det$J>0$.由Routh-Hurwitz准则知, $(u_{\ast}, v_{\ast})$是局部渐近稳定的.

(ⅱ)通过以上分析, 若$h=h_{H}$成立, 则可以得到一对纯虚数, 使得tr$_{0}=0$, det$J>0$成立.记$\lambda=p+{\rm i}q$并带入到方程(8), 分离实部与虚部得如下方程:

$ \left\{\!\! \begin{array}{l} p^{2}-q^{2}-{\rm tr}_{0}\cdot p+{\rm det}J=0, \\ 2pq-{\rm tr}_{0}\cdot q=0. \end{array} \right. $

(9)

对上述方程组的第二个方程两边关于$h$求导, 可以得到

$ \dfrac{{{\rm{d}}}p}{{\rm d}h}\Big|_{h=h_{H}}=\dfrac{u_{\ast}}{2(\rho+u_{\ast})^{2}}>0. $

因此, 陈述(ⅱ)成立.

下面考虑系统(4)仅加入自扩散时平衡点的稳定性.首先在$(u_{\ast}, v_{\ast})$处引入小扰动

$ \left\{\!\! \begin{array}{l} u(x, y, t)=u_{\ast}+U(x, y, t), |U(x, y, t)|\ll u_{\ast}, \\ v(x, y, t)=v_{\ast}+V(x, y, t), |V(x, y, t)|\ll v_{\ast}, \end{array} \right. $

(10)

其中, $U(x, y, t)=U_{0}{\rm e}^{\lambda t}{\rm e}^{{\rm i}(k_{x}x+k_{y}y)}$, $V(x, y, t)=V_{0}{\rm e}^{\lambda t}{\rm e}^{{\rm i}(k_{x}x+k_{y}y)}$; $\lambda$表示时间$t$处的扰动增长率; i表示虚数单位, $k_{x}, k_{y}$为相应的振幅; $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$表示波数; $U_{0}, V_{0}$为两个实数.得到如下的特征方程:

$ \lambda^{2}+B_{1}\lambda+B_{2}=0, $

(11)

其中$B_{1}=(d_{11}+d_{22})k^{2}-{\rm tr}_{0}$, $B_{2}=d_{11}d_{22} k^{4}-(a_{22}d_{11}+a_{11}d_{22})k^{2}+{\rm det}J$.

由上面的分析可以得到, 若$H1$-$H3$成立, 则系统(3)仅加入自扩散项时的平衡点$(u_{\ast}, v_{\ast})$是局部渐近稳定的.

注1 在系统(4)中仅加入自扩散项得到的偏微分系统, 系统的平衡点$E_{\ast}$总是局部渐近稳定的.

为了得到Turing斑图的形成条件, 接下来我们考虑交叉扩散项对系统动力学的影响.

在介绍定理1.2之前, 作如下假设:

$H5: a_{11}d_{22}+a_{22}d_{11}-a_{21}d_{12}-a_{12}d_{21}-2 \sqrt{(d_{11}d_{22}-d_{12}d_{21}){\rm det}J}>0.$

定理1.2 若$d_{11}d_{22}-d_{12}d_{21}>0$, 当$H1$-$H2$和$H4$-$H5$均成立时, 系统(3)的平衡点$(u_{\ast}, v_{\ast})$是不稳定的.

证明 记$D=\left(\begin{array}{ll} d_{11} & d_{12}\\ d_{21} & d_{22} \end{array} \right)$, $I=\left(\begin{array}{ll} 1 & 0\\ 0 & 1 \end{array} \right)$.将式(10)代入系统(3), 得到如下的特征方程:

$ \mid J-k^{2}D-\lambda I\mid=0, $

(12)

特征方程(12)的解为如下形式:

$ \lambda_{k}=\dfrac{{\rm tr}_{k}\pm\sqrt{{\rm tr}_{k}^{2}-4\delta_{k}}}{2}, $

(13)

其中

$ {\rm tr}_{k}=a_{11}+a_{22}-(d_{11}+d_{22})k^{2}, $

(14)

$ \delta_{k}=(d_{11}d_{22}-d_{12}d_{21})k^{4}-(a_{11}d_{22}+ a_{22}d_{11}-a_{12}d_{21}-a_{21}d_{12})k^{2}+{\rm det}J, $

(15)

${\rm tr}_{k}$和$\delta_{k}$分别表示矩阵$J-k^{2}D$的迹和行列式.令$\lambda_{1k}$和$\lambda_{2k}$为式(12)的解, 则${\rm tr}_{k}=\lambda_{1k}+\lambda_{2k}$和$\delta_{k}= \lambda_{1k}\lambda_{2k}$.显然, 由假设$H1$知, $\lambda_{1k}+ \lambda_{2k} < 0$.若$d_{11}d_{22}-d_{12}d_{21}>0$, 当$H1$-$H2$和$H4$-$H5$均成立时, 存在某些波数$k>0$, 使得$\delta_{k} < 0$, 则至少有一个特征值的实部大于零, 由稳定性理论知, 系统(3)的平衡点$(u_{\ast}, v_{\ast})$是不稳定的.

注2 由上述定理可知, 在$H1$-$H2$和$H4$-$H5$成立的条件下, Turing斑图的形成是由交叉扩散项引起的.

图 2讨论了当$h$变化时, 特征值的实部Re$(\lambda)$与$k$之间的关系, 即参数的变化对Turing不稳定的影响.当$h=0.034$时, 由于Re$(\lambda) < 0$, 故不可能出现Turing不稳定现象; 当$h>0.059 6$时, 即当$h$取值为$0.095, 0.128, 0.156, 0.198$时分别有$k\in(0.444 1, 1.043 9)$, $k\in(0.366 6, 1.196 4)$, $k\in(0.319 6, 1.299 8)$, $k\in(0.264 1, 1.426 2)$, 使得$\delta_{k} < 0$, 此时Re$(\lambda)>0$, 故当$h$为上述值时, 均会出现Turing不稳定现象.此时, 临界值$h_{T}=0.059 6$.

接下来, 将讨论Turing斑图的存在区域.选择$h$作为研究的分支参数.由分支理论^[15-16]可知当${\rm Im}(\lambda_{k})\neq0$和Re$(\lambda_{k})=0$在$k=0$成立时, Hopf分支发生且Hopf分支曲线为

$ \dfrac{su_{\ast}}{\beta\gamma(\beta+u_{\ast})^{2}}+ \dfrac{h}{(\rho+u_{\ast})^{2}}-\dfrac{\alpha\beta}{(\beta+u_{\ast})^{2}} =1. $

(16)

当Im$(\lambda_{k})=0$和Re$(\lambda_{k})=0$在$k=k_{T}\neq0$成立时, Turing分支发生且分支参数的临界值$h_{T}$满足如下的Turing分支曲线:

$ \Big( \dfrac{su_{\ast}^{2}}{\beta\gamma(\beta+u_{\ast})^{2}}+ \dfrac{hu_{\ast}}{(\rho+u_{\ast})^{2}}-u_{\ast}\Big)d_{22}- \dfrac{\alpha\beta u_{\ast}}{(\beta+u_{\ast})^{2}}d_{11}-\dfrac{\alpha u_{\ast}}{(\beta+u_{\ast})^{2}\gamma}d_{12}+\dfrac{su_{\ast}} {\beta+u_{\ast}}d_{21} \notag\\ -2\sqrt{(d_{11}d_{22}-d_{12}d_{21}){\rm det}J}=0. $

(17)

通过Hopf分支曲线和Turing分支曲线, 我们得到了Hopf分支区域和Turing不稳定区域, 如图 3所示. Turing空间位于Turing分支曲线上方以及Hopf分支曲线下方, 当系统(3)中的参数位于此Turing空间内时会出现Turing失稳, 从而出现Turing斑图.

2 振幅方程

本节中, 主要利用多重尺度分析法推导系统(3)的振幅方程, 获得系统(3)的Turing斑图选择结果, 从而得到Turing斑图的类型.

令$u=u_{\ast}+\widetilde{u}, v=v_{\ast}+\widetilde{v}$, 为了便于表示我们仍用原变量.系统$(3)$在平衡点$(u_{\ast}, v_{\ast})$处的Taylor展开式为:

$ \left\{\!\! \begin{array}{l} \dfrac{\partial u}{\partial t}= d_{11}\Delta u+d_{12}\Delta v+a_{11}u+a_{12}v+ \dfrac{1}{2}P_{20}u^{2}+P_{11}uv+\dfrac{1}{2}P_{02}v^{2}+ \dfrac{1}{6}P_{30}u^{3}\\[2mm] \qquad +\dfrac{1}{2}P_{21}u^{2}v+\dfrac{1}{2}P_{12}uv^{2}+ \dfrac{1}{6}P_{03}v^{3}+o(\rho^{3}), \\[2mm] \dfrac{\partial v}{\partial t}= d_{21}\Delta u+d_{22}\Delta v+a_{21}u+a_{22}v+\dfrac{1}{2}Q_{20} u^{2}+Q_{11}uv+\dfrac{1}{2}Q_{02}v^{2}+\dfrac{1}{6}Q_{30}u^{3}\\[2mm] \qquad +\dfrac{1}{2}Q_{21}u^{2}v+\dfrac{1}{2}Q_{12}uv^{2}+ \dfrac{1}{6}Q_{03}v^{3}+o(\rho^{3}), \end{array} \right. $

(18)

系统$(18)$可表示为

$ \dfrac{\partial U}{\partial t}=LU+N, $

(19)

其中

$ N=\left( \begin{array}{c}\frac{1}{2}P_{20}u^{2}+P_{11}uv+ \frac{1}{2}P_{02}v^{2}+ \frac{1}{6}P_{30}u^{3}+\frac{1}{2}P_{21}u^{2}v+\frac{1}{2} P_{12}uv^{2}+\frac{1}{6}P_{03}v^{3}\\ \frac{1}{2}Q_{20}u^{2}+Q_{11}uv+\frac{1}{2}Q_{02}v^{2}+ \frac{1}{6}Q_{30}u^{3}+\frac{1}{2}Q_{21}u^{2}v+\frac{1} {2}Q_{12}uv^{2}+\frac{1}{6}Q_{03}v^{3} \end{array} \right), $

$ L=L_{T}+(h-h_{T})M=\left( \begin{array}{cc} a^{T}_{11}+d_{11}\Delta & a^{T}_{12}+d_{12}\Delta\\[2mm] a^{T}_{21}+d_{21}\Delta & a^{T}_{22}+d_{22}\Delta \end{array} \right)+ (h-h_{T})M, $

$ M=\left( \begin{array}{cc} m_{11} & m_{12}\\ m_{21} & m_{22} \end{array} \right) =\left.\left( \begin{array}{cc} \dfrac{\partial a_{11}}{\partial h} & \dfrac{\partial a_{12}} {\partial h}\\[2mm] \dfrac{\partial a_{21}}{\partial h} & \dfrac{\partial a_{22}}{\partial h} \end{array} \right)\right|_{h=h_{T}}. $

这里, $h$取为控制参量, $h_{T}$为对应的Turing分支临界值, 临界波数$k_{T}^{\ast2}=\sqrt{\frac{a_{11}^{T}a_{22}^{T}- a_{12}^{T}a_{21}^{T}}{d_{11}d_{22}-d_{12}d_{21}}}$, 且$f=\frac{d_{22}k_{T}^{\ast2}-a_{22}^{T}}{a_{21}^{T}-d_{21}k_{T}^ {\ast2}}$, $g=\frac{d_{11}k_{T}^{\ast2}-a_{11}^{T}}{a_{21}^{T}-d_{21}k_{T}^ {\ast2}}$满足下列方程:

$ L_{T}\left( \begin{array}{l} f\\ 1 \end{array} \right) =0, \quad L^{\ast}_{T}\left( \begin{array}{l} 1\\ g \end{array} \right) =0, $

其中$L_{T}^{\ast}$为$L_{T}$的伴随算子.

在计算中, 我们只分析控制参数在临界值附近的行为, 因此将控制参数$h$展成如下的形式:

$ h-h_{T}=\varepsilon h_{1}+\varepsilon^{2}h_{2}+\varepsilon^{3}h_{3}+o(\varepsilon^{3}), $

(20)

其中$\varepsilon$为一个小参数.同样地, 将$U$和$N$也按$\varepsilon$展开:

$ U=\left( \begin{array}{l} u\\ v \end{array} \right) =\varepsilon\left( \begin{array}{l} u_{1}\\ v_{1} \end{array} \right) +\varepsilon^{2}\left( \begin{array}{l} u_{2}\\ v_{2} \end{array} \right) +\varepsilon^{3}\left( \begin{array}{l} u_{3}\\ v_{3} \end{array} \right) +o(\varepsilon^{3}), $

$ N=\varepsilon^{2}N_{1}+\varepsilon^{3}N_{2}+o(\varepsilon^{3}), $

$ N_{1}=\left(\!\! \begin{array}{c} \frac{1}{2}P_{20}u_{1}^{2}+P_{11}u_{1}v_{1}+\frac{1}{2} P_{02}v_{1}^{2}\\ \frac{1}{2}Q_{20}u_{1}^{2}+Q_{11}u_{1}v_{1}+ \frac{1}{2}Q_{02}v_{1}^{2} \end{array} \!\!\right)\Bigg|_{h=h_{T}}, $

$ N_{2}\!=\!\left(\!\!\!\! \begin{array}{c} P_{20}u_{1}u_{2}\!+\!P_{11}(u_{2}v_{1}\!+\!u_{1}v_{2})\!+\!P_{02}v_{1} v_{2}\!+\! \frac{1}{6}P_{30}u_{1}^{3}\!+\!\frac{1}{2}P_{21}u_{1}^{2}v_{1} + \frac{1}{2}P_{12}u_{1}v_{1}^{2}\!+\!\frac{1}{6}P_{03}v_{1}^{3}\\ Q_{20}u_{1}u_{2}\!+\!Q_{11}(u_{2}v_{1}\!+\!u_{1}v_{2})\!+\!Q_{02}v_{1} v_{2}\!+\! \frac{1}{6}Q_{30}u_{1}^{3}\!+\!\frac{1}{2}Q_{21}u_{1}^{2}v_{1}\!+\! \frac{1}{2}Q_{12}u_{1}v_{1}^{2} \!+\!\frac{1}{6}Q_{03}v_{1}^{3} \end{array} \!\!\!\!\right)\Bigg|_{h=h_{T}}. $

$ \left( \begin{array}{l} u_{1}\\ v_{1} \end{array} \right)= \left( \begin{array}{l} f\\ 1 \end{array} \right) \Bigg(\sum\limits_{j=1}^{3}(W_{j}\exp({\rm i}\overrightarrow{{ k}_{j}}\cdot\overrightarrow{ r}))\Bigg)+c.c., $

$ \left( \begin{array}{l} u_{2}\\ v_{2} \end{array} \right)= \left( \begin{array}{l} U_{0}\\ V_{0} \end{array} \right)+ \sum\limits_{j=1}^{3} \left( \begin{array}{l} U_{j}\\ V_{j} \end{array} \right) \exp({\rm i}\overrightarrow{{ k}_{j}}\cdot\overrightarrow{ r})+ \sum\limits_{j=1}^{3} \left( \begin{array}{l} U_{jj}\\ V_{jj} \end{array} \right) \exp(2{\rm i}\overrightarrow{{ k}_{j}}\cdot\overrightarrow{ r}) \\ + \left( \begin{array}{l} U_{12}\\ V_{12} \end{array} \right) \exp({\rm i}(\overrightarrow{{ k}_{1}}-\overrightarrow{{ k}_{2}})\cdot \overrightarrow{ r})+ \left( \begin{array}{l} U_{23}\\ V_{23} \end{array} \right) \exp({\rm i}(\overrightarrow{{ k}_{2}}-\overrightarrow{{ k}_{3}}) \cdot\overrightarrow{ r})\\ + \left( \begin{array}{l} U_{31}\\ V_{31} \end{array} \right) \exp({\rm i}(\overrightarrow{{ k}_{3}}-\overrightarrow{{ k}_{1}})\cdot\overrightarrow{ r})+c.c., $

其中

$ \left( \begin{array}{l} U_{0}\\ V_{0} \end{array} \right)= \left( \begin{array}{l} u^{0}\\ v^{0} \end{array} \right) (|W_{1}|^{2}+|W_{2}|^{2}+|W_{3}|^{2}), \quad U_{j}=fV_{j}, \\ \left( \begin{array}{l} U_{jj}\\ V_{jj} \end{array} \right)= \left( \begin{array}{l} u^{1}\\ v^{1} \end{array} \right) W_{j}^{2}, \quad \left( \begin{array}{l} U_{jk}\\ V_{jk} \end{array} \right)= \left( \begin{array}{l} u^{\star}\\ v^{\star} \end{array} \right) W_{j}\overline{W}_{k}, \\ u^{0}=\dfrac{2(l_{1}a^{T}_{22}-l_{2}a^{T}_{12})}{a^{T}_{11}a^{T} _{22}-a^{T}_{12}a^{T}_{21}}, \quad v^{0}=\dfrac{2(l_{2}a^{T}_{11}-l_{1}a^{T}_{21})}{a^{T}_{11} a^{T}_{22}-a^{T}_{12}a^{T}_{21}}, \\ u^{1}=\dfrac{l_{1}(a_{22}^{T}-4d_{22}k_{T}^{\ast2})-l_{2}(a_{12}^{T} -4d_{12}k_{T}^{\ast2})}{(a_{11}^{T}-4d_{11}k_{T}^{\ast2})(a_{22}^{T} -4d_{22}k_{T}^{\ast2})-(a_{12}^{T}-4d_{12}k_{T}^{\ast2})(a_{21}^{T}- 4d_{21}k_{T}^{\ast2})}, \\ v^{1}=\dfrac{l_{2}(a_{11}^{T}-4d_{11}k_{T}^{\ast2})-l_{1}(a_{21}^ {T}-4d_{21}k_{T}^{\ast2})}{(a_{11}^{T}-4d_{11}k_{T}^{\ast2}) (a_{22}^{T}-4d_{22}k_{T}^{\ast2})-(a_{12}^{T}-4d_{12}k_{T}^{\ast2}) (a_{21}^{T}-4d_{21}k_{T}^{\ast2})}, \\ u^{\star}=\dfrac{2l_{1}(a_{22}^{T}-3d_{22}k_{T}^{\ast2})-2l_{2} (a_{12}^{T}-3d_{12}k_{T}^{\ast2})}{(a_{11}^{T}-3d_{11}k_{T}^ {\ast2})(a_{22}^{T}-3d_{22}k_{T}^{\ast2})-(a_{12}^{T}-3d_{12} k_{T}^{\ast2})(a_{21}^{T}-3d_{21}k_{T}^{\ast2})}, \\ v^{\star}=\dfrac{2l_{2}(a_{11}^{T}-3d_{11}k_{T}^{\ast2})-2l_{1} (a_{21}^{T}-3d_{21}k_{T}^{\ast2})}{(a_{11}^{T}-3d_{11}k_{T}^{\ast2}) (a_{22}^{T}-3d_{22}k_{T}^{\ast2})-(a_{12}^{T}-3d_{12}k_{T}^{\ast2}) (a_{21}^{T}-3d_{21}k_{T}^{\ast2})}, \\ l_{1}=-\dfrac{f^{2}}{2}\cdot P_{20}-\dfrac{1}{2}P_{02}-f\cdot P_{11}, \quad l_{2}=-\dfrac{f^{2}}{2}\cdot Q_{20}-\dfrac{1}{2}Q_{02}-f\cdot Q_{11}, $

而$c.c.$表示复数共轭项.

利用中心流形理论推导得到如下的振幅方程^[15]

$ \left\{ \begin{array}{l} \tau_{0}\dfrac{\partial A_{1}}{\partial t}=\mu A_{1}+l\overline{A_{2}} \overline{A_{3}}-[g_{1}|A_{1}|^{2}+g_{2}(|A_{2}|^{2}+|A_{3}| ^{2})]A_{1}, \\[2mm] \tau_{0}\dfrac{\partial A_{2}}{\partial t}=\mu A_{2}+l \overline{A_{1}}\overline{A_{3}}-[g_{1}|A_{2}|^{2}+g_{2}(|A_{1}| ^{2}+|A_{3}|^{2})]A_{2}, \\[2mm] \tau_{0}\dfrac{\partial A_{3}}{\partial t}=\mu A_{3}+l\overline{A_{1}}\overline{A_{2}}-[g_{1}|A_{3}|^{2}+g_{2}(|A_{1}|^{2}+|A_{2}|^{2})]A_{3}, \\ \end{array} \right. $

(21)

其中, $\mu$表示到临界值的归一化距离, $\tau_{0}$表示松弛时间.通过计算得到系统$(21)$中系数$\tau_{0}$, $l$, $g_{1}$, $g_{2}$的表达式分别为

$ \begin{align*} \tau_{0}=&\dfrac{f+g}{h_{T}[fm_{11}+m_{12}+g(fm_{21}+m_{22})]}, \mu=\dfrac{h-h_{T}}{h_{T}}, \\ l=&\dfrac{L}{h_{T}[fm_{11}+m_{12}+g(fm_{21}+m_{22})]}, g_{1}=\dfrac{G_{1}}{h_{T}[fm_{11}+m_{12}+g(fm_{21}+m_{22})]}, \\ g_{2}=&\dfrac{G_{2}}{h_{T}[fm_{11}+m_{12}+g(fm_{21}+m_{22})]}, \end{align*} $

其中

$ L=-2l_{1}-2gl_{2}|_{h=h_{T}}, $

$ G_{1}=-\Big\{(f\cdot P_{20}+P_{11})(u^{0}+u^{1})+(f\cdot P_{11}+P_{02})(v^{0}+v^{1})+\dfrac{f^{3}}{2}\cdot P_{30}+ \dfrac{1}{2}P_{03}\\ +\dfrac{3f^{2}}{2}\cdot P_{21}+\dfrac{3f}{2}\cdot P_{12}+g\Big[(f\cdot Q_{20}+Q_{11})(u^{0}+u^{1})+(f\cdot Q_{11}+Q_{02})(v^{0}+v^{1})\\ +\dfrac{f^{3}}{2}\cdot Q_{30}+\dfrac{1}{2}Q_{03}+\dfrac{3f^{2}}{2}\cdot Q_{21}+\dfrac{3f}{2}\cdot Q_{12}\Big]\Big\}\Big|_{h=h_{T}}, $

$ G_{2}=-\{(f\cdot P_{20}+P_{11})(u^{0}+u^{\star})+(f\cdot P_{11}+ P_{02})(v^{0}+v^{\star})+f^{3}\cdot P_{30}+P_{03}\\ +3f^{2}\cdot P_{21}+3f\cdot P_{12}+g[(f\cdot Q_{20}+Q_{11})(u^{0}+ u^{\star})+(f\cdot Q_{11}+Q_{02})(v^{0}+v^{\star})\\ +f^{3}\cdot Q_{30}+Q_{03}+3f^{2}\cdot Q_{21}+3f\cdot Q_{12}]\}| _{h=h_{T}}. $

方程组$(21)$的每个振幅均可分解为模$\rho_{j}=|A_{j}|$与一个相应的相角$\psi_{j}$的乘积, 将$A_{j}=\rho_{j}{\rm e}^{{\rm i}\phi_{j}}$带入方程组$(21)$并且分离实部和虚部, 可得到如下由4个方程构成的方程组:

$ \left\{ \begin{array}{l} \tau_{0}\dfrac{\partial \phi}{\partial t}=-l\dfrac{\rho_{1}^{2} \rho_{2}^{2}+\rho_{1}^{2}\rho_{3}^{2}+\rho_{2}^{2}\rho_{3}^{2}} {\rho_{1}\rho_{2}\rho_{3}}\sin\phi, \\ \tau_{0}\dfrac{\partial \rho_{1}}{\partial t}=\mu\rho_{1} +l\rho_{2}\rho_{3}\cos\phi-g_{1}\rho_{1}^{3}-g_{2}(\rho_{2}^{2} +\rho_{3}^{2})\rho_{1}, \\ \tau_{0}\dfrac{\partial \rho_{2}}{\partial t}=\mu\rho_{2}+ l\rho_{1}\rho_{3}\cos\phi-g_{1}\rho_{2}^{3}-g_{2}(\rho_{1}^{2} +\rho_{3}^{2})\rho_{2}, \\ \tau_{0}\dfrac{\partial \rho_{3}}{\partial t}=\mu\rho_{3}+l\rho_{1}\rho_{2}\cos\phi-g_{1}\rho_{3}^{3}-g_{2}(\rho_{1}^{2}+\rho_{2}^{2})\rho_{3}, \\ \end{array} \right. $

(22)

其中$\phi=\phi_{1}+\phi_{2}+\phi_{3}$.系统$(22)$有$4$种类型的解^[15].

$(1)$定态解$(O)$: $\rho_{1}=\rho_{2}=\rho_{3}=0$.当$\mu < \mu_{2}=0$时, 定态解稳定; 当$\mu>\mu_{2}$时不稳定.

$(2)$条状斑图$(S)$: $\rho_{1}=\sqrt{\frac{\mu}{g_{1}}}, $ $\rho_{2}=\rho_{3}=0$.当$\mu>\mu_{3}=\frac{l^{2}g_{1}}{(g_{2}-g_{1}) ^{2}}$时, 该解稳定, 否则不稳定.

$(3)$六边形斑图$(H_{0}, H_{\pi})$: $\rho_{1}=\rho_{2}=\rho_{3}= \frac{|l|\pm\sqrt{l^{2}+4(g_{1}+2g_{2})\mu}}{2(g_{1}+2g_{2})}.$当$\mu>\mu_{1}=\frac{-l^{2}}{4(g_{1}+2g_{2})}$时, 该解存在; 当$\mu < \mu_{4}=\frac{2g_{1}+g_{2}}{(g_{2}-g_{1})^{2}}l^{2}$时, 解$\rho^{+}=\frac{|l|+\sqrt{l^{2}+4(g_{1}+2g_{2})\mu}}{2(g_{1}+ 2g_{2})}$稳定, 而解$\rho^{-}=\frac{|l|-\sqrt{l^{2}+4(g_{1}+2g_{2}) \mu}}{2(g_{1}+2g_{2})}$一直不稳定.

$(4)$混合斑图: $\rho_{1}=\frac{|l|}{g_{2}-g_{1}}$, $\rho_{2}=\rho_{3}=\sqrt{\frac{\mu-g_{1}\rho_{1}^{2}}{g_{1}+g_{2}}}$.当$g_{2}>g_{1}$, $\mu>g_{1}\rho_{1}^{2}$时, 该解存在.

通过上述讨论, 我们知道$\mu_{1} < \mu_{2} < \mu_{3} < \mu_{4}$, 并且得到类似于文献[15]中的定理.

定理2.1 若系统(3)的参数值在Turing空间中发生变化, 则

(ⅰ)当$\mu$位于区间$(\mu_{2}, \mu_{3})$时, 系统(3)仅有六边形斑图(即点状斑图);

(ⅱ)当$\mu$位于区间$(\mu_{3}, \mu_{4})$时, 系统(3)存在双稳态(即点状斑图与条状斑图共存);

(ⅲ)当$\mu>\mu_{4}$, 系统(3)仅有稳定的条状斑图.

3 数值模拟

本节中, 我们通过Matlab软件对关于系统(3)的空间斑图进行了数值模拟.所有的数值模拟都在离散格子$200\times200$内进行, 取定空间间隔$\Delta h=1$, 时间间隔$\Delta t=0.01$.初始条件是种群的初始随机分布, 所有的数值模拟都采用了齐次Neumann边界条件.由于食饵和捕食者的空间分布是类似的, 所以在接下来的数值模拟中我们仅研究食饵的空间斑图的形成.

首先, 设定参数值$\alpha=1.5$, $\beta=0.8$, $\gamma=0.8$, $\rho=0.4$, $s=1$, 变化$h$的值, 则可得到参数值$h_{T}=0.059 6$, $\tau_{0}=28.126 8$, $l=3.849 6$, $g_{1}=85.266 6$, $g_{2}=185.76$, $\mu_{1}=-0.008 1$, $\mu_{2}=0$, $\mu_{3}=0.125 1$, $\mu_{4}=0.522 8$.这里扩散系数取值为$d_{11}=1$, $d_{12}=1$, $d_{21}=14$, $d_{22}=15$. 图 4展示的是系统(3)的Turing分支图.从图 4中我们可以看出, 系统(3)存在双稳态区域$(\mu_{1}, \mu_{2})$, 即当控制参量$\mu\in(\mu_{1}, \mu_{2})$时, 六边形斑图$H_{0}$(点状斑图)和定态解均是稳定的.当$\mu\in(\mu_{2}, \mu_{3})$时, 条状斑图是不稳定的, 六边形斑图$H_{0}$(点状斑图)是稳定的.当$\mu\in(\mu_{3}, \mu_{4})$时, 系统(3)存在另一个双稳态, 即点状和条状斑图共存.当$\mu\in(\mu_{4}, +\infty)$时, 仅有稳定的条状斑图.

图 5展示了当$h=0.110$时, 食饵$u$的空间斑图演化过程. 4个子图的迭代步数分别为$0$, 20 000, 60 000, 800 000.从图 5(a)中的右边颜色条(colorbar)数值基本不变可看出, 初值选取为平衡解$E_{\ast}=(0.366 0, 0.571 8)$加上一个随机扰动, 之后出现了类似点状斑图(图 5(b)), 然后条状斑图占优(图 5(c)), 最终条状斑图(类似迷宫斑图)占据整个区域(图 5(d)), 并且系统的动力学行为不再发生变化.参数的选择满足$\mu\approx0.845 6>\mu_{4}$, 通过上述定理2.1可知系统(3)仅存在条状斑图, 即数值模拟结果验证了上述理论结果的正确性.

图 6和图 7分别展示了当$h=0.082$, $h=0.071$时, 食饵$u$的空间斑图演化过程. 图 6中四个子图的迭代步数分别为$0$, 80 000, 140 000, 800 000, 参数的选择满足$\mu_{3} < \mu\approx0.375 8 < \mu_{4}$. 图 7中4个子图的迭代步数分别为$0$, 130 000, 170 000, 800 000, 参数的选择满足$\mu_{3} < \mu\approx0.191 3 < \mu_{4}$. 图 6和图 7均出现点状斑图和条状斑图共存的现象, 这与理论结果相一致.不同的是图 6(d)中条状斑图明显占优, 而图 7(d)中点状斑图明显占优.

图 8展示了当$h=0.065$时, 食饵$u$随时间演化的空间斑图. 图 8中4个子图的迭代步数分别为$0$, 280 000, 340 000, 800 000. 图 8(a)仍选取平衡解$E_{\ast}=(0.398 8, 0.623 2)$加上一个随机扰动, 之后出现了点状和条状共存的现象(图 8(b)和图 8(c)), 最终点状斑图占据整个区域(图 8(d)), 此处$\mu_{2} < \mu\approx0.090 6 < \mu_{3}$.同样地, 数值模拟结果与理论结果相一致.

4 结论

基于收获策略对生态系统的动力学行为具有重要的影响, 本文从理论和数值两方面研究了一类具有非线性收获效应的捕食者-食饵模型的空间斑图生成与选择问题, 所得结果表明:

(1) 理论结果表明自扩散不会导致系统(3)产生Turing斑图现象.即系统(3)的Turing斑图的产生是由交叉扩散项引起的.

(2) 理论与数值结果表明非线性收获效应影响了Turing斑图的生成.即在其他参数固定的情况下, 收获能力参数$h$必须大于某个临界值$h_{T}$时, 系统才会出现Turing不稳定现象.

(3) 理论和数值结果表明非线性收获效应影响了Turing斑图的选择.即随着收获能力参数$h$的变化, 系统展现了点状, 条状以及二者共存的斑图结构.