分数阶双时滞传染病系统的捕食者_食饵模型动力学
时间：2022-12-10 12:00:04 来源：雅意学习网 本文已影响 人 
李 健, 肖 敏, 周 帅
(南京邮电大学自动化学院, 南京 210023)
近年来, 分数阶微积分发展迅速并被广泛应用于多个领域[1-2]．由于分数阶微积分的阶数可取非整数, 故基于分数阶微积分的分数阶模型能更准确地描述实时系统的性质特征和行为．运用分数阶微分方程解决现实工程问题备受关注[3-4]．例如, Rihan等[5]提出一个时滞分数阶捕食者-食饵模型, 并讨论了模型的传统响应和解的存在性问题．传染病一直是影响生态系统的重大问题之一, 传染病的爆发会导致生态系统紊乱,捕食者与食饵将失去平衡状态, 故研究传染病的动力学行为显得尤为重要[6-7]．目前, 分数阶传染病模型已成为研究热点．Singh等[8]将生态传染病学模型引入分数阶系统,研究了传染病系统平衡点的稳定性和解的性质．众所周知, 除稳定性研究外,Hopf分岔分析[9-11]是完善非线性系统动态行为的有效手段．通过分析分岔特征可捕捉动态系统的内在性质, 获知系统本质特性的关键影响因素, 从而掌握更全面的系统动态特性信息．Huang[12]、 Tao[13]、 Zhang[14]等将Hopf分岔用于非线性动力学, 通过确定Hopf分岔发生的条件,判断系统是否存在周期振荡．随着分数阶微积分的进一步发展,越来越多的学者开始聚焦于分数阶模型的Hopf分岔研究．Wang等[15]通过探讨阶次不对称的生态流行病学模型的稳定性, 指出分数阶阶次和时滞可以改变系统的稳定性并产生Hopf分岔．然而, 在现实生态系统中, 单个时滞已无法完整地刻画捕食者与食饵的动力学行为; 因此, 本文拟提出一个双时滞捕食者-食饵模型, 以捕食者的妊娠期时滞和疾病的潜伏期时滞作为分岔参数研究其动力学行为．
定义1[1]定义Caputo分数阶导数
其中分数阶次α>0,m-1≤α 考虑如下分数阶双时滞传染病系统的捕食者-食饵模型: (1) 其中S(t),I(t),Y(t)分别表示易感染食饵数量、已感染食饵数量和捕食者数量;a,k分别为食饵的内在增长率和承载能力;b为食饵之间的疾病传染率;p1,p2为捕食系数;q为食饵转化成捕食者的转化率;τ1,τ2分别为捕食者的妊娠期和疾病的潜伏期;c,d分别为已感染食饵和捕食者的死亡率．经计算得到系统(1)的唯一正平衡点为E*(S*,I*,Y*), 其中 (2) 将系统(1)在平衡点E*(S*,I*,Y*)处作线性化处理可得 (3) 系统(3)的特征矩阵为 其中s为特征根, 其余参数如下: 令det(Δ(s))=0, 故系统(1)的特征方程可改写为 A1(s)+A2(s)e-sτ1+A3(s)e-sτ2+A4(s)e-s(τ1+τ2)=0, (4) 其中 A1(s)=sα1+α2+α3+dsα1+α2+a22sα1+α3-a11sα2+α3+a22dsα1-a11dsα2-a11a22sα3-a11a22d, 将捕食者的妊娠期时滞τ1作为分岔参数寻找Hopf分岔条件, 此时疾病的潜伏期时滞τ2为正常数．由式(4)可得 A1(s)+(A2(s)+A4(s)e-sτ2)e-sτ1+A3(s)e-sτ2=0. (5) (6) 其中 G1=Re[A2]+Re[A4]cos(ωτ2)+Im[A4]sin(ωτ2), 并且 a11a22a33-a11a23a32-a12a23a31+a13a22a31, 于是, 有 (7) 由sin2(ωτ1)+cos2(ωτ1)=1, 得 (8) 根据cos(ωτ1)=f1(ω)可得 (9) 假设式(8)至少存在一个正实根, 则分岔点 (10) 当τ1=0时, 假设α1=α2=α3=α, 令sα=λ, 改写特征方程(5)为 λ3+B1λ2+B2λ+B3=0, (11) 其中 B1=-a33-u+d+a22-a11, ua11a33-a12a21a33+a13a21a32-ua13a21. 令φ1=1,φ2=B1B2-B3,φ3=B3φ2, 对式(11)作如下假设: (H1)φ1>0,φ2>0,φ3>0. 由假设(H1)易见方程(11)的根均具负实部．故当τ1=0时, 分数阶系统(1)的平衡点E*(S*,I*,Y*)是局部渐近稳定的． 当τ1>0时, 改写特征方程(5)为 E(s)+F(s)e-sτ1=0, (12) 其中E(s)=A1(s)+A2(s),F(s)=A3(s)+A4(s).取E(s)=E1+iE2,F(s)=F1+iF2, 有 E1+iE2+(F1+iF2)e-sτ1=0. (13) (14) 基于式(14)得 (15) 经计算, 方程(15)无实根．故当τ1>0时, 分数阶系统(1)的平衡点E*(S*,I*,Y*)失去稳定性, 即系统产生Hopf分岔． 引理1假设s(τ1)=ε(τ1)+iω(τ1)是方程(5)的根, 且满足ε(τ0)=0,ω(τ0)=ω0.若M1N1+M2N2>0, 则穿越条件满足以下形式: 其中 M1=Re[A′1(iω0)]+(Re[A′2(iω0)]-τ0Re[A2(iω0)])cos(ω0τ0)+(Re[A′3(iω0)]- 证明 在方程(5)两端对τ1进行求导得 于是, 有 将s=iω0和τ=τ0代入其中得 进一步可得 定理1如果假设(H1)成立,M1N1+M2N2>0且方程(15)无实根, 那么可得如下结论: i) 系统(1)的正平衡点E*(S*,I*,Y*)在τ1∈[0,τ0)处是局部渐近稳定的; ii) 当τ1=τ0时, 系统(1)的正平衡点E*(S*,I*,Y*)附近出现Hopf分岔． 首先将疾病的潜伏期时滞τ2作为分岔参数来寻找Hopf分岔条件, 此时捕食者妊娠期时滞τ1为正常数．由式(4)可得 A1(s)+(A3(s)+A4(s)e-sτ1)e-sτ2+A2(s)e-sτ1=0. (16) (17) 其中 G′1=Re[A4]cos(ωτ1)+Im[A4]sin(ωτ1)+Re[A3], 进一步可得 (18) 由sin2(ωτ2)+cos2(ωτ2)=1, 得 (19) 根据cos(ωτ2)=g1(ω)可得 (20) 假设方程(19)至少存在一个正实根, 则此时的分岔点为 (21) 当τ2=0时, 假设α1=α2=α3=α, 令sα=λ, 则特征方程(16)可改写为 λ3+C1λ2+C2λ+C3=0, (22) 其中 C1=-a33-u+d+a22-a11, 令φ1=1,φ2=C1C2-C3,φ3=C3φ2, 对式(22)作如下假设: (H2)φ1>0,φ2>0,φ3>0. 由假设(H2)易见方程(22)的根均具负实部．故当τ2=0时, 分数阶系统(1)的平衡点E*(S*,I*,Y*)是局部渐近稳定的． 当τ2>0时, 改写特征方程(18)为 P(s)+Q(s)e-sτ2=0, (23) 其中P(s)=A1(s)+A3(s),Q(s)=A2(s)+A4(s).取P(s)=P1+iP2,Q(s)=Q1+iQ2, 则有 P1+iP2+(Q1+iQ2)e-sτ2=0. (24) (25) 基于式(25)得 (26) 经计算, 方程(26)无实根．故当τ2>0时, 分数阶系统(1)的平衡点E*(S*,I*,Y*)失去稳定性, 即系统产生Hopf分岔． 引理2假设s(τ2)=ε(τ2)+iω(τ2)是方程(16)的根, 且满足ε(σ0)=0,ω(σ0)=ω0.若M′1N′1+M′2N′2>0, 则穿越条件满足以下形式: 其中 M′1=Re[A′1(iω0)]+(Re[A′3(iω0)]-σ0Re[A3(iω0)])cos(ω0σ0)+(Re[A′2(iω0)]- 证明 在方程(16)两端对τ2进行求导, 得 于是,有 将s=iω0和σ=σ0代入其中得 进一步可得 定理2如果假设(H2)和M′1N′1+M′2N′2>0成立且方程(26)无实根, 那么可得以下结论: i) 系统(1)的正平衡点E*(S*,I*,Y*)在τ2∈[0,σ0)处是局部渐近稳定的; ii) 当τ2=σ0时, 系统(1)的正平衡点E*(S*,I*,Y*)附近出现Hopf分岔． 应用MATLAB软件进行数值仿真, 验证本文所提模型的正确性和合理性．设定系统(1)的参数a=0.5,b=1,c=0.25,d=0.5,k=1,p1=0.125,p2=6,q=0.75. 当τ2一定, 仅以捕食者妊娠期时滞τ1作为分岔参数时, 易证系统(1)的参数满足假设(H1),M1N1+M2N2>0且方程(15)无实根, 故该系统在初始状态下正平衡点E*(S*,I*,Y*)保持稳定．固定时滞τ2=21, 令α1=0.8,α2=0.99,α3=0.9,初始值为(0.7,0.1,0.08)．根据式(9)计算可得分岔阈值τ0=2.55．选取τ1=1.5,3.5, 系统(1)的波形图和相位图分别如图1～2所示．由图1可见,易感食饵S(t)、已感染食饵I(t)和捕食者Y(t)均能收敛至平衡点E*(S*,I*,Y*)附近, 表明系统(1)的正平衡点E*(S*,I*,Y*)是局部渐近稳定的．由图2可见,易感食饵S(t)、已感染食饵I(t)和捕食者Y(t)产生震荡, 表明系统(1)的平衡点E*(S*,I*,Y*)失去稳定性, 并产生Hopf分岔． 图1 τ1=1.5τ0=2.55时系统(1)的波形图与相位图Fig.1 Waveform diagram and phase diagram of system (1) at τ1=1.5τ0=2.55 图2 τ1=3.5τ0=2.55时系统(1)的波形图与相位图Fig.2 Waveform diagram and phase diagram of system (1) at τ1=3.5τ0=2.55 当τ1一定, 仅以传染病潜伏期时滞τ2作为分岔参数时, 易证系统(1)的参数满足假设(H2),M′1N′1+M′2N′2>0且方程(26)无实根, 故该系统在初始状态下正平衡点E*(S*,I*,Y*)保持稳定．固定时滞τ1=3.5, 令α1=0.7,α2=0.99,α3=0.9, 初始值为(0.7,0.1,0.08)．根据式(20)计算可得分岔阈值σ0=19.15．选取τ2=17, 22, 系统(1)的波形图和相位图分别如图3～4所示．由图3可见,易感食饵S(t)、已感染食饵I(t)和捕食者Y(t)均收敛至平衡点E*(S*,I*,Y*)附近, 表明系统(1)的正平衡点E*(S*,I*,Y*)是局部渐近稳定的．由图4可见,易感食饵S(t)、已感染食饵I(t)和捕食者Y(t)产生震荡, 表明分数阶系统(1)的平衡点E*(S*,I*,Y*)失去稳定性, 并产生Hopf分岔． 图3 τ2=17σ0=19.15时系统(1)的波形图与相位图Fig.3 Waveform diagram and phase diagram of system (1) at τ2=17σ0=19.15 图4 τ2=22σ0=19.15时系统(1)的波形图与相位图Fig.4 Waveform diagram and phase diagram of system (1) at τ2=22σ0=19.15 本文构建了一个具有双时滞分数阶的三维捕食者-食饵模型, 通过调节捕食者的妊娠期时滞和传染病的潜伏期时滞给出模型的稳定性和Hopf分岔条件, 并且确定了2种不同时滞引起的分岔判据．数值仿真结果表明, 捕食者的妊娠期时滞与被捕食者的疾病潜伏期时滞可影响捕食者-食饵模型的动力学．当妊娠期时滞τ1小于临界值τ0时, 整个系统处于稳定状态, 捕食者与被捕食者的种群密度在一段时间后都会趋于平衡点E*(S*,I*,Y*)附近; 而当τ1大于τ0时, 系统失去稳定性, 捕食者与被捕食者的种群密度在一定范围内产生周期震荡, 被捕食者面临泛滥或趋于濒危的可能性．类似地, 传染病时滞τ2的存在也会导致系统失去稳定性并产生Hopf分岔, 物种失去平衡状态．今后将在模型中进一步加入控制器, 制定合理的疾病预防措施和方案, 以期控制或消除疾病．
A2(s)=-a33sα1+α2+(a23a32-a22a33)sα1+(a11a33+a13a31)sα2+a11a22a33-a11a23a32-a12a23a31+a13a22a31,
A3(s)=-usα1+α3+(ua11+a12a21)sα3-udsα1+ua11d+a12a21d,
A4(s)=ua33sα1-ua11a33-a12a21a33-a12a21a33+a13a21a32-ua13a31.3.1 捕食者的妊娠期时滞τ1
G2=Im[A2]-Re[A4]sin(ωτ2)+Im[A4]cos(ωτ2),
G3=-Re[A1]-Re[A3]cos(ωτ2)-Im[A3]sin(ωτ2),
G4=-Im[A1]+Re[A3]sin(ωτ2)-Im[A3]cos(ωτ2),
B2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
B3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+a12a21d-
Re[A3(iω0)])cos(ω0τ0)+(Re[A′4(iω0)]-(τ0+τ2)Re[A4(iω0)])cos(ω0τ0),
M2=(τ0Im[A2(iω0)]-Im[A′2(iω0)])sin(ω0τ0)+(Im[A3(iω0)]-Im[A′3(iω0)])sin(ω0τ0)+
((τ0+τ2)Im[A4(iω0)]-Im[A′4(iω0)])sin(ω0τ0),
N1=ω0Re[A2(iω0)]sin(ω0τ0)+ω0Re[A4(iω0)]sin(ω0(τ0+τ2)),
N2=ω0Im[A2(iω0)]cos(ω0τ0)+ω0Im[A4(iω0)]cos(ω0(τ0+τ2)).
3.2 疾病的潜伏期时滞τ2
G′2=Im[A4]cos(ωτ1)-Re[A4]sin(ωτ1)+Im[A3],
G′3=-Re[A2]cos(ωτ1)-Im[A2]sin(ωτ1)-Re[A1],
G′4=-Im[A2]cos(ωτ1)+Re[A2]sin(ωτ1)-Im[A1].
C2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
C3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+
a12a21d-ua11a33-a12a21a33+a13a21a32-ua13a21.
τ1Re[A2(iω0)])cos(ω0τ1)+(Re[A′4(iω0)]-τ1Re[A4(iω0)]-τ2Re[A4(iω0)])cos(ω0(σ0+τ1)),
M′2=(σ0Im[A3(iω0)]-Im[A′3(iω0)])sin(ω0σ0)+(τ1Im[A2(iω0)]-Im[A′2(iω0)])sin(ω0τ1)+
(τ1Im[A4(iω0)]-Im[A′4(iω0)]+σ0Im[A4(iω0)])sin(ω0(σ0+τ1)),
N′1=ω0Re[A3(iω0)]sin(ω0σ0)+ω0Re[A4(iω0)]sin(ω0(σ0+τ1)),
N′2=ω0Im[A3(iω0)]cos(ω0σ0)+ω0Im[A4(iω0)]cos(ω0(σ0+τ1)).
