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.
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]):with initial conditions and . 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).where , and is the standard Brownian motion increment.
From the first equation of model (2), we can obtain
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 Section 2, 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: Let be scalar Brownian motions defined on the probability space. We define to be a complete filtered probability space on a finite time horizon . We assume that is the natural filtration and , denotes the norm of an Euclidean space and C is the different parameters in the paper. For any integer . We let T to be an arbitrary positive number. , , and the time increment is . By using the similar method as Lemma 2.2 and 2.3 in [3], we have the following lemma.
Lemma 1. For any , we denote , thenwhere C is a constant.
By equation (2), we have the nonnegative invariant set of model (1) as follow:
This suggests that 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:whereare step processes. That is , , and for when
Let be the deterministic grid points of . Let . 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
Case 1. Let
Case 2. LetFrom Cases 1 and 2, we can obtainwhere , and are control functions.
3. Convergence of BIM
In this section, we shall give main results for the degrees of convergence of BIM.
For , we can define the three-step functions:
Hence, equation (9) can be rewritten in the integral formwhereand with initial values .
Lemma 2. There exists a constant C dependent on S, I, and R but independent of such thatLet s(t), i(t), r(t) be the continuous EM approximate solution, and its step processes , , and are close to each other.
Proof. For , let be the integer part of .
Then we havewhich givesUsing the same method with equation (21), we can getAdding equations (21)–(23), we obtain thatWe hence haveBy the Holder inequalityusing the Doob martingale inequality, we getwhere Substituting equation (27) with equation (26) into equation (25) yields
Theorem 1. Let S(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).where C is a constant independent of .
Proof. For any , by equations (1) and (6), we haveUsing the similarity of the method and , we can have thatAdding equations (30)–(32), we obtain thatFor any , by the Doob martingale inequality and the Holder inequality, we then computeWe give the following inequality:and substituting equation (35) into equation (34) yieldsBy Lemma 2, we havewhere the C is positive constants.
Remark 1. Theorem 1 shows that and , and , and and are close to each others.
Lemma 3. The Euler approximate (6) is bounded, i.e., there exists a constant such that
Proof. By the similarity method of Lemma 2 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
Proof. By Theorem 1, we learn that the Euler method converges with order , and in the mean square sense, we haveLet . Recalling equations (6) and (17), we getand then, in order to simplify the equations, let , , and , , ,Noticing the assumption that the control functions and are uniformly bounded, instead , , into equation (42) and , we can getThen by Lemma 3, we can deriveand by using element inequalityThe above inequality and Gronwall’s inequality implies thatWe 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:
The initial values in our numerical experiment are . The control is measurable and for any . We use BIM [14] to prove the simulation of stochastic model (13).
Figure 1 shows the solution curves that we used by the EM method; it represents the variation of population for susceptible , infected and recovered , respectively. We can see that the numerical values of , , and 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 and has not appeared negative value.
In Table 1, 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 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 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 with 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research was supported by the Ningxia Key R&D Program Key Projects (2018BFG02003).