Next Article in Journal
Modeling and Analysis of the Influence of Fear on a Harvested Food Web System
Next Article in Special Issue
Network Thermodynamics-Based Scalable Compartmental Model for Multi-Strain Epidemics
Previous Article in Journal
A Multi-Mechanism Seagull Optimization Algorithm Incorporating Generalized Opposition-Based Nonlinear Boundary Processing
Previous Article in Special Issue
A Multi-Scale Model for Cholera Outbreaks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Efficient Numerical Solutions to a SIR Epidemic Model

1
Department of Mathematics, Faculty of Science, University of Maragheh, Maragheh 55181-83111, Iran
2
Scientific Computing Group, Universidad de Salamanca, Plaza de la Merced, 37008 Salamanca, Spain
3
School of Mathematics and Information Science, Henan Polytechnic University, Jiaozuo 454000, China
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(18), 3299; https://doi.org/10.3390/math10183299
Submission received: 12 July 2022 / Revised: 7 September 2022 / Accepted: 8 September 2022 / Published: 11 September 2022
(This article belongs to the Special Issue Mathematical Methods and Models in Epidemiology)

Abstract

:
Two non-standard predictor-corrector type finite difference methods for a SIR epidemic model are proposed. The methods have useful and significant features, such as positivity, basic stability, boundedness and preservation of the conservation laws. The proposed schemes are compared with classical fourth order Runge–Kutta and non-standard difference methods (NSFD). The stability analysis is studied and numerical simulations are provided.
MSC:
65L06; 03H05; 03H10

1. Introduction

Ordinary differential conditions are extensively used in demonstrating numerous natural and physical applications. Mathematical strategies dependent on finite differences [1,2], Taylor series [3], interpolation, such as Runge–Kutta, Euler, and multistep techniques [4,5], and some different strategies [6,7] are broadly utilized. A large number of problems in mathematical epidemiology are modeled by autonomous systems of nonlinear ordinary differential equations, which implies that the boundaries of the model are autonomous, regarding time. In these models, the factors address subpopulations of susceptibles, infected, recovered, etc. Consequently, the solutions of the ODE framework portray the advancement of the various classes of subpopulations in the model over the long haul. The use of different schemes brings up such issues as what the truncation error is or how large the region of steadiness is. Numerous standard strategies, such as forward Euler, Runge–Kutta and others, fail to demonstrate the non-actual motions, bifurcations and chaos (see, for example, [8]).
One method for forestalling these kinds of mathematical problems is the development of mathematical methods dependent on non-standard finite-difference techniques. This sort of technique was initially devised by Mickens [9,10]. Piyanwong et al. [11] and Jansen and Twizell [1] have planned positive and genuinely stable plans for the SIR and SEIR epidemic models, individually. In any case, in their created solutions, they have not applied the preservation law unequivocally, which can prompt impossible or unreasonable arrangements.
The best situation to accurately resolve an ODE based-model is the point at which a careful difference solution can be developed in [12]. Although the non-standard finite difference technique preserve the standard properties of the approximation solution, such as consistency and convergence, it can also preserve the qualitative properties of the solution, such as boundedness, monotonicity, positivity, and so on [13,14,15,16,17,18,19,20,21,22]. In this paper, we have developed two predictor–corrector types of NSFD techniques to obtain mathematical results for the SIR epidemic model. The new mathematical strategy has the heuristic properties of the ODE framework solution and is basically stable. These strategies are also very practical due to the huge time step used.
The sections of the paper are as follows: In Section 2, we introduce the SIR epidemic model and its steady-state point numerical model. In Section 3, we present the new methods. In Section 4 and Section 5, we explore positivity and basic stability results. In Section 6, other NSFD techniques are constructed to be used as correctors in the event of larger stepsizes. In Section 7, we introduce the results of the new methods compared to different schemes used for inspection, such as the the classical Runge–Kutta (RK4) strategy, ode45, etc. We finish the paper with the conclusions and discussion.

2. The SIR Model

We consider the SIR epidemic model as given in [11]. This model incorporates three independent variables to represent the individuals:
s ( t ) : susceptible i ( t ) : infected r ( t ) : recovered
A variant of this model was utilized in [23] to depict the elements of beating hack epidemics in London involving occasional varieties in susceptibility. The model is depicted by a nonlinear ODEs framework as follows:
s ( t ) = μ N μ s ( t ) β N i ( t ) s ( t ) , i ( t ) = ( μ + ν ) i ( t ) + β N s ( t ) i ( t ) , r ( t ) = μ r ( t ) + ν i ( t ) , t > 0 , i ( 0 ) = i 0 , s ( 0 ) = s 0 , r ( 0 ) = r 0 ,
in which
  • β is the coefficient of transmission;
  • μ is the death rate or the birth rate, which are equal to each other;
  • ν is the recovery rate;
  • N is the total number of individuals, N ( t ) = i ( t ) + s ( t ) + r ( t ) .
All are assumed to be positive.
The equations of the continuous system, if added together, satisfy the conservation law
d N d t = 0 .
Since N is constant, system (1) can be written as follows by referring to [24] with the assumption S ( t ) = s ( t ) N , I = i ( t ) N and R ( t ) = r ( t ) N
S ( t ) = μ μ S ( t ) β I ( t ) S ( t ) , I ( t ) = ( μ + ν ) I ( t ) + β S ( t ) I ( t ) , R ( t ) = μ R ( t ) + ν I ( t ) , t > 0 , I ( 0 ) = I 0 , S ( 0 ) = S 0 , R ( 0 ) = R 0 ,
where
s ( t ) : susceptible i ( t ) : infected r ( t ) : recovered
The reproduction number of (2) is
R 0 = β μ + ν .
If 1 > R 0 , model (2) has steady-state (disease-free) E 0 * = ( 1 , 0 , 0 ) , and it has a unique endemic steady-state, then
E * = 1 R 0 , μ μ + ν ( 1 1 R 0 ) , ν μ + ν ( 1 1 R 0 ) ,
if R 0 > 1 .
Definition 1.
A discrete version of (2) is called the nonstandard method provided that one of the conditions below is satisfied:
 (i)
In the discrete derivatives of d S d t , d I d t and d R d t , a non-negative function φ ( h ) substitutes the step size h, such that
φ ( h ) = h + O ( h 2 ) a s 0 < h 0 ;
 (ii)
Nonlinear terms in the right hand side of (2) are approximated in a nonlocal way, that is to say, by an appropriate function of some points in the mesh.
For instance,
S 2 S n S n + 1 , S 2 S n S n + 1 , S 3 α S n 2 S n + 1 + ( 1 α ) S n S n + 1 S n 1 .

3. Construction of New Schemes

We express system (2) as
S n + 1 S n φ ( h ) = μ μ ( δ S n + γ S n + 1 ) β ( η S n + 1 + θ S n ) I n , I n + 1 I n φ ( h ) = β ( η S n + 1 + θ S n ) I n μ ( γ I n + 1 + δ I n ) ν ( θ I n + η I n + 1 ) , R n + 1 R n φ ( h ) = ν ( θ I n + η I n + 1 ) μ ( γ R n + 1 + δ R n ) ,
with 1 = η + θ = γ + δ ,
φ ( h ) = e x p ( h ) 1 .
We are free to choose the parameters keeping the condition 1 = η + θ = γ + δ , and here we chose theta to achieve a set of stable explicit methods.
Remark 1.
The condition 0 < φ ( h ) < 1 for h > 0 is essential, and there exist several functions satisfying this condition, e.g., φ ( h ) = e x p ( h ) 1 , φ ( h ) = 1 e x p ( h ) or φ ( h ) = s i n ( h ) .
It should be noted that by writing the Taylor expansion of S n + 1 and φ ( h ) = h + O ( h 2 ) in the first relation of Equation (4), we have
S n + h S n + h 2 2 ! S n + S n h + O ( h 2 ) = μ μ ( δ S n + γ ( S n + h S n + h 2 2 ! S n + ) ) β ( η ( S n + h S n + h 2 2 ! S n + ) + θ S n ) I n
with h 0 , from which we have
S n = μ μ ( δ S n + γ S n ) β ( η ( S n + θ S n ) ) I n = μ μ ( δ + γ ) S n β ( ( η + θ ) S n ) I n
On the other hand, δ + γ = 1 and η + θ = 1 , then
S n = μ μ S n β S n I n
The second and third relations of the equation can be written in the same way. Therefore, we can move from the discrete form to the continuous form with h 0 .
Lemma 1.
The new NSFD family in (4) preserves the conservation law.
Proof. 
We use induction to prove it. Since S + I + R = 1 , for the initial values, it is S 0 + I 0 + R 0 = 1 . Hence, for n = 0 , S 1 + I 1 + R 1 1 = μ φ ( h ) γ ( 1 ( S 1 + I 1 + R 1 ) ) , and thus S 1 + I 1 + R 1 = 1 , which implies S n + 1 + I n + 1 + R n + 1 = 1 ; as a result, the new family preserves the conservation law. □
After investigating the characteristics of the new family, we present the following methods:

3.1. Scheme 1

Letting γ = 3 2 , θ = 0 , δ = 1 2 , η = 1 leads to
S n + 1 = μ φ ( h ) + S n ( 1 + 1 2 φ ( h ) μ ) 3 2 μ φ ( h ) + φ ( h ) β I n + 1 = F ( I n , S n ) ,
I n + 1 = I n 1 2 μ φ ( h ) + 1 + β φ ( h ) F ( S n , I n ) 3 2 μ φ ( h ) + ν φ ( h ) + 1 = G ( S n , I n )
R n + 1 = 1 I n + 1 S n + 1 .

3.2. Scheme 2

Letting θ = 0 , γ = 2 , η = 1 , δ = 1 gives
S n + 1 = φ ( h ) μ + ( 1 + φ ( h ) μ ) S n 1 + 2 φ ( h ) μ + φ ( h ) I n β = F ¯ ( S n , I n ) ,
I n + 1 = I n μ φ ( h ) + 1 + β φ ( h ) F ¯ ( I n , S n ) 1 + 2 φ ( h ) μ + ν φ ( h ) = G ¯ ( I n , S n )
R n + 1 = 1 I n + 1 S n + 1 .

4. Positivity

In this part, we analyze the positivity features of the proposed strategies. With positivity, we mean that the component-wise non-negativity of the initial vector is preserved in time for the numerical solution. It must be mentioned that the positivity property of a mathematical technique is significant when it is employed to address differential models emerging in population science because these state factors address subpopulations, which never take negative qualities. Many papers have been written about the positivity property (see for instance [25]).
Definition 2.
A finite difference technique is said to preserve the positivity property, if, for h and y 0 R + n , y k R + n for every k N .
Theorem 1.
NSFD schemes (5) and (6) are positivity preserving.
Proof. 
Because ν , β , μ > 0 , the (5) and (6) positive NSFD schemes for any φ ( h ) > 0 if 0 < S n < 1 , 0 < I n < 1 and 0 < R n < 1 for all n 0 . □

5. Elementary Stability

We give sufficient conditions for schemes (5) and (6) to maintain the stability properties of the steady points of the model (2). A difference scheme with this stability property is called an elementary stable scheme [26].
Definition 3.
If the linear stability features of the discrete and differential models are the same, the finite-difference scheme is said to be elementary stable.
The following result can be readily obtained by using the well-known Jury conditions [27].
Lemma 2.
Consider λ 2 A λ + B = 0 . Both roots satisfy | λ i | < 1 , if
  • B < 1 ,
  • 1 + A + B > 0 ,
  • 1 A + B > 0 .
Theorem 2.
Schemes in (5) and (6) are elementary stable.
Proof. 
Consider only Equations (5a), (5b), (6a) and (6b): Steady-state points of (5) are E 0 * and E * of (2).
The Jacobian of (5) is J ( S n , I n ) = F S F I G S G I , where
F S = 1 + 1 2 μ φ ( h ) 3 2 φ ( h ) μ + φ ( h ) β I n + 1 , F I = β φ ( h ) μ φ ( h ) + S n ( 1 + 1 2 μ φ ( h ) ) ( 1 + 3 2 μ φ ( h ) + β φ ( h ) I n ) 2 , G S = β φ ( h ) I n F S 1 + 3 2 μ φ ( h ) + ν φ ( h ) = β φ ( h ) I n ( 1 + 1 2 μ φ ( h ) ) ( 1 + 3 2 μ φ ( h ) + φ ( h ) β I n ) ( 1 + 3 2 μ φ ( h ) + ν φ ( h ) ) ,
G I = 1 + 1 2 μ φ ( h ) + β φ ( h ) F ( I n , S n ) + β I n φ ( h ) F I 3 2 μ φ ( h ) + 1 + φ ( h ) ν = 1 + 1 2 μ φ ( h ) + β φ ( h ) μ φ ( h ) + S n ( 1 + 1 2 φ ( h ) ) 1 + 3 2 μ φ ( h ) + φ ( h ) β I n β 2 φ 2 ( h ) I n μ φ ( h ) + S n ( 1 + 1 2 μ φ ( h ) ) ( 1 + 3 2 μ φ ( h ) + β φ ( h ) I n ) 2 1 + 3 2 μ φ ( h ) + ν φ ( h ) .
Writing E 0 * into J ( S n , I n ) ,
J E 0 * = 2 + φ ( h ) μ 2 + 3 μ φ ( h ) 2 β φ ( h ) 2 + 3 μ φ ( h ) 0 2 + μ φ ( h ) + 2 β φ ( h ) 2 + 3 μ φ ( h ) + 2 φ ( h ) ν ,
from which we obtain the two eigenvalues
λ 1 = 2 + μ φ ( h ) 2 + 3 μ φ ( h ) , λ 2 = 2 + μ φ ( h ) + 2 β φ ( h ) 2 + 3 μ φ ( h ) + 2 ν φ ( h ) .
Hence, we have | λ 1 | < 1 . On the other hand, if R 0 < 1 , i.e., β < μ + ν , then | λ 2 | < 1 . However, if R 0 < 1 , the disease-free steady-state point is asymptotically stable; otherwise, it is unstable.
If 1 < R 0 , the system (2) has an endemic steady state. We have
J E * = a b c b 2 a f b d 1 d ( a + c b c f b 2 ) ,
where
a = 1 + 1 2 μ φ ( h ) > 1 , b = 1 + 3 2 μ φ ( h ) + β φ ( h ) I * = 1 + 3 2 μ φ ( h ) + μ φ ( h ) ( R 0 1 ) = 1 + 1 2 μ φ ( h ) + μ φ ( h ) R 0 > 1 , c = β φ ( h ) φ ( h ) μ + S * ( 1 + 1 2 μ φ ( h ) ) = β φ ( h ) μ φ ( h ) + 1 R 0 ( 1 + 1 2 μ φ ( h ) ) = ( μ + ν ) μ φ 2 ( h ) R 0 + ( μ + ν ) μ φ ( h ) ( 1 + 1 2 μ φ ( h ) ) > 0 , f = β φ ( h ) I * = μ φ ( h ) ( R 0 1 ) > 0 , d = 1 + 3 2 μ φ ( h ) + ν φ ( h ) > 0 .
The characteristic equation obtained from J ( E * ) is λ 2 A λ + B = 0 in which
A = T r a c e J ( E * ) = a b 2 + b d a + c ( b f ) b 2 d ,
B = D e t J ( E * ) = a 2 b 2 + a c b c a f b 3 d + c a f b 3 d = a 2 b 2 + a c b b 3 d = a 2 b + a c b 2 d .
Thus, E * is stable if Lemma 2 holds. Clearly B > 0 since b > f , it is also A > 0 . Since A and B are positive, we have
1 + A + B > 0 .
Additionally, since a b < 1 and
1 d ( a + c b c f b 2 ) < 1 d ( a + c b ) = a b + c b d = 1 ,
we have
1 A + B = ( 1 a b ) 1 1 d ( a + c b c f b 2 ) + a c f b 3 d > 0 .
Finally, we have
B = a 2 b d + a c b 2 d < a 2 b d + c b d = a 2 + c b d = ( 1 + 1 2 μ φ ( h ) ) 2 + μ β φ ( h ) φ ( h ) + 1 R 0 ( 1 + 1 2 μ φ ( h ) ) ( 1 + 1 2 φ ( h ) μ + μ φ ( h ) R 0 ) ( 1 + 3 2 μ φ ( h ) + ν φ ( h ) ) < 1 ,
From (7), (8) and (9), we observe that Lemma 2 holds. Hence, J ( E * ) has eigenvalues which are < 1 in modulus, regardless of the size of h, as long as R 0 > 1 . Therefore, we have verified the dynamical consistency between system (2) and scheme (5) around every steady state that results in the elementary stability of (5).
Similarly, we have J ( S n , I n ) = F ¯ S F ¯ I G ¯ S G ¯ I , in which
F ¯ S = 1 + μ φ ( h ) 1 + 2 μ φ ( h ) + φ ( h ) β I n , F ¯ I = β φ ( h ) μ φ ( h ) + S n ( 1 + μ φ ( h ) ) ( 1 + 2 μ φ ( h ) + β φ ( h ) I n ) 2 , G ¯ S = β φ ( h ) I n F ¯ S 1 + 2 μ φ ( h ) + ν φ ( h ) = β φ ( h ) I n ( 1 + μ φ ( h ) ) ( 1 + 2 μ φ ( h ) + φ ( h ) β I n ) ( 1 + 2 μ φ ( h ) + ν φ ( h ) ) , G ¯ I = 1 + φ ( h ) μ + β φ ( h ) F ¯ ( S n , I n ) + β φ ( h ) I n F ¯ I 1 + 2 μ φ ( h ) + ν φ ( h ) , = 1 + φ ( h ) μ + β φ ( h ) μ φ ( h ) + S n ( 1 + φ ( h ) μ ) 1 + 2 μ φ ( h ) + φ ( h ) β I n β 2 φ 2 ( h ) I n μ φ ( h ) + S n ( 1 + φ ( h ) μ ) ( 1 + 2 μ φ ( h ) + β φ ( h ) I n ) 2 1 + 2 μ φ ( h ) + ν φ ( h ) .
Writing E 0 * in J ( S n , I n ) , we obtain
J E 0 * = 1 + μ φ ( h ) 1 + 2 μ φ ( h ) β φ ( h ) 1 + 2 μ φ ( h ) 0 1 + μ φ ( h ) + β φ ( h ) 1 + 2 μ φ ( h ) + ν φ ( h ) ,
hence, we obtain the two eigenvalues
λ 1 = 1 + μ φ ( h ) 1 + 2 μ φ ( h ) , λ 2 = 1 + μ φ ( h ) + β φ ( h ) 1 + 2 μ φ ( h ) + ν φ ( h ) .
It is clear that | λ 1 | < 1 always holds. On the other hand, if R 0 < 1 , i.e., β < μ + ν , then | λ 2 | < 1 . Therefore, if R 0 < 1 , the disease-free equilibrium is asymptotically stable; otherwise, it is unstable. If R 0 > 1 , (2) has an endemic steady state.
The Jacobian at the endemic steady state is
J E * = a b c b 2 a f b d 1 d ( a + c b c f b 2 ) ,
where
a = 1 + φ ( h ) μ > 1 , b = 1 + 2 μ φ ( h ) + β φ ( h ) I * = 1 + 2 μ φ ( h ) + μ φ ( h ) ( R 0 1 ) = 1 + μ φ ( h ) + μ φ ( h ) R 0 > 1 , c = β φ ( h ) μ φ ( h ) + S * ( 1 + μ φ ( h ) ) = β φ ( h ) μ φ ( h ) + 1 R 0 ( 1 + μ φ ( h ) ) > 0 , f = 1 β φ ( h ) I * = φ ( h ) μ ( R 0 1 ) > 0 , d = 1 + 2 μ φ ( h ) + φ ( h ) ν > 0 .
Hence, from λ 2 A λ + B = 0 , we obtain the eigenvalues as
A = T r a c e J ( E * ) = a b 2 + b d a + c ( b f ) b 2 d ,
B = D e t J ( E * ) = a 2 b 2 + a c b c a f b 3 d + c a f b 3 d = a 2 b 2 + a c b b 3 d = a 2 b + a c b 2 d .
Thus, E * is stable if Lemma 2 holds. Clearly B > 0 , and, since b > f , it is also A > 0 . Since A and B are positive, we have
1 + A + B > 0 .
Again, since a b < 1 and
1 d ( a + c b c f b 2 ) < 1 d ( a + c b ) = a b + c b d = 1
we have
1 A + B = ( 1 a b ) 1 1 d ( a + c b c f b 2 ) + a c f b 3 d > 0 .
Finally, we have
B = a 2 b d + a c b 2 d < a 2 b d + c b d = a 2 + c b d = ( 1 + μ φ ( h ) ) 2 + β φ ( h ) μ φ ( h ) + 1 R 0 ( 1 + μ 1 2 φ ( h ) ) ( 1 + μ φ ( h ) + μ φ ( h ) R 0 ) ( 1 + 2 μ φ ( h ) + ν φ ( h ) ) < 1 .
From (10), (11) and (12) we observe that Lemma 2 holds. Accordingly, J ( E * ) has eigenvalues which are less than 1 in modulus, regardless of h, provided that R 0 > 1 . We have therefore verified the dynamical consistence between system (2) and scheme (6) around each steady state that gives us the elementary stability of (6). □

6. The New Nonstandard Discretizations of SIR Model

To improve the performance of the NSFD schemes (5) and (6), we formulate them by a predictor–corrector type method. In the NSFD methods, if step size h is relatively large, some of their features, including convergence, conservation law and positivity, may be lost. Therefore, it is useful to build robust computational methods that can cope with these drawbacks. In general, explicit methods produce larger errors than implicit methods, but implicit methods require solving nonlinear models. Next we will introduce two NSFD schemes of predictor-corrector type schemes (5) and (6), where, for larger step sizes, h maintains the aforementioned significant features. To develop scheme (5), we employ (5) as a predictor scheme, that is
S p n + 1 = μ φ ( h ) + S n ( 1 + 1 2 μ φ ( h ) ) 1 + 3 2 μ φ ( h ) + φ ( h ) β I n ,
I p n + 1 = 1 + 1 2 μ φ ( h ) + β φ ( h ) S p n + 1 I n 1 + 3 2 μ φ ( h ) + ν φ ( h ) ,
Now, we introduce an implicit NSFD method for solving the system in (2)
S n + 1 S n φ ( h ) = μ μ ( 3 S n + 1 S n 2 ) β S n + 1 I n + 1 S n + 1 φ ( h ) + S n + 1 φ ( h ) , I n + 1 I n φ ( h ) = β S n + 1 I n + 1 μ ( 3 I n + 1 I n 2 ) ν I n + 1 , R n + 1 R n φ ( h ) = ν I n + 1 μ ( 3 R n + 1 R n 2 ) .
Thus,
S c n + 1 = μ φ ( h ) + S n ( 1 + 1 2 μ φ ( h ) ) + S p n + 1 2 + 3 2 μ φ ( h ) + φ ( h ) β I p n + 1 ,
I c n + 1 = 1 + 1 2 μ φ ( h ) I n + β φ ( h ) S p n + 1 I p n + 1 1 + 3 2 μ φ ( h ) + ν φ ( h ) ,
R n + 1 = 1 S c n + 1 I c n + 1 .
Similarly, to develop method (6),
S p n + 1 = μ φ ( h ) + S n ( 1 + φ ( h ) μ ) 1 + 2 μ φ ( h ) + φ ( h ) β I n ,
I p n + 1 = 1 + μ φ ( h ) + β φ ( h ) S p n + 1 I n 1 + 2 μ φ ( h ) + ν φ ( h )
Applying the implicit NSFD method to solve (2), we obtain
S n + 1 S n φ ( h ) = μ μ ( 2 S n + 1 S n ) β S n + 1 I n + 1 S n + 1 φ ( h ) + S n + 1 φ ( h ) , I n + 1 I n φ ( h ) = β S n + 1 I n + 1 μ ( 2 I n + 1 I n ) ν I n + 1 , R n + 1 R n φ ( h ) = ν I n + 1 μ ( 2 R n + 1 R n ) .
By the conservation law (since the population is constant), we have
S c n + 1 = μ φ ( h ) + S n ( 1 + μ φ ( h ) ) + S p n + 1 2 + 2 μ φ ( h ) + φ ( h ) β I p n + 1 ,
I c n + 1 = 1 + μ φ ( h ) I n + β φ ( h ) S p n + 1 I p n + 1 1 + 2 μ φ ( h ) + ν φ ( h )
R n + 1 = 1 I c n + 1 S c n + 1 .
The algorithm to get the computational solution may be written as follows:
  • Step 1 Choose 0 < ϵ < < 1 , and I 0 , S 0 , R 0 such that S 0 + I 0 + R 0 = 1 .
  • Step 2 For n = 0 , 1 , do
  • Step 3 Evaluate S p n + 1 .
  • Step 4 Via S p n + 1 and I n , evaluate I p n + 1 .
  • Step 5 Correct the value S c n + 1 , using S n , S p n + 1 , I p n + 1 .
  • Step 6 Correct the value I c n + 1 , using I n , S p n + 1 , I p n + 1 .
  • Step 7 If S c n + 1 S n < ϵ and I c n + 1 I n < ϵ then
  • Step 8 Calculate R n + 1 , else S n = S c n + 1 , I n = I c n + 1 and go to step 5.

7. Numerical Results

In this section, we address some useful simulations to affirm the hypothetical outcomes and illustrate the upsides of the developed NSFD methods. Consider the model in (2) by constants
β = 123 , ν = 24 , μ = 0.04 .
If R 0 < 1 , (2) has an asymptotically stable disease-free point and an endemic equilibrium point if R 0 > 1 . Figure 1 illustrates that for R 0 > 1 , the present method will approach the correct epidemic point and generate positive quantities for all time t. At the same time, the conventional Matlab programming procedures do not meet or sometimes produce unreasonable negative qualities for the contaminated population. It should be noted that the fourth request graph of Runge–Kutta with a period step of h = 0.005 and the new plan with a cycle step of h = 0.01 are combined with the correct popularity points. The comparisons of the absolute errors provided by different methods for different values of h show that the proposed schemes are more precise than the other techniques. The absolute value of the errors with the new methods, the method presented in [28], and Runge-Kutta method as a reference method are illustrated in Figure 2. Figure 3 shows that the behavior of the P-C schemes with h = 0.01 is similar to the Runge–Kutta method with h = 0.005 , while the NSFD schemes do not work like them. This means that P-C methods perform better than NSFD methods. We can also see in Figure 4 that P-C 1 has lower dissipation than P-C 2.

8. Conclusions

In this article, we introduced two competing non-standard finite differences (NSFD) of the of predictor-correctdor type for the classic SIR epidemic model. The model has two biological equilibrium points: one is disease-free equilibrium point E 0 * , if and only if R 0 < 1 is an asymptotically stable node, and the other is endemic equilibrium point E * . The positivity, boundedness and stability of the proposed scheme are shown. The numerical comparison between the NSFD numerical scheme proposed in this paper and the Runge-Kutta type scheme shows that the NSFD numerical scheme satisfying the conservation law is unconditionally stable for R 0 < 1 , and converges to the time step of the long disease-free equilibrium point. In addition, the same behavior is obtained for R 0 > 1 . Furthermore, we showed that the well-known methods in the Matlab software package did not converge to the popular equilibrium point. We conclude that the developed non-standard scheme retains the basic characteristics of the continuous SIR model. In the new scheme, large time steps can be used, making them more economical to reach a steady state over a long period of time.

Author Contributions

Funding acquisition, S.-W.Y.; methodology, M.M.K., A.S., H.R., S.-W.Y. and M.M. All authors planned the scheme, developed the mathematical modeling and examined the theory validation. The manuscript was written through the contribution by all authors. All authors have read and agreed to the published version of the manuscript.

Funding

National Natural Science Foundation of China (No. 71601072), the Fundamental Research Funds for the Universities of Henan Province (No. NSFRF210314) and Innovative Research Team of Henan Polytechnic University (No. T2022-7). This research received funding support from the National Science, Research and Innovation Fund (NSRF), Thailand.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ashyralyev, A.; Agirseven, D.; Agarwal, R.P. Stability estimates for delay parabolic differential and difference equations. Appl. Comput. Math. 2020, 19, 175–204. [Google Scholar]
  2. Ashyralyev, A.; Erdogan, A.S.; Tekalan, S.N. An investigation on finite difference method for the first order partial differential equation with the nonlocal boundary condition. Appl. Comput. Math. 2019, 18, 247–260. [Google Scholar]
  3. Odibat, Z. Fractional power series solutions of fractional differential equations by using generalized Taylor series. Appl. Comput. Math. 2020, 19, 47–58. [Google Scholar]
  4. Khalsaraei, M.M.; Shokri, A. The new classes of high order implicit six-step P-stable multiderivative methods for the numerical solution of schrödinger equation. Appl. Comput. Math. 2020, 19, 59–86. [Google Scholar]
  5. Khalsaraei, M.M.; Shokri, A. A new explicit singularly P-stable four-step method for the numerical solution of second order IVPs. Iranian J. Math. Chem. 2020, 11, 17–31. [Google Scholar]
  6. Ramos, H.; Popescu, P. How many k-step linear block methods exist and which of them is the most efficient and simplest one? Appl. Math. Comput. 2018, 316, 296–309. [Google Scholar] [CrossRef]
  7. Ramos, H.; Rufai, M.A. Numerical solution of boundary value problems by using an optimized two-step block method, Numer. Algorithms 2020, 84, 229–251. [Google Scholar] [CrossRef]
  8. Lambert, J.D. Computational Methods in Ordinary Differential Equations; Wiley: New York, NY, USA, 1973. [Google Scholar]
  9. Mickens, R.E. Nonstandard Finite Difference Models of Differential Equations; World Scientific: Singapore, 1994. [Google Scholar]
  10. Mickens, R.E. Nonstandard finite difference schemes for differential equations. J. Differ. Equations Appl. 2002, 8, 823–847. [Google Scholar] [CrossRef]
  11. Piyawong, W.; Twizell, E.H.; Gumel, A.B. An unconditionally convergent finite difference scheme for the SIR model. Appl. Math. Comput. 2003, 146, 611–625. [Google Scholar] [CrossRef]
  12. Ramos, H. Contributions to the development of differential systems exactly solved by multistep finite-difference schemes. Appl. Math. Comput. 2010, 217, 639–649. [Google Scholar] [CrossRef]
  13. Shokri, A.; Khalsaraei, M.M.; Molayi, M. Nonstandard Dynamically Consistent Numerical Methods for MSEIR Model. J. Appl. Comput. Mech. 2021, 8, 196–205. [Google Scholar]
  14. Shokri, A.; Khalsaraei, M.M.; Molayi, M. Dynamically Consistent NSFD Methods for Predator-prey System. J. Appl. Comput. Mech. 2021, 7, 1565–1574. [Google Scholar] [CrossRef]
  15. Anguelov, R.; Lubuma, J.M.S. Contributions to the mathematics of the nonstandard finite difference method and applications. Numer. Meth. Par. Diff. Eq. 2001, 17, 518–543. [Google Scholar] [CrossRef]
  16. Chen, C.; Wang, W.; Wang, X.; Wise, S.M. Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential. J. Comput. Phys. 2019, X3, 100031. [Google Scholar] [CrossRef]
  17. Dong, L.; Wang, C.; Zhang, H.A.; Zhang, Z. A positivity-preserving, energy stable and convergent numerical scheme for the cahn-hilliard equation with a flory-huggins-degennes energy. Commun. Math. Sci. 2019, 17, 921–939. [Google Scholar] [CrossRef]
  18. Iskenderov, N.S.; Allahverdiyeva, S.I. An inverse boundary value problem for the boussinesq-love equation with nonlocal integral condition. TWMS J. Pure Appl. Math. 2020, 11, 226–237. [Google Scholar]
  19. Qalandarov, A.A.; Khaldjigitov, A.A. Mathematical and numerical modeling of the coupled dynamic thermoelastic problems for isotropic bodies. TWMS J. Pure Appl. Math. 2020, 11, 119–126. [Google Scholar]
  20. Hale, J.K. Ordinary Differential Equations; Wiley-Interscience: New York, NY, USA, 1969. [Google Scholar]
  21. Roeger, L.W.; Barnard, R.W. Preservation of local dynamics when applying central difference methods: Application to SIR model. J. Differ. Equations Appl. 2007, 13, 333. [Google Scholar] [CrossRef]
  22. Khalsaraei, M.M.; Shokri, A.; Ramos, H.; Heydari, S. A positive and elementary stable nonstandard explicit scheme for a mathematical model of the influenza disease. Math. Comput. Simul. 2021, 182, 397–410. [Google Scholar] [CrossRef]
  23. Duncan, C.J.; Duncan, S.R.; Scott, S. Whooping cough epidemic in London, 1701–1812: Infection dynamics seasonal forcing and the effects of malnutrition. Proc. R. Soc. Lond. B 1996, 263, 445–450. [Google Scholar]
  24. Hethcote, H.W. The Mathematics of Infectious Diseases. SIAM Rev. 2000, 42, 599–653. [Google Scholar] [CrossRef]
  25. Mickens, R.E.; Jordan, P.M. A positivity-preserving nonstandard finite difference scheme for the Damped Wave Equation. Numer. Meth. Partial Diff. Eq. 2004, 20, 639–649. [Google Scholar] [CrossRef]
  26. Bacaer, N.; Ouifki, R.; Pretorius, C.; Wood, R.; Williams, B. Modeling the joint epidemics of TB and HIV in a south African township. J. Math. Biol. 2008, 57, 557–593. [Google Scholar] [CrossRef] [PubMed]
  27. Brauer, F.; Castillo-Chavez, C. Mathematical Models in Population Biology and Epidemiology; Texts in Applied Mathematics; Springer: New York, NY, USA, 2012; Volume 40. [Google Scholar] [CrossRef]
  28. Arenasa, A.J.; González-Parrab, G.; Chen-Charpentier, B.M. A nonstandard numerical scheme of predictor—Corrector type for epidemic models. Comput. Math. Appl. 2010, 59, 3740–3749. [Google Scholar] [CrossRef] [Green Version]
Figure 1. The numerical results of system (2) using P-C schemes with h = 0.01 , ode45 and RK4 scheme with h = 0.005 taking parameter values ( R 0 > 1 ) and I.C. S 0 = 0.9 , I 0 = 0.05 , ν = 24 , μ = 0.04 , β = 123 , R 0 = 0.05 .
Figure 1. The numerical results of system (2) using P-C schemes with h = 0.01 , ode45 and RK4 scheme with h = 0.005 taking parameter values ( R 0 > 1 ) and I.C. S 0 = 0.9 , I 0 = 0.05 , ν = 24 , μ = 0.04 , β = 123 , R 0 = 0.05 .
Mathematics 10 03299 g001
Figure 2. Absolute errors for the system (2) with h= 0.01 by the new schemes, the proposed method in [28] and I.C. ( S 0 , I 0 , R 0 ) = ( 0.9 , 0.05 , 0.05 ) with β = 123 , μ = 0.04 , ν = 24 , ( R 0 > 1 ) using Runge–Kutta as a reference solution.
Figure 2. Absolute errors for the system (2) with h= 0.01 by the new schemes, the proposed method in [28] and I.C. ( S 0 , I 0 , R 0 ) = ( 0.9 , 0.05 , 0.05 ) with β = 123 , μ = 0.04 , ν = 24 , ( R 0 > 1 ) using Runge–Kutta as a reference solution.
Mathematics 10 03299 g002aMathematics 10 03299 g002b
Figure 3. The numerical results of system (2) using the P-C schemes, NSFD schemes with h = 0.01 and RK4 with h = 0.005 following parameter values ( R 0 > 1 ) and I.C. S 0 = 0.24 , I 0 = 0.007 , R 0 = 0.753 , β = 123 , μ = 0.04 , ν = 24 .
Figure 3. The numerical results of system (2) using the P-C schemes, NSFD schemes with h = 0.01 and RK4 with h = 0.005 following parameter values ( R 0 > 1 ) and I.C. S 0 = 0.24 , I 0 = 0.007 , R 0 = 0.753 , β = 123 , μ = 0.04 , ν = 24 .
Mathematics 10 03299 g003aMathematics 10 03299 g003b
Figure 4. The numerical results of system (2) using the P-C 1 schemes and P-C 2 schemes with h = 0.01 following parameter values ( R 0 > 1 ) and I.C. S 0 = 0.9 , I 0 = 0.05 , R 0 = 0.05 , β = 123 , μ = 0.04 , ν = 24 .
Figure 4. The numerical results of system (2) using the P-C 1 schemes and P-C 2 schemes with h = 0.01 following parameter values ( R 0 > 1 ) and I.C. S 0 = 0.9 , I 0 = 0.05 , R 0 = 0.05 , β = 123 , μ = 0.04 , ν = 24 .
Mathematics 10 03299 g004
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mehdizadeh Khalsaraei, M.; Shokri, A.; Ramos, H.; Yao, S.-W.; Molayi, M. Efficient Numerical Solutions to a SIR Epidemic Model. Mathematics 2022, 10, 3299. https://doi.org/10.3390/math10183299

AMA Style

Mehdizadeh Khalsaraei M, Shokri A, Ramos H, Yao S-W, Molayi M. Efficient Numerical Solutions to a SIR Epidemic Model. Mathematics. 2022; 10(18):3299. https://doi.org/10.3390/math10183299

Chicago/Turabian Style

Mehdizadeh Khalsaraei, Mohammad, Ali Shokri, Higinio Ramos, Shao-Wen Yao, and Maryam Molayi. 2022. "Efficient Numerical Solutions to a SIR Epidemic Model" Mathematics 10, no. 18: 3299. https://doi.org/10.3390/math10183299

APA Style

Mehdizadeh Khalsaraei, M., Shokri, A., Ramos, H., Yao, S. -W., & Molayi, M. (2022). Efficient Numerical Solutions to a SIR Epidemic Model. Mathematics, 10(18), 3299. https://doi.org/10.3390/math10183299

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop