一類考慮復(fù)吸的時滯戒煙模型*

張子振，張偉詩

(安徽財經(jīng)大學(xué)管理科學(xué)與工程學(xué)院，安徽 蚌埠 233030)

1 問題的提出

生物數(shù)學(xué)領(lǐng)域?qū)W者認為吸煙可以通過社交活動接觸傳播，因此，近年來他們一直努力基于傳染病動力學(xué)原理，通過數(shù)學(xué)建模的方式向人類宣傳控?zé)熀徒錈煹谋匾?尤其是Garsow等[1]提出基本的戒煙模型以后，學(xué)者構(gòu)建了各類戒煙模型，如考慮臨時戒煙群體的戒煙模型[2]和考慮偶爾吸煙群體的戒煙模型[3].由于現(xiàn)實世界中多數(shù)問題具有一定的不確定性，因此學(xué)者進一步構(gòu)建了具不確定性的戒煙模型，如隨機戒煙模型[4]和具有非線性發(fā)生率的戒煙模型[5-6].Ullah等[7]在文獻[3-4]的基礎(chǔ)上，提出了如下考慮復(fù)吸的戒煙模型：

(1)

其中：P(t)，O(t)，S(t)，Qt(t)和Qp(t)分別為潛在吸煙者、偶爾吸煙者、重度吸煙者、暫時戒煙者和永久戒煙者在時刻t的數(shù)量；Λ，β，μ，α1，α2，γ，δ為模型(1)的參數(shù)，Λ為潛在吸煙者的輸入率，β為潛在吸煙者和重度吸煙者之間的有效接觸率，μ為所有個體的自然死亡率，α1為偶爾吸煙者變?yōu)橹囟任鼰熣叩谋壤?，?為暫時戒煙者復(fù)吸的比例，γ為重度吸煙者戒煙的比例，δ為永久戒煙者在所有戒煙者中的占比，1-δ為暫時戒煙者在所有戒煙者中的占比.

顯然，模型(1)假設(shè)偶爾吸煙者變?yōu)橹囟任鼰熣呤撬矔r的，這與現(xiàn)實不相符，因為偶爾吸煙者需要經(jīng)過一定的時間周期才會對尼古丁產(chǎn)生依賴而變?yōu)橹囟任鼰熣?因此，筆者在模型(1)中引入偶爾吸煙者變?yōu)橹囟任鼰熣咝枰?jīng)歷的時間周期時滯τ：

(2)

2 基本再生數(shù)和吸煙平衡點

于是

進而有

這里S*是方程

(3)

的正根.這里：

l0=μ((μ+γ)(μ+α1)μ-α1βΛ)；

l1=μ(μ+γ)(μ+α1)(α2+β)-α1α2βΛ-α2γμ(1-δ)(μ+α1)；

l2=α2β(μ+α1)(μ+γδ).

3 局部漸近穩(wěn)定性和Hopf分岔的存在性

模型(2)在吸煙平衡點E*(P*,O*,S*,Qt*,Qp*)的雅克比矩陣

其中：

m11=-(βS*+μ)；m13=-βP*；m21=βS*；m22=-μ；m23=βP*；

n22=-α1；m33=α2Qt*-(μ+γ)；m34=α2S*；n32=α1；

m43=γ(1-δ)-α2Qt*；m44=-(α2S*+μ)；m53=γδ；m55=-μ.

相應(yīng)的特征方程為

λ5+Θ4λ4+Θ3λ3+Θ2λ2+Θ1λ+Θ0+(Ξ4λ4+Ξ3λ3+Ξ2λ2+Ξ1λ+Ξ0)e-λτ=0.

(4)

其中：

Θ0=-m11m22m55(m33m44+m34m43)；

Θ1=m11m22(m33m44+m33m55+m44m55)+m33m44m55(m11+m22)+

m34m43(m11m22+m11m55+m22m55)；

Θ2=-(m34m43(m11+m22+m55)+m11m22(m33+m44+m55))-

((m11+m22)(m33m44+m33m55+m44m55)+m33m44m55)；

Θ3=(m11+m22)(m33+m44+m55)+m11m22+m33m44+m33m55+m44m55+m34m43；

Θ4=-(m11+m22+m33+m44+m55)；

Ξ0=m11m44m55(m23n32-m33n22)-m55(m13m21m44n32+m11m34m43n22)；

Ξ1=m34m43n22(m11+m55)+(m44+m55)(m11m33n22+m13m21n32)+

m44m55n22(m11+m33)-m23n32(m11m44+m11m55+m44m55)；

Ξ2=m23n32(m11+m44+m55)-m13m21n32-(m11m33+m34m43)n22-

n22(m44m55+(m11+m33)(m44+m55))；

Ξ3=n22(m11+m33+m44+m55)-m23n32；

Ξ4=-n22.

當τ=0時，方程(4)變?yōu)?/p>

λ5+Λ4λ4+Δ3λ3+Δ2λ2+Δ1λ+Δ0=0.

其中：

Δ0=Θ0+Ξ0；Δ1=Θ1+Ξ1；Δ2=Θ2+Ξ2；Δ3=Θ3+Ξ3；Δ4=Θ4+Ξ4.

根據(jù)Hurwitz穩(wěn)定性判據(jù)，當如下方程成立且τ=0時，模型(2)是局部漸近穩(wěn)定的：

L1=Δ4>0，

當τ>0時，假設(shè)λ=i?(?>0)是方程(4)的根，那么

(5)

進而得到關(guān)于?的代數(shù)方程

?10+Κ4?8+Κ3?6+Κ2?4+Κ1?2+Κ0=0.

(6)

其中：

令?2=χ，則方程(6)變?yōu)?/p>

χ5+Κ4χ4+Κ3χ3+Κ2χ2+Κ1χ+Κ0=0.

(7)

其中

對方程(4)左右兩邊同時求λ關(guān)于τ的導(dǎo)數(shù)，得到

于是

根據(jù)文獻[10]中關(guān)于動力系統(tǒng)產(chǎn)生Hopf分岔的定理，可得如下結(jié)果：

定理1如果R0>1，那么當τ∈(0,τ0)時，模型(2)局部漸近穩(wěn)定；當τ=τ0時，模型(2)在吸煙平衡點E*(P*,O*,S*,Qt*,Qp*)處產(chǎn)生Hopf分岔，并產(chǎn)生一簇分岔周期解.

4 仿真示例

參考文獻[7]中的參數(shù)取值，同時考慮到模型(2)產(chǎn)生Hopf分岔的充分條件，定義Λ=1，β=0.14，μ=0.001，α1=0.05，α2=0.002 5，γ=0.8，δ=0.48，代入模型(2)可得

(8)

于是R0=1.713 5>1，方程(3)變?yōu)?/p>

(9)

求解方程(9)得到唯一正根S*=2.173 1，由此可得模型(2)的唯一吸煙平衡點E*(3.276 2,19.543 6,2.173 1,140.532 4,834.470 4)，進而得到?0=1.605 6，τ0=25.270 9.當τ=23.379 2∈(0,τ0)時，模型(8)是局部漸近穩(wěn)定的，此時模型(8)的狀態(tài)軌跡和相圖分別如圖1和圖2所示；當τ=29.903 9>τ0時，模型(8)失去穩(wěn)定性，在E*(3.276 2,19.543 6,2.173 1,140.532 4,834.470 4)附近產(chǎn)生Hopf分岔和一簇分岔周期解，此時模型(8)的狀態(tài)軌跡和相圖分別如圖3和圖4所示.

圖1 當τ=23.379 2時，模型(8)的狀態(tài)軌跡Fig. 1 State Trajectories of Model (8) When τ=23.379 2

圖2 當τ=23.379 2時，模型(8)的相圖Fig. 2 Phase Plots of Model (8) When τ=23.379 2

圖3 當τ=29.903 9時，模型(8)的狀態(tài)軌跡Fig. 3 State Trajectories of Model (8) When τ=29.903 9

圖4 當τ=29.903 9時，模型(8)的相圖Fig. 4 Phase Plots of Model (8) When τ=29.903 9

5 結(jié)語

在文獻[7]的基礎(chǔ)上引入偶爾吸煙者變?yōu)橹囟任鼰熣咝枰?jīng)歷的時間周期時滯，研究了一類考慮復(fù)吸的時滯戒煙模型.首先計算出模型的基本再生數(shù)和吸煙平衡點，然后以偶爾吸煙者變?yōu)橹囟任鼰熣咝枰?jīng)歷的時間周期時滯為分岔參數(shù)，推導(dǎo)出模型的局部漸近穩(wěn)定和產(chǎn)生Hopf分岔的充分條件，進而計算出模型產(chǎn)生Hopf分岔時滯的臨界點.研究結(jié)果表明，時滯取值足夠小時，模型處于理想的穩(wěn)定狀態(tài)，此時有利于控制吸煙的流行傳播；時滯取值越過臨界點τ0時，模型失去穩(wěn)定，產(chǎn)生Hopf分岔，此時不利于控制吸煙的流行傳播.值得注意的是，筆者只研究了模型(2)的局部穩(wěn)定性和Hopf分岔的存在性，并未對Hopf分岔的方向及穩(wěn)定性等展開討論.下一步，筆者將深入研究Hopf分岔的性質(zhì).