Application of operational matrices to numerical solution of stochastic SIR model

The goal of this paper is to give useful method for solving a problem in biologic system that is formulated by stochastic Volterra integral equations. Here, we consider triangular functions, block pulse functions and their operational matrices of integration. Illustrative example is included to demonstrate the validity and applicability of the operational matrices.

The numerical study and simulation of stochastic Volterra integral equations have been an active field of research for the last years [9,12]. Most SVIEs do not have analytic solutions and hence it is of great importance to provide numerical schemes. The paper [12] solved stochastic Volterra integral equations by stochastic operational matrix based on block pulse functions (BPFs). Also, [9] considered stochastic operational matrix of triangular functions.
In the next section, stochastic concepts of Brownian motion and Ito integral are presented. In Sect. 3, we review block pulse functions, triangular functions and their operational matrices of integration both in the deterministic and stochastic case. In Sect. 4, mathematical model of SIR is presented, then these operational matrices are applied to solve stochastic SIR model and numerical results are shown. Finally, Sect. 5 provides a brief conclusion of this work.  [15].
Proof See [15]. [15] Let f ∈ ν(S, T ). Then the Itô integral of f (from S to T) is defined by

Definition 2.4 (The Itô integral)
where φ n is a sequence of elementary functions such that Proof See [15].
Definition 2.6 (one-dimensional Itô processes) [15] Let B(t) be one-dimensional Brownian motion on ( , F, P). A one-dimensional Itô process (stochastic integral) is a stochastic process X (t) on ( , F, P) of the form where P t 0 v 2 (s, ω)ds < ∞, for all t ≥ 0 = 1, Theorem 2.7 (The one-dimensional Itô formula) Let X (t) be an Itô process given by (1) and Then is again an Itô process, and is computed according to the rules Proof See [15].
Some properties of the Itô integral are listed as follows. Let f, g ∈ ν(0, T ) and let 0 ≤ S < U < T . Then For more details see [10,15].

Operational matrix of integration for some typical functions
The goal of this section is to recall notations and definitions of some typical functions that are used in the next sections. Their properties and operational matrices are reviewed.
The BPFs are defined on the time interval [0, T ) by where i = 1, . . . , m and for convenience we put h = T m . The BPFs on [0, T ) have disjointness, orthogonality and completeness property. The set of BPFs can be written as a vector From the above representation and disjointness property, it follows that where F is a m-dimensional vector andF = diag(F). The expansion of a function where Now, operational matrix of integration is considered where Q is operational matrix of integration that is given by So, The Itô integral of BPFs is defined by where stochastic operational matrix of integration is given by So, the Itô integral of every function f (t) can be approximated by

Triangular functions (TFs)
Two m-sets of triangular functions are defined [5] over the interval [0, T ) as e l s e w h e r e (17) where i = 0, . . . , m − 1 and m is the number of elementary functions and h = T m . Moreover, From the definition of TFs, it is clear that TFs are disjoint, orthogonal, and complete [5].
We can write TFs in vector form and where T (t) is called the TF vector. Put where diag(T (t)) is a 2m × 2m diagonal matrix. Furthermore, The expansion of f (t) over [0, T ) with respect to 1D-TFs, is given by where C1 i and C2 i are samples of f , for example, The vector C is called the TF coefficient vector. Consider where where M is operational matrix of integration base on TFs. Therefore, we can use the below approximation Now, consider the following definitions: Therefore, where is stochastic operational matrix of integration in the TF domain. So, the Itô integral of f (t) is approximated with

Stochastic SIR model
Most models for the transmission of infectious diseases descend from the classical SIR model of Kermack and McKendrick established in 1927, [8]. One of the most basic SIR models is ⎧ ⎨ where the parameters , β, , μ, γ are positive constants, and S(t), I (t), R(t) denote the number of the individuals susceptible to the disease, number of infected members and number of members who have been removed from the possibility of infection through full immunity, respectively. Some notable features of the model: The influx of individuals into the susceptibles is given by a constant , it is assumed that the natural death rates are assumed to be equal (denoted by constant μ) and individuals in I (t) suffer an additional death due to disease with rate constant , β and γ represent the disease transmission coefficient and the rate of recovery from infection respectively. The basic reproduction number R 0 = β μ(μ+ +γ ) is the threshold of the system for an epidemic to occur. If R 0 ≤ 1, model (38) has only the disease-free equilibrium E 0 = ( μ , 0, 0) which is globally asymptotical stable. This means the disease will disappear and the entire population will become susceptible. If R 0 > 1, E 0 becomes unstable and there exists a global asymptotically stable endemic equilibrium E * = ( μ+ +γ β , μ+ +γ − μ β , γ μ(μ+ +γ ) − γ β ) which implies the disease always remains [7]. However, in the real world, epidemic models are inevitably affected by environmental noise, which is an important component in realism. Environmental noise can provide an additional degree of realism compared to their deterministic counterparts. Stochastic models depend on the chance variations in risk of exposure, disease and other illness dynamics. In this paper, our approach to include stochastic perturbation is analogous to that of Imhof and Walcher [6]. Here we assume that stochastic perturbations are of a white noise type which are directly proportional to S(t), I (t), R(t), influenced on theṠ(t),İ (t),Ṙ(t) in the model (38). By this way, the model (38) will be deduced to the form: where B i (t) are independent standard Brownian motions and σ 2 i ≥ 0 represent the intensities of B i (t), i = 1, 2, 3.
It is proved that the stochastic model led to extinction even though the deterministic counterpart predicts persistence. In paper [7], unique global positive solution of this model was considered. Also, the asymptotic behavior of this solution was investigated.

Application of operational matrices to solve SIR model
By integration of Eq. (39), this method is applied for the following system We approximate function S(t), S 0 (t), I (t), I 0 (t), R(t), R 0 (t) and t by BPFs and TFs.
For convenient, we choose (t), as a vector form of these functions: ( (t), m-dimensional vector form of BPFs, T (t), 2m-dimensional vector form of TFs and P as a operational matrix of them, also, we put P S as a stochastic operational matrix of the discussed functions.) where the vectors S, S 0 , I, I 0 , R, R 0 are function coefficient of S(t), S 0 (t), We replace S T (t) T (t)I by S TĨ (t). So, After replacing with =, By solving nonlinear system (51) we find S, I, R and finally S(t), I (t), R(t) of (41) are approximated.

Numerical example
In practice, we typically observed only a single realization of this process, a single path, out of a multitude of possible paths. The simulated sample path can then be analyzed by usual statistical method to determine how good the approximation is and in what sense it is closed to the exact solution.
In this section, we solve Eq.  It is easy to compute that R 0 = 0.8 ≤ 1. The computing result is in good agreement with Theorem (4.1). It is obvious that the fluctuations of the curves increase as the noise level increases. Also, we can believe the solution is stochastically asymptotically stable in the large.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Conclusion
Because we cannot solve some stochastic problems analytically, in this article we presented a new technique for solving systems of SVIEs numerically. We considered some operational matrices of integration. The benefits of this method are lower cost of setting up the system of equations without any integration, coupled with the low computational operation cost. Finally, convergence of TFs is faster than BPFs [9] and order of convergence is O(h 2 ). These advantages make the method very simple.