譚偉 劉茂省

DOI：10.16783/j.cnki.nwnuz.2024.01.006

收稿日期：20221108;修改稿收到日期：20230522

基金項目：國家自然科學(xué)基金資助項目（12071445，12001501）

作者簡介：譚偉（1996—），男，河北張家口人，碩士研究生.主要研究方向為傳染病動力學(xué).

Email：1848231722@qq.com

*通信聯(lián)系人，男，教授，博士.主要研究方向為傳染病動力系統(tǒng)及復(fù)雜網(wǎng)絡(luò).Email：liumaoxing@126.com

摘要：研究一個具有飽和發(fā)生率和疫苗接種率的隨機離散SIR傳染病模型的穩(wěn)定性.首先，引入一個具有飽和發(fā)生率和疫苗接種率的確定性SIR模型，考慮到隨機噪聲對疾病傳播有很大影響，應(yīng)用非標(biāo)準(zhǔn)有限差分（NSFD）方法將模型離散化，最終得到一個隨機離散的SIR傳染病模型.這種離散方法是對系統(tǒng)的右側(cè)進行局部離散，得出離散模型，然后系統(tǒng)左側(cè)用廣義前向差分法對一階導(dǎo)數(shù)進行逼近，并且要選取恰當(dāng)?shù)姆帜负瘮?shù).其次，應(yīng)用Lyapunov函數(shù)和矩陣方法給出了系統(tǒng)平衡解穩(wěn)定的充分條件，并且得到了非線性差分方程概率穩(wěn)定的充分條件和線性差分方程漸近均方穩(wěn)定的充分條件.最后，通過數(shù)值仿真對理論分析結(jié)論進行驗證.

關(guān)鍵詞：飽和發(fā)生率；隨機離散模型；非標(biāo)準(zhǔn)有限差分方法；Lyapunov函數(shù)；漸近均方穩(wěn)定

中圖分類號：O 175.7??? 文獻標(biāo)志碼：A??? 文章編號：1001-988Ⅹ（2024）01-0030-09

Nonstandard numerical discretization of

a stochastic SIR epidemic model

TAN Wei1，LIU Mao-xing1，2

（1.College of Mathematics，North University of China，Taiyuan 030051，Shanxi，China;

（2.College of Science，Beijing University of Civil Engineering and Architecture，Beijing 102616，China）

Abstract：This paper studies the stability of a stochastic discrete SIR epidemic model with saturated incidence and vaccination rate.A deterministic SIR model with saturated incidence and vaccination rate is introduced.Considering the great influence of stochastic noise on disease transmission，the model is discretized by nonstandard finite difference（NSFD） method，and finally a stochastic discrete SIR epidemic model is obtained.This discretization method is to locally discretize the right side of the system to obtain the discrete model，and then use the generalized forward difference method to approximate the first derivative on the left side of the system，and select the appropriate denominator function.The sufficient conditions for the stability of the equilibrium solution of the system are given by using Lyapunov function method and matrix method，and the sufficient conditions for the probabilistic stability of nonlinear difference equations and the sufficient conditions for the asymptotic mean square stability of linear difference equations are also proposed.Finally，the conclusion is verified by numerical simulation.

Key words：saturation incidence;stochastic discrete model;nonstandard finite difference method;Lyapunov function;asymptotic mean square stability

0? 引言

近幾十年來，數(shù)學(xué)模型被認為是理解傳染病傳播機制和評估其控制策略的關(guān)鍵工具［1-3］.在這方面，Yusuf等［4］引入并分析了一個具有疫苗接種和治療的SIR流行模型，其中假設(shè)發(fā)病率是雙線性的.但經(jīng)過長期的驗證發(fā)現(xiàn)，雙線性發(fā)生率對于較大的種群規(guī)模來說并不符合實際，為了應(yīng)對易感個體產(chǎn)生的心理或抑制效應(yīng)，Capasso等［5］提出了飽和發(fā)病率βSI/（1+αI），其中β是疾病的傳染率，α是衡量由擁擠效應(yīng)或社會行為變化引起的抑制效應(yīng)的飽和因子.當(dāng)前，飽和發(fā)病率已廣泛應(yīng)用于確定性傳染病模型［6-9］和隨機傳染病模型［10-12］.本文研究具有飽和發(fā)生率和疫苗接種率的SIR確定性傳染病模型［13］：

dSdt=Λ-βSI1+αI-（μ+ν）S，

dIdt=βSI1+αI-（μ+γ+δ）I，

dRdt=νS+γI-μR.（1）

其中S≡S（t）為t時刻易感人群數(shù)量，I≡I（t）為t時刻已感染人群數(shù)量，R≡R（t）為t時刻染病后免疫人群的數(shù)量，且S（0）≥0，I（0）≥0，R（0）≥0；Λ為從外部進入種群的凈遷移率，μ為人口的自然死亡率，ν為接種疫苗的人口比率，γ為患病種群恢復(fù)的概率，δ為因為感染疾病而死亡的人口概率，且所有參數(shù)都是正的.

注意到系統(tǒng)（1）的第一和第二個方程中不包含R，系統(tǒng)（1）的第三個式子不影響其動力學(xué)行為，因此，我們只需研究以下系統(tǒng)：

dSdt=Λ-βSI1+αI-（μ+ν）S，

dIdt=βSI1+αI-（μ+γ+δ）I.（2）

通過簡單計算可知，系統(tǒng)（2）中存在一個無病平衡點E0和一個地方病平衡點E+，并且無論在什么條件下，無病平衡點

E0=（S0，I0）=Λμ+ν，0

都存在.當(dāng)R0>1時地方病平衡點E+=（S+，I+）存在，且

S+=αΛ+（μ+γ+δ）α（μ+ν）+β，

I+=（μ+ν）（R0-1）α（μ+ν）+β，

基本再生數(shù)

R0=Λβ（μ+ν）（μ+γ+β）.

由于白噪聲對疾病的傳播有很大的影響，所以與確定性模型相比，隨機模型更加符合實際.本文在模型（2）的基礎(chǔ)上加入隨機噪聲，得到

dSdt=Λ-βSI1+αI-（μ+ν）S+

σ1（S-S*）dB1（t），

dIdt=βSI1+αI-（μ+γ+δ）I+

σ2（I-I*）dB2（t），（3）

其中Bi（t）（i=1，2）為白噪聲，也就是常說的布朗運動，σ2i（i=1，2）為白噪聲Bi（t）的強度［14］.假設(shè)隨機擾動強度與狀態(tài)S（t）和I（t）偏離平衡點的偏差成正比，即當(dāng)系統(tǒng)狀態(tài)偏離平衡點的偏差增大時，隨機擾動強度也隨之增加［15］.

傳染病模型通常是一個非線性微分方程組，因此很難得到其解析解，通常需要數(shù)值近似來得到可靠的結(jié)果.然而，許多標(biāo)準(zhǔn)的數(shù)值方法，如Euler方法、Runge-Kutta方法可能無法保持連續(xù)模型的動力學(xué)特性，例如收斂到錯誤的平衡點或數(shù)值不穩(wěn)定［16-18］.為克服這些方法的缺點，Mickens［19］提出了非標(biāo)準(zhǔn)有限差分（NSFD）方法.

如果以下兩個條件之一被滿足，那么一階微分方程組的這種數(shù)值格式即為NSFD格式：

（i）用廣義前向差分近似系統(tǒng)中的一階導(dǎo)數(shù)

dundt≈u（n+1）-u（n）ψ，

其中un=u（tn），分母函數(shù)ψ=ψ（h）滿足ψ（h）=h+O（h2），h為時間步長.

（ii）非線性項是非局部近似的，例如

u2（tn）≈unun+1.

然而，注意到NSFD格式對一個微分方程來說并不是唯一的，所以為了確保NSFD格式與原始連續(xù)模型具有相同的動力學(xué)性態(tài)，我們需要對NSFD格式進行動力學(xué)特性的分析.

假定條件（i）成立，我們根據(jù)系統(tǒng)（3）提出以下NSFD格式：

S（n+1）-S（n）ψ=

Λ-βS（n）I（n）1+αI（n）-（μ+ν）S（n）+

σ11ψ（S（n）-S*）ω1（n+1），

I（n+1）-I（n）ψ=

βS（n）I（n）1+αI（n）-（μ+γ+δ）I（n）+

σ21ψ（I（n）-I*）ω2（n+1），（4）

其中分母函數(shù)ψ是與h相關(guān)的函數(shù)，初始條件為S（0）=φ1（0），I（0）=φ2（0）.

系統(tǒng)（4）是對系統(tǒng)（3）進行廣義前向差分和右側(cè)局部近似得到的，通過一些簡單的代數(shù)計算，我們可以很容易地驗證，具有NSFD格式的模型（4）與初始連續(xù)模型（2）具有完全相同的平衡點.數(shù)學(xué)期望我們通常由E來表示，ωi（n+1）（i=1，2）（n∈Z）是一個相互獨立且服從標(biāo)準(zhǔn)正態(tài)分布N（0，1）的隨機序列，并且對于ωi（n+1）（i=1，2，n∈Z），滿足

E［ωi（n）］=0， E［ω2i（n）］=1，

E［ωi（n）］≠E［ωj（n）］， i≠j.（5）

通過直接計算，容易驗證分母函數(shù)［20］

ψ（h）=1-e-μhμ.（6）

對離散NSFD系統(tǒng)（4）可重新排列，得到其顯式形式：

S（n+1）=S（n）+

Λ-βS（n）I（n）1+αI（n）-（μ+ν）S（n）ψ+

σ1ψ（S（n）-S*）ω1（n+1），

I（n+1）=I（n）+

βS（n）I（n）1+αI（n）-（μ+γ+δ）I（n）ψ+

σ2ψ（I（n）-I*）ω2（n+1）.（7）

近幾十年來，具有噪聲擾動的確定性模型的隨機流行病模型［21-23］已經(jīng)被廣泛研究，并且越來越多的學(xué)者將差分方程模型［24，25］運用到傳染病學(xué)的建模中.然而，多數(shù)學(xué)者都是將隨機模型和差分格式分別建模來研究的［21］，但在傳染病學(xué)中對差分方程所描述的隨機離散模型的穩(wěn)定性的研究很少，所以本文主要對隨機離散的傳染病模型進行研究.

1? 預(yù)備知識

對系統(tǒng)（7）的地方病平衡點進行平移變換，令u（n）=S（n）-S+，m（n）=I（n）-I+可得

u（n+1）=u（n）+

Λ-β（u（n）+S+）（m（n）+I+）1+α（m（n）+I+）-

（u+v）（u（n）+S+）ψ+

σ1ψu（n）ω1（n+1），??????? （8）

m（n+1）=m（n）+

β（u（n）+S+）（m（n）+I+）1+α（m（n）+I+）-

（u+γ+δ）（m（n）+I+）ψ+

σ2ψm（n）ω2（n+1）.

系統(tǒng)（8）的零解等價于系統(tǒng)（7）的正解E+=（S+，I+）.

在平衡點E+=（S+，I+）處對系統(tǒng)（7）進行線性化，得到線性近似系統(tǒng)

x（n+1）=1-βI+ψ1+αI+-（μ+v）ψx（n）-

βS+ψ（1+αI+）2y（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=βI+ψ1+αI+x（n）+????? （9）

1+βS+ψ（1+αI+）2-

（μ+γ+δ）ψy（n）+

σ2ψy（n）ω2（n+1）.

類似地進行平移變換.對無病平衡點，令u（n）=S（n）-Λ/（μ+v），v（n）=I（n），則

u（n+1）=u（n）+

-βu（n）+Λμ+vm（n）1+αm（n）-

（u+v）u（n）ψ+

σ1ψu（n）ω1（n+1），??????? （10）

m（n+1）=m（n）+

βu（n）+Λμ+vm（n）1+αm（n）-

（u+γ+δ）m（n）ψ+σ2ψm（n）ω2（n+1）.

在無病平衡點處對系統(tǒng)（7）線性化，得到相應(yīng)的線性近似系統(tǒng)為

x（n+1）=（1-（μ+v）ψ）x（n）-Λβψμy（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=1+Λβψμ-（μ+γ+δ）ψy（n）+

σ2ψy（n）ω2（n+1）.（11）

令ρ（n）=（u（n），m（n））T，z（n）=（x（n），y（n））T，φ（n）=（φ1（n），φ2（n））T，T代表轉(zhuǎn)置.

下面引入了文獻［26］的一些定義和定理來研究系統(tǒng)的動態(tài)行為.

定義1（文獻［26］定義7.1）? 對初始函數(shù)S（0）=φ1（0），I（0）=φ2（0），如果對ε1>0和ε2>0，存在一個δ>0使得系統(tǒng)（8）或（10）的解ρ（n）=ρ（n，φ）滿足不等式Psupn∈Zρ（n）>ε1<ε2和P{φ<δ}=1，那么稱系統(tǒng)（8）或（10）是依概率穩(wěn)定的.

定義2（文獻［26］定義1.2）? 如果對每個ε>0，存在δ（ε）>0，使得E［z（n）2］<ε，n∈Z，并且E［φ（n）2］<δ對任意初始條件都被滿足，那么線性系統(tǒng)（9）或（11）的零解稱為均方穩(wěn)定的；如果系統(tǒng)（9）或（11）是均方穩(wěn)定的且其解滿足limn→∞E［z（n）2］=0，則稱它們是漸近均方穩(wěn)定的.

對非負函數(shù)Vi=V（i，z（0），z（1），…，z（i）），i∈Z，定義

ΔVi=V（i+1，z（0），z（1），…，z（i+1））-

V（i，z（0），z（1），…，z（i））.

定理A（文獻［26］定理1.1）? 對于線性系統(tǒng)（9）或（11），若存在一個滿足條件E［V（0，φ）］≤c1φ2和E［ΔVi］≤-c2E［z（i）2］（i∈Z）的非負函數(shù)Vi=V（i，z（0），z（1），…，z（i）），其中c1和c2是正常數(shù)，那么系統(tǒng)（9）或（11）的零解是漸近均方穩(wěn)定的.

考慮以下系統(tǒng)：

x（n+1）=a11x（n）+a12y（n）+

σ1ψx（n）ω1（n+1），

y（n+1）=a21x（n）+a22y（n）+

σ2ψx（n）ω2（n+1），（12）

其中σ1，σ2是隨機項中的常數(shù)，ψ表示分母函數(shù)，ωi（n+1）（i=1，2）是相互獨立的隨機序列且滿足條件（5）.可以看出，系統(tǒng)（12）為系統(tǒng)（9）和（11）的一般形式，因而，要探究系統(tǒng)（9）和（11）解的穩(wěn)定性，只需研究系統(tǒng)（12）解的穩(wěn)定性.

令

A=a11a12a21a22，

D=d11d12d12d22，

P=p11p12p12p22，（13）

其中D和P是對稱矩陣.對于實對稱矩陣P和Q，如果P-Q是一個正定矩陣，那么P>Q.

定理B（文獻［26］定理5.1）? 假設(shè)存在正定矩陣P，使得形為（13）的半正定矩陣D滿足矩陣方程ATDA-D=-P，同時P要滿足條件

P>d11σ21ψ00d22σ22ψ，（14）

那么系統(tǒng)（12）的零解是漸近均方穩(wěn)定的.

證明? 由（13）式，系統(tǒng)（12）可表示為

z（n+1）=（A+θ（ω（n+1）））z（n），

其中

z（n）=x（n）y（n），

θ（ω（n+1））=

σ1ψω1（n+1）00σ2ψω2（n+1）.

考慮Lyapunov函數(shù)V（n）=zT（n）Dz（n），則

ΔV（n）=V（n+1）-V（n）=

zT（n+1）Dz（n+1）-zT（n）Dz（n）.

然后計算ΔV的期望，得到

E［ΔV（n）］=E［zT（n+1）Dz（n+1）-

zT（n）Dz（n）］=

E［zT（n）］·［（A+θ（ω（n+1）））T×

D（A+θ（ω（n+1）））-D］z（n）=

E［zT（n）］［ATDA-D+θT（ω（n+1））×

Dθ（ω（n+1））］z（n）=

EzT（n）［-P+θT（ω（n+1））×

Dθ（ω（n+1））］z（n）.

由關(guān)系式（5），可以求得

E［θT（ωω（+1）Dθ）Dθ（+1）］=

Eσ1ψw1（n+1）00σ2ψw2（n+1）×

d11d12d12d22×

σ1ψw1（n+1）00σ2ψw2（n+1）=

d11σ21ψ00d22σ22ψ，

然后根據(jù)（14）式，可以得到

E［ΔV（n）］=EzT（n）-P+

d11σ21ψ00d22σ22ψz（n）≤

-cE［z（n）2］，

其中c是正數(shù).依據(jù)定理A可知系統(tǒng)（12）的零解是漸近均方穩(wěn)定的.? 】

2? 地方病平衡點的穩(wěn)定性

因為系統(tǒng)（9）是系統(tǒng)（12）的特殊形式，可以推導(dǎo)出系統(tǒng)（9）的系數(shù)矩陣

A=A11A12A21A22，（15）

其中

A11=1-βI+ψ1+αI+-（μ+v）ψ，

A12=-βS+ψ（1+αI+）2，

A21=βI+ψ1+αI+，

A22=1+βS+ψ（1+αI+）2-（μ+γ+δ）ψ.

從矩陣方程ATDA-D=-P可解出

-p11=（a211-1）d11+2a11a21d12+a221d22，

-p12=a11a12d11+（a21a12+a11a22-1）d12+

a21a22d22，

-p22=a212d11+2a12a22d12+（a222-1）d22.（16）

將矩陣A的元素帶入（16）式可得

p11=-1-βI+ψ1+αI+-（μ+v）ψ2-1d11-

2βI+h1-βI+ψ1+αI+-（μ+v）ψ1+αI+ψd12-

βI+ψ1+αI+2d22，

p12=1-βI+ψ1+αI+-（μ+v）ψ×

βS+ψ（1+αI+）2d11+

βI+ψ1+αI+·βS+ψ（1+αI+）2d12-

-1-βI+ψ1+αI+-（μ+v）ψ×

1-βS+ψ（1+αI+）2-（μ+γ+δ）ψ-1d12-

βI+ψ1+αI+1+βS+ψ（1+αI+）2-

（μ+γ+δ）ψd22，

p22=-βS+ψ（1+αI+）22d11+2βS+ψ（1+αI+）2×

1+βS+ψ（1+αI+）2-（μ+γ+δ）ψd12-

1+βS+ψ（1+αI+）2-（μ+γ+δ）ψ2-1d22.（17）

再由（14）式可導(dǎo)出

p11-d11σ21ψ>0， p11p22-p212>0，

（p11-d11σ21ψ）（p22-d22σ22ψ）-p212>0.（18）

根據(jù)以上的分析我們得出以下推論.

推論1? 若存在一個正定矩陣P，使得條件（18）對于矩陣方程（17）的解（d11，d12，d22）被滿足，并且矩陣D是半正定的，那么系統(tǒng)（9）的零解是漸近均方穩(wěn)定的.因此，系統(tǒng)（8）的零解是依概率穩(wěn)定的，這等同于系統(tǒng)（7）的正平衡點Ee是依概率穩(wěn)定的.

當(dāng)R0>1時，系統(tǒng)（7）存在地方病平衡點E+=（S+，I+）.對于數(shù)值仿真，假設(shè)初值

S（0）=5， I（0）=2.（19）

其它參數(shù)設(shè)置為

Λ=4，α=0.01，β=0.5，μ=0.5，

γ=0.1， ν=0.2，δ=0.1，h=0.1.（20）

通過這些參數(shù)能夠求出分母函數(shù)為ψ=0.0976，基本再生數(shù)R0=4.0816>1，從而地方病平衡點存在，且E+=（1.4597，4.2629）.通過對（15）求解得到矩陣

A=0.7322-0.06550.19950.9972；

假設(shè)P為單位矩陣，即

P=1001，

求解（17）式得到半正定矩陣

D=7.86907.38447.384412.3981.

由推論1可知，若σ1，σ2滿足條件

σ1<1d11ψ=1.1410，

σ2<1d22ψ=0.9091，（21）

那么，系統(tǒng)（9）的零解是漸近均方穩(wěn)定的.因此，系統(tǒng)（8）的零解是依概率穩(wěn)定的，這等同于系統(tǒng)（7）的地方病平衡點E+是依概率穩(wěn)定的.

下面借助數(shù)值模擬證明之前的分析.（19）和（20）式給出了初值和一些參數(shù)值.圖1顯示了噪聲擾動強度σ1=σ2=0的情況，系統(tǒng)（9）和線性系統(tǒng)（7）的解曲線分別由圖1（a）和圖1（b）表示.從中可以看出，線性系統(tǒng)（9）的解曲線最終收斂到0，系統(tǒng)（7）的解曲線最終收斂于地方病平衡點E+=

（1.4597，4.2629）.因此，系統(tǒng)（9）的零解是漸近均方穩(wěn)定的，系統(tǒng)（7）的地方病平衡點E+是依概率穩(wěn)定的.

對具有隨機擾動的系統(tǒng)（7）的解進行模擬.（19）和（20）式已給出初值和其它的參數(shù)值.假設(shè)σ1，σ2的值滿足條件（21），我們?nèi)ˇ?=σ2=0.8，得到具有噪聲擾動的系統(tǒng)（7）的解軌跡，見圖2.

從圖2可以看出，隨機噪聲并不影響曲線的最終走向，系統(tǒng)（7）的解軌跡最終還是收斂到地方病平衡點E+，這就表明系統(tǒng)（7）在正平衡點處是依概率穩(wěn)定的.

圖2? σ1=σ2=0.8時系統(tǒng)（7）的解軌跡

Fig 2Solution trajectories of system （7） when σ1=σ2=0.8

3? 無病平衡點的穩(wěn)定性

由系統(tǒng)（11）可以寫出它的系數(shù)矩陣

A=1-（μ+v）ψ-Λβψμ

01+Λβψμ-（μ+γ+δ）ψ.（22）

同理，由方程（16）解得

p11=-（（1-（μ+ν）ψ）2-1）d11，

p12=（1-（μ+ν）ψ）·Λβψμd11-

（1-（μ+ν）ψ）1+Λβψμ-

（μ+γ+δ）ψ-1d12，

p22=-Λβψμ2d11+

2Λβψμ1+Λβψμ-（μ+γ+δ）ψd12-

1+Λβψμ-（μ+γ+δ）ψ2-1d22.（23）

再由（18）式，我們同樣能推導(dǎo)出以下結(jié)論.

推論2? 對于（22）式中的矩陣A，如果存在一個正定矩陣P，使條件（18）被方程（23）的解（d11，d12，d22）所滿足，且矩陣D為半正定的，那么稱系統(tǒng)（11）的零解是漸近均方穩(wěn)定的.因此，系統(tǒng)（10）的零解是依概率穩(wěn)定的，這等同于是依概率穩(wěn)定的系統(tǒng)（7）的邊界平衡點E0.

當(dāng)R0<1時，系統(tǒng)（7）不存在地方病平衡點，只有一個無病平衡點.取參數(shù)值

Λ=4，μ=0.01，γ=0.9，β=0.01，

δ=0.2，ν=0.1，h=0.1.（24）

此時R0=0.167<1，所以只有一個無病平衡點存在，且E0=（20，0）.借助（6）式求得分母函數(shù)ψ=0.0995，給定初值為S（0）=10，I（0）=5.由以上參數(shù)能夠求出（22）式中的矩陣

A=0.9801-0.039800.9204.

仍取P為單位矩陣，借助（23）式求得半正定矩陣

D=25.3807-10.1124-10.112411.4208.

并由推論1可知，若σ1，σ2滿足條件：

σ1<1d11ψ=0.6293，

σ2<1d22ψ=0.9381，（25）

那么系統(tǒng)（11）的零解是漸近均方穩(wěn)定的，因此，是依概率穩(wěn)定的系統(tǒng)（10）的零解，這等同于是依概率穩(wěn)定的系統(tǒng)（7）的無病平衡點.

根據(jù)（24）中的參數(shù)值和給定的初值，當(dāng)隨機擾動強度為0時，我們得到系統(tǒng)（7）和（11）的解曲線，見圖3.從中可以看出，線性系統(tǒng)（11）和系統(tǒng)（7）的解曲線最終分別收斂于它們的平衡點.因此，通過圖像，我們也可以得到系統(tǒng)（11）的零解是漸近均方穩(wěn)定的，而系統(tǒng)（7）的無病平衡點E0是依概率穩(wěn)定的.

當(dāng)隨機噪聲強度不為0，且σ1和σ2的值滿足條件（25）時，我們?nèi)ˇ?=σ2=0.5可以得到具有隨機噪聲擾動的系統(tǒng)（7）的解軌跡，見圖4.從圖4可以直觀地看到，白噪聲只在開始的一段時間內(nèi)對疾病傳播具有較小的影響，但系統(tǒng)（7）的最終走向仍舊收斂到邊界平衡點E0=（20，0），這也表明系統(tǒng)（7）的邊界平衡點在噪聲擾動下依然是依概率穩(wěn)定的.

when σ1=σ2=0

when σ1=σ2=0.5

當(dāng)R0>1時，考慮系統(tǒng)（7）中地方病平衡點存在時無病平衡點的穩(wěn)定性.仍?。?9）和（20）式中的初值和參數(shù)，此時求得（22）式中的矩陣

A=0.9317-0.390401.3221.

P仍為單位矩陣，推導(dǎo)出矩陣

D=7.511511.786911.786913.4013，

通過計算可知，矩陣D不是半正定的，所以根據(jù)推論1可知，系統(tǒng)（7）的無病平衡點不穩(wěn)定.

當(dāng)隨機噪聲強度為0時，圖5中代表染病者的曲線I趨向于無窮，所以此時系統(tǒng)（11）的零解是不穩(wěn)定的.當(dāng)加入一個較小的噪聲強度后，系統(tǒng)的解軌跡是非常不規(guī)則的（圖6）.因此，當(dāng)一個正平衡點存在時，無論是否存在隨機擾動，系統(tǒng)（7）的邊界平衡點都是不穩(wěn)定的.

4? 結(jié)論

本文討論了一個具有飽和發(fā)生率和疫苗接種率的連續(xù)SIR傳染病模型，鑒于疾病在傳播過程中會受到各種外部因素的影響，例如隨機白噪聲的影響，我們在系統(tǒng)中加入隨機擾動項.由于離散模型的動力學(xué)性態(tài)比連續(xù)模型表現(xiàn)的更加豐富和準(zhǔn)確，我們還構(gòu)造了保持連續(xù)模型基本性質(zhì)的非標(biāo)準(zhǔn)有限差分格式，并選擇了一個合適的分母函數(shù).隨后用Lyapunov函數(shù)方法探究了所構(gòu)建的隨機差分模型的動力學(xué)性質(zhì)，建立了模型在平衡點處穩(wěn)定的充分條件.數(shù)值模擬表明，當(dāng)所選參數(shù)和噪聲強度滿足給定條件時，適當(dāng)?shù)脑肼晱姸炔粫绊懴到y(tǒng)在平衡點處的穩(wěn)定性；并且當(dāng)正平衡點存在時，無論隨機擾動強度是多少，邊界平衡點的穩(wěn)定性都不會被改變，在這種情況下，邊界平衡點是不穩(wěn)定的.

參考文獻：

［1］? CASTILLOCHAVEZ C.Mathematical Models in Population Biology and Epidemiology［M］.New York：Springer，2012.

［2］? BRAUER F.Mathematical epidemiology：past，present，and future［J］.Infectious Disease Modelling，2017（2）：113.

［3］? REZAPOUR S，MOHAMMADI H，SAMEI M E.SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order［J］.Advances in Difference Equations，2020，2020（1）：490.

［4］? YUSUF T T，BENYAH F.Optimal control of vaccination and treatment for an SIR epidemiological model［J］.World Journal of Modelling & Simulation，2012，8（3）：194.

［5］? CAPASSO V，SERIO G.A generalization of the Kermack-Mc Kendrick deterministic epidemic model［J］.Mathematical Biosciences，1978，42（1/2）：43.

［6］? XU R，MA Z，WANG Z.Global stability of a delayed SIRS epidemic model with saturation incidence and temporary immunity［J］.Computers & Mathematics with Applications，2010，59（9）：3211.

［7］? JANA S，NANDI S K，KAR T K.Complex dynamics of an SIR epidemic model with saturated incidence rate and treatment［J］.Acta Biotheoretica，2016，64：65.

［8］? MATHUR K S，NARAYAN P.Dynamics of an epidemic model with vaccination and saturated incidence rate［J］.Nature Public Health Emergency Collection，2018，4（5）：1.

［9］? LIU J.Bifurcation analysis for a delayed SEIR epidemic model with saturated incidence and saturated treatment function［J］.Journal of Biological Dynamics，2019，13（1）：461.

［10］? RAJASEKAR S P，PITCHAIMANI M，ZHU Q.Dynamic threshold probe of stochastic SIR model with saturated incidence rate and saturated treatment function［J］.Physica A：Statistical Mechanics and its Applications，2019，535：122300.

［11］? RAJASEKAR S P，PITCHAIMANI M，ZHU Q.Progressive dynamics of a stochastic epidemic model with logistic growth and saturated treatment［J］.Physica A：Statistical Mechanics and its Applications，2019，536：122649.

［12］? LI D，ZHAO Y，SONG S.Dynamic analysis of stochastic virus infection model with delay effect［J］.Physica A：Statistical Mechanics and its Applications，2019，528：121463.

［13］? SURYANTO A，DARTI I.On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination［J］.AIMS Mathematics，2020，6（1）：141.

［14］? MAO X.Stochastic Differential Equations and Applications［M］.New York：Academic Press，2006.

［15］? BERETTA E，KOLMANOVSKII V，SHAIKHET L.Stability of epidemic model with time delays influenced by stochastic perturbations［J］.Mathematics and Computers in Simulation，1998，45（3/4）：269.

［16］? SURYANTO A.A Dynamically consistent nonstandard numerical scheme for epidemic model with saturated incidence rate［J］.Int J Math Comput，2011，13：112.

［17］? SURYANTO A.Stability and bifurcation of a discrete SIS epidemic model with delay［C］//Proceedings of the 2nd International Conference on Basic Sciences，Malang，Indonesia，2012：1.

［18］? HU Z，TENG Z，JIANG H.Stability analysis in a class of discrete SIRS epidemic models［J］.Nonlinear Analysis：Real World Applications，2012，13（5）：2017.

［19］? MICKENS R E.Nonstandard Finite Difference Models of Differential Equations［M］.Singapore：World Scientific，1993.

［20］? MICKENS R E.Calculation of denominator functions for nonstandard finite difference schemes for differential equations satisfying a positivity condition［J］.Numerical Methods for Partial Differential Equations，2010，23（3）：672.

［21］? LIN Y，JIANG D.Long-time behaviour of a perturbed SIR model by white noise［J］.Discrete & Continuous Dynamical Systems（Series B），2013，18（7）：1873.

［22］? GRAY A，GREENHALGH D，HU L，et al.A stochastic differential equation SIS epidemic model［J］.SIAM Journal on Applied Mathematics，2011，71（3）：876.

［23］? DANG H N，TRAN K，DIEU N T，et al.Extinction and permanence in a stochastic SIRS model in regime-switching with general incidence rate［J］.Nonlinear Analysis：Hybrid Systems，2019，34：121.

［24］? SHAIKHET L.Stability of equilibrium states for a stochastically perturbed exponential type system of difference equations［J］.Journal of Computational & Applied Mathematics，2015，290：92.

［25］? PALMER P.Application of a discrete It formula to determine stability（instability） of the equilibrium of a scalar linear stochastic difference equation［J］.Computers & Mathematics with Applications，2012，64：2302.

［26］? SHAIKHET L.Lyapunov Functionals and Stability of Stochastic Difference Equations［M］.London：Springer Science & Business Media，2011.

（責(zé)任編輯? 馬宇鴻）