狄根虎 許 勇 徐 偉 顧仁財

(西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，西安 710129)

一類復(fù)雜流行病學(xué)模型的混沌研究*

狄根虎 許 勇徐 偉 顧仁財

(西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，西安 710129)

(2010年5月15日收到;2010年9月17日收到修改稿)

研究了一類周期變化的非線性復(fù)雜發(fā)病率的廣義流行病學(xué)模型(SIR(susceptible，infected，recovered)模型).通過一系列坐標(biāo)變換將原模型轉(zhuǎn)化為Hamilton系統(tǒng)，運(yùn)用Melnikov方法證明了該系統(tǒng)存在混沌運(yùn)動，給出了發(fā)生同宿分岔的條件，并用數(shù)值模擬驗(yàn)證了上述結(jié)果.

SIR(susceptible，infected，recovered)模型，混沌運(yùn)動，Melnikov方法，同宿分岔

PACS:05.45.-a，05.45.Pq，02.30.Qz

1.引 言

數(shù)學(xué)模型和方法已經(jīng)用于傳染病的傳播原理研究.Kermack和 McKendrick［1］在研究1665—1666年黑死病的流行規(guī)律時，構(gòu)造了著名的((SIR)倉室模型，現(xiàn)今最常用的是Anderson和May［2］所建立的傳染病數(shù)學(xué)模型，即假設(shè)個體對傳染病是免疫情形下的標(biāo)準(zhǔn)SIR模型:

其中，模型(1)假定人口是不變的，將人口分為三類:S(t)代表t時刻易感人群;I(t)代表t時刻的感染人群;R(t)代表t時刻的恢復(fù)且不再感染這類疾病的人群.參數(shù)μ表示人口出生率或死亡率(此模型假設(shè)出生率和死亡率相等)，γ表示感染人群的恢復(fù)率，而B(I，S，t)代表隨時間周期變化的發(fā)病率函數(shù).

在以往研究中，通常取B(I，S，t)=β(t)IS，其中傳染率β(t)或者是常數(shù)或者周期變化，如文獻(xiàn)［3—5］中，取β(t)=β0(1+β1sin(ωt))，并數(shù)值證明了當(dāng) β1充分大時，該模型存在混沌運(yùn)動.文獻(xiàn)［6—10］研究了復(fù)雜發(fā)病率函數(shù)，得到了具有不同于雙線性發(fā)病率函數(shù)的動力學(xué)行為，包括霍普夫分岔、鞍結(jié)分岔和同宿分岔.這些分岔意味著傳染病的爆發(fā)和滅絕，由參數(shù) β，p，q決定.文獻(xiàn)［11—17］基于復(fù)雜網(wǎng)絡(luò)對一類傳染病模型進(jìn)行研究.因此，具有復(fù)雜發(fā)病率函數(shù)的SIR模型(1)的研究是重要和必要的.然而，具有復(fù)雜發(fā)病率函數(shù)的 SIR模型(1)的混沌行為還研究較少，本文取復(fù)雜發(fā)病率函數(shù)為B(I，S，t)=β(t)IpSq，傳染率函數(shù)取為

其中，b1，b2，b3，Ω以及 p，q都為常數(shù)且 p>1，使得2pb2-(p-1)b21>0，ε>0是小參數(shù).這里我們所取的傳染率函數(shù)在 p=2，q=1時就退化到文獻(xiàn)［18］中所考慮的傳染率情形，它是在文中第二部分的模型變換時，為了將模型變換為Hamilton系統(tǒng)的標(biāo)準(zhǔn)形式，對傳染率函數(shù)進(jìn)行了適當(dāng)?shù)男拚紤]了參數(shù)p，q對傳染率的影響，這也是合理的.同時假定模型(1)中的恢復(fù)率函數(shù)為γ=μ(1+b1ε2)/(p-1).

本文將在復(fù)雜發(fā)病率函數(shù)(2)下研究模型(1)的混沌動力學(xué)行為.首先通過一系列坐標(biāo)變換將模型(1)化為Hamilton系統(tǒng)，然后應(yīng)用Melnikov方法計(jì)算了混沌產(chǎn)生的臨界值，得到同宿分岔的條件，并通過數(shù)值實(shí)驗(yàn)說明了解析結(jié)果的正確性.

2.模型變換

假定模型(1)中的總?cè)丝谑浅?shù)，且可表示為S +R+I=1.

將S=1-I-R代入模型(1)，模型簡化為

由于μ>0，定義新時間變量 t′=μt，消去 μ，方程(3)為

其中，β=μβ′，γ=μγ′.

令y=γ′I-R，即 I=γ′(-1)(y+R)，代入(4)式有

其中，(5)式中的導(dǎo)數(shù)是對新的時間變量t′而言.現(xiàn)在計(jì)算(5)式穩(wěn)態(tài)點(diǎn).當(dāng)y=0時，有

顯然R=0是方程(6)的一個解，即具有穩(wěn)態(tài)點(diǎn)(0，0)，此時稱(S，I，R)=(1，0，0)為無病平衡點(diǎn).進(jìn)一步，方程(6)簡化為

記 H=γ′/(1+γ′)，σ=(1+γ′)/β′γ′(1-p)，則Rp-1(1-R/H)q=σ，R∈［0，H］

定義函數(shù)f(R)=Rp-1(1-R/H)q，其中，R∈［0，H］.由極值原理，可知f(R)在［0，H］上的最大值為 σ*，其中，σ*=(p-1)p-1qqHp-1/(p+q-1)p+q-1.

當(dāng)σ>σ*時，系統(tǒng)(7)的正平衡點(diǎn)不存在，當(dāng)σ=σ*時，存在一個正平衡點(diǎn)，當(dāng) σ<σ*時，存在兩個正平衡點(diǎn).這里我們給出如下定理:

定理 如果β′<(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)，系統(tǒng)(7)除無病平衡點(diǎn)外不存在正平衡點(diǎn);如果 β′=(p+q-1)p+q-1(1+γ′)p/ ((p-1) q)，系統(tǒng)將除無病平衡點(diǎn)外還存在一個正平衡點(diǎn);如果β′>(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)，系統(tǒng)除無病平衡點(diǎn)外還存在兩個平衡點(diǎn).

由上述定理可得，系統(tǒng)(7)在β′=βs′n=(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)處發(fā)生鞍結(jié)分岔，在此處的平衡點(diǎn)為

系統(tǒng)(7)在穩(wěn)態(tài)(R，0)處的Jacobi矩陣為

其中，

注意到發(fā)生鞍結(jié)分岔時，J21=0.如果J22=0，則穩(wěn)態(tài)點(diǎn)就處在更加退化的情形，此時將穩(wěn)態(tài)點(diǎn)稱為Takens-Bogdanov點(diǎn)，由［19］知，在 Takens-Bogdanov點(diǎn)附近可以將方程化為近似 Hamilton系統(tǒng)的標(biāo)準(zhǔn)形式.

在鞍結(jié)分岔時，如果γ′=1/(p-1)，則矩陣A的跡消失，即J22=0.

取上述參數(shù)值，可得退化平衡點(diǎn)為

由于霍普夫分岔的曲線是在此退化平衡點(diǎn)終止的，曲線是由 J21<0，J22=0給定的.為擾動退化點(diǎn)，定義新變量和新參數(shù)

其中，x，g，b表示正的小擾動量.

將其帶入(5)式，忽略x，g，b的三次項(xiàng)及高階項(xiàng)，有

其中，

引進(jìn)小參數(shù)ε>0(與(2)式中的小參數(shù)ε相同)，令

將(13)式代入(12)式，忽略高次項(xiàng)得

定義新的坐標(biāo)變量和時間變量x=ε2u，y=ε3v，τ=εt′，系統(tǒng)(11)為

令z=u-((p-1)2/(p2(p+q-1)))b1，則系統(tǒng)(15)化為

定義變量k=(p(p+q-1)/(2q(1-p)))z，w=(p(p+q-1)/(2q(1-p)))v，方程(16)為

其中，c2=(p2(p+q-1)/(2(p-1)q))b2- (p(p+q-1)/(4q)，由前面假設(shè)可知 c2是嚴(yán)格大于0的常數(shù).以后取c>0.

當(dāng)b3=0，方程(17)是Takens-Bogdanov分岔的標(biāo)準(zhǔn)形式.假定 b1>0，b3=0以及ε是充分小的正常數(shù)，則方程有兩平衡點(diǎn)S±=(±c，0)，S+是一鞍點(diǎn)，S-有如下三種情形:

3.SIR模型中的混沌

本節(jié)運(yùn)用 Melnikov方法［20—24］研究 SIR模型經(jīng)過變換后的Hamilton系統(tǒng)(17)，證明系統(tǒng)的混沌存在性以及同宿分岔的條件.

如果ε=0，系統(tǒng)(17)的未擾形式為

其中，導(dǎo)數(shù)是對新時間變量τ而言，且c>0.

系統(tǒng)(18)的Hamilton函數(shù)和勢函數(shù)分別為

方程(18)有兩平衡點(diǎn)，(c，0)是鞍點(diǎn)，(-c，0)是中心.中心是被連續(xù)的一族周期軌道(kT(τ)，wT(τ))所圍繞，而且這些周期軌道是有界限的，是被同宿軌道(kh(τ)，wh(τ))所控制，這是因?yàn)?/p>

系統(tǒng)的Melnikov函數(shù)為

未擾系統(tǒng)的同宿解為［25］

代入(19)式并積分得

令b1=0，給出一簡單表達(dá)形式，其中

令b3=0，由(21)式給出系統(tǒng)無周期激勵的同宿分岔?xiàng)l件為b1=10(p-1)c/(7p)，即同宿分岔是由參數(shù)空間(b1，b2)上的曲線

確定的.

4.數(shù)值結(jié)果

本節(jié)運(yùn)用數(shù)值實(shí)驗(yàn)來驗(yàn)證 Melnikov方法所得到的解析結(jié)果，即證實(shí)系統(tǒng)存在混沌運(yùn)動.針對系統(tǒng)(17)，取如下參數(shù):

圖1給出了上述參數(shù)情形下的Poincaré截面，而圖1中點(diǎn)的雜亂分布，由 Poincaré理論知系統(tǒng)具有混沌運(yùn)動，與解析結(jié)果是一致的.圖2給出了系統(tǒng)(17)的相圖，我們發(fā)現(xiàn)系統(tǒng)也是混沌的.很明顯，數(shù)值實(shí)驗(yàn)與解析結(jié)果是一致的.

圖1 系統(tǒng)(17)的Poincaré截面 (a)bM=2.7919，(b)bM=6.6050

圖2 系統(tǒng)(17)的相圖 (a)bM=2.7919，(b)bM=6.6050

5.結(jié) 論

本文通過解析方法即 Melnikov方法證明復(fù)雜周期激勵的傳染病模型在 Takens-Bogdanov分岔點(diǎn)的很小范圍內(nèi)具有混沌運(yùn)動，給出出現(xiàn)混沌的臨界表達(dá)式(22)式.數(shù)值模擬驗(yàn)證了解析結(jié)果的正確性.從圖1，2可知，參數(shù)達(dá)到臨界值時，傳染病會處于爆發(fā)狀態(tài)，因此在預(yù)防傳染病的傳播時，要對參數(shù)進(jìn)行合理控制，使得傳染病在可控范圍內(nèi).

［1］Kermack W O，McKendrick A G 1927 Proc.Roy.Soc.A 115 700

［2］Anderson R M，May R M 1991 Infectious Diseases of Humans (Oxford:Oxford University Press)

［3］Aron J L，Schwartz I B 1984 Theor.Biol.110 665

［4］Kuznetsov，Yu A，Piccardi C 1994 Math.Biol.32 109

［5］Olsen L F，Shaffer W M 1990 Science 249 499

［6］Hethcote H W，Driessche P V D 1991 Math.Biol.29 271

［7］Liu W M，Hethcote H W，Levin S A 1987 Math.Biol.25 359

［8］Liu W M，Levin S A，Iwasa Y 1986 Math.Biol.23 187

［9］Yi H T 2009 MS dissertation(Hunan:Hunan university)(in Chinese)［尹海濤2009碩士學(xué)位論文(湖南:湖南大學(xué))］

［10］Li J 2005 MS dissertation(Nanjing:Nanjing university of science and technology)(in Chinese)［李 靜 2005碩士學(xué)位論文(南京:南京理工大學(xué))］

［11］Ni S J，Weng W G，F(xiàn)an W C 2009 Acta Phys.Sin.58 3707(in Chinese)［倪順江、翁文國、范維澄2009物理學(xué)報58 3707］

［12］Liu M X，Ruan J 2009 Chin.Phys.B 18 2115

［13］Liu M X，Ruan J 2009 Chin.Phys.B 18 5111

［14］Shi H J，Duan Z S，Chen G R，Li R 2009 Chin.Phys.B 18 3309

［15］Pei W D，Chen Z Q，Yuan Z Z 2008 Chin.Phys.B 17 0373

［16］Zhong L，Weng J Q 2000 Acta Phys.Sin.49 0626(in Chinese)［鐘 玲、翁甲強(qiáng)2000物理學(xué)報49 0626］

［17］Xu D，Li X，Wang X F 2007 Acta Phys.Sin.56 1313(in Chinese)［許 丹、李 翔、汪小帆2007物理學(xué)報56 1313］

［18］Glendinning P，Perry L P 1997 Math.Biol.34 359

［19］Chow S N，Li C Z，Wang D 1994 Normal forms and bifurcation of plannarvectorfields(Cambridge:Cambridge University Press)

［20］Yang X L，Xu W，Sun Z K 2006 Acta Phys.Sin.55 1678(in Chinese)［楊曉麗、徐 偉、孫中奎2006物理學(xué)報55 1678］

［21］Liu Z R 2002 The analytical method of chaos(Shanghai: Shanghai university press)p56(in Chinese)［劉增榮 2002混沌研究中的解析方法(上海:上海大學(xué)出版社)第56頁］

［22］Lei Y M，Xu W 2007 Acta Phys.Sin.56 5103(in Chinese)［雷佑銘、徐 偉2007物理學(xué)報56 5103］

［23］GuckenheimerJ， Holmes P 1983 NonlinearOscillations，Dynamical Systems and Bifurcations of Vectors Fields(New York: New York Springer)

［24］Yang X L，Xu W 2009 Acta Phys Sin.58 3722(in Chinese)［楊曉麗、徐 偉2009物理學(xué)報58 3722］

［25］Baesens N C，Nicolis G 1983 Z.Phys.B 52 345

PACS:05.45.-a，05.45.Pq，02.30.Qz

Chaos for a class of complex epidemiological models*

Di Gen-Hu Xu YongXu Wei Gu Ren-Cai
(Department of Applied Mathematics，Northwestern Polytechnical University，Xi'an 710129，China)

15 May 2010;revised manuscript

17 September 2010)

We study the well-known SIR(susceptible，infected，recoverd)model with nonlinear complex incidence rates. Firstly，a series of coordinate transformations are carried out to change the equations as the amenable Hamiltonian systems. Secondly the Melnikov′s method is used to establish the conditions of existence of chaotic motion and find the analytically critical values of homoclinic bifurcation.Good agreement can be found between numerical results and analytical results.

SIR(susceptible，infected，recoverd)model，chaotic motion，Melnikov's method，homoclinic bifurcation

*國家自然科學(xué)基金(批準(zhǔn)號:10972181，11002001，10872165，10932009)、西北工業(yè)大學(xué)基礎(chǔ)研究基金和翱翔之星資助的課題.

*Project supported by the National Natural Science Foundation of China(Grant Nos.10972181，11002001，10872165，10932009)，the NPU Foundation for Fundamental Research and the“Aoxiang Star”Plant of Northwestern Polytechnical University，China.