An explicit unconditionally stable scheme: application to diffusive Covid-19 epidemic model

An explicit unconditionally stable scheme is proposed for solving time-dependent partial differential equations. The application of the proposed scheme is given to solve the COVID-19 epidemic model. This scheme is first-order accurate in time and second-order accurate in space and provides the conditions to get a positive solution for the considered type of epidemic model. Furthermore, the scheme’s stability for the general type of parabolic equation with source term is proved by employing von Neumann stability analysis. Furthermore, the consistency of the scheme is verified for the category of susceptible individuals. In addition to this, the convergence of the proposed scheme is discussed for the considered mathematical model.


Introduction
Mathematical modeling of epidemic diseases is one of the branches of modeling concerned with somehow estimating and predicting some insight into actual disease. In the literature, the constructed mathematical models for epidemic diseases were the first-order differential equations system that might have been constructed on some assumptions. SIR models belong to the constructed mathematical models of epidemic diseases that can describe some relationships between susceptible, infected, and recovered individuals in COVID-19 epidemic disease. In [1] presented the SIR model that contained health medication factor, initial infected, transmits factor, and human contact factor. One of the concluded results was decreasing COVID-19 spreading by choosing a low contact factor and high medication factor. Reference [2] has consisted of the SEIR model of COVID-19 that contained isolation factors and vaccination as model parameters. The basic reproduction number is found by using the generation matrix method, and the global stability of the given model has also been discussed. The modification of the classical SIR model has been shown in [2] by proposing a susceptible-infected-removed-sick (SIRSi) computational model. The proposed model in [2] considered the level of immunity within the population and asymptomatic cases. The SEIR epidemic model given in [3] used a convex incidence rate. The simulations were obtained by applying the nonstandard finite difference method. An SL 1 L 2 I 1 I 2 A 1 A 2 R epidemic model has been formulated in [4] for spreading an epidemic within the population. The model used an Erlang distribution of time of sojourn in considered compartments. In [5], a simple compartmental Kermack-McKendrick-type epidemic model was introduced with homogeneous and heterogeneously mixed populations. For the dynamics of COVID-19, a primer for analyzing, formulating, and simulating mathematical models were also given.
Some existing models can be used for the COVID-19 epidemic to see some insight into the epidemic disease. The mathematical model considered in [6] has consisted of susceptible, exposed, asymptomatic, infected, and recovered individuals. Exposed were those individuals that had pathogen but cannot transmit it to other individuals. At the same time, asymptomatic individuals could transmit the pathogen. However, they do not know about it and are infected. They knew that they had the disease and can transmit it since quarantine is another category of individuals considered in COID-19 epidemic disease. However, in the present modeling, quarantine individuals are not considered, although they can be regarded as infected individuals. Because if someone is infected, then it means that the individual knows about the disease. This infected individual can be considered one of the quarantine individuals, but quarantine can be regarded as the category of infected and under treatment people. So, for COVID-19 disease, infected and quarantine individuals can be considered to be the same. It is also assumed that the recovered individuals are not shifted to exposed or asymptomatic or infected individuals. The recovered individuals can be assumed to have an ignorable chance of being infected again. For the present modeling of COVID-19, it is also assumed that exposed people cannot be shifted into the category of recovered people.
Some of the numerical solutions for epidemic models included the diffusion effects that can be found in [7][8][9]. [10] has presented a predictor-corrector system to find a solution to large time values for obtaining insight into an epidemic to limit behavior. Variational iteration method and successive approximation methods have been applied in [11] to solve the SIR epidemic model with a constant vaccination strategy. The existing variational iteration method was shown to be inaccurate for the large domain. The existing variational iteration method was improved and identical to the successive approximation method. The modified method was more accurate than the existing one, given in [11]. The susceptible, exposed, infected, diagnosed, recovered (SEIJR) epidemic model was considered in [12] with effects of net inflow of people into a region. Different initial population distributions were considered with the considered model, and it was solved by the numerical method for analyzing the transmission of disease. A diffusive epidemic model [13] has been investigated for describing the transmission of influenza as an epidemic. The spread of the disease showed that diffusion and initial population distribution played an important role.
The COVID 19 pandemic is a worldwide destructive disease that raised severe health issues. It is considered one of the most devastating crises after World War II as it increased the death toll by 1,458,000, which is still rampant. This pandemic surged social issues and the economic recession and environmental disability that led to the destruction of habits, trade, economic relations, forms of work, and political organizations. Reportedly, this disease imparts curb on social movement more than 4 billion people.
Globally, all government and private healthcare departments were unprepared for this trauma which was simply a matter of time that arrived now. Pathogenic disorders wreak havoc in society in the past few months, for which mathematical modeling is the best way to investigate and control them once they enter the community. Nowadays, coronavirus is an essential topic for researchers as regards finding its treatment as the effectiveness and deaths are unbridled.
The virus was first reported in December 2019 in a Chinese city, Wuhan, as an infective agent named coronavirus [14,15]. The viral disease COVID-19 is primarily transferred using droplets produced by infected persons. The infection is transmitted through droplets that are so dense than air particles and immediately fall on the ground, created due to sneezing and coughing of the infected person. COVID-19 confirmed cases reached 4 million earlier in more than 180 countries, and approximately 1,458,000 people have become victims of this dreadful virus [16].
A retrofitted state SIR system to task the overall number of sick circumstances and the specialized obligations on hemodialysis units and hospitals are presented [17]. Nesteruk observed the coronavirus epidemic trying to spread numbers based on assumptions throughout mainland China regrettably. The majority of casualties of COVID are predicted to become much higher than that forecast on February 2020; two days later, 12289 confirmed cases were added. Additional research focuses on updating predictions using up-to-date data and applying more convoluted mathematical representations. There does not exist any approved vaccine or diagnostic drugs for the avoidance and cure of coronavirus. However, research studies on potential antiviral drugs and vaccine candidates are under way in several countries. Vaccine evaluations, growth, and allocation are usually a big task than clinical trials. It is unlikely that the COVID-19 flu shot will be mentally prepared by 2021 within the shortest time possible. The dreadful germ can spread quickly in closely packed locations. Social detachment or low contact rate refers to increasing the disk environment among both people to delay infection spread [18][19][20][21]. They have studied the SIR model to guesstimate the adult location of the coronavirus infectious disease [22].

Diffusive epidemic models
In the literature, some diffusive epidemic models have existed. From these, [23] investigated a diffusive model for the transmission of influenza as an epidemic. The equations have been tackled with the splitting method using different initial conditions of population density. Another diffusive epidemic model of H1N1 has been formulated [24] for gaining a basic understanding of virus behavior. It was assumed that all newborns were susceptible, and also, it was assumed that the mortality rate of individuals is greater than the natural mortality rate. Among the diffusive epidemic models, a vaccinated diffusive epidemic model has also been developed [25] for exploring the impact of diffusion and vaccination on the transmission of dynamics of influenza. In this work [25], a basic reproduction number was found for both vaccinated and non-vaccinated cases. Based on parameters in the system, sensitivity analysis of the reproduction number has been investigated. HIV/AIDS is incurable for human beings mentioned in [26], and a diffusive compartmental model of HIV/AIDS has been proposed with a delay process. The proposed scheme had the advantage of producing a positive solution. The stability and consistency have been given. A nonlinear model for Immunodeficiency Virus (HIV ) has been proposed in [27]. For boundedness and wellposedness, theorems and propositions have been constructed, and the model was solved by employing the evolutionary Padé-approximation technique. This is some literature given on diffusive epidemic models. The reader can find more work on diffusive epidemic models by referring to [23][24][25][26][27]. In [28] an SEIR epidemic model for COVID-19 is constructed using several common control strategies, including hospitalization, quarantine, and external input. The particle swarm optimization (PSO) algorithm is used to estimate the system's parameters using data from Hubei province.
In this contribution, a numerical scheme is proposed to solve the COVID-19 epidemic model with diffusion. The scheme is shown to be unconditionally stable for the considered type of epidemic models. The scheme is first-order accurate because it is constructed to provide the first-order accuracy in time and second-order accuracy in space for diffusion contained epidemic models. The scheme provides a conditionally positive solution. The conditions of finding positive solutions are found in the construction of the proposed scheme. The scheme is constructed on three-time levels, and it is an explicit scheme. The convergence of the modified epidemic model's scheme is also discussed by applying the condition of convergence of infinite geometric series. Since the scheme is constructed on three-time levels, so it requires another scheme to be implemented at the first time level. The proposed scheme can be useful in those mathematical models where the positive solution is required to be found. Other than epidemic models, it can also be applied to solving any time-dependent partial differential equations which contain a diffusion term.
The model given in [6] is modified with diffusion effects, and the modified diffusive COVID-19 model is expressed as The boundary conditions corresponding to Eq. (1)-(5) are expressed as: and the initial conditions are expressed as Moreover, Table 1 provides a summary of the physical meaning of the parameters used in this model. For (1)-(5) become ordinary differential equations, and linear stability of the system is determined by the Jacobean at the equilibrium point (1, 0, 0, 0, 0) given in [1] as In [1], two eigenvalues of the Jacobean (8) are zero, and the remaining eigenvalues can be found from the following polynomial: Using the Routh Hurwitz criteria for linear stability of the system (1)-(5), real parts of the eigenvalues must be negative, and this gives the condition for stability [1].

Numerical scheme
The proposed numerical scheme can solve systems of Eqs. (1)- (5). At this stage, the construction of the numerical scheme is given. The scheme is constructed for Eq. (1) and the remaining Eqs. (2)-(5) will be discretized later on. Consider the following difference equation with unknown parameter a: We expand S n+1 i and S n-1 i using Taylor series, as follows: Substituting the Taylor series expansions (11) and (12) into Eq. (10) yields Equation (13) is further simplified as Comparing coefficients of t( ∂S ∂t ) n i on both sides of Eq. (14) gives Solving Eq. (15) gives the expression for a, Thus Eqs. (1)-(5) are discretized as The expressions for a 1 , a 2 , a 3 , a 4 and a 5 are , Equations (17) a = a 1 , a 2 , a 3 , a 4 , a 5 ≥ 0. Then Eqs. (17)- (21) can be expressed as So, the explicit relationships (22a)-(22e) show that the scheme will produce a positive solution at each time level with the first-order accuracy in time and second-order accuracy in space. The two positive initial conditions can be provided to apply the proposed scheme instead of employing any other scheme on the first time level, which will be constructed on two time levels.

Stability analysis
The stability of the proposed scheme is checked by employing the von Neumann stability criteria. Before starting the procedure of von Neumann stability criteria, consider the general form of the epidemic model given by where u can be considered as one of the susceptible, exposed, infectious, asymptomatic, or recovered individuals and α 1 , β 1 are some rates. Employing the proposed scheme in Eq. (23) yields where a = 2 1-2dtβ 1 . By following the von Neumann stability criteria, the dependent components in the scheme (24) are expressed as where I = √ -1. Applying transformation (25) to Eq. (24), one obtains Simplification of (26) leads to where d = t ( x) 2 . Equation (27) can be expressed as where b 1 = 2adα 1 cos θ 1+(2α 1 d+β 1 t)a , b 2 = 1 1+a(2α 1 d+β 1 t) . Consider one more equation of the form The matrix-vector equation can be constructed as For this case, the amplification factor is a matrix, and the condition of stability can be imposed on the eigenvalues of this matrix, which are expressed as where λ 1 = b 1 2 -1 2 b 2 1 + 4b 2 , λ 2 = b 1 2 + 1 2 b 2 1 + 4b 2 . Let d = α 1 d and β = β 1 t, then b 1 and b 2 can be expressed as Since the eigenvalues λ 1 and λ 2 contain a fractional power, before finding stability conditions, one can first find the region that corresponds to real eigenvalues. For this reason, the expression b 2 1 + 4b 2 should be non-negative. So, For cos θ = 0 So real eigenvalues correspond to the region d ≤ 1-β 2 so in this region condition on the eigenvalue λ 1 can be expressed as Since | cos θ | ≤ 1, consider first cos θ = -1 so inequality (35) can be expressed as 0 ≤ 8β 2 + 8β + 16βd which holds for every value of β and d.
So all inequalities for -1 ≤ λ 2 and λ 2 ≤ 1 are satisfied with the extreme value of cos θ . Thus |λ 1 | ≤ 1 and |λ 2 | ≤ 1 are true for every β and d when d ≤ 1-β 2 . The complex eigenvalues are obtained for the region d ≥ 1-β 2 . For this case, the stability conditions |λ 1 | ≤ 1 and |λ 2 | ≤ 1 are expressed as Hence, -b 2 ≤ 1 and so which is valid for every value of β and d. Thus the proposed scheme is unconditionally stable for Eq. (23), which is considered for the general type of epidemic disease model or parabolic partial differential equations having some source term(s).

Consistency of scheme
Taylor series expansions prove the consistency of Eq. (1). For this reason, consider the Taylor series expansions for S n+1 i and S n-1 i given in (11) and (12) and the following Taylor series expansions: So, using expansions (36)-(37), the following equation can be obtained: Substituting (11), (12), and (38) into the discretized form (10) yields

∂S ∂t
Canceling the same terms on both sides of Eq. (41) yields

∂S ∂t
By incorporating consistency conditions t → 0, x → 0 in (42) yields the original Eq. (1) evaluated at the ith grid point and at time level n.
For n = 1 in Eq. (54) Since e 0 = 0 due to the initial condition, Eq. (55) can be expressed as For n = 2 in Eq. (54) Let the error at the first time obtained be e 1 ≤ M then (57) is expressed as For n = 3 in Eq. (54) For n = 4 in Eq. (54) For n = 5 in Eq. (54) For n = 6 in Eq. (54) If this is continued for a finite number of n then for even n e 2n ≤ δ n 3 M + (δ n-1 Equation (63) is obtained by considering even exponent terms. For odd n e 2n-1 ≤ δ n-1 For large n the series 1 + δ 3 + · · · + δ n-1 3 will converge if |δ 3 | ≤ 1. Similarly, convergence can be found for the cases when max(e n-1 , e n ) = e n . This gives convergence of the proposed scheme for the first two Eqs. (1) and (2). Similarly, convergence can be found for the remaining Eqs. (3)-(5).

Application
The proposed numerical scheme has been constructed and employed for solving the diffusive epidemic model. Initially, the model comprised ordinary differential equations constructed in [6], but it is modified with diffusion effects in this contribution. Due to the fact that diffusion is dependent on a spatial variable, it also makes use of information to transmit individuals from one location to another. The ODEs model [6] only uses the information in time variable, but here time and space both have been utilized, and a diffusive epidemic model has been presented. The main concern is the numerical scheme. Among existing numerical schemes, the present attempt is made to construct a numerical scheme The other advantage of this scheme is first-order accuracy for solving partial differential equations. Nonstandard finite difference method is not even first-order accurate, so it may produce some doubt full results but present strategy of constructing scheme is based on applying Taylor series, so theoretically it is first-order accurate which has the advantage for consumption of less time than one of the nonstandard schemes and this can be seen by drawing the graphs on the spatial variable when solving partial differential equations.
Since the von Neumann type boundary conditions are employed on the boundary, it makes the proposed explicit scheme implicit. An additional iterative approach of the Gauss-Seidel iterative method is also employed to solve the resulting difference equations. The iterative approach tackles the von Neumann type boundary condition on the left boundary. The von Neumann type boundary condition can be incorporated explicitly if it is employed on the last grid point. In this case, the backward difference formula can be considered to find each dependent variable's value on the last grid point. But, for the first gird point, the first-order forward difference formula using Gauss-Seidel iterative method and this can be expressed in the following manner: which implies Similarly, this formula can be employed on all boundary conditions imposed on exposed, asymptomatic, infected and recovered individuals at the left endpoint. The boundary condition on the right endpoint can be tackled explicitly. Using the Gauss-Seidel iterative method, it reads where a 1 = a is given in (16). Figure 1 shows the behaviors of susceptible and exposed individuals over time when the rate γ varies. Figure 1 also indicates that both susceptible and exposed individuals are de- creasing with an increase in the conversion rate γ . The exposed people will decrease due to their transmission from the exposed category to the asymptomatic category, and since both asymptomatic and infected individuals have to increase. Decreasing behaviors so susceptible people decrease mostly, but these people also have increasing behavior which is very small and can only be seen on a very small scale. Figure 2 shows the asymptomatic and infected individuals over time. Both categories have increasing and decreasing behavior due to the increase in conversion rate γ . Figure 3 presents asymptomatic and infected individuals over time. Figure 3 clearly shows the enhancement and decay of infected and asymptomatic people by enhancing the growth rate σ because asymptomatic people will shift the category of asymptomatic individuals to infected individuals. Figure 4 shows the Two dimensional phase portrait of susceptible and exposed individuals behavior of susceptible and exposed individuals. Figure 4 shows that susceptible people are increasing and exposed people are decreasing by enhancing the recovery rate parameter μ. Figure 5 shows the behavior of asymptomatic and infected individuals over time.
It is seen clearly from    Three dimensional phase portrait of susceptible, exposed and asymptomatic individuals Figure 10 Three dimensional phase portrait of exposed, asymptomatic and infected individuals Figures 12-14 are drawn to elucidate the comparison of three schemes. The solution obtained by the proposed scheme and the nonstandard finite difference method has been displayed in Figs. 13 and 14, respectively. The values on the y-axis can be seen for clear comparison between three scheme. The proposed scheme produced the solution near to first-order method but nonstandard finite difference method produced the solution a little away from the solution obtained by first-order method (forward Euler method). These figures show the advantage of using proposed scheme. Since the Euler method does not converge on those time levels which are used by proposed scheme so larger number of time levels (N t ) are used to get converged solution. Figures 15 and 16 are surface plots for exposed and infective individuals. These figures are three dimensional views of exposed and infective individuals over t and x-axis. From Figs. 15 and 16, it can be observed how  plotted individuals behave on t-and x-axis. In the captions of these figures N x and N t denote the number of grid points and number of time levels, respectively. The parameter a and b are used for choosing a particular scheme.

Conclusion
A first-order in time and second-order in explicit space scheme has been proposed. The scheme is unconditionally stable according to von Neumann's stability criteria. The pro- can be considered to solve related epidemic models and other parabolic partial differential equations.