key: cord-1055398-51axv011
authors: Nawaz, Yasir; Arif, Muhammad Shoaib; Ashraf, Muhammad Usman
title: Development of Explicit Schemes for Diffusive SEAIR COVID-19 Epidemic Spreading Model: An Application to Computational Biology
date: 2021-09-13
journal: Iran J Sci Technol Trans A Sci
DOI: 10.1007/s40995-021-01214-0
sha: a985985ce61be5ae745a1c7bbc66650f9289b3e6
doc_id: 1055398
cord_uid: 51axv011

In this contribution, a first-order time scheme is proposed for finding solutions to partial differential equations (PDEs). A mathematical model of the COVID-19 epidemic is modified where the recovery rate of exposed individuals is also considered. The linear stability of the equilibrium states for the modified COVID-19 model is given by finding its Jacobian and applying Routh–Hurwitz criteria on characteristic polynomial. The proposed scheme provides the first-order accuracy in time and second-order accuracy in space. The stability of the proposed scheme is given using the von Neumann stability criterion for standard parabolic PDEs. The consistency for the proposed scheme is also given by expanding the involved terms in it using the Taylor series. The scheme can be used to obtain the condition of getting a positive solution. The stability region of the scheme can be enlarged by choosing suitable values of the contained parameter. Finally, a comparison of the proposed scheme is made with the existing non-standard finite difference method. The results indicate that the non-standard classical technique is incapable of preserving the unique characteristics of the model’s epidemiologically significant solutions, whereas the proposed approaches are capable of doing so. A computational code for the proposed discrete model scheme may be made available to readers upon request for convenience.

Mathematical models in physical sciences are useful for numerical testing. Since verifying some scientific phenomena consumes a significant amount of time, so that numerical procedures can be adopted. The mathematical models in epidemic diseases describe some physical happening for fractions of individuals. SIR is one of the models containing three categories of existing mathematical models: susceptible, infected, and recovered. This SIR model can be extended by adding more categories of people. These models may estimate the real situation, and some are based on the assumption(s). The real model must reflect the physical happening by doing experiments. Once the model is constructed, it can estimate the spreading of disease, and the consumption of time is reduced by solving the model rather than doing experiments .

In (Basnarkov 2021) , author proposed a mathematical model of COVID-19. The model has consisted of five equations. The condition of linear stability was found by applying the Routh-Hurwitz criterion on characteristic polynomial. Routh-Hurwitz criterion is an alternative way for proving negative eigenvalues. Once the mathematical model is constructed, numerical schemes are used to solve the considered models because certain numerical approaches to epidemic models have been described in the literature. But most of the existing numerical schemes do not provide the guarantee to get positive solutions. Even highorder schemes exist for solving epidemic models, but stability and positivity are the issues to use in particular situations. The non-standard finite difference method is unconditionally stable among numerical schemes and guarantees a positive solution without applying any restriction and condition. The non-standard explicit finite difference method uses functions that are replaced with fixed step sizes. But these methods do not provide accurate results, which can be seen by comparing the results obtained by non-standard finite difference method and some accurate scheme. In (Sweilam et al. 2017) , the mathematical tuberculosis (TB) model was studied by applying non-standard finite difference methods. The comparison was also made with Euler implicit and Runge-Kutta methods.

Anwar Zeb and Abdullah Alzahrani (Zeb and Alzahrani 2021) formulated a mathematical model of smoking. The reversion class was added, and techniques for finding local and global stability related to smoking-free and positive smoking equilibrium points were also given. In (Anguelov et al. 2014) , two mathematical models of biological systems were presented. Richardson's extrapolation was applied to the solution obtained by the non-standard finite difference method. It was concluded that Richardson's extrapolation improved accuracy. In , two schemes are given for finding numerical solutions to the continuous-time model. One of these methods is the implicitly derived-explicit Euler method, and the second method is the non-standard finite difference method. Still, the former method did not satisfy some essential features. In (Garbab et al. 2011) , finite difference methods are constructed dynamically for the considered model. One method did not ensure some important features of the model, including spurious behavior and numerical instability. Development of non-standard finite difference has been given in Arenas et al. (2010) . The scheme was the predictor-corrector type, but it did not satisfy some features of the model. Some more work on the non-standard finite difference method and compact numerical schemes can be seen in Yanga (2016); Raza et al. 2021; Martin-Vaquero et al. 2018; Arenas and Gilberto Gonz'alez-Parra, Benito M. Chen-Charpentier, 2017 ; Muhammad Rafiq et al. 2021; Arif et al. 2019; Nawaz and Arif 2021 Nawaz et al. 2021a) . A Caputo fractional derivative-based COVID-19 pandemic model has been formulated in Mohammad et al. (2021) . The total reported cases were simulated based on Riesz wavelets. An analysis approach was used to explore and solve the mathematical model for coronavirus Bats-Hosts-Reservoir-People Coronavirus (Gao 2020) . This method is based on the construction of various functional and finds the solution in iterations forms. One of the essential ingredients of this method is the Lagrange multiplier that can be found using the chosen form of a linear operator. The disadvantage of the method is producing the so-called noise terms that can be responsible for time consumption. Another type of analytical method named as q-homotopy analysis transform method was implemented for finding solutions of the fractional epidemic model of ordinary differential equations (Wei Gao et al. 2020) . The method found the solution in series form, and it has been shown that the considered scheme was highly emphatic. The applicability and effects of the fractional derivatives were confirmed to real-world problems. For studying the dynamics of a fractional COVID-19 model, a computational scheme was proposed (Harendra Singh and Srivastava 2020) based on the memory principle and discretization of the domain. The stability of the proposed scheme was discussed, and its efficiency was shown by listing the CPU time. The proposed scheme involved low computational cost, and it was shown that it would work for long-time behavior. An infection system of the coronavirus with a non-local operator has been investigated (Gao and 1, 2, P. Veeresha 3, D. G. 2020) by employing the fractional natural decomposition method, the considered method based on natural transform and Adomian decomposition methods. The reported model was composed of the data taken from the city of Wuhan, China. The Adomian decomposition method finds the solution in different components form. The sum of these components is the required solution of the considered model, and the so-called Adomian polynomials handle nonlinear terms. The method requires a differential operator, which should be invertible. This can be considered as one of the deficiencies that arise in Adomian's method. A COVID-19 model was focused on nearly pharmaceutical intervention (Zamir et al. 2021) . The most active/sensitive parameters were presented after finding the basic reproduction number R 0 . To control the sensitive parameters, non-pharmaceutical interventions have been applied. More work on epidemic models can be seen in Danane et al. (2021) ; P. Veeresha et al. 2021) and references therein.

In numerical schemes, the first question can be considered to be the accuracy of the considered or proposed scheme. A scheme that is not even first-order accurate may consume more time to converge. The high-order schemes converge faster than the lower-order schemes. Therefore, if a scheme is not even first-order accurate, how can it provide faster convergence of the acquired solution? Therefore, faster convergence cannot be achieved by implementing the majority of non-standard finite difference methods that provide solution positivity. This is the main drawback of using most of the existing non-standard finite difference methods. But, the question may arise that the time step size in non-standard methods can vary to achieve first-order accuracy. Again the first-order accuracy will not be guaranteed, and the choice of variable step size must depend on the use of the Taylor series so that the resulting difference equation will provide the first-order accuracy.

On the other hand, the simple basic forward Euler scheme is first-order accurate and explicit. Still, it is conditionally stable, so a larger value of diffusion number in solving time-dependent partial equations cannot be chosen. Therefore, it has the drawback of choosing a limited diffusion number. The Dufort-Frankel method is another numerical scheme that provides the first-order accuracy. It is explicit and unconditionally stable, but the conditions of getting a positive solution cannot be obtained by employing the Dufort-Frankel method in epidemic disease models. One of the advantages of explicit schemes is their faster convergence compared to implicit schemes. Implicit schemes, on the other hand, can provide unconditional stability but may take longer to converge than explicit schemes.

In this work, a time and space discretization scheme provide first-order and second-order accuracy in space. All the presented schemes are parameter dependent. So the stability region can be different for different values of the parameter. Since a particular non-standard scheme does not provide the first-order accuracy, the presented scheme provides the first-order accuracy using Taylor series expansions. It has the limitation of providing a conditionally positive solution instead of getting the unconditional positive solution. So, the positivity of the solution depends on the found conditions. These conditions are determined during the construction of the proposed scheme for the considered type of model. So, if the solution satisfies the conditions and the initial conditions are chosen to be positive, the proposed scheme will provide a positive solution.

Consider the epidemic model discussed in Basnarkov (2021) . In this model, the recovery rate of exposed individuals was not given, which is modified by incorporating the healing rate of exposed individuals. Let S denotes fractions of susceptible individuals, E denotes fractions of exposed people, A denotes fractions of asymptomatic individuals, I denotes infectious individuals, and R denotes fractions of recovered individuals. Let a be the infecting rate of asymptomatic individuals while b be the infecting rate of infectious ones. Let c be the rate at which exposed become asymptomatic. Let r be the enhancement of the fractions of infectious hosts. Let l be the same recovery rate of exposed, asymptomatic, and infectious individuals. The mathematical model with effects of diffusion can be expressed as

From Eqs.

(1)- (5), it can be verified that the total number of individuals in all states is constant, i:e:

(1)-(5) is 1; 0; 0; 0; 0; 0 ð Þwhich corresponds to the states for the absence of the pathogen.

For finding linear stability of the equilibrium states, the system of Eqs. (1)-(5) is linearized first so that the Jacobean can be expressed as

Substituting the equilibrium 1; 0; 0; 0; 0 ð Þinto Jacobian (6), it yields

Since Jacobian (7) contains two columns having all zero entries, so two eigenvalues of Jacobian are zero, and the remaining three can be found from the characteristic polynomial, which is given as

where

According to the Routh-Hurwitz criterion, the real parts of nonzero eigenvalues of Jacobian (7) are negative if and only if a 2 [ 0 and a 1 a 2 [ a 0 [ 0. Since a 2 [ 0, so the remaining task is to prove a 1 a 2 [ a 0 and a 0 [ 0:

First, assume that a 0 [ 0; i:e:;

this implies

The expression for a 1 a 2 À a 0 becomes 8cl 2 À ac 2 þ 2c 2 l þ cr 2 þ c 2 r þ 2lr 2 þ 8l 2 r þ 8l 3 À 2acl À acr þ bcr þ 6clr

this contains three negative terms, which can be replaced with the terms from the following three inequalities

clr þ lr 2 þ l 2 r þ cr 2 [ acr ð13Þ

These inequalities (11)-(13) are obtained from (9). Substituting (11)-(13) into (10) and simplifying yields positive terms, which proves the inequality a 1 a 2 [ a 0 . Thus, the assumption (8a) is the condition for linear stability of equilibrium state.

To solve Eqs. (1)-(5), a numerical scheme is proposed. The construction of the numerical scheme is started with the following procedure. Let the difference equation for Eq. (1) can be expressed as

Expanding S nþ1 i and S nÀ1 i using Taylor series

Substituting (15) and (16) into (14) gives 

Comparing coefficient of S n i and Dt oS ot À Á n i on both side of Eq. (17) yields

Solving Eq. (19) for c gives

Similarly, the difference equations for Eqs.

(2)-(5) can be constructed by using the proposed scheme.

To check the stability of the proposed scheme, the von Neumann stability criterion is employed. Before this, consider the time-dependent partial differential equation with only diffusion effects.

For solving Eq. (21), a scheme can be constructed in the following manner

where a þ b ¼ 1, c ¼ 1þb 1À2d . According to von Neumann stability criterion, the following transformations are considered

By using transformation (23) into Eq. (22) and simplifying yields

where a 1 ¼ 2cdcosh

Since the proposed scheme is constructed on a threetime level, so one more equation is required to find

Combing eqs. (24) and (25) gives matrix-vector equation

The amplification matrix is a 1 b 1 1 0 ! which have eigenvalues.

k 1 ¼ a 1 2 À 1 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 1 þ 4b 1 p and k 2 ¼ a 1 2 þ 1 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 1 þ 4b 1 p . The stability condition can be imposed on k 1 and k 2 as k 1 j j 1; and; k 2 j j 1 ð27Þ

Two conditions (27) can be split into four cases/conditions. From the graphical illustration, it can be verified that for a [ 0; the region of stability is d:5; and for a; the stability region is d [ 0:5. Therefore, by combining a [ 0 and a, the stability region is the whole positive set of real numbers except d ¼ 0:5 and d ¼ 0:5 is not included because the denominator for the expression of c becomes zero.

The consistency of the scheme is checked by employing the Taylor series

15 By using (15) (1) evaluated at a point t n ; x 2 ð Þ. Thus, proposed scheme is consistent for Eq. (1), and similarly, the consistency for remaining Eqs. (2)-(5) can be proved.

Theorem: The proposed scheme for Eqs. (1)-(5) converges if the following conditions are satisfied,

Proof. By combining eqs.

(1)-(2) in a single equation in the form

Discretize Eq. (32), it is obtained 

Comparing coefficients of oE ot À Á n i on both sides of Eq. (34), it yields 1 ¼ Àb þ c 1 À 2d À cDt À lDt ð Þ ,

The discretization of Eq. (33) becomes 

Let e n 1;i be the error between exact and numerical scheme at x i th grid point and at t n th time level, i.e., S n i À S n i ¼ e n 1;i . 

1 þ 2cdd 2 Dt þ ccDt þ clDt j j e nþ1 ae n þ be nÀ1 þ cDt

where e nþ1 ¼ max max .

Equation (38) can be expressed as

where

First, consider max e n ; e nÀ1 ð Þ¼e nÀ1 , then Eq. (39) becomes

Let d 1 þ d 2 ¼ d 3 , then Eq. (41) becomes

Since fixed initial conditions are considered to apply, so e 0 ¼ 0. Also, let the error at the first-time step is bounded, i.e., e 1 M.

Substituting n ¼ 1 in Eq. (42) gives

Substituting n ¼ 2 in Eq. (42) gives

Similarly, substituting n ¼ 4 in Eq. (42) gives

If this continues, then for odd n

when d 3 j j 1 and n ! 1 in Eq. (44), the finite geometric series becomes infinite expressed as 1 1Àd 3 ¼ :: þ d nÀ1 3 þ . . . þ d 3 þ 1 and the geometric series

The cases when e n ¼ max e n ; e nÀ1 ð Þ can be discussed similarly.

Also, eqs.

(3)-(5) can be discussed similarly by combining eqs. (3)-(5), and then convergence conditions can be found.

The proposed numerical scheme is applied to the epidemic model of diffusive COVID-19. The stability of the scheme is based on a parameter. So, a class of schemes can be constructed by choosing different values of the parameter. Two different signs of parameters provide different stability regions. The negative values of a parameter provide schemes that converge for large values of diffusion number, whereas the scheme can be constructed for small values of the diffusion number. So, in this manner, coupled unconditionally stable scheme can be achieved by choosing two different values of a parameter a. The results obtained by the scheme are closer to the first-order forward Euler method, but the results obtained by the non-standard finite difference method are a little far away from the first-order Euler method. The comparison can be seen by looking at Figs. 1, 3 . These results are compared over space coordinate instead of the time coordinate. The first-order finite difference discretizations are chosen at each end of the boundaries due to Neumann boundary conditions. This type of boundary condition requires any additional iterative method, so the Gauss-Seidel iterative method is employed. Due to the larger stability region, a large diffusion number can be chosen, reducing the computational time and providing more accurate results than the non-standard finite difference approach. These schemes may consume small time than the standard implicit schemes for solving nonlinear differential equations. Figures 1, 2 and 3 show the solution obtained by applying the forward Euler method, proposed scheme, and non-standard finite difference method over the spatial variable x and time, respectively. Due to the fact that the Euler method can be used in circumstances where the solution is positive and stable for specified parameter values and step sizes in space and time. So, in these circumstances, the forward Euler method's solution can be considered a standard first-order solution. Therefore, methods can be compared with this method to check the first-order accuracy. In Figs. 1 and 2 , the more accurate first-order solution will be considered that which remains near Euler's solution. The non-standard finite difference method is not accurate enough to provide the appropriate solution at a reduced value of diffusion number diffusive epidemic models. The solution obtained by the non-standard finite difference method will approach the first-order 4 , 5, 6 and 7 show three-dimensional surface plots for absolute error with underneath contours. The following data are used to draw these plots. a ¼ 0:9; b ¼ 0:7; c ¼ 0:9; r ¼ 0:95; l ¼ 0:97, d i ¼ 0:9 for i ¼ 12; 3; 4; 5, N x ¼ 40; N t ¼ 400; t f ¼ 10; L ¼ 10,where N x , N t ; t f andL denote the number of grid points, the number of time levels, final time, and domain length, respectively, with a ¼ 0:1 is used to apply a particular scheme in Figs. 1, 2, 3 , 4, 5, 6 and 7. From these Figs. 1, 2, 3, 4, 5, 6 and 7, it can be observed that the solutions obtained by the non-standard finite difference method remain away from the first-order accurate solution in most of the considered cases. The solution obtained by the non-standard finite difference method is not accurate enough over spatial variable x, as well as overtime variable t. 

Explicit schemes containing two parameters have been proposed to solve time-dependent partial differential equations, which was better in terms of accuracy than most existing non-standard finite difference methods. The scheme could find conditions of getting a positive solution using the positive choice of initial conditions. The proposed numerical approach was used to solve a modified COVID-19 model. The proposed scheme is first-order accurate in time, and due to the use of second-order space discretization, it was second-order accurate in space. From the comparison of schemes, it can be concluded that the proposed scheme provided a more accurate solution for the considered set of partial differential equations than that obtained by the existing non-standard difference method.

Furthermore, results indicate that the non-standard classical technique cannot find an accurate solution to the model, whereas the proposed approaches can do so. Over the spatial variable, susceptible were seen to be in decreasing behavior. In contrast, the remaining four categories of individuals had increasing behavior, and these increasing and decreasing behaviors were so small that they can be considered constants. Due to this reason, the comparison plots over spatial variable x were shown to be in straight lines. Following the completion of this work, it Iran J Sci Technol Trans Sci is possible to propose additional applications for the current methodology 2021a; 2021b) .

Dynamically consistent non-standard finite differences schemes for epidemiological models

A non-standard numerical scheme of predic-tor_corrector type for epidemic models

Construction of non-standard finite difference schemes for the SI and SIR epidemic models of fractional order

A Reliable Numerical analysis for stochastic hepatitis b virus epidemic model with the migration effect

New investigation of bats-hosts-reservoir-people coronavirus model and application to 2019-nCoV system

Dynamicallyconsistent non-standard finite difference method for an epidemic model

Combination of non-standard schemes and Richardson's extrapolation to improve the numerical solution of population models

Numerical simulation and stability analysis for the fractional-order dynamics of COVID-19

Variable step length algorithms with high-order extrapolated non-standard finite difference schemes for a SEIR model

Kottakkaran Sooppy Nisar, Mathematical analysis and simulation of a stochastic COVID-19 Levy jump model with isolation strategy

SEAIR Epidemic spreading model of COVID-19

The dynamics of COVID-19 in the UAE based on fractional derivative modeling using Riesz wavelets simulation

A new class of a-stable numerical techniques for ordinary differential equations: application to boundary-layer flow

Modified class of explicit and enhanced stability region schemes: application to mixed convection flow in a square cavity with a convective wall

Modified class of explicit and enhanced stability region schemes: application to mixed convection flow in a square cavity with a convective wall

A finite difference method and effective modification of gradient descent optimization algorithm for MHD fluid flow over a linearly stretching surface

An explicit fourth-order compact numerical scheme for heat transfer of boundary layer flow

A class of second-order schemes with application to chemically reactive radiative natural convection flow in a rectangular enclosure

New Dynamical Behaviour of the Coronavirus (2019-Ncov) Infection System with Non-Local Operator from Reservoirs to People

Nisar Kottakkaran Sooppy, Iqbal Zafar (2021) Mathematical analysis and design of the non-standard computational method for an epidemic model of computer virus with delay Effect: applicationof mathematical biology in computer science

Non-standard finite difference method for solving the multi-strain TB model

A new study of unreported cases of 2019-nCOV epidemic outbreaks

Prakasha 4 and Haci Mehmet Baskonus, Novel dynamic structures of 2019-nCoV with Nonlocal Operator via Powerful Computational Technique

Non-standard finite difference scheme for a diffusive within-host virus dynamics model with bothvirus-to-cell and cell-to-cell transmissions

Hammouch Zakia (2021) Threshold condition and non pharmaceutical interventions's control strategies for elimination of COVID-19

Non-standard finite difference scheme and analysis of smoking model with reversion class

Acknowledgements The authors are grateful to the reviewers for their valuable comments and suggestions in advance. We are also very thankful for providing an excellent study atmosphere and infrastructure for the Vice-Chancellor of Air University Islamabad.Author Contributions All authors contributed equally to this research work.Funding The author(s) received no financial support for the research, authorship, and/or publication of this article.Code availability A computational code for the proposed discrete model scheme may be made available to readers upon request for convenience.

Conflict of interests The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

The manuscript included all required data and information for its implementation.