高 鹤,李秀玲

(吉林财经大学 应用数学学院,长春 130117)

0 引 言

在自然界中,两个生物种群之间的相互关系有捕食者与被捕食者、寄生关系、竞争关系、互惠共存关系等.如果两个物种之间的相互作用对这两个物种都有利,则两个物种之间的共生关系称为互惠共生的.互惠共生可通过两种方式帮助猎物: 一是提高猎物的生长速度;二是阻止捕食者捕食猎物.在第二种情形中,关于时滞对种群动力系统影响的研究目前已有很多结果,例如: Ojha等[1]给出了一个考虑休眠时滞的浮游生物动力学模型;Djilali等[2]分析了一个捕食-食饵扩散模型的复杂动力学行为;Lavanya等[3]提出了具有食肉动物和捕食者孕育时滞的Holling-Ⅱ型捕食者-食饵模型,并研究了所有平衡点的存在性,以及时滞模型的局部稳定性和Hopf分支.目前,关于捕食者食饵共生系统的研究已有一些结果[4-12],但对于具有时滞的捕食者食饵共生模型系统的研究文献报道较少.本文在文献[13-14]的基础上,讨论棉蚜、蚂蚁、瓢虫之间的种间关系.在棉田中,棉蚜是主要害虫之一,它们通过吸食汁液、传染病菌使棉花滋生霉菌,植株矮小,叶片变小,从而导致棉花大量减产,造成重大经济损失.随着农药的大量应用,棉蚜的抗药性逐渐增强,此时需通过引入天敌进行防治.瓢虫以棉蚜为食,蚂蚁会驱赶瓢虫,从而获得更多棉蚜分泌的营养物质,因此棉蚜、蚂蚁、瓢虫三者相互制约,构成了捕食-食饵共生系统.通过控制三者的数量关系可解决棉田虫害问题,从而更好地保护生态环境平衡.以食饵的发育时间和共生者的发育时间作为时滞,考虑如下具双时滞的捕食-食饵共生模型:

(1)

其中x(t),y(t),z(t)分别表示食饵的种群密度、共生者的种群密度和捕食者的种群密度,r表示食饵的内禀增长率,K表示环境能维持的食饵最大数量,β表示捕食者对食饵的转化率,k1表示食饵对共生者的共生效应,d1表示共生者的自然死亡率,v表示共生者因拥挤而导致的死亡率,k2表示共生者对捕食者的驱赶效应,d2表示捕食者的自然死亡率,δ表示捕食者因种群内部竞争导致的死亡率,τ1,τ2分别表示食饵的发育时间和共生者的发育时间.所有参数都是正常数,考虑到其实际生物学意义,发育时间只对其中一些变量有影响,死亡率不受时滞影响.本文考虑双时滞对模型(1)的影响,首先,给出系统平衡点的稳定性及Hopf分支的存在性;其次,利用规范型理论和中心流形定理给出Hopf分支方向和分支周期解的稳定性;最后,利用Mathematica软件进行数值模拟,验证所得结论的正确性.通过对模型结果的分析和讨论,给出模型物种之间的基本关系,从而对棉田的生物防治提供一定的理论基础.

1 系统内部平衡点的稳定性和Hopf分支的存在性

令系统(1)各等式右边为零,可得系统(1)的内部正平衡点E=(x*,y*,z*),其中

系统(1)在平衡点E=(x*,y*,z*)处的Jacobi矩阵为

其中

a22=k1x*-d1-2vy*,a33=βx*-k2y*-d2-2δz*,

系统(1)在平衡点E=(x*,y*,z*)处对应的特征方程为

λ3+P02λ2+P01λ+P00+(P12λ2+P11λ+P10)e-λτ1+P20e-λ(τ1+τ2)=0,

(2)

其中

P00=-a11a22a33,P01=a11a22+a11a33+a22a33,

P02=-a11-a22-a33,P10=a22a13b31-a22a33b11,

P11=b11a22+b11a33-a13b31,P12=-b11,P20=-a13b21c32.

根据时滞τ1和τ2取值的不同,下面分4种情形讨论.

情形1)τ1=τ2=0.

此时,方程(2)可转化为如下形式:

λ3+J12λ2+J11λ+J10=0,

(3)

其中

J10=P00+P10+P20,J11=P01+P11,J12=P02+P12.

假设:

(H1)J10>0,J11>0,J12>0,J12J11>J10.

根据Routh-Hurwitz稳定性判据,可得以下结论:

定理1当且仅当条件(H1)成立时,方程(3)的所有根都具有负实部,即系统(1)在正平衡点E=(x*,y*,z*)处是渐近稳定的.

下面考虑时滞对系统(1)正平衡点E=(x*,y*,z*)稳定性的影响.

情形2)τ1>0,τ2=0.

此时,方程(2)可转化为

λ3+P02λ2+P01λ+P00+(P12λ2+P11λ+P10+P20)e-λτ1=0.

(4)

令λ=iω1(ω1>0)为方程(4)的根,分离实部和虚部可得

其中

于是可得

(5)

其中

θ10=-P00P10-P00P20,θ11=-P00P11+P01P11+P01P20,

θ12=-P01P11+P02P10+P02P20+P12P00,

θ13=P02P11-P10-P20-P01P12,θ14=P11-P02P12,θ15=P12,

将式(5)中两式两端取平方后再相加可得

(6)

其中

假设:

(H2)φ0<0或φ0≥0且Δ>0,s1>0,h(s1)≤0.

引理1当条件(H2)成立时,方程(5)存在正根.

s3+φ2s2+φ1s+φ0=0.

(7)

再令

h(s)=s3+φ2s2+φ1s+φ0,

(8)

由式(8)可知

因此当s1>0,h(s1)≤0时,方程(7)有正根.

这里

令λ(τ1)=α(τ1)+iω(τ1)是方程(4)在时滞τ1=τ10处满足α(τ1)=0,ω(τ1)=ω10的根.

假设:

(H3) 3s2+2φ2s+φ1>0.

证明: 将方程(4)两边对τ1求导,可得

由方程(4)可得

(P12λ2+P11λ+P10+P20)e-λτ1=-λ3-P02λ2-P01λ-P00,

于是有

(9)

将λ=iω10代入式(9)可得

综上,可得以下结论:

定理2对于系统(1),当τ1>0,τ2=0时,如果条件(H2),(H3)成立,则当τ1∈[0,τ10)时,系统(1)的平衡点E=(x*,y*,z*)是渐近稳定的;当τ1>τ10时,系统(1)的平衡点E=(x*,y*,z*)是不稳定的;当τ1=τ10时,系统(1)在正平衡点E=(x*,y*,z*)处产生Hopf分支.

情形3)τ1=0,τ2>0.

此时,方程(2)可转化为如下形式:

λ3+F02λ2+F01λ+F00+F10e-λτ2=0,

(10)

其中

F00=P00+P10,F01=P01+P11,F02=P02+P12,F10=P20.

令λ=iω2(ω2>0)为方程(10)的根,分离实部和虚部可得

(11)

将式(11)中两式两端取平方后再相加可得

(12)

假设:

(H4)φ5<0或φ5≥0且Δ2>0,r1>0,h(r1)≤0.

引理3若条件(H4)成立,则方程(12)存在正根.

证明: 令

则方程(12)变为

r3+φ3r2+φ4r+φ5=0.

(13)

再令

g(r)=r3+φ3r2+φ4r+φ5.

由g(r)可知,

因此当r1>0,g(r1)≤0时,方程(10)有正根.

令λ(τ2)=α(τ2)+iω(τ2)是方程(12)在时滞τ2=τ20处满足α(τ2)=0,ω(τ2)=ω20的根.

假设:

(H5) 3r2+2φ3r+φ4>0.

证明: 将式(10)两边对τ2求导,得

(14)

由于F10e-λτ2=-λ3-F02λ2-F01λ-F00,将λ=iω20代入式(14)可得

综上,可得以下结论:

定理3对于系统(1),当τ1=0,τ2>0时,若条件(H4),(H5)成立,则当τ2∈[0,τ20)时,系统(1)的平衡点E=(x*,y*,z*)是渐近稳定的;当τ2>τ20时,系统(1)的平衡点E=(x*,y*,z*)是不稳定的;当τ2=τ20时,系统(1)在平衡点E=(x*,y*,z*)处产生Hopf分支.

情形4)τ1>0,τ2∈(0,τ20),即把τ2限制在其稳定区间内,将τ1作为参数.

系统(1)在平衡点E=(x*,y*,z*)处对应的特征方程为

λ3+P02λ2+P01λ+P00+(P12λ2+P11λ+P10)e-λτ1+P20e-λ(τ1+τ2)=0.

(15)

令λ=iω3(ω3>0)为方程(15)的根,则

其中

进一步可得

(16)

其中

e2=2P12P20,e1=2P11P20,e0=-2P10P20.

假设:

引理5若条件(H6)成立,则方程(16)至少存在一个正根.

证明: 令

令λ(τ1)=α(τ1)+iω(τ1)是方程(15)在时滞τ1=τ30处满足α(τ1)=0,ω(τ1)=ω30的根.

假设:

(H7)QRPR+QIPI>0.

证明: 将式(15)对τ1求导,得

(17)

将λ=iω30代入式(17)可得

其中

综上,可得以下结论:

定理4对于系统(1),当τ1>0,τ2∈(0,τ20)时,如果条件(H6),(H7)成立,则当τ1∈[0,τ30)时,系统(1)的平衡点E=(x*,y*,z*)是渐近稳定的;当τ1>τ30时,系统(1)的平衡点E=(x*,y*,z*)是不稳定的;当τ1=τ30时,系统(1)在平衡点E=(x*,y*,z*)处产生Hopf分支.

2 Hopf分支的方向和稳定性

(18)

其中

ht=(h1,h2,h3)T∈C([-1,0],3),

由Riesz表示定理可知,存在3×3的矩阵函数η(θ,u)和变量θ∈[-1,0],满足

事实上,可以选择

对于φ∈C([-1,0],3),定义

则式(1)可转化为

(19)

定义A的伴随算子A*为

直接计算可得如下结论:

根据文献[15],计算相关系数如下:

其中

这里E1,E2满足

于是,经计算可得

因此,系数μ2,β2,τ3决定了Hopf分支的性质.其中μ2决定Hopf分支的方向,当μ2>0,τ1>τ30时,系统在平衡点附近发生超临界Hopf分支;当μ2<0,τ1<τ30时,系统在平衡点附近发生亚临界Hopf分支.且在上述相应时滞范围内存在周期解,当T2>0(或T2<0)时,周期解的周期增加(或减少);当β2<0(β2>0)时,周期解是渐近稳定的(或不稳定的).

3 数值模拟

对于系统(1),取如下参数值:r=0.1,K=40,k1=0.03,k2=0.02,d1=0.02,d2=0.04,β=0.02,δ=0.006,v=0.05,则系统(1)可变为

图1 当τ1=τ2=0时,系统(1)的时间序列曲线Fig.1 Time series curves of system (1) when τ1=τ2=0

经计算可得该系统的正平衡点为E=(7.087,3.514,4.114 3).显然有J10=0.000 982,J11=0.020 3,J12=0.235,J12J11=0.004 8>J10.根据Routh-Hurwitz稳定性判据可知,当τ1=τ2=0时,该系统的正平衡点E是渐近稳定的,如图1所示.

当τ1>0,τ2=0时,经计算可得τ10=10.291 4.当τ1<τ10时,该系统的平衡点E是渐近稳定的,如图2所示.当τ1>τ10时,该系统的平衡点E是不稳定的,此时在平衡点E处出现Hopf分支,如图3所示.

当τ1>0,τ2>0时,通过计算可知条件(H6)和(H7)均满足.选取τ1=6<τ30=8.574 7,τ2=5∈(0,τ20)=(0,21.501 7),则系统的平衡点E是渐近稳定的,如图4所示.选取τ1=8.63>τ30=8.574 7,τ2=5,则系统的平衡点E是不稳定的,此时在平衡点E处发生Hopf分支,如图5所示.

图2 当τ1=5,τ2=0时,系统(1)的时间序列曲线(A),(B),(C)及相图(D)Fig.2 Time series curves (A),(B),(C) and phase diagram (D) of system (1) when τ1=5,τ2=0

图3 当τ1=10.292,τ2=0时,系统(1)的时间序列曲线(A),(B),(C)及相图(D)Fig.3 Time series curves (A),(B),(C) and phase diagram (D) of system (1) when τ1=10.292,τ2=0

图4 当τ1=6,τ2=5时,系统(1)的时间序列曲线(A),(B),(C)及相图(D)Fig.4 Time series curves (A),(B),(C) and phase diagram (D) of system (1) when τ1=6,τ2=5

图5 当τ1=8.63,τ2=5时,系统(1)的时间序列曲线(A),(B),(C)及相图(D)Fig.3 Time series curves (A),(B),(C) and phase diagram (D) of system (1) when τ1=8.63,τ2=5

综上所述,本文对一类捕食-食饵共生模型引入了食饵和共生者的发育时间两个时滞.首先,得到了模型平衡点的稳定性及Hopf分支存在性的充分条件,结果表明,在一定的条件下,系统的平衡点是渐近稳定的,且当时滞超过一定阈值后,系统在平衡点处发生Hopf分支,并得到了一系列的周期解;其次,利用规范型理论和中心流形定理,给出了上述模型的Hopf分支方向和分支周期解的稳定性;最后,利用Mathematica软件对该模型进行数值模拟,验证了所得结论的正确性.