李 健, 肖 敏, 周 帥

(南京郵電大學自動化學院, 南京 210023)

近年來, 分數(shù)階微積分發(fā)展迅速并被廣泛應用于多個領域[1-2]．由于分數(shù)階微積分的階數(shù)可取非整數(shù), 故基于分數(shù)階微積分的分數(shù)階模型能更準確地描述實時系統(tǒng)的性質特征和行為．運用分數(shù)階微分方程解決現(xiàn)實工程問題備受關注[3-4]．例如, Rihan等[5]提出一個時滯分數(shù)階捕食者-食餌模型, 并討論了模型的傳統(tǒng)響應和解的存在性問題．傳染病一直是影響生態(tài)系統(tǒng)的重大問題之一, 傳染病的爆發(fā)會導致生態(tài)系統(tǒng)紊亂,捕食者與食餌將失去平衡狀態(tài), 故研究傳染病的動力學行為顯得尤為重要[6-7]．目前, 分數(shù)階傳染病模型已成為研究熱點．Singh等[8]將生態(tài)傳染病學模型引入分數(shù)階系統(tǒng),研究了傳染病系統(tǒng)平衡點的穩(wěn)定性和解的性質．眾所周知, 除穩(wěn)定性研究外,Hopf分岔分析[9-11]是完善非線性系統(tǒng)動態(tài)行為的有效手段．通過分析分岔特征可捕捉動態(tài)系統(tǒng)的內在性質, 獲知系統(tǒng)本質特性的關鍵影響因素, 從而掌握更全面的系統(tǒng)動態(tài)特性信息．Huang[12]、 Tao[13]、 Zhang[14]等將Hopf分岔用于非線性動力學, 通過確定Hopf分岔發(fā)生的條件,判斷系統(tǒng)是否存在周期振蕩．隨著分數(shù)階微積分的進一步發(fā)展,越來越多的學者開始聚焦于分數(shù)階模型的Hopf分岔研究．Wang等[15]通過探討階次不對稱的生態(tài)流行病學模型的穩(wěn)定性, 指出分數(shù)階階次和時滯可以改變系統(tǒng)的穩(wěn)定性并產生Hopf分岔．然而, 在現(xiàn)實生態(tài)系統(tǒng)中, 單個時滯已無法完整地刻畫捕食者與食餌的動力學行為; 因此, 本文擬提出一個雙時滯捕食者-食餌模型, 以捕食者的妊娠期時滯和疾病的潛伏期時滯作為分岔參數(shù)研究其動力學行為．

1 預備知識

定義1[1]定義Caputo分數(shù)階導數(shù)

其中分數(shù)階次α>0,m-1≤α

2 模型的建立

考慮如下分數(shù)階雙時滯傳染病系統(tǒng)的捕食者-食餌模型:

(1)

其中S(t),I(t),Y(t)分別表示易感染食餌數(shù)量、已感染食餌數(shù)量和捕食者數(shù)量;a,k分別為食餌的內在增長率和承載能力;b為食餌之間的疾病傳染率;p1,p2為捕食系數(shù);q為食餌轉化成捕食者的轉化率;τ1,τ2分別為捕食者的妊娠期和疾病的潛伏期;c,d分別為已感染食餌和捕食者的死亡率．經計算得到系統(tǒng)(1)的唯一正平衡點為E*(S*,I*,Y*), 其中

(2)

將系統(tǒng)(1)在平衡點E*(S*,I*,Y*)處作線性化處理可得

(3)

系統(tǒng)(3)的特征矩陣為

其中s為特征根, 其余參數(shù)如下:

令det(Δ(s))=0, 故系統(tǒng)(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,
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 雙時滯影響下系統(tǒng)的穩(wěn)定性及Hopf分岔分析

3.1 捕食者的妊娠期時滯τ1

將捕食者的妊娠期時滯τ1作為分岔參數(shù)尋找Hopf分岔條件, 此時疾病的潛伏期時滯τ2為正常數(shù)．由式(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),
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),

并且

a11a22a33-a11a23a32-a12a23a31+a13a22a31,

于是, 有

(7)

由sin2(ωτ1)+cos2(ωτ1)=1, 得

(8)

根據(jù)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,
B2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
B3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+a12a21d-

ua11a33-a12a21a33+a13a21a32-ua13a21.

令φ1=1,φ2=B1B2-B3,φ3=B3φ2, 對式(11)作如下假設:

(H1)φ1>0,φ2>0,φ3>0.

由假設(H1)易見方程(11)的根均具負實部．故當τ1=0時, 分數(shù)階系統(tǒng)(1)的平衡點E*(S*,I*,Y*)是局部漸近穩(wěn)定的．

當τ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時, 分數(shù)階系統(tǒng)(1)的平衡點E*(S*,I*,Y*)失去穩(wěn)定性, 即系統(tǒng)產生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)]-
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)).

證明 在方程(5)兩端對τ1進行求導得

于是, 有

將s=iω0和τ=τ0代入其中得

進一步可得

定理1如果假設(H1)成立,M1N1+M2N2>0且方程(15)無實根, 那么可得如下結論: i) 系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)在τ1∈[0,τ0)處是局部漸近穩(wěn)定的; ii) 當τ1=τ0時, 系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)附近出現(xiàn)Hopf分岔．

3.2 疾病的潛伏期時滯τ2

首先將疾病的潛伏期時滯τ2作為分岔參數(shù)來尋找Hopf分岔條件, 此時捕食者妊娠期時滯τ1為正常數(shù)．由式(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],
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].

進一步可得

(18)

由sin2(ωτ2)+cos2(ωτ2)=1, 得

(19)

根據(jù)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,
C2=a22a32-a22a33+a11a33+a13a31-ud+ua11+a12a21+ua33+a22d-a11d,
C3=a11a22a33-a11a23a32-a12a23a31+a13a22a31+ua11d+
a12a21d-ua11a33-a12a21a33+a13a21a32-ua13a21.

令φ1=1,φ2=C1C2-C3,φ3=C3φ2, 對式(22)作如下假設:

(H2)φ1>0,φ2>0,φ3>0.

由假設(H2)易見方程(22)的根均具負實部．故當τ2=0時, 分數(shù)階系統(tǒng)(1)的平衡點E*(S*,I*,Y*)是局部漸近穩(wěn)定的．

當τ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時, 分數(shù)階系統(tǒng)(1)的平衡點E*(S*,I*,Y*)失去穩(wěn)定性, 即系統(tǒng)產生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)]-
τ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)).

證明 在方程(16)兩端對τ2進行求導, 得

于是,有

將s=iω0和σ=σ0代入其中得

進一步可得

定理2如果假設(H2)和M′1N′1+M′2N′2>0成立且方程(26)無實根, 那么可得以下結論: i) 系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)在τ2∈[0,σ0)處是局部漸近穩(wěn)定的; ii) 當τ2=σ0時, 系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)附近出現(xiàn)Hopf分岔．

4 數(shù)值模擬與仿真

應用MATLAB軟件進行數(shù)值仿真, 驗證本文所提模型的正確性和合理性．設定系統(tǒng)(1)的參數(shù)a=0.5,b=1,c=0.25,d=0.5,k=1,p1=0.125,p2=6,q=0.75.

當τ2一定, 僅以捕食者妊娠期時滯τ1作為分岔參數(shù)時, 易證系統(tǒng)(1)的參數(shù)滿足假設(H1),M1N1+M2N2>0且方程(15)無實根, 故該系統(tǒng)在初始狀態(tài)下正平衡點E*(S*,I*,Y*)保持穩(wěn)定．固定時滯τ2=21, 令α1=0.8,α2=0.99,α3=0.9,初始值為(0.7,0.1,0.08)．根據(jù)式(9)計算可得分岔閾值τ0=2.55．選取τ1=1.5,3.5, 系統(tǒng)(1)的波形圖和相位圖分別如圖1～2所示．由圖1可見,易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)均能收斂至平衡點E*(S*,I*,Y*)附近, 表明系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)是局部漸近穩(wěn)定的．由圖2可見,易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)產生震蕩, 表明系統(tǒng)(1)的平衡點E*(S*,I*,Y*)失去穩(wěn)定性, 并產生Hopf分岔．

圖1 τ1=1.5τ0=2.55時系統(tǒng)(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時系統(tǒng)(1)的波形圖與相位圖Fig.2 Waveform diagram and phase diagram of system (1) at τ1=3.5τ0=2.55

當τ1一定, 僅以傳染病潛伏期時滯τ2作為分岔參數(shù)時, 易證系統(tǒng)(1)的參數(shù)滿足假設(H2),M′1N′1+M′2N′2>0且方程(26)無實根, 故該系統(tǒng)在初始狀態(tài)下正平衡點E*(S*,I*,Y*)保持穩(wěn)定．固定時滯τ1=3.5, 令α1=0.7,α2=0.99,α3=0.9, 初始值為(0.7,0.1,0.08)．根據(jù)式(20)計算可得分岔閾值σ0=19.15．選取τ2=17, 22, 系統(tǒng)(1)的波形圖和相位圖分別如圖3～4所示．由圖3可見,易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)均收斂至平衡點E*(S*,I*,Y*)附近, 表明系統(tǒng)(1)的正平衡點E*(S*,I*,Y*)是局部漸近穩(wěn)定的．由圖4可見,易感食餌S(t)、已感染食餌I(t)和捕食者Y(t)產生震蕩, 表明分數(shù)階系統(tǒng)(1)的平衡點E*(S*,I*,Y*)失去穩(wěn)定性, 并產生Hopf分岔．

圖3 τ2=17σ0=19.15時系統(tǒng)(1)的波形圖與相位圖Fig.3 Waveform diagram and phase diagram of system (1) at τ2=17σ0=19.15

圖4 τ2=22σ0=19.15時系統(tǒng)(1)的波形圖與相位圖Fig.4 Waveform diagram and phase diagram of system (1) at τ2=22σ0=19.15

5 結論

本文構建了一個具有雙時滯分數(shù)階的三維捕食者-食餌模型, 通過調節(jié)捕食者的妊娠期時滯和傳染病的潛伏期時滯給出模型的穩(wěn)定性和Hopf分岔條件, 并且確定了2種不同時滯引起的分岔判據(jù)．數(shù)值仿真結果表明, 捕食者的妊娠期時滯與被捕食者的疾病潛伏期時滯可影響捕食者-食餌模型的動力學．當妊娠期時滯τ1小于臨界值τ0時, 整個系統(tǒng)處于穩(wěn)定狀態(tài), 捕食者與被捕食者的種群密度在一段時間后都會趨于平衡點E*(S*,I*,Y*)附近; 而當τ1大于τ0時, 系統(tǒng)失去穩(wěn)定性, 捕食者與被捕食者的種群密度在一定范圍內產生周期震蕩, 被捕食者面臨泛濫或趨于瀕危的可能性．類似地, 傳染病時滯τ2的存在也會導致系統(tǒng)失去穩(wěn)定性并產生Hopf分岔, 物種失去平衡狀態(tài)．今后將在模型中進一步加入控制器, 制定合理的疾病預防措施和方案, 以期控制或消除疾?。?/p>