Development of Explicit Schemes for Diffusive SEAIR COVID-19 Epidemic Spreading Model: An Application to Computational Biology.

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. © Shiraz University 2021.

Introduction

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 high-order 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 (González-Parra et al. 2010), 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 -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 . 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 denotes fractions of susceptible individuals, denotes fractions of exposed people, denotes fractions of asymptomatic individuals, denotes infectious individuals, and denotes fractions of recovered individuals. Let be the infecting rate of asymptomatic individuals while be the infecting rate of infectious ones. Let be the rate at which exposed become asymptomatic. Let be the enhancement of the fractions of infectious hosts. Let 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, for . One trivial solution of Eqs. (1)–(5) is 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 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 . Since , so the remaining task is to prove and

First, assume that this implies

The expression for becomes this contains three negative terms, which can be replaced with the terms from the following three inequalities

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

Construction of Numerical Schemes

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 and using Taylor series

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

Comparing coefficient of and on both side of Eq. (17) yields

Solving Eq. (19) for gives

Similarly, the difference equations for Eqs. (2)–(5) can be constructed by using the proposed scheme.

Stability of Numerical Schemes

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 , .

According to von Neumann stability criterion, the following transformations are considered

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

where and where , .

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

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

The amplification matrix is which have eigenvalues.

and .

The stability condition can be imposed on and as

Two conditions (27) can be split into four cases/conditions.

From the graphical illustration, it can be verified that for the region of stability is and for the stability region is . Therefore, by combining and , the stability region is the whole positive set of real numbers except and is not included because the denominator for the expression of becomes zero.

Consistency of Numerical Schemes

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

By using (15), (16), and (30) into difference Eq. (14) yields

Consider the relationship , and dividing on both sides of Eq. (31) gives

This implies

Applying the limits as in Eq. (33) results from the original Eq. (1) evaluated at a point . Thus, proposed scheme is consistent for Eq. (1), and similarly, the consistency for remaining Eqs. (2)–(5) can be proved.

Convergence of Numerical Schemes

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

Error at first-time level is bounded

Proof. By combining eqs. (1)– (2) in a single equation in the form

Discretize Eq. (32), it is obtained

Applying Taylor series for in Eq. (33), it gives

Comparing coefficients of on both sides of Eq. (34), it yields

,

The discretization of Eq. (33) becomes

Let be the error between exact and numerical scheme at grid point and at time level, i.e., . The error notation equation for Eq. (36) can be expressed aswhere .

Equation (38) can be expressed aswhere

First, consider , then Eq. (39) becomes

Let , then Eq. (41) becomes

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

Substituting in Eq. (42) gives

Substituting in Eq. (42) gives

.

Substituting in Eq. (42) gives

Similarly, substituting in Eq. (42) gives

If this continues, then for odd when and in Eq. (44), the finite geometric series becomes infinite expressed as

and the geometric series converges if and .

For even

The infinite series converges if and .

The cases when can be discussed similarly.

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

Results and Discussion

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 . 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.

Fig. 1

Comparison of solutions over spatial variable

Fig. 3

Comparison of absolute error for susceptible individuals

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 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 solution if spatial step size decreases and temporal step size approaches zero. But this may consume more time to get a solution near to first-order method’s solution. The step size in space and time for the non-standard finite difference method, forward Euler method, and proposed scheme are chosen to be the same. The values of the parameters used in these simulations are given as

Fig. 2

Comparison of solutions over time variable

, and the value of .

Figures 3, 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.

Fig. 4

Comparison of absolute error for exposed individuals

Fig. 5

Comparison of absolute error for asymptomatic individuals

Fig. 6

Comparison of absolute error for infected individuals

Fig. 7

Comparison of absolute error for recovered individuals

, , ,where , denote the number of grid points, the number of time levels, final time, and domain length, respectively, with 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 , as well as overtime variable .

Conclusion

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 were shown to be in straight lines. Following the completion of this work, it is possible to propose additional applications for the current methodology (Nawaz and Arif 2021; Nawaz et al. 2020; 2021a; 2021b).