key: cord-0900419-onbqiy3j
authors: Nawaz, Yasir; Arif, Muhammad Shoaib; Abodayeh, Kamaleldin; Shatanawi, Wasfi
title: An explicit unconditionally stable scheme: application to diffusive Covid-19 epidemic model
date: 2021-08-03
journal: Adv Differ Equ
DOI: 10.1186/s13662-021-03513-7
sha: 3f9819ab4ea7fbaed79e83ec84fbda09f83ed340
doc_id: 900419
cord_uid: onbqiy3j

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.

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

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

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.

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.

Consider the second case when cos θ = 1 and inequality -1 ≤ λ 1 yields

This implies 0 ≤ 4(β + 4d) 2 + 4β 2 + 16d which is also true for every value of β and d.

Consider the case now when λ 1 ≤ 1 and cos θ = -1,

and hence -2 -2β -8d ≤ -4β 2 -16d + 4. This is also true because the negative number is always less than the positive number for d ≤ 1-β 2 . The fourth case, when λ 1 ≤ 1 and cos θ = 1, yields -2 -2β ≤ -4β 2 -16d + 4, This is also true for every β and d.

Therefore four cases for eigenvalues λ 1 have been discussed, and inequality |λ 1 | ≤ 1 holds for every value of β and α when d ≤ 1-β 2 . Four cases for the second eigenvalue |λ 2 | ≤ 1 can be discussed at extreme values of cos θ , so four cases are given as when -1 ≤ λ 2 and cos θ = -1, and the following inequality can be obtained:

when -1 ≤ λ 2 and cos θ = 1, the inequality can be obtained in the form -2 -2β -8d ≤ -4β 2 -16d + 4 when λ 2 ≤ 1 and cos θ = -1, an inequality can be expressed in the form of 0 ≤ 4β 2 + 16βd + (2β + 8d) 2 + 4(2β + 8d) and, in the last case when λ 2 ≤ 1 and cos θ = 1, an inequality can be obtained in the form 0 ≤ 2β 2 + 2β + 4d 2 .

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

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 

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

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.

The stability has been proved for the linear differential equation representing any susceptible, exposed, asymptomatic, infectious, and recovered individuals. The convergence is given for the first two Eqs. (1) and (2), which are merged and form a single linear equation expressed as

Discretize Eq. (43) using a proposed scheme with unknown parameter in the following manner: 

This implies the expression for the unknown parameteř

Thus Eq. (44) is expressed in the form of

Let e n+1 2,i and e n+1 1,i be the difference between the exact solution and numerical solution for exposed and susceptible individuals computed at grid point i and at time level n. Thus the corresponding error equation for Eq. (47) is expressed as

whereǎ is given in Eq. (46) . Combining e n+1 2,i on the left side of Eq. (48) gives

(1 + 4dd 2ǎ + 2 tǎγ )e n+1 2,i = e n-1 2,i + 2 tǎ

Applying absolute values on both sides of Eq. (49) yields

Let e n = max{|e n 1,i |, |e n 2,i |}, then the inequality (50) can be expressed as

Equation (51) can be expressed as

where δ = |1 + 4dd 2ǎ + 2 tǎγ | -2 |ǎ|( |1+2dd 1 |d 2 + 2d|d 1 | ( x) 2 ), δ 1 = 1 + |ǎ||1 + 2dd 1 | and δ 2 = 2|ǎ||d 2 |d + 4|ǎ||d 1 |d.

Let δ > 0, then inequality (52) yields

Let e n-1 = max(e n-1 , e n ) and δ 1 +δ 2 δ = δ 3 , then Eq. (53) is given by

where C 1 is the coefficient of the leading remainder term.

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

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.

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.

Corona Covid-19 spread -a nonlinear modeling and simulation

SIRSi compartmental model for COVID-19 pandemic with immunity loss

Study of global dynamics of Covid-19 via a new mathematical model

A simple model for Covid-19

A primer on using mathematics to understand COVID-19 dynamics: modeling, analysis and simulations

SEAIR epidemic spreading model of COVID-19

Numerical analysis of diffusive susceptible-infected-recovered epidemic model in three space dimension

Structure Preserving Numerical Analysis of HIV and CD4+T-Cells Reaction Diffusion Model in Two Space Dimensions

Computational Modelling and Bifurcation Analysis of Reaction Diffusion Epidemic System with Modified Nonlinear Incidence Rate

Numerical modelling of an SIR epidemic model with diffusion

A numerical method for approximating' solutions to the functional equations arising in the epidemic model of Hoppensteadt and Waltman

Variational iteration and successive approximation methods for a SIR epidemic model with constant vaccination strategy

Numerical study of Sars epidemic model with the inclusion of diffusion in the system

SEIR modeling of the COVID-19 and its dynamics

Strong social distancing measures in the United States reduced the Covid-19 growth rate: study evaluates the impact of social distancing measures on the growth rate of confirmed Covid-19 cases across the United States

Clinical Features, and Short-term Outcomes of 102 Patients with Corona Virus Disease

Breaking down of healthcare system: mathematical modeling for controlling the novel coronavirus (2019-nCoV) outbreak in Wuhan

Statistics-based predictions of coronavirus epidemic spreading in mainland China

The global macroeconomic impacts of COVID-19: seven scenarios

Viability of intertwined supply networks: extending the supply chain resilience angles towards survivability. A position paper motivated by the Covid-19 outbreak

A minimal model of hospital patients' dynamics in Covid-19

Estimation of the final size of the coronavirus epidemic by the SIR model

Numerical study of an influenza epidemic model with diffusion

Numerical study of a diffusive epidemic model of influenza with variable transmission coefficient

A numerical study on an influenza epidemic model with vaccination and diffusion

Numerical and bifurcation analysis of spatio-temporal delay epidemic model

Treatment of HIV/AIDS epidemic model with vertical transmission by using evolutionary Padé-approximation

SEIR modeling of the Covid-19 and its dynamics

We thank the anonymous referees warmly. The first two authors are grateful to Vice-Chancellor, Air University, Islamabad, for providing an excellent research environment and facilities. The third and fourth author also thanks Prince Sultan University for funding this research work.

This research project is funded by the Research and initiative Centre, Prince Sultan University.

The manuscript included all required data and implementing information.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The authors have equally made contributions. All authors read and approved the final manuscript. 

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 17 February 2021 Accepted: 14 July 2021