# An Accurate Predictor-Corrector-Type Nonstandard Finite Difference Scheme for an SEIR Epidemic Model

Authors: Asma Farooqi; Riaz Ahmad; Rashada Farooqi; Sayer O. Alharbi; Dumitru Baleanu; Muhammad Rafiq; Ilyas Khan; M. O. Ahmad
Journal: Journal of Mathematics
 (2020)

---

## Abstract

The present work deals with the construction, development, and analysis of a viable normalized predictor-corrector-type nonstandard finite difference scheme for the SEIR model concerning the transmission dynamics of measles. The proposed numerical scheme double refines the solution and gives realistic results even for large step sizes, thus making it economical when integrating over long time periods. Moreover, it is dynamically consistent with a continuous system and unconditionally convergent and preserves the positive behavior of the state variables involved in the system. Simulations are performed to guarantee the results, and its effectiveness is compared with well-known numerical methods such as Runge–Kutta (RK) and Euler method of a predictor-corrector type.


---

## Body

## 1. Introduction

Mathematicians and biologists have been working for a long time on the biological process of life science. They succeeded in evaluating some remarkable results from their work. An important mark in mathematical biology is the mathematical modeling of infectious diseases [1]. Measles is one of such a highly infectious childhood disease, caused by respiratory infection by a Morbilli virus, measles virus. Arthur Ransom first observed the irregular cyclic behavior of measles that is considered as a mainly conspicuous facet of measles. The age structure of the population, contact, immigration rate, and the school seasons were known as a crucial phase for the swell of measles [2–4]. William Hammer in 1906 published a discrete numerical model for the transmission of measles epidemic. Later assumption of “Mass Action” is applied to that model which is the basic rule to the current theory of deterministic modeling of infectious diseases [5].Society has a keen concern in knowing the major evolution for the spread of diseases. Analytical results give a solution to these problems but for limited cases and causes many problems. The homotopy perturbation method and variational iteration method can be used for the solution of the epidemic models [6]. However, the first choice to solve these laws of nature is the numerical method based on a difference scheme for good approximations [7–9]. In general, already developed numerical schemes such as Euler, Runge–Kutta, and others at times stop working by generating nonphysical results. These unnecessary oscillations contrived chaos and false fixed points [10]. Moreover, some methods are unsuccessful if we check them on larger step sizes [11]. To avoid such discrepancies, the numerical schemes based on the “nonstandard finite difference method (NSFD)” are established. These techniques were first developed by R. E. Mickens [7, 12, 13]. The created numerical schemes preserve the essential properties such as dynamical consistency, stability, and equilibrium points [14–21]. Researchers have developed competitive NSFD schemes for epidemic diseases. Many of these NSFD schemes are consistent for small step sizes with the continuous model, but for large step sizes, the unwanted oscillations have been observed. Piyanwong, Jansen, and Twizel have constructed a positive and unconditionally stable scheme for SIR and SEIR models, respectively [22, 23]. Nevertheless, the lack of application of conservation law in their developed schemes explicitly caused impracticable and unrealistic solutions, while Abraham and Gilberto have developed NSFD schemes of the SIR epidemic model to obtain the physically realistic solutions for all step sizes, where they apply the conservation law in addition to nonlocal approximation [24].In this paper, we have developed a normalized NSFDM of predictor-corrector- (PC-) type inspired by the previous work discussed to double refine the numerical solution of a nonlinear dynamics regarding the transmission of measles. To keep the method explicit, we will use the forward difference approximations for the first derivative terms. The nonlocal approximations are used to tackle the nonlinear terms withφh as a nonclassical denominator function. By using this idea, the measles model will converge to equilibrium points, even for the larger step sizes.


## 2. The Mathematical Frame of Work

In this work, the dynamics of the measles epidemic described by the SEIR mathematical model suggested by Jansen and E. H. Twizel is considered [23]. In the SEIR model, the total human population is categorized as susceptible, exposed, infectious, and recovered subpopulations denoted by S, E, I, and R, respectively. Consider the flow of the SEIR model for the measles epidemic, as shown in Figure 1.HereS= susceptible individualsE=exposed individualsI= infected individualR=recovered individualμ= birth rate and death rateβ= the rate at which susceptible individuals are infected by those who are infectiousσ= the rate at which exposed individuals become infectedγ= the rate at which infected individuals recover.Figure 1
Flowchart of the SEIR Model for measles epidemic.Hereμ,β,σ,andγ are considered as positive parameters. Furthermore, we assume that E0,I0≠0. Consider N is the constant size of the population so that the number of recovered individuals R=Rt at time t is defined by R=N−S−E−I. A susceptible is a move to the exposed model, where the individual is infected but yet not infectious. After sometimes the individual becomes infectious and enters into the infected compartment, and in this way, the disease spreads into the population.The mathematical model is written as(1)dSdt=μN−μS−βIS,t>0,S0=S0,dEdt=βIS−μE−σE,t>0,E0=E0,dIdt=σE−μI−γI,t>0,I0=I0,dRdt=γI−μR,t>0,R0=R0,where(2)S+E+I+R=N.Equation (2) shows that the total size of the population remains constant, which is connected with the continuous system (1). The following two points give the equilibrium points of (1):The disease-free equilibrium (DFE) pointN,0,0The endemic equilibrium (EE) pointN/R0,μN/μ+σ1−1/R0,μ/βR0−1Here,R0=σβN/μ+σμ+γ is the basic reproductive number associated with the measles model.For the sake of brevity, we are not mentioning RK-4 and Euler-PC scheme.


## 3. Numerical Modeling

To the construction of NSFD scheme, the continuous system (1) is discretized using the forward difference approximation for the first-order time derivatives. Thus, if ft is differentiable, the f′t can be approximated by(3)dftdt=ft+h−ftφh+Ohash⟶0,where φh is a real-valued function satisfying the condition φh⟶0 as h⟶0. We have(4)φh=1−e−h.Thus, the NSFD scheme for system (1) takes the form(5)Sn+1−Snφh=μN−μSn+1−βSn+1In,En+1−Enφh=βSn+1In−μEn+1−σEn+1,In+1−Inφh=σEn+1−μIn+1−γIn+1,Rn+1−Rnφh=γIn+1−μRn+1.From equation (5), we have(6)1+φhμSn+1+En+1+In+1+Rn+1=φhμN−φhμSn+En+In+Rn.Thus, ifSn+En+In+Rn=N for all n≥0, then Sn+1+En+1+In+1+Rn+1 = N for all n≥0.Thus, the conservation law property is satisfied. System (5) is written as(7)Sn+1=Sn+μφhN1+φhμ+βIn,En+1=En+φhβSn+1In1+φhμ+σ,In+1=In+φhσEn+11+φhμ+γ.Now, the positivity ofSn+1,En+1,In+1,andRn+1 is guaranteed if 0<Sn<1, 0<En<1,0<In<1, and 0<Rn<1 for all n≥0. Thus, the constructed scheme has the following properties:The conservation law is satisfiedThe positivity and boundedness of the solution is satisfied: for system (7), we have that if 0<Sn<1, 0<En<1, 0<In<1,and0<Rn<1, then 0<Sn+1<1, 0<En+1<1, 0<In+1<1,and0<Rn+1<1 for all n≥0.


## 4. NSFD Predictor-Corrector Scheme

This section is improved by the approach of a predictor-corrector-type NSFD scheme (7) to obtain the benefits of both methods. For the development of this scheme, firstly, system (7) is taken as a predictor scheme, i.e.,(8)Spn+1=Sn+μφhN1+φhμ+βIn,Epn+1=En+φhβSpn+1In1+φhμ+σ,Ipn+1=In+φhσEpn+11+φhμ+γ.Now, we evaluate system (1) at time t+φh and introduce the term ϵ−1Sn+1/φh (where ϵ is just like the accelerating factor and its range is 0<ϵ<1). Thus, the expression will be NSFDCL scheme that preserves conservation law represented as(9)Sn+1−Snφh=μN−μSn+1−βSn+1In+1−ϵ−1Sn+1φh+ϵ−1Sn+1φh,En+1−Enφh=βSn+1In−μ+σEn+1,In+1−Inφh=σEn+1−μ+γIn,Rn+1−Rnφh=γIn+1−μRn+1.Thus, the corrector scheme is obtained as(10)Scn+1=Sn+μφhN+ϵ−1Spn+11+ϵ−1+μφh+βφhIpn+1,Ecn+1=En+φhβScn+1Ipn+11+φhμ+σ,Icn+1=In+φhσEcn+11+φhμ+γ.where Rn+1=N−Scn+1−Ecn+1−Icn+1.


## 5. Convergence Analysis

In this section, the unconditional convergence of the numerical solution is presented by the proposed method. Taking the values as(11)F1S,E,I=S+μNϕh1+ϕhμ+βI,F2S,E,I=E+βϕhIF11+ϕhμ+σ,F3S,E,I=I+σϕhF21+ϕhμ+γ,G1S,E,I=S+μNϕh+ϵ−1F11+ϵ−1+μϕh+βϕlF3,G2S,E,I=E+βϕhG1F31+ϕhμ+σ,G3S,E,I=I+σϕhG21+ϕhμ+γ,where F1,F2,andF3 are given as in equation (8).The convergence and stability analysis of scheme (10) is carried out by calculating the eigenvalues of the Jacobian of the linearized scheme and studying its behavior and evaluated at the fixed points. If (S∗,E∗,I∗) is the fixed point of system (1), then the Jacobian matrix JG is given by(12)JGS∗,E∗,I∗=∂G1S∗,E∗,I∗∂S∂G1S∗,E∗,I∗∂E∂G1S∗,E∗,I∗∂I∂G2S∗,E∗,I∗∂S∂G2S∗,E∗,I∗∂E∂G2S∗,E∗,I∗∂I∂G3S∗,E∗,I∗∂S∂G3S∗,E∗,I∗∂E∂G3S∗,E∗,I∗∂I,where(13)∂F1∂S=11+μϕh+βϕhI∂F1∂E,∂F1∂I=−S+μNϕh.βϕh1+μϕh+βϕhI2,∂F2∂S=βϕhI∂F1/∂S1+ϕhμ+σ,∂F2∂E=1+βϕhI∂F1/∂E1+ϕhμ+σ,∂F2∂I=βϕhF1+I∂F1/∂I1+ϕhμ+σ,∂F3∂E=σϕh∂F2/∂E1+ϕhμ+γ,∂F3∂I=1+σϕh∂F2/∂S1+ϕhμ+γ,∂G1∂S=1+ϵ−1∂F1/∂S1+ϵ−1+μϕh+βϕhF3−S+μNϕh+ϵ−1F1.βϕh∂F3/∂S1+ϵ−1+μϕh+βϕlF32,∂G1∂E=ϵ−1∂F1/∂E1+ϵ−1+μϕh+βϕhF3−S+μNϕh+ϵ−1F1.βϕh∂F3/∂E1+ϵ−1+μϕh+βϕhF32,∂G1∂I=ϵ−1∂F1/∂I1+ϵ−1+μϕh+βϕhF3−S+μNϕh+ϵ−1F1.βϕh∂F3/∂I1+ϵ−1+μϕh+βϕhF32,∂G2∂S=βϕh1+ϕhμ+σ∂G1∂S.F3+G1∂F3∂S,∂G2∂E=1+βϕh∂G1/∂E.F3+G1∂F3/∂E1+ϕhμ+σ,∂G2∂I=βϕh1+ϕhμ+σ∂G1∂I.F3+G1∂F3∂I,∂G3∂S=σϕh∂G2/∂S1+ϕhμ+γ,∂G3∂E=σϕh∂G2/∂E1+ϕhμ+γ,∂G3∂I=1+σϕh∂G2/∂I1+ϕhμ+γ.


### 5.1. Disease-Free EquilibriumSo,Eo,Io=N,0,0

The values of functions and its derivatives at DFE pointN,0,0 are as(14)F1N,0,0=N,F2N,0,0=0,F3N,0,0=0,G1N,0,0=N,G2N,0,0=0,G3N,0,0=0,∂F1N,0,0∂S=11+μϕh,∂F1N,0,0∂E=0,∂F1N,0,0∂I=−Nβϕh1+μϕh,∂F2N,0,0∂S=0,∂F2N,0,0∂E=11+ϕhμ+σ,∂F2N,0,0∂I=βϕh1+ϕhμ+σ,∂F3N,0,0∂S=0,∂F3N,0,0∂E=σϕh1+ϕhμ+σ.1+ϕhμ+γ,∂F3N,0,0∂I=1+σϕhβNϕh/1+ϕhμ+γ1+ϕhμ+γ,∂G1N,0,0∂S=11+μϕh,∂G1N,0,0∂E=−Nσϕhβϕh1+μϕh+ϵ−11+ϕhμ+σ1+ϕhμ+γ,∂G1N,0,0∂I=−Nβϕh1+ϵ−1+μϕhϵ−11+μϕh+11+ϕhμ+γ+σNβϕ2h1+ϕhμ+σ1+ϕhμ+γ,∂G2N,0,0∂S=0,∂G2N,0,0∂E=11+ϕhμ+σ1+Nβσϕ2h1+ϕhμ+σ1+ϕhμ+γ,∂G2N,0,0∂I=βϕh1+ϕhμ+σ1+ϕhμ+σ+Nσβϕ2h1+ϕhμ+σ1+ϕhμ+γ∂G3N,0,0∂S=0,∂G3N,0,0∂E=σϕh1+ϕhμ+σ1+ϕhμ+γ1+σβNϕ2h1+ϕhμ+σ1+ϕhμ+γ∂G3N,0,0∂I=11+ϕhμ+σ+Nσβϕ2h1+ϕhμ+σ1+ϕhμ+γ2+N2β2σ2ϕ4h1+ϕhμ+σ21+ϕhμ+σ2.For ease in the calculation, we use the following conventions:(15)α=μNϕh>0,η1=μϕh>0,η2=βϕh>0,η3=σϕh>0,θ=1+ϵ−1+μϕh=1+ϵ−1+η1>1,δ1=1+ϕhμ+σ>1,δ2=1+ϕhμ+γ>1.It is clear thatθ>η1,1+η1>δ1,1+η1>δ2,andδ2=θ−ϵ−1+γϕh. If R0=σβN/μ+σμ+γ<1, the Jacobian calculated at disease-free points S0∗,E0∗,I0∗=N,0,0 is given as(16)JGN,0,0=∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I.Since the partial derivatives using the conventions are(17)∂G1N,0,0∂S=11+η1,∂G1N,0,0∂E=−Nη2η3θδ1δ2,∂G1N,0,0∂I=−Nη2θϵ−11+η1+1δ2+Nη2η3δ1δ2,∂G2N,0,0∂S=0,∂G2N,0,0∂E=1δ1+Nη2η3δ12δ2,∂G2N,0,0∂I=Nη2δ1δ2+N2η22η3δ12δ2,∂G3N,0,0∂S=0,∂G3N,0,0∂E=η3δ1δ2+Nη2η32δ12δ22,∂G3N,0,0∂I=1δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.So,(18)JGN,0,0=11+η1−Nη2η3θδ1δ2−Nη2θϵ−11+η1+1δ2+Nη2η3δ1δ201δ1+Nη2η3δ12δ2Nη2δ1δ2+N2η22η3δ12δ20η3δ1δ2+Nη2η32δ12δ221δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.To calculate the eigenvalues, putdetJ−λI=0,λ1=1/1+η1=1/1+μϕh<1. And(19)λ2−λδ1δ2+Nη2η3δ12δ2+δ12δ2+δ1Nη2η3+N2η22η32δ12δ22+δ1δ2+Nη2η3δ12δ2δ12δ2+δ1Nη2η3+N2η22η32δ12δ22−Nη2δ1+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ22=0,λ2−Aλ+B=0,where(20)A=δ1δ2+Nη2η3/δ12δ2+δ12δ2+δ1Nη2η3+N2η22η32/δ12δ22=traceJ∗,B=δ1δ2+Nη2η3δ12δ2δ12δ2+δ1Nη2η3+N2η22η32δ12δ22−Nη2δ1+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ22=detJ∗,J∗=1δ1+Nη2η3δ12δ2Nη2δ1δ2+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ221δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.To calculate the eigenvalues ofJ∗, the following lemma is proved.Lemma 1.
For the quadratic equation,λ2−λA+B=0, both roots satisfy λi<1,i=2,3 if and only if the following conditions are satisfied [25]:(21)11−A+B>0,21+A+B>0,3B<1.Let us defineA=TraceJ∗,B=DetJ∗. Therefore,(22)A=δ1+δ2Nη2η3+δ1δ2+N2η22η32δ12δ22>0,B=1δ1δ2+Nη2η3δ12δ22<1,f−1=1+A+B>0,f0=B=1δ1δ2+Nη2η3δ12δ22<1δ1δ2+Nη2η3δ1δ2=1+Nη2η3δ1δ2=1+R0μ+σμ+γϕ2h1+ϕ2hμ+σμ+γ+ϕhμ+σμ+γ<1sinceR0<1,⇒B<1.Now, for(23)f1=1−A+B=δ12δ22+δ1+δ2Nη2η3+δ1δ2+N2η22η32+Nη2η3+δ1δ2δ12δ22>0.As all the conditions of the lemma hold forR0<1, therefore the absolute value of both eigenvalues of J∗ is less than one, whenever R0<1. Thus, the numerical scheme in equations (8) and (10) will converge unconditionally to disease-free equilibrium N,0,0 for any value of time step h whenever R0<1 which is also computationally verified in Figure 2(a).Figure 2
Spectral radius of the Jacobian at equilibrium points.
(a)(b)


### 5.2. Endemic EquilibriumSe,Ee,Ie=N/R0,μN/μ+σ1−1/R0,μ/βR0−1

The value of functions and its derivatives at an endemic point is as(24)F1Se,Ee,Ie=NR0,F2Se,Ee,Ie=μN1−1/R0μ+σ,F3Se,Ee,Ie=μβR0−1,G1Se,Ee,Ie=NR0,G2Se,Ee,Ie=μβR0−1,G3Se,Ee,Ie=μN1−1/R0μ+σ,(25)∂F1Se,Ee,Ie∂S=11+μϕhR0,∂F1Se,Ee,Ie∂E=0,∂F1Se,Ee,Ie∂I=−Nβϕh1+μϕhR0,∂F2Se,Ee,Ie∂S=μϕhR0−11+μϕhR01+ϕhμ+σ,∂F2Se,Ee,Ie∂E=11+ϕhμ+σ,∂F2Se,Ee,Ie∂I=−NβϕhμϕlR02−2μϕhR0−1R01+ϕhμ+σ1+ϕhμ+γ,∂F3Se,Ee,Ie∂S=μσϕ2hR0−11+ϕhμ+σ1+ϕhμ+γ1+μϕhR0,∂F3Se,Ee,Ie∂E=σϕh1+ϕhμ+σ1+ϕhμ+γ,∂F3Se,Ee,Ie∂I=R01+ϕhμ+σ1+μϕhR0−σNβϕ2hμϕhR02−2μϕh−1R01+ϕhμ+σ1+μϕhR01+ϕhμ+γ,(26)∂G1Se,Ee,Ie∂S=11+μϕhR0−Nμβσϕ3hR0−1R01+ϵ−1+μϕhR01+ϕhμ+σ1+μϕhR01+ϕhμ+γ,∂G1Se,Ee,Ie∂E=−Nβσϕ2hR01+ϵ−1+μϕhR01+ϕhμ+σ1+ϕhμ+γ,∂G1Se,Ee,Ie∂I=−Nβϕh1+ϵ−1+μϕhR01+μϕhR0ϵ−1+R01+ϕhμ+σ1+μϕhR0−σNβϕ2hμϕhR02−2μϕhR0−1R021+ϕhμ+σ1+μϕhR01+ϕhμ+γ,∂G2Se,Ee,Ie∂S=μβϕhR0−1β1+ϕhμ+σ1+μϕhR0−Nμβ2σϕ4hμR0−12βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+Nμβσϕ3hR0−1R01+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G2Se,Ee,Ie∂E=11+ϕhμ+σ+Nβσϕ2hR01+ϕhμ+σ21+ϕhμ+γ−Nμβ2σϕ3hR0−1βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR0,∂G2Se,Ee,Ie∂I=−Nμβ2ϕ2hϵ−1R0−1β1+ϕhμ+σ1+μϕhR01+ϵ−1+μϕhR0+Nμβ2ϕ2hR0−1R01+ϕhμ+σ1+μϕhR0−Nβσϕ2hμϕhR02−2μϕhR0−1βR021+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+NβϕhR01+ϕhμ+σ1+ϕhμ+γ−N2β2σϕ3hμϕhR02−2μϕhR0−1R021+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G3Se,Ee,Ie∂S=σϕh1+ϕhμ+γμβϕhR0−1β1+ϕhμ+σ1+μϕhR0−Nμβ2σϕ4hμR0−12βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+Nμβσϕ3hR0−1R01+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G3Se,Ee,Ie∂E=σϕh1+ϕhμ+γ11+ϕhμ+σ−Nμβ2σϕ3hR0−1βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR0+Nβσϕ2hR01+ϕhμ+σ21+ϕhμ+γ,∂G3Se,Ee,Ie∂I=σϕh1+ϕhμ+γ−Nμβ2ϕ2hϵ−1R0−1β1+ϕhμ+σ1+μϕhR01+ϵ−1+μϕhR0+Nμβ2ϕ2hR0−1R01+ϕhμ+σ1+μϕhR0−Nβσϕ2hμϕhR02−2μϕhR0−1βR021+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+NβϕhR01+ϕhμ+σ1+ϕhμ+γ−N2β2σϕ3hμϕhR02−2μϕhR0−1R021+ϕhμ+σ21+ϕhμ+γ1+μϕhR0.The Jacobian evaluated at an endemic pointSe∗,Ee∗,Ie∗=Se,Ee,Ie=N/R0,μN/μ+σ1−1/R0,μ/βR0−1 is evaluated as(27)JGSe∗,Ee∗,Ie∗=∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I.It is tedious to calculate the eigenvalues of the Jacobian of endemic equilibrium analytically, so we have plotted the largest eigenvalue against each step size, and Figure2(b) shows that, for all step sizes, the spectral radius of the Jacobian of EE remains less than one, if R0>1, which implies that the numerical scheme in equations (8) and (10) is unconditionally convergent if R0>1, for all step sizes.


## 5.1. Disease-Free EquilibriumSo,Eo,Io=N,0,0

The values of functions and its derivatives at DFE pointN,0,0 are as(14)F1N,0,0=N,F2N,0,0=0,F3N,0,0=0,G1N,0,0=N,G2N,0,0=0,G3N,0,0=0,∂F1N,0,0∂S=11+μϕh,∂F1N,0,0∂E=0,∂F1N,0,0∂I=−Nβϕh1+μϕh,∂F2N,0,0∂S=0,∂F2N,0,0∂E=11+ϕhμ+σ,∂F2N,0,0∂I=βϕh1+ϕhμ+σ,∂F3N,0,0∂S=0,∂F3N,0,0∂E=σϕh1+ϕhμ+σ.1+ϕhμ+γ,∂F3N,0,0∂I=1+σϕhβNϕh/1+ϕhμ+γ1+ϕhμ+γ,∂G1N,0,0∂S=11+μϕh,∂G1N,0,0∂E=−Nσϕhβϕh1+μϕh+ϵ−11+ϕhμ+σ1+ϕhμ+γ,∂G1N,0,0∂I=−Nβϕh1+ϵ−1+μϕhϵ−11+μϕh+11+ϕhμ+γ+σNβϕ2h1+ϕhμ+σ1+ϕhμ+γ,∂G2N,0,0∂S=0,∂G2N,0,0∂E=11+ϕhμ+σ1+Nβσϕ2h1+ϕhμ+σ1+ϕhμ+γ,∂G2N,0,0∂I=βϕh1+ϕhμ+σ1+ϕhμ+σ+Nσβϕ2h1+ϕhμ+σ1+ϕhμ+γ∂G3N,0,0∂S=0,∂G3N,0,0∂E=σϕh1+ϕhμ+σ1+ϕhμ+γ1+σβNϕ2h1+ϕhμ+σ1+ϕhμ+γ∂G3N,0,0∂I=11+ϕhμ+σ+Nσβϕ2h1+ϕhμ+σ1+ϕhμ+γ2+N2β2σ2ϕ4h1+ϕhμ+σ21+ϕhμ+σ2.For ease in the calculation, we use the following conventions:(15)α=μNϕh>0,η1=μϕh>0,η2=βϕh>0,η3=σϕh>0,θ=1+ϵ−1+μϕh=1+ϵ−1+η1>1,δ1=1+ϕhμ+σ>1,δ2=1+ϕhμ+γ>1.It is clear thatθ>η1,1+η1>δ1,1+η1>δ2,andδ2=θ−ϵ−1+γϕh. If R0=σβN/μ+σμ+γ<1, the Jacobian calculated at disease-free points S0∗,E0∗,I0∗=N,0,0 is given as(16)JGN,0,0=∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I∂G1N,0,0∂S∂G1N,0,0∂E∂G1N,0,0∂I.Since the partial derivatives using the conventions are(17)∂G1N,0,0∂S=11+η1,∂G1N,0,0∂E=−Nη2η3θδ1δ2,∂G1N,0,0∂I=−Nη2θϵ−11+η1+1δ2+Nη2η3δ1δ2,∂G2N,0,0∂S=0,∂G2N,0,0∂E=1δ1+Nη2η3δ12δ2,∂G2N,0,0∂I=Nη2δ1δ2+N2η22η3δ12δ2,∂G3N,0,0∂S=0,∂G3N,0,0∂E=η3δ1δ2+Nη2η32δ12δ22,∂G3N,0,0∂I=1δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.So,(18)JGN,0,0=11+η1−Nη2η3θδ1δ2−Nη2θϵ−11+η1+1δ2+Nη2η3δ1δ201δ1+Nη2η3δ12δ2Nη2δ1δ2+N2η22η3δ12δ20η3δ1δ2+Nη2η32δ12δ221δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.To calculate the eigenvalues, putdetJ−λI=0,λ1=1/1+η1=1/1+μϕh<1. And(19)λ2−λδ1δ2+Nη2η3δ12δ2+δ12δ2+δ1Nη2η3+N2η22η32δ12δ22+δ1δ2+Nη2η3δ12δ2δ12δ2+δ1Nη2η3+N2η22η32δ12δ22−Nη2δ1+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ22=0,λ2−Aλ+B=0,where(20)A=δ1δ2+Nη2η3/δ12δ2+δ12δ2+δ1Nη2η3+N2η22η32/δ12δ22=traceJ∗,B=δ1δ2+Nη2η3δ12δ2δ12δ2+δ1Nη2η3+N2η22η32δ12δ22−Nη2δ1+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ22=detJ∗,J∗=1δ1+Nη2η3δ12δ2Nη2δ1δ2+N2η22η3δ12δ2η3δ1δ2+Nη2η32δ12δ221δ2+Nη2η3δ1δ22+N2η22η32δ12δ22.To calculate the eigenvalues ofJ∗, the following lemma is proved.Lemma 1.
For the quadratic equation,λ2−λA+B=0, both roots satisfy λi<1,i=2,3 if and only if the following conditions are satisfied [25]:(21)11−A+B>0,21+A+B>0,3B<1.Let us defineA=TraceJ∗,B=DetJ∗. Therefore,(22)A=δ1+δ2Nη2η3+δ1δ2+N2η22η32δ12δ22>0,B=1δ1δ2+Nη2η3δ12δ22<1,f−1=1+A+B>0,f0=B=1δ1δ2+Nη2η3δ12δ22<1δ1δ2+Nη2η3δ1δ2=1+Nη2η3δ1δ2=1+R0μ+σμ+γϕ2h1+ϕ2hμ+σμ+γ+ϕhμ+σμ+γ<1sinceR0<1,⇒B<1.Now, for(23)f1=1−A+B=δ12δ22+δ1+δ2Nη2η3+δ1δ2+N2η22η32+Nη2η3+δ1δ2δ12δ22>0.As all the conditions of the lemma hold forR0<1, therefore the absolute value of both eigenvalues of J∗ is less than one, whenever R0<1. Thus, the numerical scheme in equations (8) and (10) will converge unconditionally to disease-free equilibrium N,0,0 for any value of time step h whenever R0<1 which is also computationally verified in Figure 2(a).Figure 2
Spectral radius of the Jacobian at equilibrium points.
(a)(b)


## 5.2. Endemic EquilibriumSe,Ee,Ie=N/R0,μN/μ+σ1−1/R0,μ/βR0−1

The value of functions and its derivatives at an endemic point is as(24)F1Se,Ee,Ie=NR0,F2Se,Ee,Ie=μN1−1/R0μ+σ,F3Se,Ee,Ie=μβR0−1,G1Se,Ee,Ie=NR0,G2Se,Ee,Ie=μβR0−1,G3Se,Ee,Ie=μN1−1/R0μ+σ,(25)∂F1Se,Ee,Ie∂S=11+μϕhR0,∂F1Se,Ee,Ie∂E=0,∂F1Se,Ee,Ie∂I=−Nβϕh1+μϕhR0,∂F2Se,Ee,Ie∂S=μϕhR0−11+μϕhR01+ϕhμ+σ,∂F2Se,Ee,Ie∂E=11+ϕhμ+σ,∂F2Se,Ee,Ie∂I=−NβϕhμϕlR02−2μϕhR0−1R01+ϕhμ+σ1+ϕhμ+γ,∂F3Se,Ee,Ie∂S=μσϕ2hR0−11+ϕhμ+σ1+ϕhμ+γ1+μϕhR0,∂F3Se,Ee,Ie∂E=σϕh1+ϕhμ+σ1+ϕhμ+γ,∂F3Se,Ee,Ie∂I=R01+ϕhμ+σ1+μϕhR0−σNβϕ2hμϕhR02−2μϕh−1R01+ϕhμ+σ1+μϕhR01+ϕhμ+γ,(26)∂G1Se,Ee,Ie∂S=11+μϕhR0−Nμβσϕ3hR0−1R01+ϵ−1+μϕhR01+ϕhμ+σ1+μϕhR01+ϕhμ+γ,∂G1Se,Ee,Ie∂E=−Nβσϕ2hR01+ϵ−1+μϕhR01+ϕhμ+σ1+ϕhμ+γ,∂G1Se,Ee,Ie∂I=−Nβϕh1+ϵ−1+μϕhR01+μϕhR0ϵ−1+R01+ϕhμ+σ1+μϕhR0−σNβϕ2hμϕhR02−2μϕhR0−1R021+ϕhμ+σ1+μϕhR01+ϕhμ+γ,∂G2Se,Ee,Ie∂S=μβϕhR0−1β1+ϕhμ+σ1+μϕhR0−Nμβ2σϕ4hμR0−12βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+Nμβσϕ3hR0−1R01+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G2Se,Ee,Ie∂E=11+ϕhμ+σ+Nβσϕ2hR01+ϕhμ+σ21+ϕhμ+γ−Nμβ2σϕ3hR0−1βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR0,∂G2Se,Ee,Ie∂I=−Nμβ2ϕ2hϵ−1R0−1β1+ϕhμ+σ1+μϕhR01+ϵ−1+μϕhR0+Nμβ2ϕ2hR0−1R01+ϕhμ+σ1+μϕhR0−Nβσϕ2hμϕhR02−2μϕhR0−1βR021+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+NβϕhR01+ϕhμ+σ1+ϕhμ+γ−N2β2σϕ3hμϕhR02−2μϕhR0−1R021+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G3Se,Ee,Ie∂S=σϕh1+ϕhμ+γμβϕhR0−1β1+ϕhμ+σ1+μϕhR0−Nμβ2σϕ4hμR0−12βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+Nμβσϕ3hR0−1R01+ϕhμ+σ21+ϕhμ+γ1+μϕhR0,∂G3Se,Ee,Ie∂E=σϕh1+ϕhμ+γ11+ϕhμ+σ−Nμβ2σϕ3hR0−1βR01+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR0+Nβσϕ2hR01+ϕhμ+σ21+ϕhμ+γ,∂G3Se,Ee,Ie∂I=σϕh1+ϕhμ+γ−Nμβ2ϕ2hϵ−1R0−1β1+ϕhμ+σ1+μϕhR01+ϵ−1+μϕhR0+Nμβ2ϕ2hR0−1R01+ϕhμ+σ1+μϕhR0−Nβσϕ2hμϕhR02−2μϕhR0−1βR021+ϕhμ+σ21+ϕhμ+γ1+ϵ−1+μϕhR01+μϕhR0+NβϕhR01+ϕhμ+σ1+ϕhμ+γ−N2β2σϕ3hμϕhR02−2μϕhR0−1R021+ϕhμ+σ21+ϕhμ+γ1+μϕhR0.The Jacobian evaluated at an endemic pointSe∗,Ee∗,Ie∗=Se,Ee,Ie=N/R0,μN/μ+σ1−1/R0,μ/βR0−1 is evaluated as(27)JGSe∗,Ee∗,Ie∗=∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I∂G1Se∗,Ee∗,Ie∗∂S∂G1Se∗,Ee∗,Ie∗∂E∂G1Se∗,Ee∗,Ie∗∂I.It is tedious to calculate the eigenvalues of the Jacobian of endemic equilibrium analytically, so we have plotted the largest eigenvalue against each step size, and Figure2(b) shows that, for all step sizes, the spectral radius of the Jacobian of EE remains less than one, if R0>1, which implies that the numerical scheme in equations (8) and (10) is unconditionally convergent if R0>1, for all step sizes.


## 6. Numerical Results and Discussion

All the methods showed the convergence for small step sizes in Table1. However, for large values of step sizes, only the normalized NSFDPC converge to the correct disease-free point for β=0.1×10−5 as well as the correct endemic point for β=0.3×10−5. This means that, for the value of the basic reproductive number (threshold parameter), R0<1; then, DFE is stable, i.e., solution to the system S,E,I,R are converging to it. This can be seen from Figures 3(a)–3(c) for h=0.01. Also, for the value of the basic reproductive number (threshold parameter) R0>1, then EE is stable, i.e., solution to the system S,E,I,R is converging to it. This can be seen from Figures 3(d)–3(f) for h=0.01.Table 1
Comparison of NSFD scheme with Euler and RK-4 schemes.
l0.010.11001000Euler PCConvergenceDivergenceDivergenceDivergenceRK-4ConvergenceDivergenceDivergenceDivergenceNSFDPCConvergenceConvergenceConvergenceConvergenceFigure 3
Time plots of susceptible, exposed, and infected population using the NSFDPC method withh=0.01.
(a)(b)(c)(d)(e)(f)Figures4(a)–4(c) show how the NSFDPC method of DFE converges to the equilibrium point for different step sizes. Figures 4(d)–4(f) show how the NSFDPC method of EE converges to the equilibrium point for different step sizes.Figure 4
Time plots of susceptible and infected population using the NSFDPC method with different step sizes: (a)h=0.1, (b) h=10, and (c) h=1000.
(a)(b)(c)(d)(e)(f)For comparison with the well-known RK-4 method to system (1), using the parameter values given in Table 2, it is found that, in Figures 5(a)–5(f), the numerical solution converges for h=0.01 for both equilibrium points, and Figure 6 shows overflow for a critical value of step size, i.e., h=0.1 of DFE and EE, respectively.Table 2
Parameters value of the measles model.
ParametersNμσγ5×107year−10.0245.6year−173year−1DFES0,E0,I0=N,0,0,β=0.1×10−5EESe,Ee,Ie=2.435×107,1.125×104,7022.671,β=0.3×10−5Figure 5
Time plots of susceptible, exposed, and infected population using the RK-4 method withh=0.01.
(a)(b)(c)(d)(e)(f)Figure 6
Divergence graphs of the population using the RK-4 method with step sizeh=0.1.
(a)(b)Another popular numerical method is the Euler predictor-corrector (E-PC) technique, which employs the explicit Euler method as a predictor and the trapezoidal rule as the corrector. This combination is second-order accurate but has similar stability properties in PECE mode to the explicit Euler method alone. Figures7(a)–7(f) show convergence results for h=0.01 in both cases that are DFE and EE. The value of h=0.1 which produces overflow when solving the system with the parameter values of Table 2 for the PECE combination can be seen in Figure 8.Figure 7
Time plots of population using the Euler-PC method withh=0.01.
(a)(b)(c)(d)(e)(f)Figure 8
Divergence graphs of the population using the Euler-PC method with step sizeh=0.1.
(a)(b)Thus, the presented numerical results demonstrated that NSFDPC has better convergence property following the Euler PC and RK-4, as shown in Figure9. It is proved that the approximations made by other standard numerical methods experience difficulties in preserving either the stability or the positivity of the solutions or both but NSFDPC is unconditionally convergent.Figure 9
Comparative convergence graphs of the NSFDPC method.
(a)(b)


## 7. Conclusions

In this paper, a normalized NSFDPC for the SEIR model concerning the transmission dynamics of measles is constructed and analysed. This proposed numerical scheme is very competitive. It is qualitatively stable, that is, it double refines the solution and gives realistic results even for large step sizes. It is dynamically consistent with a continuous system and unconditionally convergent and satisfies the positivity of the state variables involved in the system. Simulations are carried out and its usefulness is compared with a well-known numerical method of standard difference schemes such as RK4 and Euler predictor-corrector method. The standard finite difference schemes are highly dependent on step sizes and initial value problems, but NSFDPC is independent of these two features which make it more practical. Also, this method saves computation time and memory.



---

*Source: 8830829-2020-12-15.xml*