具有避难所和双时滞的捕食系统的Hopf分支

于智博, 吕堂红, 周林华

(长春理工大学 数学与统计学院， 吉林 长春 130022)

研究Lotka-Volterra模型的目的是描述和分析捕食者和猎物之间的相互作用，反映真实的自然现象，揭示种群动力学的本质。生态学中的概念，如物种迁移、功能响应函数、种内竞争、Allee效应、食饵避难所、恐惧效应和时间延迟等因素被添加到Lotka-Volterra方程中，以获得对种群数量变化和动力学行为更准确的描述。在现实世界中，食饵寻找避难所是一种可以提高存活率的策略。避难所的存在对捕食者种群和猎物种群的共存产生明显的影响，被视为影响种群动态的关键性因素。研究表明，随着躲入避难所的猎物数量增加，捕食者种群的数量会先增加后减小，直到灭绝[1-2]。早期的研究工作主要使用线性猎物避难所，即有一定比例的猎物被避难所保护。此时系统平衡点的稳定性与存在性取决于庇护强度[3-5]。近年来，有学者考虑更为实际的情况，即躲入避难所的猎物数量受到猎物种群数量和捕食者种群数量的影响，因此提出非线性避难所，且以非线性避难所作为控制项研究系统的动力学性质，并得出丰富的结论[6-7]。

食饵种群数量的减少，除了受到捕食者种群的直接捕杀外，还受到恐惧效应的间接影响。2016年，Wang Xiao-ying等[8]研究了如下形式的具有恐惧效应的捕食系统:

(1)

将追捕猎物的时间、食饵的妊娠期、捕食者的进食时间纳入经典的Lotka-Volterra模型更符合自然界的实际情况。在生物模型中引入时滞比常微分方程具有更复杂的动力学特性，因为它可以使平衡点的稳定性发生改变，并产生极限环[10-11]。受到上述研究内容的启发，本文在模型(1)的基础上，提出了一个包括非线性避难所与双时滞的捕食-食饵模型，并研究其动力学性态。改进后的系统模型如下：

(2)

1 正平衡点的存在性

系统(2)正平衡点E*记做

如果系统(2)存在平衡点，则应满足

(3)

由式(3)中的第二个方程可得

(4)

将式(4)代入式(3)中的第一个方程，有

η1(x*)4+η2(x*)3+η3(x*)2+η4x*+η5=0,

(5)

其中

η1=-c3p3δ3d2-d2c3p3δ2k，

η3=3c2p2dδ2r-3c2p2δ2dd1-3cpd2d2δ-2d1c2p2kδd-d2d2cpk-2c3p3δd-c2p3dk，

η4=3rcpd2δ-3cpd2d1δ-d2d3-d1d2cpk-cp2d2，η5=(r-d1)d3。

根据Descarte’s符号准则，当下列条件成立时，至少存在一个正平衡点E*：

(H1)η1<0,η2<0,η3<0,η4>0,η5>0或η1<0,η2<0,η3<0,η4<0,η5>0。

2 局部稳定性与Hopf分支

系统(2)在E*处的Jacobi矩阵为

(6)

则有系统(2)的特征方程：

λ2+A1λ+A2e-λτ1+(A3λ+A4)e-λτ2+A5=0,

(7)

其中A1=-(a11+a22)，A2=-a12a21′，A3=-a22′，A4=a11a22′，A5=a11a22。

情形一τ1=τ2=0。

特征方程(7)变为

λ2+(A1+A3)λ+A2+A4+A5=0。

(8)

根据Routh-Hurwitz准则，若

(H2)A1+A3>0，A2+A4+A5>0,

则系统(2)的正平衡点E*在没有时滞的影响下是稳定的。

情形二τ1=0,τ2>0。

特征方程(7)变为

λ2+A1λ+A2+A5+(A3λ+A4)e-λτ2=0。

(9)

令λ=iω1(ω1>0)是该方程的根，代入式(9),可得

(10)

简化可得

(11)

(12)

于是有如下结论。

定理1 对于系统(2)，当τ1=0,τ2>0与(H1)、(H2)成立，存在一个τ10，使得当τ2∈(0,τ10)时，E*是稳定的；当τ2>τ10时，E*不稳定；当τ2=τ10时，系统(2)发生Hopf分支。

情形三τ1>0,τ2=0。

证明过程类似于情形二。有

(13)

(14)

对于系统(2)，当τ1>0、τ2=0与(H1)成立，当τ1未超过τ20时，平衡点E*是稳定的；当τ1>τ20时，平衡点E*不稳定；当τ1=τ20时，系统(2)发生Hopf分支。

情形四τ1=τ2=τ>0。

特征方程(7)简化为

λ2+A1λ+A5+(A3λ+A2+A4)e-λτ=0。

(15)

令λ=iω3(ω3>0)是该方程的根，代入式(15)，可得

(16)

简化可得

(17)

(18)

剩余证明过程与情形三类似，则对于系统(2)，当τ1=τ2=τ与(H1)成立，当τ未超过临界值τ30时，平衡点E*是稳定的；当τ>τ30时，平衡点E*不稳定；当τ=τ30时，系统(2)发生Hopf分支。

情形五τ1>0,τ2>0。

考虑式(2)中τ1在稳定的区间，τ2作为参数。设λ=iω4为式(7)的根，代入式(7)有

化简得

(19)

若(H6)式(19)至少具有一个非零解成立，则将式(19)的非零解表示为ω40。显然

k=0,1,2,…。

对式(7)关于τ2求导，令τ2=τ40，经计算有

(20)

假设(H7)R≠0，则

定理2 对于系统(2)，固定τ1∈(0,τ20)，若(H1)、(H6)和(H7)成立，则当τ2∈(0,τ40)时，平衡点E*是渐近稳定的；若τ2>τ40，平衡点E*不稳定；当τ2=τ40，系统(2)发生Hopf分支。

3 Hopf分支方向及其稳定性

本节在τ1=τ2=τ=τ30条件下，运用Hassard[12]理论，得到确定E*(x*,y*)附近Hopf分支性质的显示公式。

令u(t)=(u1(t),u2(t))T∈R2，其中u1(t)=x(τt)，u2(t)=y(τt)，τ=τ30+μ，μ∈R,则当μ=0时是系统(2)的Hopf分支值，在C=C([-1,0],R2)上记做

(21)

其中Lμ:C→R2,F:R×C→R2为

(22)

F(μ,φ)=(τ30+μ)(F1(μ,φ),F2(μ,φ))T,

(23)

其中

φ=(φ1,φ2)∈C([-1,0],R2)，

c11=cp-2cpδy*，c12=-d，c13=-4cpδx*，

由Riesz表示定理，可得

η(θ,μ):[-1,0]→R2,

使得

(24)

这里

其中δ(θ)是Dirac-delta函数。

对于φ∈C1([-1,0],R2)，定义

于是，系统(21)可改写为

(25)

这里u=(u1,u2)，ut(θ)=u(t+θ),θ∈[-1,0]。

对于Ψ∈C1([-1,0],(R2)*)，定义A=A(0)的伴随算子A*：

对于φ∈C1([-1,0],R2)，Ψ∈C1([-1,0],(R2)*)，为了使算子A和伴随算子A*的特征向量规范化，可以定义以下的双线性内积:

设q(θ)和q*(s)是对应于特征根iω3τ30与-iω3τ30的特征向量。于是

A(0)q(θ)=iω3τ30q(θ),A*(0)q*(s)=-iω3τ30q*(s),

计算得

这里

令Xt为μ=0时方程(25)的解。定义

z(t)=〈q*,Xt〉,

(26)

在中心流型C0上，有

(27)

(28)

在由式(27)和(28)得

Xt(θ)=W(t,θ)+2Re{z(t)q(θ)},

综合式(23)得

比较系数可得

为了确定g21,下面计算W20(θ)和W11(θ)。由式(25)和(26)得到

即

(29)

其中

(30)

结合式(26)和(27)得到

其中

(31)

通过比较系数得到

(2iω0-A)W20(θ)=H20(θ), -AW11(θ)=H11(θ)。

当θ∈[-1,0)时

(32)

比较式(31)和(32)有

其中W20(θ)、W11(θ)的结果如下：

可得

(33)

(34)

这里，C1(0)由式(33)给出，易得出μ2、β2、T2的值。因此有了定理3。

定理3 当τ=τ30时，由式(34)可得以下结论：

(1)μ2确定Hopf分支的方向，如果μ2>0，则Hopf分支是超临界的，反之是次临界的；

(2)β2确定分支周期解的稳定性，如果β2<0，则周期解是稳定的，反之是不稳定的；

(3)T2确定分支周期解的周期，如果T2>0，则周期解的周期是增加的，反之是减少的。

4 数值模拟

为验证情形三、情形五的结论(情形二与情形三类似)，考虑参数如下：

(35)

根据Descarte’s符号准则，系统(35)有唯一平衡点E*(x*,y*)=E*(17.902 3,3.919 1)。当τ1=τ2=0时，两种群密度的数量变化趋于稳定。

当τ1>0,τ2=0时，τ20=4.67。当τ1=4.4<τ20=4.67时，种群密度保持稳定，见图1；当τ1=6>τ20=4.67时，系统出现周期解，见图2。

当τ1>0,τ2>0时，取τ1=3，τ*=3.849 6。当τ2=2.0<τ*=3.849 6时，平衡点E*是稳定的，两种群数量最终趋向于平衡点E*,如图3所示；当τ2=5>τ*=3.849 6时，平衡点E*附近出现周期解，如图4所示。

(a) x的变化 (b) y的变化 (c) 系统的变化图1 当τ1=4.4<τ20=4.67,τ2=0时，系统(35)的动力学行为

(a) x的变化 (b) y的变化 (c) 系统的变化图2 当τ1=6>τ20=4.67,τ2=0时，系统(35)的动力学行为

(a) x的变化 (b) y的变化 (c) 系统的变化图3 当τ1=3,τ2=2.0<τ*=3.849 6时，系统(35)的动力学行为

(a) x的变化 (b) y的变化 (c) 系统的变化图4 当τ1=3,τ2=5>τ*=3.849 6时，系统(35)的动力学行为

为验证情形四的结论，选择参数：r=0.4,d1=0.05,d2=0.01,p=0.1,δ=0.04,c=0.7,d=0.07,k=0.03。此时有η1<0,η2<0,η3<0,η4<0,η5>0。根据Descarte’s符号准则，至少存在一个正平衡点E*。

当τ1=τ2=τ时，ω3=0.284 5，τ30=4.15，C1(0)=-1.160 5-2.147 7i，μ2=-71.24，β2=-2.53，T2=-7.54。

当τ=4.0<τ30=4.15时，平衡点E*是稳定的,见图5；当τ=5.1>τ30=4.15时，平衡点E*失去稳定性，见图6。

(a) x的变化 (b) y的变化 (c) 系统的变化图5 当τ1=τ2=τ=4.0<τ30=4.15时，系统(35)的动力学行为

(a) x的变化 (b) y的变化 (c) 系统的变化图6 当τ1=τ2=τ=5.1>τ30=4.15时，系统(35)的动力学行为

表1是在情形四的条件下避难所δ的变化对Hopf分支临界值的影响。在情形四的条件下，只改变δ的值，而其余参数不变。研究发现通过增大δ的值，时滞也会逐渐增大。意味着庇护所的存在，对种群数量的稳定起到积极作用。

表1 在情形四的条件下避难所δ的变化对Hopf分支临界值的影响

5 结论

本文研究一类具有避难所和双时滞的捕食系统Hopf分支问题，分析了时滞在不同情况下系统正平衡点的稳定性状况。理论分析表明，时滞的存在是决定系统发生Hopf分支的充分条件。其次选择恰当的参数，验证了理论分析的正确性。从生物学的角度来看，在种群的初始密度固定时，当反馈时滞很小(小于临界值时)，系统将维持稳定的状态，两种群的数量最终趋向一个固定的数值；当反馈时滞(大于临界值时)，捕食者和猎物种群之间的动态平衡是呈现出一种相互限制和周期性循环的形式。此外，数值分析表明，避难所系数的增大，能够改变时滞的临界值，从而增强系统的稳定性。然而，本文并没有考虑恐惧效应对系统稳定性的影响，这一方面的研究可以在未来的工作中进行补充。