高燕鑫， 石劍平

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

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

近幾十年，分?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)定性以及分岔等動(dòng)力學(xué)性質(zhì)[11-13]。本文基于前人的研究結(jié)果，對(duì)一類時(shí)滯計(jì)算機(jī)病毒傳播模型[14]引入分?jǐn)?shù)階進(jìn)行描述，通過增加系統(tǒng)的時(shí)滯項(xiàng)，獲得一個(gè)分?jǐn)?shù)階時(shí)滯計(jì)算機(jī)病毒模型。并以時(shí)滯項(xiàng)作為控制參數(shù)，研究系統(tǒng)在正平衡點(diǎn)處的穩(wěn)定性以及出現(xiàn)Hopf分岔的條件。

1 模型描述

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

(1)

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

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

(2)

0

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

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

2 Hopf分岔分析

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

計(jì)算系統(tǒng)(1)的平衡點(diǎn)為

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

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

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

(3)

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

(4)

其中

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

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

(5)

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

其中I是單位矩陣。

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

(6)

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

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

(7)

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

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

在本節(jié)中，以時(shí)滯作為參數(shù)分析系統(tǒng)(1)發(fā)生Hopf分岔的條件。首先，在正平衡點(diǎn)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)在原點(diǎn)處的線性化系統(tǒng)為

(9)

其中

a11=-(μ+k),

a31=k,a12=-βS*,

a32=γ,a13=-βI*,

a33=-μ。a21=βI*,

(10)

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

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

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

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

(11)

計(jì)算得

為了分析τ=0時(shí)(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ù)實(shí)部。

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

對(duì)系統(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時(shí)，線性化系統(tǒng)(9)的特征方程有一對(duì)純虛根。不妨假設(shè)s=ωi=ω(cos(π/2)+isin(π/2))(ω>0，i為虛數(shù)單位)是(13)式的根，代入分離實(shí)部和虛部并化簡，將cos(ωτ)和sin(ωτ)看作未知項(xiàng)，可得到以下方程組

(14)

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

求解方程組(14)可得

(15)

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

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

(17)

條件3 橫截條件分析。

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

(18)

其中

通過計(jì)算，可以得到

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

給出以下假設(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時(shí)，系統(tǒng)(1)在正平衡點(diǎn)E*處產(chǎn)生Hopf分岔，τ0為(17)式定義的時(shí)滯臨界點(diǎn)。

3 數(shù)值模擬

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

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

(20)

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

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

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

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

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

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

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

4 結(jié)論

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

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