# The Positivity of Numerical Method for Susceptible-Infected-Recovered-Susceptible Epidemic Model

Authors: Feng Feng; Zong Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Sudden environmental perturbations may affect the positivity of the solution of the susceptible-infected-recovered-susceptible (SIRS) model. Most of the SIRS epidemic models have no analytical solution. Thus, in order to find the appropriate solution, the numerical technique becomes more essential for us to solve the dynamic behavior of epidemics. In this paper, we are concerned with the positivity of the numerical solution of a stochastic SIRS epidemic model. A new numerical method that is the balanced implicit method (BIM) is set, which preserves the positivity under given conditions. The BIM method can maintain positive numerical solution. An illustrative numerical instance is presented for the numerical BIM of the stochastic SIRS model.


---

## Body

## 1. Introduction

In recent years, mathematical models have been an important tool in analyzing the dynamic behaviors of infectious diseases, such as susceptible-infected-recovered (SIR) models, susceptible-infected-vaccinated-susceptible (SIVS) epidemic models, and SIRS models (see, e.g., [1–7]). Especially SIRS mathematical models attract the attention of many scholars. The positivity is a basic characteristic in plenty of areas, such as economy and ecology areas. For example, the susceptible number, infected number, and recovered number are positive in SIRS models, it means to find the positive numerical solution of SIRS model.In order to predict and control the epidemics, SIRS epidemic models have drawn much attention in recent years. For example, Lahrouz et al. [8] investigated the dynamics of a stochastic SIR epidemic model with saturated incidence and verified that a random effect may lead the disease to extinction under a simple condition. Cai et al. [9] gave the conditions of persistence and extinction for the dynamic properties around the positive equilibrium of the SIRS model. Guo et al. [10] used the Markov semigroup theory to prove the global behavior of a stochastic SIRS model with media convergence. An SIRS model with a nonlinear incidence rate and vaccination can be described by the following set of differential equations (see, e.g., [11]):(1)dSt=Λ−μSt−βStIt1+αIt2+γRtdt−σStIt1+αIt2dBt,dIt=βStIt1+αIt2−μ+ν+δItdt+σStIt1+αIt2dBt,dRt=νIt−μ+γRtdt,with initial conditions S0≥0,I0≥0 and R0≥0. Λ is the recruitment rate of the susceptible population, μ denotes the natural death rate, β presents the transmission coefficient. γ is the rate of recovered individuals who lose immunity and returned to the susceptible class, ν denotes the natural recovery rate of the infected individuals, δ is the disease-induced death rate, α presents the parameter in the incidence rate that accounts for the psychological or inhibitory effect, and σ is the intensity of the noise. Guo et al. gave the Euler–Mayuyama (EM) numerical method of SIRS [11] (see Section 5).(2)Sk+1=Sk+Λ−μSk−βSkIk1+αIk2+γRkΔt−σSkIk1+αIk2ΔBt,Ik+1=Ik+βSkIk1+αIk2−μ+ν+δIkΔt+σSkIk1+αIk2ΔBt,Rk+1=Rk+νIk−μ+γRkΔt,where ΔBt=Btk+1−Btk, and ΔBt is the standard Brownian motion increment.From the first equation of model (2), we can obtain(3)PSk+1Sk=PSk−βSkIk1+αIk2+γRkΔ−σSkIk1+αIk2ΔB<0>0.By equation (3), we can see that the EM approximation cannot guarantee the positivity of the SIRS model. So a question raises, how it is possible to prevent a numerical integration scheme from becoming negative.The idea of balancing was introduced by Ruan and Wang [12], where BIM turned out to be efficient for preserving positivity of stochastic differential equations (SDEs). Motivated by the balancing idea in [12], we have formulated an improved algorithm which can paint the numerical solution of the SIRS model to remain positive and converged to analytical solution. Our contributions are as follows:(i)
We constructed a new BIM numerical method for the stochastic SIRS model with local Lipschitz condition, and the numerical scheme will keep positive solution for the stochastic system.(ii)
The maximum error estimates for the stochastic SIRS model have been derived, which demonstrated that our schemes have stability.The outline of the paper is arrayed as follows: in Section2, we have presented some assumptions and preparations for studying the numerical method for stochastic SIRS equation (1). In Section 3, we have constructed a new BIM numerical method preserving positivity for the stochastic SIRS system and proved the numerical solutions to have strong convergence in probability. In Section 4, we have given the numerical simulations of the EM and BIM methods. Finally, some conclusions are drawn in Section 5.


## 2. Preliminaries and Model Formulation

Throughout this paper, we have used the following notation: LetBt be scalar Brownian motions defined on the probability space. We define Ω,ℱ,ℱt0≤t≤T,ℙ to be a complete filtered probability space on a finite time horizon 0,T. We assume that ℱt0≤t≤T is the natural filtration and ℱ=ℱT, ⋅ denotes the norm of an Euclidean space and C is the different parameters in the paper. For any integer p≥2. We let T to be an arbitrary positive number. T=N⋅Δt, tn=n⋅Δt, and the time increment is Δ=T/N≪1. By using the similar method as Lemma 2.2 and 2.3 in [3], we have the following lemma.Lemma 1.
For anyη>0, we denote Xi=St,It,RtT, then(4)∑i=13Esup0≤t≤TXitη≤C,where C is a constant.By equation (2), we have the nonnegative invariant set of model (1) as follow:(5)F=xt∈R+3:Λμ+δ<X1t+X2t+X3t≤Λμ⊂R+3.This suggests thatXiti=1,2,3 is bounded by constant Λ/μ, which means that the inequality (5) holds. By using Euler—Maruyama (EM) method, we obtain the corresponding discretization equation of the SIRS model as follows:(6)st=s0+∫0tΛ−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu∫0tσs¯ui¯u1+αi¯2udBu,it=i0+∫0tβs¯ui¯u1+αi¯2u−μ+ν+δi¯udu+∫0tσs¯ui¯u1+αi¯2udBu,rt=r0+∫0tνi¯u−μ+γr¯udu,where(7)s¯t=∑k=0T/Δsk1tk,tk+1t,i¯t=∑k=0T/Δik1tk,tk+1t,r¯t=∑k=0T/Δxk1tk,tk+1t,are step processes. That is s¯t=sk, i¯t=ik, and r¯t=rk for t∈tk,tk+1 when(8)k=0,1,2,3,…,TΔ.Lett1,t2,…,tm be the deterministic grid points of 0,T. Let Δt=tn+1−tn,ti=i⋅Δt. ΔBn=Btn+1−Btn denotes the increments of the time and Brownian motion, respectively. For system (1), the balanced implicit method (BIM) approximate solution is defined by the scheme(9)Sn+1=Sn+Λ−μSn−βSnIn1+αIn2+γRnΔt−σSnIn1+αIn2ΔBn+CnSn−Sn+1,In+1=In+βSnIn1+αIn2−μ+ν+δInΔt+σSnIn1+αIn2ΔBn+DnIn−In+1,Rn+1=Rn+νIn−μ+γRnΔt+EnRn−Rn+1.Case 1.
Let(10)F=Λ−μSn−βSnIn1+αIn2+γRn,(11)C1t,Sn,In=μSn+βSnIn1+αIn2−1,ifF≥0,Λ+γRn−1,ifF<0.Case 2.
Let(12)G=βSnIn1+αIn2−μ+ν+δIn,(13)D1t,In=μ+ν+δIn−1,ifG≥0,βSnIn1+αIn2−1,ifG<0.
From Cases1 and 2, we can obtain(14)Cn=C1t,Sn,InΔ+μαInσΛΔBn,Dn=D1t,InΔ+αΛσSnμΔBn,En=μμνIn−Λμ−γΛΔ,where Cn, Dn, and En are control functions.


## 3. Convergence of BIM

In this section, we shall give main results for the degrees of convergence of BIM.Fort∈tn,tn+1, we can define the three-step functions:(15)S¯t≔S¯t=Sn,t∈tn,tn+1,I¯t≔I¯t=In,t∈tn,tn+1,R¯t≔R¯t=Rn,t∈tn,tn+1,(16)C¯t≔Cn,D¯t≔Dn,E¯t≔En,t∈tn,tn+1,B¯t≔Btn+1−Btn,t∈tn,tn+1.Hence, equation (9) can be rewritten in the integral form(17)stB=s0+∫0t1+Cs−1Λ−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu−∫0t1+Cs−1σs¯ui¯u1+αi¯2udBu,itB=i0+∫0t1+Ds−1βs¯ui¯u1+αi¯2u−μ+ν+δi¯udu+∫0t1+Ds−1σs¯ui¯u1+αi¯2udBu,rtB=r0+∫0t1+Es−1νi¯u−μ+γr¯udu,where(18)Cs=c0S¯t,I¯t,R¯tΔ+c1S¯t,I¯tB¯ts,Ds=d0S¯t,I¯t,R¯tΔ+d1S¯t,I¯B¯ts,Es=e0I¯t,R¯tΔ,and with initial values S0>0,I0>0.Lemma 2.
There exists a constant C dependent onS, I, and R but independent of Δ such that(19)sup0≤t≤Tst−s¯t2+it−i¯t2+rt−r¯t2≤CΔ1−1/p.
Let s(t), i(t), r(t) be the continuous EM approximate solution, and its step processes s¯t, i¯t, and r¯t are close to each other.Proof.
Fort∈0,T, let T/Δ be the integer part of t/Δ.
Then we have(20)st−s¯t=∫T/ΔΔt−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu+∫T/ΔΔtσs¯ui¯u1+αi2¯udBu,which gives(21)st−s¯t2≤2−μs¯t−βs¯ti¯t1+αi¯2t+γr¯t2Δ2+2σs¯tit1+αi¯2t2Bt−BtΔ2≤6μ2s¯2t+β2s¯2ti¯2t1+αi¯2t2+γ2r¯2tΔ2+2σs¯tit1+αi¯2t2Bt−BtΔ2.
Using the same method with equation (21), we can get(22)it−i¯t2≤4β2s¯2ti¯2t1+αi¯2t2+4μ+ν+δ2i¯2tΔ2+2σ2s¯2ti¯2t1+αi¯2t2Bt−BtΔ2,(23)rt−r¯t2≤4ν2i¯2t+4μ+γ2r¯2tΔ2.
Adding equations (21)–(23), we obtain that(24)st−s¯t2+it−i¯t2+rt−r¯t2≤10β2s¯2ti¯2t1+αi¯2t2+6μ2s¯2t+4μ+ν+δ2+4ν2i¯2t+4μ+γ2+6γ2r¯2tΔ2+4σ2s¯2ti¯2t1+αi¯2t2Bt−BtΔ2≤10β2α2+Ms¯t+i¯t+r¯t2Δ2+4σ2s¯2ti¯2t1+αi¯2t2Bt−BtΔ2≤10β2α2+MΛμ2Δ2+4σ2s¯2ti¯2t1+αi¯2t2Bt−BtΔ2.
We hence have(25)sup0≤t≤Tst−s¯t2+it−i¯t2+rt−r¯t2≤10β2α2+MΛμ2Δ2+4σ2s¯2ti¯2t1+αi¯2t2sup0≤t≤T∧qBt−BΔtΔ2≤10β2α2+MΛμ2Δ2+4σ2α2sup0≤t≤TBt−BΔtΔ2.
By the Holder inequality(26)Esup0≤t≤TBt−BΔtΔ2≤Esup0≤t≤TBt−BΔtΔ2p1/p,using the Doob martingale inequality, we get(27)sup0≤t≤TBt−BΔtΔ2p≤Esup0≤k≤T/ΔsupkΔ≤r≤Δk+1Br−BΔk2p≤∑k=0T/ΔEsupkΔ≤r≤Δk+1Br−BΔk2p≤2p2p−12p∑k=0T/ΔEBΔk+1−BΔk2p≤2p2p−12p∑k=0T/Δ2p−1!!Δp≤2p2p−12p2p−1!!Δp−1T+1,where 2p−1!!=2p−1×2p−3×⋯×3×1. Substituting equation (27) with equation (26) into equation (25) yields(28)Esup0≤t≤Tst−s¯t2+it−i¯t2+rt−r¯t2≤10β2α2+MΛμ2Δ2+4σ2α22p2p−12p2p−1!!Δp−1T+11/p≤CΔ1−1/p.Theorem 1.
LetS(t), I(t), R(t) be the solution and s(t), i(t), r(t) be the continuous EM approximate solution to the stochastic SIRS model of (1).(29)Esup0≤t≤TSt−st2+It−it2+Rt−rt2≤CΔ1−1/p,where C is a constant independent of Δ.Proof.
For any0≤t≤T, by equations (1) and (6), we have(30)St−st2≤2μ2∫0tSu−s¯udu2+2β2∫0tSuIu1+αI2u−s¯ui¯u1+αi¯2udu2+2γ2∫0tRu−r¯udu2−2σ2∫0tSuIu1+αI2u−s¯ui¯u1+αi¯2udBu2.
Using the similarity of the method andξ=μ+ν+δ, we can have that(31)It−it2≤2β2∫0tSuIu1+αI2u−s¯ui¯u1+αi¯2udu2+2ξ2∫0tIu−i¯udu2+2σ2∫0tSuIu1+αI2u−s¯ui¯u1+αi¯2udBu2,(32)Rt−rt2≤2ν2∫0tIu−i¯udu2+2μ+γ2∫0tRu−r¯udu2.
Adding equations (30)–(32), we obtain that(33)St−st2+It−it2+Rt−rt2≤4β2∫0tSuIu1+αI2u−s¯ui¯u1+αi¯2udu2+2μ2∫0tSu−s¯udu2+2ξ2+2ν2∫0tIu−i¯udu2+2γ2+2μ+γ2∫0tRu−r¯udu2.
For anyt1∈0,T, by the Doob martingale inequality and the Holder inequality, we then compute(34)sup0≤t≤TSt−st2+It−it2+Rt−rt2≤4β2TE∫0t1SuIu1+αI2u−s¯ui¯u1+αi¯2u2du+2μ2TE∫0t1Su−s¯u2du+2ξ2+2ν2TE∫0tIu−i¯udu2+2γ2+2μ+γ2TE∫0tRu−r¯udu2≤8β2TE∫0t1SuIu1+αI2u−suiu1+αi2u2du+8β2TE∫0t1SuIu1+αI2u−s¯ui¯u1+αi¯2u2du+4μ2TE∫0t1Su−su2du+4μ2TE∫0t1Su−s¯u2du+22ξ2+2ν2TE∫​0tIu−iudu2+22ξ2+2ν2TE∫​0t∧qIu−i¯udu2+22γ2+2μ+γ2TE∫0tRu−rudu2+22γ2+2μ+γ2TE∫0tRu−r¯udu2.
We give the following inequality:(35)∫0t1ESuIu1+αI2u−suiu1+αi2u2du+∫0t1Esuiu1+αi2u−s¯ui¯u1+αi¯2u2du≤I1+αI2∫0t1Su−su2du−S1+αI2∫0t1Iu−iu2du+i1+αi2∫0t1su−s¯u2du−s1+αi2∫0t1iu−i¯u2du,and substituting equation (35) into equation (34) yields(36)sup0≤t≤TSt−st2+It−it2+Rt−rt2≤4μ2+8β2I1+αI2∫0t1Su−su2du+4ξ2+4ν2−8β2S1+αI2∫0tIu−iu2du+4γ2+4μ+γ2∫0tRu−ru2du+4μ2+8β2i1+αi2∫0t1su−s¯u2du+4ξ2+4ν2−8β2s1+αi2∫0tiu−i¯u2du+4γ2+4μ+γ2∫0tru−r¯u2du.
By Lemma2, we have(37)sup0≤t≤TSt−st2+It−it2+Rt−rt2≤C∫0tSt−st2+It−it2+Rt−rt2du≤CΔ1−1/p,where the C is positive constants.Remark 1.
Theorem1 shows that St and st, It and it, and Rt and rt are close to each others.Lemma 3.
The Euler approximate (6) is bounded, i.e., there exists a constant C>0 such that(38)sup0≤t≤Tst2+it2+rt2≤C.Proof.
By the similarity method of Lemma2 in [13].Theorem 2.
The BIM (17) numerical solution convergences to the true solution of system (1) with strong order 4, i.e., there is constant A such that(39)sup0≤t≤TSt−stB2+It−itB2+Rt−rtB2≤CΔ1−1/p.Proof.
By Theorem1, we learn that the Euler method converges with order 1−1/p, and in the mean square sense, we have(40)sup0≤t≤TSt−stB2+It−itB2+Rt−rtB2=sup0≤t≤TSt−st+st−stB2+It−it+it−itB2+Rt−rt+rt−rtB2≤2sup0≤t≤TSt−st2+2sup0≤t≤Tst−stB2+2sup0≤t≤TIt−it2+2sup0≤t≤Tit−itB2+2sup0≤t≤TRt−rt2+2sup0≤t≤Trt−rtB2≤2CΔ1−1/p+2sup0≤t≤Tst−stB2+it−itB2+rt−rtB2.
LetV=st−stB2+it−itB2+rt−rtB2. Recalling equations (6) and (17), we get(41)dV=2st−stBCs1+Cs−1Λ−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu−σs¯ui¯u1+αi¯2udBu+12Cs2+Csσs¯ui¯u1+αi¯2u2du+2it−itBDs1+Ds−1βs¯ui¯u1+αi¯2u−μ+ν+δi¯udu+σs¯ui¯u1+αi¯2udBu+12Ds2+Dsσs¯ui¯u1+αi¯2u2du+2rt−rtBEs1+Es−1νi¯u−μ+γr¯udu,and then, in order to simplify the equations, let SE≔st, IE≔it, RE≔rt and SB≔stB, IB≔itB, RB≔rtB,(42)dV≤2SE−SBCsΛ−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu−σs¯ui¯u1+αi¯2udBu+12Cs2+Csσs¯ui¯u1+αi¯2u2du+2IE−IBDsβs¯ui¯u1+αi¯2u−μ+ν+δi¯udu+σs¯ui¯u1+αi¯2udBu+12Ds2+Dsσs¯ui¯u1+αi¯2u2du+2RE−RBEsνi¯u−μ+γr¯udu.
Noticing the assumption that the control functionsC1t,Sn,In and D1t,In are uniformly bounded, instead Cs, Ds, Es into equation (42) and dB2=dt, we can get(43)V≤2SE−SBC1t,Sn,InΔt+μαIσΛΔB∫0tΛ−μs¯u−βs¯ui¯u1+αi¯2u+γr¯udu−∫0tσs¯ui¯u1+αi¯2udBu+12Cs2∫0tσs¯ui¯u1+αi¯2u2du+2IE−IBD1t,InΔt+αΛσSμΔB∫0tβs¯ui¯u1+αi¯2u−μ+ν+δi¯udu+∫0tσs¯ui¯u1+αi¯2udBu+12Ds2∫0tσs¯ui¯u1+αi¯2u2du+2RE−RBμμνIn−Λμ−γΛΔt∫0tνi¯u−μ+γr¯udu.
Then by Lemma3, we can derive(44)dV≤2SE−SBΔ2−Δ+12Δ2+2IE−IBΔ2+Δ+Δ3/2+12Δ2+Δ2≤2SE−SBΔ2−Δ+2IE−IBΔ2+Δ+Δ3/2+2Δ2,and by using element inequality(45)st−stB2+it−itB2+rt−rtB2≤A∫0Tst−stB2+it−itB2+rt−rtB2du+∫0TΔ2−Δ2+Δ2+Δ+Δ3/22+2Δ4dt.
The above inequality and Gronwall’s inequality implies that(46)sup0≤t≤TSt−stB2+It−itB2+Rt−rtB2≤CΔ1−1/p.
We have proof of the BIM (17) numerical solution of SIRS model convergences to the true solution of system (1).


## 4. Numerical Simulation

In this section, we perform numerical simulation to verify the validity of our theoretical results.We use the EM method [14] to conduct the simulation of stochastic model (2), and the parameter values are chosen as follows:(47)Λ=0.06,β=0.04,μ=0.03,m=0.021,δ=0.01,γ=0.1,a=1,σ=0.05,α=0.1,v=0.1;The initial values in our numerical experiment areS0=2,I0=1.5,R0=0.4. The control Ck,Dk,Ek is measurable and for any t∈0,T. We use BIM [14] to prove the simulation of stochastic model (13).Figure1 shows the solution curves that we used by the EM method; it represents the variation of population for susceptible St, infected It and recovered St, respectively. We can see that the numerical values of St, It, and Rt emerge negative. Figure 2 shows the solution curves that are used by the BIM method for the stochastic SIRS model, and we can see that the numerical value of St and It has not appeared negative value.Figure 1
The path of variables ofS(t), I(t), R(t) with the EM numerical method.Figure 2
The path of variables ofS(t), I(t), R(t) with the numerical BIM.In Table1, we present the percentage of negative paths of susceptible S(t), infected I(t), and recovered R(t) with SIRS system (1), respectively. The results are obtained by averaging all of the 1000 samples. The percentage of the negative numerical solution of susceptible St with the EM method governed by equation (2) condition (27.75%; 25.67%; 15.42%) under response stepsize (0.125; 0.025; 0.0125).Table 1
The experience of two numerical methods.
Number
S
I
R
Stepsize
Δ
=
0.25
Δ
=
0.025
Δ
=
0.0125
Δ
=
0.25
Δ
=
0.025
Δ
=
0.0125
Δ
=
0.25
Δ
=
0.025
Δ
=
0.0125
EM (%)
27.75
25.67
15.42
24.52
22.07
15.22
14.68
13.42
12.31
BIM (%)
0
0
0
0
0
0
0
0
0Table1 shows that the numerical BIM can preserve positivity for all step sizes Δ, independent of the length of integration interval [0, T]. Just as the theoretical results, the numerical solution of the SIRS model with the BIM method can remain positive.


## 5. Concluding Remarks

In this paper, we proposed a novel numerical method based on the EM numerical method to keep positivity of the SIRS model. Preserving positivity of the proposed method is proved under given conditions. Numerical results indicate that BIM is efficient for maintaining positivity of the SIRS epidemic model. Then, we also give the prove that the BIM numerical solutions have convergence order with1−1/p in probability.However, in this model the numerical solutions of the SIRS model have low convergence precision. Therefore, it is interesting to find a high-precision numerical solution of the epidemic model. We leave this for further consideration.



---

*Source: 6825284-2020-03-19.xml*