Nonstandard Finite Difference Schemes for an SIR Epidemic Model

: This paper aims to present two nonstandard ﬁnite difference (NFSD) methods to solve an SIR epidemic model. The proposed methods have important properties such as positivity and boundedness and they also preserve conservation law. Numerical comparisons conﬁrm that the accuracy of our method is better than that of other existing standard methods such as the second-order Runge–Kutta (RK2) method, the Euler method and some ready-made MATLAB codes.


Introduction
The modeling of physical, chemical, and biological systems is commonly based on nonlinear ordinary differential equations (ODEs). Mathematical models and their simulations are important for understanding the quality and quantity of these systems. Numerical methods, such as Euler and Runge-Kutta methods for solving the mentioned equations, are typically used [1][2][3][4][5][6][7][8][9][10][11][12]. A downside of these methods is that the qualitative characteristics of a true solution are not transferred to the numerical solution. To overcome such a deficiency, the NSFD methods were introduced by Mickens [13,14]. Indeed, these methods were proposed in order to preserve the main characteristics of systems of differential equations, which preserve positivity, boundedness, stability of equilibrium points and conservation law . For a method to be designated as an NSFD, at least one of the following criteria should be fulfilled:

•
In the first-order discrete derivatives, step-size h in the denominator is replaced by a function φ(h), such that • Nonlinear terms in f (y) can be approximated non-locally (by a suitable function of several of the mesh points), for instance: The main goal of this article is to introduce and present two NSFD methods to solve an SIR epidemic model. The constructed methods preserve the characteristics of the mathematical model such as positivity, boundedness, stability of equilibrium points and conservation law. The rest of the paper is structured as follows: In Section 2, we introduce the mathematical model for the SIR epidemic model with its equilibrium points. We also investigate the positivity properties of the model. In Section 3, the new schemes are proposed. The positivity property is investigated in Section 4. In Section 5, we investigate the elementary stability of the new schemes. In Section 6, we present the results obtained by the new schemes and some schemes such as the Euler method, ode45 (a MATLAB's standard solver for ODEs ), and so on.

Mathematical Model
Here, the mathematical model of the SIR epidemic introduced in [36] is presented. This model consists of the following three characteristics: The mathematical model of the SIR epidemic presented in [42] is related to rubella disease in London. The model is introduced by a nonlinear differential equation system as follows: and they are all positive.
The consecutive system of equations, if added together, clearly satisfies conservation law, which implies that the population is constant. The population can be normalized to one, since it is assumed to be constant: The above equation must be valid for any numerical method. The reproduction number of system (2) is as follows: the system (2) has two equilibrium points as follows:

Construction of New Schemes
In this section, we propose a family of NSFD schemes for the system (2) in the following form: with γ + δ = 1, η + θ = 1 and Note that this scheme satisfies R (t) + S (t) + I (t) = 0, as n → ∞. In addition, it follows from (3) that Thus, if S n + R n + I n = 1 for all n ≥ 1 and, since, γ = 1 − δ then S n+1 + R n+1 + I n+1 = 1 for all n ≥ 1. Therefore, the conservation law property holds.
We have chosen two different formulas of (3), given as follows:

Positivity
In this section, the positive property of the proposed methods in the previous section is investigated. Positivity means that the non-negativity of initial data is preserved over time in the approximate solution. It should be noted that the positive property of a numerical method used to solve differential models related to population biology is very important because the variables used in the subpopulation models can never be negative values (see, for example, [43]).

Definition 1. A finite difference method is positive if
∀h > 0(step size), and y 0 ∈ R + (initial condition), the other solutions are also positive, that is, we have y k ∈ R + for any k ∈ N. Proof. Since the parameters B, N, ν, a and µ are positive, then the numerical schemes (4) and (5) represent positive NSFD methods for any ψ(h) > 0 when the initial values S(0), I(0) and R(0) are non-negative.

Elementary Stability
In this section, the stability properties of the proposed methods are investigated and sufficient conditions to preserve this property of the equilibrium points of model (2) are proposed. The finitedifference methods that have this stability property are called elementary stable methods [44].

Definition 2.
The finite difference method is called elementary stable, if for any value of the step size h, its only equilibrium pointsĒ are those of the original differential system, the linear stability properties of eachĒ being the same for both the differential system and the discrete method.

Theorem 2. Equations
Proof. In order to examine the convergence of schemes (4) and (5), it is enough to only consider the first and second equations of (4) and (5). The equilibrium points of (4) are exactly the equilibria E * 0 and E * of system (2). The Jacobian J of the Equation (4) from which therefore, it is obvious that |λ 1 | < 1 always holds and if R 0 < 1, that is, NB < µ + ν, then |λ 2 | < 1. Hence, the disease equilibrium is asymptotically stable if R 0 < 1 and otherwise it will be unstable. If R 0 > 1 system (2) has an endemic equilibrium point. The Jacobian matrix of Equation (4) evaluated at the endemic equilibrium is Eigenvalues of J(E * ) are roots of the quadratic equation since c > d and b > e and, also, Therefore, the equilibrium E * is stable if and only if all conditions of Lemma 1 hold. Since Λ and ω are positive, we have Since 1 b < 1 and also, From (6)- (8) we see that the conditions of Lemma 1 hold. Thus, when R 0 > 1, the eigenvalues of J(E * ) are less than the unity in the modulus, regardless of the size of h. Therefore dynamic consistence between the system (2) and Equation (4) is held around all equilibria, which leads to the elementary stability of (4).
Similarly, the equilibrium points of (5) are the same equilibria of the original system.
The Jacobian J of (5) has the form J(S n , I n ) = Substituting E * 0 in J(S n , I n ), we obtain Therefore, it is clear that |λ 1 | < 1 always holds and if R 0 < 1, that is, NB < µ + ν, then |λ 2 | < 1. Hence, the disease equilibrium is asymptotically stable if R 0 < 1; otherwise, it will be unstable. If R 0 > 1, system (2) has an endemic equilibrium point. The Jacobian matrix of Equation (5) evaluated at the endemic equilibrium is Eigenvalues of J(E * ) are roots of the quadratic equation since c − d > 0 and b − e > 0, and also, Therefore, the equilibrium E * is stable if and only if all conditions of Lemma 1 hold. Since Λ and ω are positive we have Since e b < 1 and From (9)-(11), we see that the conditions of Lemma 1 holds. Then, the eigenvalues of J(E * ) are less than the unity in the modulus, regardless of the size of h, when R 0 > 1. Therefore dynamic consistence between the system (2) and Equation (5) is held around all equilibria, which leads to the elementary stability of (5).

Numerical Results
In this section, we present some numerical results obtained by applying the new methods. These results indicate the efficiency of the methods. Consider model (2) with the parameters µ = 0.04, ν = 24, NB = 123.
The epidemic model (2) has a disease free equilibrium point for R 0 < 1 and an endemic equilibrium point for R 0 > 1. Figure 1 shows that, for R 0 > 1, the new methods converge to the endemic equilibrium point and produce positive results for all times, while the ready-made MATLAB codes do not converge or sometimes produce negative values.  It is shown that all numerical solutions obtained by the proposed NSFD schemes with a = 0.05, 1, 5, 10 are convergent to the EE point. It can be seen that NSFD schemes with a great inhibitory constant have better convergence to an endemic equilibrium point. In Figures 5 and 6, we used the same parameters as above for µ and ν, but we choose NB = 23 so that R 0 = 0.9567. As can be seen for each h, new schemes retain the positivity property, while the Euler method and the RK2 method give negative values when using the SIR model. It should be noted that, due to the very large numerical outputs at intervals larger than the range given in Figure 6, the system is unable to display those values.

Conclusions
In this work, two NSFD methods were proposed to solve the SIR pandemic model. The proposed model has two equilibrium points, which are dubbed disease free equilibrium points, that is, E * 0 , in which case R 0 < 1, and the endemic equilibrium point E * , wherein R 0 > 1. It is shown that the proposed methods preserve features such as the positivity property, elementary stability and conservation law. By comparing the numerical results obtained from the new methods and the Runge-Kutta methods, it was shown that the NSFD numerical methods are made unconditionally stable by preserving conservation law for R 0 < 1. They are also convergent to the disease-free equilibrium point for any step size. This behavior is also satisfied for R 0 > 1. Moreover, the same behavior is obtained for R 0 > 1. It is also confirmed that the ready-made MATLAB codes do not converge to the endemic equilibrium point. Therefore, it can be concluded that the new NSFD methods preserve features of the SIR epidemic model.