一類隨機(jī)SIR傳染病模型的非標(biāo)準(zhǔn)數(shù)值離散化

譚偉 劉茂省

DOI：10.16783/j.cnki.nwnuz.2024.01.006

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

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

作者簡介：譚偉（1996—），男，河北張家口人，碩士研究生.主要研究方向?yàn)閭魅静?dòng)力學(xué).

Email：1848231722@qq.com

*通信聯(lián)系人，男，教授，博士.主要研究方向?yàn)閭魅静?dòng)力系統(tǒng)及復(fù)雜網(wǎng)絡(luò).Email：liumaoxing@126.com

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

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

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

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

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

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

其中S≡S（t）為t時(shí)刻易感人群數(shù)量，I≡I（t）為t時(shí)刻已感染人群數(shù)量，R≡R（t）為t時(shí)刻染病后免疫人群的數(shù)量，且S（0）≥0，I（0）≥0，R（0）≥0；Λ為從外部進(jìn)入種群的凈遷移率，μ為人口的自然死亡率，ν為接種疫苗的人口比率，γ為患病種群恢復(fù)的概率，δ為因?yàn)楦腥炯膊《劳龅娜丝诟怕剩宜袇?shù)都是正的.

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

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

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

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

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

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

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

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

基本再生數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

（ii）非線性項(xiàng)是非局部近似的，例如

u2（tn）≈unun+1.

然而，注意到NSFD格式對(duì)一個(gè)微分方程來說并不是唯一的，所以為了確保NSFD格式與原始連續(xù)模型具有相同的動(dòng)力學(xué)性態(tài)，我們需要對(duì)NSFD格式進(jìn)行動(dòng)力學(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）是對(duì)系統(tǒng)（3）進(jìn)行廣義前向差分和右側(cè)局部近似得到的，通過一些簡單的代數(shù)計(jì)算，我們可以很容易地驗(yàn)證，具有NSFD格式的模型（4）與初始連續(xù)模型（2）具有完全相同的平衡點(diǎn).數(shù)學(xué)期望我們通常由E來表示，ωi（n+1）（i=1，2）（n∈Z）是一個(gè)相互獨(dú)立且服從標(biāo)準(zhǔn)正態(tài)分布N（0，1）的隨機(jī)序列，并且對(duì)于ω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）

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

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

對(duì)離散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）

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

1? 預(yù)備知識(shí)

對(duì)系統(tǒng)（7）的地方病平衡點(diǎn)進(jìn)行平移變換，令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(yù)（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）的零解等價(jià)于系統(tǒng)（7）的正解E+=（S+，I+）.

在平衡點(diǎn)E+=（S+，I+）處對(duì)系統(tǒng)（7）進(jìn)行線性化，得到線性近似系統(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）.

類似地進(jìn)行平移變換.對(duì)無病平衡點(diǎn)，令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(yù)（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）.

在無病平衡點(diǎn)處對(duì)系統(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)置.

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

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

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

對(duì)非負(fù)函數(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（文獻(xiàn)［26］定理1.1）? 對(duì)于線性系統(tǒng)（9）或（11），若存在一個(gè)滿足條件E［V（0，φ）］≤c1φ2和E［ΔVi］≤-c2E［z（i）2］（i∈Z）的非負(fù)函數(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是隨機(jī)項(xiàng)中的常數(shù)，ψ表示分母函數(shù)，ωi（n+1）（i=1，2）是相互獨(dú)立的隨機(jī)序列且滿足條件（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是對(duì)稱矩陣.對(duì)于實(shí)對(duì)稱矩陣P和Q，如果P-Q是一個(gè)正定矩陣，那么P>Q.

定理B（文獻(xiàn)［26］定理5.1）? 假設(shè)存在正定矩陣P，使得形為（13）的半正定矩陣D滿足矩陣方程ATDA-D=-P，同時(shí)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）.

然后計(jì)算Δ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? 地方病平衡點(diǎn)的穩(wěn)定性

因?yàn)橄到y(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? 若存在一個(gè)正定矩陣P，使得條件（18）對(duì)于矩陣方程（17）的解（d11，d12，d22）被滿足，并且矩陣D是半正定的，那么系統(tǒng)（9）的零解是漸近均方穩(wěn)定的.因此，系統(tǒng)（8）的零解是依概率穩(wěn)定的，這等同于系統(tǒng)（7）的正平衡點(diǎn)Ee是依概率穩(wěn)定的.

當(dāng)R0>1時(shí)，系統(tǒng)（7）存在地方病平衡點(diǎn)E+=（S+，I+）.對(duì)于數(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，從而地方病平衡點(diǎn)存在，且E+=（1.4597，4.2629）.通過對(duì)（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）的地方病平衡點(diǎn)E+是依概率穩(wěn)定的.

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

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

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

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

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

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

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

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

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

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

此時(shí)R0=0.167<1，所以只有一個(gè)無病平衡點(diǎn)存在，且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）的無病平衡點(diǎn).

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

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

when σ1=σ2=0

when σ1=σ2=0.5

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

A=0.9317-0.390401.3221.

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

D=7.511511.786911.786913.4013，

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

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

4? 結(jié)論

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

參考文獻(xià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é)任編輯? 馬宇鴻）