Positivity preserving operator splitting nonstandard finite difference methods for SEIR reaction diffusion model

Abstract In this work, we will introduce two novel positivity preserving operator splitting nonstandard finite difference (NSFD) schemes for the numerical solution of SEIR reaction diffusion epidemic model. In epidemic model of infection diseases, positivity is an important property of the continuous system because negative value of a subpopulation is meaningless. The proposed operator splitting NSFD schemes are dynamically consistent with the solution of the continuous model. First scheme is conditionally stable while second operator splitting scheme is unconditionally stable. The stability of the diffusive SEIR model is also verified numerically with the help of Routh-Hurwitz stability condition. Bifurcation value of transmission coefficient is also carried out with and without diffusion. The proposed operator splitting NSFD schemes are compared with the well-known operator splitting finite difference (FD) schemes.


Introduction
In childhood epidemic diseases, Measles is considered as highly infectious disease, spread due to respiratory infection by a traveling virus. In 19th century, Arthur Ransom observed the unevenly recurring nature of Measles. Since 1896, age structures, contact rates and school seasons are known as critical agents for the transmission of Measles in a population. Hamer [1] in 1906 presented a model for transmission dynamics of Measles. Later, the principle of "mass action" was introduced which became a fundamental statute in present theory of infectious disease modeling [2][3][4][5]. Infectious disease dynamics is an important application of dynamical systems. Fixed points, steady states or the equilibrium points of a continuous dynamical systems are the values of variables that do not change over time. If a system starts at a nearby state and converges to the equilibrium point, then this equilibrium point is attractive. Many dynamical systems in di erent elds of science and engineering show the chaotic behavior.
If a numerical method produces chaos and it is not the feature of continuous system, then this is called contrived chaos. Many existing numerical techniques show the contrived chaos. If a continuous dynamical system demonstrates a property P, then it is necessary for the numerical method to preserve that property P. This is called dynamical consistency, in the literature, Twizzell et al and various authors [6][7][8][9][10][11][12][13] discussed about the contrived chaos and loss of the dynamical consistency in the nite di erence scheme. Similarly, Fernández [14] discussed about the same drawbacks exhibited in di erent series solution techniques. The aim of this work is to develop and investigate such numerical schemes whose equilibrium points coincide with the equilibrium points of the continuous system and preserve all the essential properties demonstrated by the continuous system. In this work, we consider SEIR epidemic model with di usion presented by Al-Shawaikh and Twizell [15]. This model describes the transmission dynamics of measles. In this model, the total population is taken as 'N' and it is categorized in four classes i.e. susceptible population, exposed population, infected population and recovered population denoted by S, E, I and R respectively. It is assumed that the total population 'N' is constant i.e. birth and death rates in the population are equal. Mathematical model of SEIR epidemic model with di usion is Since R is not present in rst three equations so we can write the system as with the initial conditions and homogeneous Neumann boundary conditions are Ix As discussed above, the purpose of this work is to construct positivity preserving operator splitting nonstandard nite di erence schemes for the system (5)- (7). Since SEIR reaction di usion system (5)-(7) is a population model therefore negative values of susceptible, exposed and infected population is meaningless. Nonstandard nite di erence method is a structural preserving numerical method which is developed by Mickens [16]. NSFD method is an important tool to solve ordinary and partial di erential equations [9-11, 17, 18] as it preserves the properties possessed by continuous system.

Equilibrium Points
There are two steady states [15] of the system (5)-(7), disease-free and endemic steady state. Disease-free state is (S , E , I ) = (N, , ) and the endemic steady state is where R is the reproductive and if R < then the system possesses disease free state and if R > then system possesses endemic state.

Stability of Equilibrium of the System
The system (5)- (7) is linearized about the equilibrium point (S * , E * , I * ) to have small perturbations S(x, t),E(x, t) and I(x, t) as describe in [19] Suppose that the system (14)- (16) possesses Fourier series solution as mentioned in [19], then where k represents the wave number for the node n and k = nπ/ , (n = , , , ...). Replacing the values of S (x, t), E (x, t) and I (x, t) in the system, it is transformed into The variational matrix V is for the equation (20) where The characteristics equation for matrix V is The Routh-Hurwitz stability criterion gives η > , η > , and η η − η > .

Bifurcation Value of Transmission Coe cient β
The values of S * , E * and I * are substituted in a , a , · · ·, to obtain The Routh-Hurwitz criterion for stability gives The equation f (β) = provides the bifurcation value of β, which shifts the stability of the equilibrium point (i.e. from stable to unstable). The solution for f (β) gives .
. The point of equilibrium is unstable for any value of β less than .
. Similarly, the bifurcation value of β for system without di usion is .
. The equilibrium point of the system without di usion is not stable for any value of β less than .
. It is concluded that bifurcation value of β with di usion is less than bifurcation value of β without di usion. In order to nd the bifurcation value of β we consider case 1 in Table 1 and d S = . , d E = . , d I = . .

Numerical Methods
In order to handle the complexity and nonlinearity of the reaction-di usion system (5)-(7), we implement two classical and two proposed nonstandard operator splitting nite di erence methods. Operator splitting FD schemes are used by various authors to solve complex problems containing di erential equations [19,[21][22][23][24][25][26][27][28][29][30]. Generally, the SEIR reaction-di usion equations are split into two systems of equations. The nonlinear reaction steps which are used for the rst half step for time and the linear di usion equation which are implemented for the second half step of time. Firstly, we present forward Euler FD and Backward Euler schemes to solve the equations.
The procedure for the backward Euler FD scheme is Now we construct nonstandard FD schemes. By applying rules de ned by Mickens [16], both nonstandard FD schemes at rst half time step areS Positivity of the solution requires that if The procedure is di erent for both NSFD operator splitting schemes at second half-time step. For explicit NSFD operator splitting scheme, we use The procedure for the implicit NSFD operator splitting scheme at second half time step is where

. Accuracy and Stability of the Numerical Methods
In operator splitting methods, the stability and consistency depends on their split solutions [21,22]. In all the above splitting methods, the reaction step is solved exactly and time derivative has O(τ) accuracy [21,22]. Similarly, the di usion step has O(h ) accuracy and the cumulative accuracy of all numerical schemes is rst order in time and second order in space [21,22]. For the stability, the reaction step is unconditionally stable in all cases as it is solved exactly [21,22]. In forward Euler method operator splitting FD method and explicit NSFD operator splitting method (32)-(34), the stability region of di usion step is In the backward Euler method operator splitting FD method and implicit NSFD operator splitting method, the di usion step is unconditionally stable.

. Positivity of NSFD methods
It can be observed that NSFD methods in reaction step (39) This implies that which is the stability condition (49) of NSFD method (39)-(41). This proves that the rst NSFD method preserves the positivity property in its stability region. In order to prove the positivity of unconditionally stable implicit NSFD method (46)-(48), we use M-matrix theory [33].

Proof
The system (46)-(48) can be written as where A, B and C are the square matrices · · · · · · · · · · · · a a and If we consider that S m > ,E m > and I m > , then from the property of M-matrix and (56)-(58) we get S m+ > ,E m+ > and I m+ > . So by the induction, the theorem is proved. It follows from Theorem 5.2.1 that implicit NSFD operator splitting scheme preserve positivity property unconditionally.

Numerical Experiment
In this section, we perform a numerical experiment for all the methods at both equilibrium points. For the numerical test, the following values of parameters are used for disease-free equilibrium [15,20] N = × , µ = . , σ = . , γ = , β = . × − .
In this experiment we take d S = . , d E = . , d I = .
. The initial conditions for the system (5)- (7) are Now the graphs of all the schemes at disease free equilibrium point is presented. First of all we show the graphs of Forward Euler operator splitting FD scheme. The gures 1-3 indicate the graphs of disease free equilibrium using forward Euler operator splitting FD scheme at h = . and λ = .
. Figure 2 and gure 3 show that forward Euler operator splitting FD scheme fails to preserve the positivity property. The gures 4-6 represent the mesh graphs of susceptible, exposed and infected individuals using forward Euler explicit operator splitting FD scheme at h = . and λ = . , λ = . and λ = .
. It can be seen that all the graphs in gures 4-6 show over ow and divergence. Remember that these are the graphs at endemic equilibrium point. Forward Euler explicit operator splitting FD scheme fails to converge to the endemic equilibrium point. Secondly we present the simulation of backward Euler implicit operator splitting FD scheme. The gures 7-9 re ect the graphs of disease free equilibrium using backward Euler implicit operator splitting FD scheme at h = . and λ = .
. Figure 8 and gure 9 show that backward Euler implicit operator splitting FD scheme also fails to preserve the positivity property. The gures 10-12 represent the mesh graphs of susceptible, exposed and infected individuals using backward Euler implicit operator splitting FD scheme at h = . and λ = . , λ = . and λ = .
. It can be seen that all the graphs in gures 10-12 show over ow and divergence. Remember that these are the graphs at endemic equilibrium point. Backward Euler implicit operator splitting FD scheme fails to converge to the endemic equilibrium point. Thirdly, the graphs of susceptible, exposed and infected population at disease free and endemic equilibrium point using nonstandard operator splitting explicit FD method are presented. The gures 13-15 re ect the graphs of disease free equilibrium using NSFD explicit operator splitting scheme at h = . and λ = .
. Graphs clearly show that NSFD operator splitting explicit scheme preserves the positivity property and converges to the disease free equilibrium point (S , E , I ). The gures 16-18 represent the graphs of endemic equilibrium point for susceptible, exposed and infected individuals using NSFD operator splitting explicit scheme at h = . and λ = . , λ = . and λ = .
. Graphs show that NSFD operator splitting explicit scheme preserves the positivity property and converges to the endemic equilibrium point (S * , E * , I * ). At the end, we present the simulations using NSFD operator splitting implicit method. , λ = . and λ = .
. Graphs clearly show that NSFD operator splitting implicit scheme preserves the positivity property and converges to the disease free equilibrium point (S , E , I ). The gures 22-24 represent the graphs of endemic equilibrium point for susceptible, exposed and infected individuals using NSFD operator splitting implicit scheme at h = . and λ = . , λ = . and λ = .
. Graphs show that NSFD operator splitting implicit method preserves the positivity property and converges to the endemic equilibrium point (S * , E * , I * ).

Conclusion
In this work, we developed two positivity preserving nonstandard nite di erence schemes which did not bring contrived chaos. The proposed NSFD schemes are operator splitting FD schemes. We showed that classical operator splitting FD schemes bring the contrived chaos which lead to the numerical instabilities. We also give the stability analysis of SEIR reaction-di usion system and nd out the bifurcation value of transmission coe cient with the help of Routh-Hurwitz criteria. The proposed NSFD operator splitting schemes should prove to be of value in the solution of nonlinear continuous dynamical systems. Our future plane includes the construction of NSFD operator splitting techniques for hidden attractors in continuous dynamical systems [34][35][36][37].