具扩散SIR传染病模型的平衡态的稳定性分析

李遵先,苟长义

(天津理工大学理学院,天津 300384)

0 引言

传染病对人类社会的发展具有重要影响.传染病所引起的人类死亡数远超过去所有战争所导致的人口死亡数[1].用数学建模方法来研究传染病的传染规律是研究传染病并对其进行控制的有效方法之一.自1927年文献[2-3]提出SIR仓室模型以来,大量的描述传染病的数学模型被建立,可见文献[1,4].文献[5]给出了一个经典的SIR传染病模型：

模型(1)将人口分为易感者、感染者以及恢复者三类.其中：s=s(t),i=i(t)和r=r(t)分别代表t时刻三类人口所占总人口的比例,正参数μ,β,γ分别代表出生(死亡)率、感染率和恢复率.方程(1)等价于：

考虑模型(2)的初值问题,文献[5]有如下结论.

命题1设 (s(t),i(t))为方程组(2)在 T={(s,i)s≥0,i≥0,s+i≤1}中的解.记若σ≤1或i(0)=0,则所有从T中出发的轨道趋于平衡点(1,0).若σ＞1,则存在正平衡点(se,ie)=,使得当i(0)＞0时,所有从T中出发的轨道均趋于此正平衡点.

基于模型(2)的众多传染病模型被建立,近期工作如文献[6-11].由命题1知,模型(2)具有阈值称为基本再生数.其生物意义为单位时间内每个感染者传染易感者个数.以下考虑具扩散的模型(2)是否仍有此阈值性质?一般地,无扩散时平衡态的稳定性与具扩散时不同,即Turing不稳定性[12].如文献[13]考虑一类具扩散的传染病模型,得到当基本再生数大于1时,无论是否存在扩散项,无疾平衡点(态)均不稳定.但文献[14]在研究一类具年龄结构的捕食-被捕食模型时,得到其无扩散时的平衡点在考虑扩散效应后从不稳定变成了局部渐近稳定,即产生了Turing不稳定性.为此,考虑以下系统：

其中：Δ为Laplace算子,Ω⊂Rn是一个具光滑边界∂Ω的有界区域.齐次Neumann边界条件表明模型(3)是封闭的,即在边界上没有人口的流入和流出.模型(3)对应的平衡态方程如下：

易见,当σ≤1时,系统(4)有一个非负常数平衡态(1,0)；当σ＞1时,系统(4)有两个非负常数平衡态(1,0)和(se,ie)=.以下分析此模型(3)的非负常数平衡态的(局部)稳定性.

1 平衡态的稳定性分析

参照文献[15],令0=λ0＜λ1＜λ2＜λ3＜…是 -Δ在区域Ω上具Neumann边界条件的特征值,并令X={(s,i)∈[C1(Ω)]2：∂νs=∂νi=0,x∈∂Ω},E(λ)={φ： - Δφ =λφ,x∈Ω, ∂νφ =0,x∈∂Ω}(λ∈R).设是 E(λi)的一组标准正交基,以及 Xij={cφij：c∈ R2},那么.其中：

当σ＜1时,对非负平衡态(1,0),有如下结论.

定理1若σ＜1,模型(3)的平衡态(1,0)是局部渐近稳定的.

证明 将模型(3)的平衡态(1,0)平移到(0,0)得：

方程(5)的线性化方程为：

相应地,当σ＞1时,对非负平衡态(1,0)和(se,ie),有如下两个结论.

定理2若σ＞1,模型(3)的平衡态(1,0)是不稳定的.

证明 与定理1证明类似,因考虑平衡态(1,0)的稳定性问题需考虑算子L1的特征值分布情况.由于此时σ ＞1,即β-γ-μ＞0.从而无法得到对任意j={0,1,…},都有det A1j＞0及Tr A1j＜0成立.事实上,取j=0,可知：

由此可知A10存在具正实部的特征值,故(1,0)不稳定.

定理3若σ＞1,模型(3)的平衡态(se,ie)是局部渐近稳定的.

证明 与定理1证明类似,将模型(3)的平衡态(se,ie)平移到(0,0)得：

方程(7)的线性化方程为：

2 数值模拟

如下对模型(3)进行数值模拟.参考文献[5,13],分别取以下两组方程参数：

图1 s(x,t)趋于1Fig.1 s(x,t)tends to 1

由图1～2可见,平衡态(1,0)局部渐近稳定.当σ=3时,模型(3)数值模拟如图3～4所示：

图3 s(x,t)趋于0.333 3Fig.3 s(x,t)tends to 0.333 3

图4 i(x,t)趋于0.031 7Fig.4 i(x,t)tends to 0.031 7

由图3～4可知,此时平衡态(1,0)不稳定且(se,ie)局部渐近稳定.

3 结果与讨论

1)由定理1～3知,具扩散的SIR传染病模型具有与命题1类似的阈值性质,即对于模型(3),决定平衡态稳定性的仍然是基本再生数,而与扩散无关.

2)命题1考虑了σ=1的情况,但定理2～3没有考虑此情形.作为此临界问题,需要进一步讨论方程高阶项的情况,这将在后续工作中考虑.

3)由图4可知,当σ＞1时,传染病在传染初期出现局部的峰值即小的爆发,这一结论对模型(3)也存在[5].

4 附录

附录 1aj≜ det A1j=d1d2+{μd2-(β -γ -μ)d1}λj-μ(β -γ -μ) ＞ 0,bj≜Tr A1j=-(d1+d2)λj+(β-γ-μ)-μ ＜0.设为A1j的两个特征值,故=aj＞0,+=bj＜0.对任意j={0,1,…},考虑两种情况：

附录2已知aj≜det A2j=d1d2+ μσd2λj+ μ(σ -1)(γ + μ) ＞ 0,bj≜TrA2j=-(d1+d2)λjμσ ＜0.设为A2j的两个特征值,故-=aj＞0,+=bj＜0.∀j={0,1,…},考虑两种情况：