高燕鑫， 石劍平

(昆明理工大學(xué) 理學(xué)院， 云南 昆明 650500)

現(xiàn)實問題里很多系統(tǒng)的變化過程不僅依賴當(dāng)前時刻的狀態(tài)，還依賴于過去某個時刻或某段時刻的狀態(tài)，這種特性稱為時滯。由于物質(zhì)和能量在變化的時候往往不能瞬間傳遞，時滯發(fā)生在幾乎所有類型的自然和社會系統(tǒng)中[1-2]。計算機病毒模型[3-4]是一種以時間為狀態(tài)變量的微分方程系統(tǒng)，引入時滯來描述和分析其動力學(xué)特征早就成為研究者的共識[5-6]，并取得了很多有價值的成果。

近幾十年，分?jǐn)?shù)階微積分[7]在眾多科學(xué)與工程領(lǐng)域得到了成功的應(yīng)用。由于能夠描述存在于醫(yī)學(xué)、物理學(xué)、生物學(xué)、工程應(yīng)用等系統(tǒng)中所固有的記憶和遺傳特性[8-10]，相比整數(shù)階微積分，分?jǐn)?shù)階微積分為系統(tǒng)性能的描述提供了更為豐富的自由度，有助于更加科學(xué)有效地研究系統(tǒng)的穩(wěn)定性以及分岔等動力學(xué)性質(zhì)[11-13]。本文基于前人的研究結(jié)果，對一類時滯計算機病毒傳播模型[14]引入分?jǐn)?shù)階進行描述，通過增加系統(tǒng)的時滯項，獲得一個分?jǐn)?shù)階時滯計算機病毒模型。并以時滯項作為控制參數(shù)，研究系統(tǒng)在正平衡點處的穩(wěn)定性以及出現(xiàn)Hopf分岔的條件。

1 模型描述

由于狀態(tài)變量在分?jǐn)?shù)階的情況下可以考慮不同的變化率，因此本文對文獻[14]研究的一類時滯計算機病毒傳播模型引入不同的分?jǐn)?shù)階，并對最后一個方程也加入時滯項，獲得下述的改進系統(tǒng)模型

(1)

研究分?jǐn)?shù)階微分方程系統(tǒng)的穩(wěn)定性和分岔問題，目前并沒有完善的理論體系。文獻[15]利用線性化近似系統(tǒng)在平衡點處的穩(wěn)定性來研究原系統(tǒng)的穩(wěn)定性和分岔問題，并給出了明確的定義，說明系統(tǒng)時滯項τ滿足一定的條件，在平衡點處會出現(xiàn)Hopf分岔。此研究思想的核心是以時滯τ作為控制參數(shù)分析分?jǐn)?shù)階系統(tǒng)平衡點的穩(wěn)定性。

定義1[15]對于以下n維時滯分?jǐn)?shù)階系統(tǒng)

(2)

0

注1：定義1中，條件(1)說明當(dāng)時滯τ=0時，分?jǐn)?shù)階微分方程系統(tǒng)(2)在0

下文將根據(jù)這個定義，結(jié)合系統(tǒng)(1)研究其在正平衡點處的穩(wěn)定性和Hopf分岔。

2 Hopf分岔分析

2.1 平衡點及基本再生數(shù)

計算系統(tǒng)(1)的平衡點為

由于系統(tǒng)參數(shù)均大于零，E0是系統(tǒng)(2)的一個無病平衡點。但是當(dāng)I(t)=0時，系統(tǒng)中已沒有感染節(jié)點。本文主要研究存在感染節(jié)點的平衡狀態(tài)下系統(tǒng)的穩(wěn)定性，為此，需要考慮系統(tǒng)參數(shù)滿足什么條件時，E1為正平衡點。由于模型(1)的參數(shù)較多，下面通過引入基本再生數(shù)來分析系統(tǒng)存在正平衡點的條件。

基本再生數(shù)是傳染病模型的重要概念，描述已感染的病人在平均患病時間內(nèi)感染易感染者的人數(shù)，記作R0。在計算機病毒傳播模型中，可用基本再生數(shù)描述病毒傳播的情況：若R0>1時，病毒侵入易感節(jié)點；若R0<1，經(jīng)過一定的時間，病毒將被徹底清除。下面采用再生矩陣方法[16]計算 。

系統(tǒng)(1)對應(yīng)的整數(shù)階系統(tǒng)為

(3)

令U(t)=(I(t),S(t),R(t))T，則當(dāng)τ=0時，(3)式可以表示為

(4)

其中

f(U(t))描述新感染個體的比率，ν(U(t))則表示轉(zhuǎn)移比率。

計算f(U(t))和ν(U(t))在無病平衡點E0處的Jacobian矩陣

(5)

由文獻[16]可知再生矩陣為FV-1，基本再生數(shù)R0是再生矩陣的譜半徑ρ(FV-1)，即再生矩陣FV-1的特征值λ0模的最大值。由(5)式可得

其中I是單位矩陣。

故系統(tǒng)(3)的基本再生數(shù)為

(6)

定理1 若R0>1，則系統(tǒng)(1)存在唯一的正平衡點E*=E1，若R0≤1，則系統(tǒng)(1)不存在正平衡點。

證明由(6)式，系統(tǒng)(1)的平衡點E1的各分量可寫為

(7)

顯然恒有S*>0,當(dāng)R0>1時，I*>0，R*>0，即系統(tǒng)(1)存在唯一的正平衡點E*=E1(S*,I*,R*)，當(dāng)R0≤1時，有I*≤0，即系統(tǒng)不存在正平衡點。

2.2 正平衡點的穩(wěn)定性及 Hopf分岔分析

在本節(jié)中，以時滯作為參數(shù)分析系統(tǒng)(1)發(fā)生Hopf分岔的條件。首先，在正平衡點E*做變換

x(t)=S(t)-S*,y(t)=I(t)-I*,z(t)=R(t)-R*,

則系統(tǒng)(1)可以轉(zhuǎn)化為

(8)

系統(tǒng)(8)在原點處的線性化系統(tǒng)為

(9)

其中

a11=-(μ+k),

a31=k,a12=-βS*,

a32=γ,a13=-βI*,

a33=-μ。a21=βI*,

(10)

下面根據(jù)定義1，分析系統(tǒng)(1)在正平衡點發(fā)生Hopf分岔的條件。

條件1τ=0時,線性化系統(tǒng)(9)系數(shù)矩陣的特征根分析。

當(dāng)τ=0時，系統(tǒng)(9)的系數(shù)矩陣對應(yīng)的特征方程為

λ3-(a11+a13+a33)λ2+(a11a33+a13a33-a12a21)λ+a12a21a33=0,

(11)

計算得

為了分析τ=0時(11)式的根λj(j=1,2,3)的取值情況，不妨作如下假設(shè):

(H1) -(a11+a13+a33) ·(a11a33+a13a33-a12a21) -a12a21a33>0。

引理1 如果假設(shè)(H1)成立，那么線性化系統(tǒng)(9)系數(shù)矩陣的所有特征根λj(j=1,2,3)均具有負(fù)實部。

條件2τ>0時系統(tǒng)(1)產(chǎn)生Hopf分岔的臨界值分析。

對系統(tǒng)(9)的兩側(cè)分別作Laplace變換[17]，得

(12)

其中

Δ(s)被看作系統(tǒng)(9)的特征矩陣，其中系數(shù)a11、a12、a13、a21、a31、a32、a33由(10)式?jīng)Q定，故系統(tǒng)(9)的特征方程為

P1(s)+P2(s)·e-sτ=0,

(13)

其中

P1(s)=sq1+q2+q3-a11sq2+q3-a33sq1+q2+a11a33sq2,

P2(s)=a13a33sq2-a13sq2+q3-a12a21sq3+a12a21a33。

由定義1可知，系統(tǒng)(1)發(fā)生Hopf分岔的條件之一是當(dāng)τ=τ0時，線性化系統(tǒng)(9)的特征方程有一對純虛根。不妨假設(shè)s=ωi=ω(cos(π/2)+isin(π/2))(ω>0，i為虛數(shù)單位)是(13)式的根，代入分離實部和虛部并化簡，將cos(ωτ)和sin(ωτ)看作未知項，可得到以下方程組

(14)

其中A2和B2分別是(13)式實部和虛部經(jīng)過整理后，cos(ωτ)和sin(ωτ)項的系數(shù)部分，A1和B1則是剩余的常數(shù)項。

求解方程組(14)可得

(15)

當(dāng)系統(tǒng)(1)中的所有參數(shù)給定時，聯(lián)立公式sin2(ωτ)+cos2(ωτ)=1，可以計算出ω的值。代入(15)式求解得

根據(jù)時滯的實際意義，主要關(guān)注出現(xiàn)Hopf分岔的最小正值，因此，將該分岔點描述為

(17)

條件3 橫截條件分析。

根據(jù)隱函數(shù)求導(dǎo)法則，在(13)式兩邊分別對τ求導(dǎo)，可得到

(18)

其中

通過計算，可以得到

其中Mj和Nj(j=1,2)分別是M(s)、N(s)的實部和虛部。其解析過程略。

給出以下假設(shè)：

(H2)M1N1+M2N2>0,

則得到橫截條件成立的引理。

引理2 若假設(shè)(H2)成立，令s(τ)=γ(τ)+iω(τ)是(13)式在τ=τ0附近滿足γ(τj)=0，ω(τj)=ω0的根，下面的橫截條件成立

(19)

綜上所述，由定義1得到以下結(jié)論：

定理2 假設(shè)(H1)和(H2)成立，給定參數(shù)組(p,b,μ,β,γ,k)，則當(dāng)τ=τ0時，系統(tǒng)(1)在正平衡點E*處產(chǎn)生Hopf分岔，τ0為(17)式定義的時滯臨界點。

3 數(shù)值模擬

在本節(jié)中，基于文獻[18]介紹的Adama-Bashforth-Moulton預(yù)估-校正方法，給出數(shù)值實例來驗證前述理論分析方法的可行性和結(jié)果的正確性，其中步長取h=0.01。

為了更具有對比性，模擬所用的系統(tǒng)參數(shù)均來自文獻[14]，p=0.9，b=1，β=0.4，γ=0.1，μ=0.1，k=0.1，則系統(tǒng)(1)為

(20)

為了驗證此結(jié)果的正確性，本文模擬了兩種情況：選取初值為(0.2,4.5,5.4)，分?jǐn)?shù)階為q1=0.92,q2=0.95,q3=0.98，取時滯τ=0.98<τ0=1.024 8時，系統(tǒng)(20)在平衡點(S*,I*,R*)=(0.5,4.0,5.5)處是漸近穩(wěn)定的(見圖1)。此結(jié)果說明隨著時間的推移，易感節(jié)點、感染節(jié)點和恢復(fù)節(jié)點均趨于穩(wěn)定值，雖然系統(tǒng)仍然存在感染節(jié)點，但是數(shù)量穩(wěn)定，有利于采取恰當(dāng)措施對系統(tǒng)進行干預(yù)，徹底清除病毒，恢復(fù)健康的網(wǎng)絡(luò)狀態(tài)。當(dāng)取時滯τ=1.07>τ0=1.024 8時，系統(tǒng)(20)在平衡點(0.5,4.0,5.5)處發(fā)生Hopf分岔(見圖2)，說明此時系統(tǒng)在平衡點處是不穩(wěn)定的，易感節(jié)點、感染節(jié)點和恢復(fù)節(jié)點的數(shù)量產(chǎn)生了隨著時間t的推移而出現(xiàn)的周期振蕩，這對于清除病毒，調(diào)節(jié)和控制網(wǎng)絡(luò)系統(tǒng)恢復(fù)到健康狀態(tài)是極為不利的。

圖1 τ=0.98<τ0時，系統(tǒng)(20)的相圖和各分量的波形圖

圖2 τ=1.07>τ0時，系統(tǒng)(20)的相圖和各分量的波形圖

下面討論分?jǐn)?shù)階的變化對系統(tǒng)(20)的時滯分岔臨界點τ0的影響，具體的做法是保持其中兩個階不變，考察另一階變化對于分岔點的影響。由于時滯τ表示的是因計算機病毒的潛伏期造成的延遲，故處于恢復(fù)狀態(tài)的計算機節(jié)點上的分?jǐn)?shù)階q3的變化對于分岔臨界點基本不產(chǎn)生影響(見表1)，但是分?jǐn)?shù)階q1、q2的變化都對分岔臨界點有較大的影響。q1從0.5到0.6的變化過程中，時滯τ0隨著q1的增大而增大，在0.6到0.7之間出現(xiàn)轉(zhuǎn)折，之后隨著q1的不斷增大，時滯τ0越來越小(見表1，圖3)。而對于分?jǐn)?shù)階q2，在0.5到1之間，分岔點一直隨著q2的增大而增大(見表1，圖4)。顯然，可以通過調(diào)節(jié)各變量分?jǐn)?shù)階的取值來改變系統(tǒng)分岔臨界值的大小，從而調(diào)節(jié)系統(tǒng)的穩(wěn)定域。

表1 qi變化對于系統(tǒng)(20)Hopf分岔臨界值(ω0,τ0)的影響

圖3 q2=0.95,q3=0.98時, 系統(tǒng)(20)中τ0隨q1的變化 圖4 q1=0.92,q3=0.98時, 系統(tǒng)(20)中τ0隨q2的變化

上述模擬結(jié)果顯示，與文獻[15]研究的整數(shù)階系統(tǒng)比較，分?jǐn)?shù)階的引入延遲了系統(tǒng) Hopf分岔的發(fā)生，放大了穩(wěn)定區(qū)間。通過調(diào)節(jié)分?jǐn)?shù)階的大小，可以在一定范圍內(nèi)有效控制系統(tǒng)正平衡點的穩(wěn)定域。分?jǐn)?shù)階系統(tǒng)獲得了比整數(shù)階系統(tǒng)更為靈活的控制方式。

4 結(jié)論

本文通過研究一類具有不同分?jǐn)?shù)階的時滯計算機病毒傳播模型的Hopf分岔，討論了系統(tǒng)正平衡點的穩(wěn)定性問題。由于該模型參數(shù)較多，首先引入流行病學(xué)中基本再生數(shù)的概念討論了系統(tǒng)存在正平衡點的條件。繼而以時滯為分岔參數(shù)研究系統(tǒng)在正平衡點的穩(wěn)定性，分析了該模型發(fā)生Hopf分岔的3個顯式條件，結(jié)果表明時滯是造成系統(tǒng)不穩(wěn)定的主要因素之一，而分?jǐn)?shù)階的引入不僅以更多的自由度豐富了系統(tǒng)的性能，也影響時滯的變化，繼而影響系統(tǒng)正平衡點的穩(wěn)定域。為了驗證理論分析的正確性，選擇了恰當(dāng)?shù)膮?shù)做數(shù)值模擬，結(jié)果說明時滯臨界值確實是分?jǐn)?shù)階計算機病毒模型出現(xiàn)Hopf分岔的一個分水嶺。在給定參數(shù)值的情況下，時滯的大小是決定系統(tǒng)穩(wěn)定性的重要因素之一，對于調(diào)節(jié)和控制系統(tǒng)從病毒侵害狀態(tài)恢復(fù)到健康狀態(tài)有直接的影響。此外，數(shù)值模擬也說明了分?jǐn)?shù)階變化引起系統(tǒng)Hopf分岔臨界值的變化，進一步解釋了分?jǐn)?shù)階對于模型穩(wěn)定域控制的有效性。

需要指出的是，分?jǐn)?shù)階模型線性化以后，通過Laplace變換得到對應(yīng)的特征方程，其純虛根的求解是基于假設(shè)，反代入方程后通過計算求出的，并沒有采用理論的方法獲得嚴(yán)密的證明，這個具有理論意義的問題將在后續(xù)的工作中加以研究。