A new unconditionally stable implicit numerical scheme for fractional diffusive epidemic model

1.   Introduction

When it comes to understanding the impact of infectious disease in a community, epidemiological studies are critical. Mathematical modeling involves checking models, estimating parameters, measuring sensitivity by varying parameters, and computing numerical simulations using models that have been constructed. This type of research is utilized to manage disease parameters and the ratio of illness spread in a population ^[1]. Infectious diseases are used to describe diseased models (i.e., the disease transferred from one person to another person). Infectious diseases include measles, rubella, chickenpox, mumps, aids, and gonorrhea syphilis ^[2].

Coronaviruses are responsible for transmitting severe acute respiratory syndrome (SARS). SARS is a valuable tool in the analysis of coronaviruses ^[3]. Coronaviruses are a group of viruses that can infect both humans and animals. There are several different types of coronaviruses. Coronaviruses in humans are primarily responsible for respiratory infections, ranging from mild colds to serious lung illnesses. They can also be accompanied by digestive diseases such as gastroenteritis, a type of diarrhea. To limit this pandemic, which is most likely still in its early stages, and prevent the collapse of healthcare systems, governments must take immediate and effective actions in a time of great uncertainty. A group of researchers studying the coronavirus family for over 30 years found that SARS and coronavirus have similar traits, such as biology and pathology ^[4]. Coronaviruses, RNA-enveloped viruses that infect people, mammals, and birds, spread worldwide. Coronaviruses are responsible for various respiratory, gastrointestinal, hepatic, and neurological disorders ^[5]. Coronaviruses are responsible for six different types of sickness in humans ^[6]. Chinese authorities declared a major outbreak of coronavirus disease 2019 (COVID-19) in the country in 2019, and this outbreak has the potential to spread around the world. Now underway, interventions and real-time data are required to control the coronavirus outbreak ^[7].

The COVID-19 pandemic is bringing quantitative mathematical modeling to the forefront of public debate more than ever before. Policymakers and the general public are turning to science, and modeling, in particular, to gain insight into the complex dynamics of the epidemic from both a local and global perspective and predict the consequences of possible interventions on the number of cases, hospitalizations, and deaths. The number of cases, hospitalizations, and deaths has increased dramatically. Bernoulli's study to analyze the efficiency of immunization against smallpox was the beginning of mathematical modeling in epidemiology, which began in 1760. Epidemiologists are building, testing and fine-tuning models to simulate the development of the COVID-19 infectious disease to gain a better knowledge of the disease and optimize measures to manage it.

Prior investigations have employed real-time data to explain the transmission of the virus from one person to another, the severity of the virus, and the pathogen's history during the first week of an outbreak. In December 2019, a group of Wuhan residents was taken to the hospital, where they were all diagnosed with idiopathic pneumonia. The majority of individuals believed that consuming wet market meat and shellfish was the root cause of pneumonia in children. With the assistance of the Wuhan local health authorities, the Chinese Center for Disease Control and Prevention (China CDC) launched an investigation into the etiology and epidemiology of the disease on December 31, 2019. Time-delay distributions, which included the admission date to the hospital and death, were used to assess whether the epidemiology was changing. According to the results of a clinical investigation conducted on the COVID-19 virus, signs of coronavirus infection develop seven days after the onset of illness ^[8]. The time interval between hospital admission and death is particularly crucial to minimize underestimating case fatality risk when calculating case fatality risk ^[9]. The fractional-order, which comprises integration and transects differentiation, is utilized better to understand real-world problems than the standard integer order and simulate genuine phenomena related to memory characterization and hereditary features using a fractional calculus ^[10,11]. Riemann Liouville is credited for introducing the concept of fractional derivative, which is based on power law. Caputo and Fabrizio present a novel fractional derivative that uses the exponential kernel, which is based on the work of Caputo and Fabrizio ^[12]. Previous results have proposed a novel fractional derivative with a non-singular kernel incorporating exponential and trigonometric functions ^{[13,14,15,16]}. This paper illustrates several new ways for epidemic models connected to this new derivative ^{[17,18,19,20,21,22,23,24,25]}. In Baleanu et al., significant results relating to this novel operator are established, and several examples are offered to illustrate the point ^[26].

The standard integer-order models do not benefit from the subsequent memory effects observed in many biological simulations. Instead, the ABC operator is used to teach the concept of heredity, which is an important property of many biological processes. Thus, when it comes to simulating biological processes, fractional operators are gaining more and more attention from various viewpoints.

For example, in ^[27], the authors proposed and studied the dynamics behavior of COVID-19, which was accomplished by analyzing a SEIHDR model based on the ABC operator. After analyzing the current results, they demonstrated the transmission dynamics of COVID-19. See ^{[28,29,30,31,32,33,34,35,36,37,38,39]} and the citations therein for further information on the current results.

It would be appropriate to include a population of COVID-19 disease-exposed individuals (E) in the model for the disease. In addition, it is possible to refine the description of specific epidemics by including assumptions about transmission from one class to another by including more compartments; see, for example, ^[40,41,42].

By introducing a generalized version of fractional models, D. Baleanu et al. ^[43] studied the effects of isolation and quarantine in the COVID-19 pandemic. The study considered the real clinical observations from China, and this was used to compute the basic reproduction number and determine the parameters. Asymptotic behavior of immunogenic tumor dynamics by using a new fractional model was studied in ^[44]. The equilibrium points and stability of the new model were discussed, the predictor-corrector numerical scheme was modified, and the results were compared with some real experimental data. Moreover, a tracking control method was employed to decrease the development of the tumor cell population. The fractional calculus theory was adopted in ^[45] to study the motion of a beam on an internally bent nanowire. The generalized fractional Lagrangian was introduced, and the fractional Euler-Lagrange equation was provided for the motion of the beam on the nanowire.

Solving the model with a numerical scheme pointed out that the fractional Euler-Lagrange equation provided a model that possessed more information than the classical description. A fractional COVID-19, SEIR epidemic model was solved in ^[46] using a predictor-corrector method. For the given time delay fractional differential equations, unique global solutions were derived under a mild Lipschitz condition using Mittag-Leffler functions, properties of a weighted norm, and Banach fixed point theorem. Real numerical data based on a case study of Wuhan, China, was used for graphical simulations. A new vaccinated model of COVID-19 ware provided ^[47] in the sense of new generalized Caputo-type and Caputo-Fabrizio fractional derivatives. The study of model dynamics such as boundedness, positivity, and local stability analysis were presented. Using fixed point theory, the existence of a unique solution to the fractional model was discussed. The numerical schemes for solving fractional differential equations have an advantage over exact analytical methods. Since the exact solution of all models in fractional calculus cannot be found analytically, some numerical approaches can be considered for solving those models. A fractional numerical algorithm has been proposed in ^[48] to get the exact solutions of generalized fractional-order differential equations based on the Caputo fractional derivative. Some simulations were also provided, including the solution of a co9mputer virus model to illustrate the applications of results. The existing treatments for the dangerous disease of COVID-19 are the use of the vaccine. From a mathematics aspect, an SEIRS dynamical model was proposed in ^[49] that comprised vaccine rate. A fractional model based on the Atangana-Baleanu derivative was formulated. For solving the model, a predictor-corrector scheme has been employed to solve the model. Stability analysis of the applied numerical method was also provided.

Since the non-standard finite difference method (NSFD) is one of the numerical methods used to find a positive solution, the method is unconditionally stable. But, its main drawback is its accuracy, which can be proved by applying Taylor series expansions. So, the solutions obtained by NSFD are not accurate when these are computed on a smaller number of grid points. NSFD is not first-order accurate in time for diffusive epidemic models but second-order accurate in space. Similarly, its fractional version can be constructed by considering the fractional Taylor series. But again, it will have an accuracy problem in time. It can be verified that the NSFD is also not suitable for finding the solution to the space variable. This contribution proposes a conditionally positivity preserving, unconditionally stable fractional numerical scheme, second-order in space and time. The scheme can be used to obtain conditions for getting a positive solution. These conditions depend on the considered mathematical model and the number of grid points and time levels. The stability of the numerical scheme can be checked by applying Fourier series analysis of Von Neumann stability criteria. This criterion is used for linear differential equations. Even though the equation is non-linear, its linearised form can be evaluated to estimate the stability criteria. So, this stability criterion is considered to check the stability of the proposed scheme. Due to the non-linear nature of the investigated diffusive epidemic model, the proposed fractional numerical approach is stable for the generalized kind of linear diffusive epidemic model.

The paper is organised in the following manner. A model of ordinary differential equations is established, and its local stability is given. After this, a numerical scheme is constructed for the general type fractional differential equation. Later on, the scheme is constructed for the fractional diffusive partial differential equation of susceptible individuals. Then the stability and convergence of the scheme are provided, and a comparison of the existing scheme for the classical parabolic equation is also provided. After this, the application of the scheme is given with simulations.

The SEAIR epidemic spreading model has been proposed in ^[50], but modification in that model is given here. Let $S$ be susceptible individuals, $E$ be the exposed individuals, $I$ be infected individuals, $Q$ be quarantined individuals, and $R$ be recovery people. Let $\bar{\alpha }$ be the rate at which susceptible become exposed, $\beta$ be the rate at which susceptible become infected, $\gamma$ be the rate at which exposed become infected, $\sigma$ be the rate at which infected become quarantined, $\bar{\sigma }$ be the rate at which infected individuals become healed or recovered, and $\mu$ be the rate at which quarantined become recovered. The mathematical model is given as:

Subject to the initial conditions

where ${S}_{0}, {E}_{0}, {I}_{0}, {Q}_{0}$ and ${R}_{0}$ can be any non-negative constants. The transmission diagram of the model (1)–(5) is given in Figure 1.

2.   Stability of the epidemic system

One of the equilibrium points of the system (1)–(5) is $(\mathrm{1, 0}, \mathrm{0, 0}, 0)$. The next linear stability procedure is given about the mentioned equilibrium point. The linear stability procedure starts with the Jacobian matrix for the system of Eqs (1)–(5) in a form

Using the equilibrium point $\left(\mathrm{1, 0}, \mathrm{0, 0}, 0\right)$ Jacobian (7) is expressed as

To find the characteristic polynomial, consider the equation

where $I.D$ is an identity matrix of order $5\times 5$. Equation (9) can be expressed as

This implies $\lambda = \mathrm{0, 0}, -\mu$ are three Eigenvalues, and the remaining two Eigenvalues can be found from the following characteristic polynomial

Now according to Ruth Hurwitz's criteria, the second-degree polynomials have both roots in the open left half-plane if and only if

So the system will be stable if conditions (14) are met.

Since the ordinary differential equations model only describes the changes in individuals over time, which does not discuss the variations of susceptible, exposed, infected, quarantined and recovered individuals over space. For space variation, both convection and diffusion effects can be included. Still, only diffusive effects are considered with each category of individuals in this study so that the resulting differential equations become partial differential equations (PDEs) of parabolic type. These partial differential equations are more general than the corresponding model of ordinary differential equations (ODEs). If the coefficients of diffusions are chosen to be zero, the model of PDEs becomes the model of epidemic disease comprised of ODEs. Also, a time-fractional differential equations model is considered instead of the classical model of the differential equation.

A numerical scheme has been proposed ^[51] for a fractional diffusive SEAIR model with a non-linear incidence rate. A numerical scheme will be proposed for the fraction diffusive SEIQR model. For constructing a fractional diffusive model, the following fractional partial differential equations model is presented by incorporating fractional Caputo time derivative and diffusion in the system of Eqs (1)–(5).

The fractional diffusive epidemic model for COVID-19 is given as,

In the first stage, a fractional numerical scheme will be constructed for the general type of parabolic equation and one of the equations in the model (15)–(19). Later on, stability and convergence will be provided. Further, it is assumed that the solution of model (15)–(19) exists and it is unique.

3.   Construction of a fractional numerical scheme

Before constructing the proposed scheme for the considered model, it is constructed for the general form of the fractional time-dependent parabolic equation of the form,

where $\stackrel{ˇ}{d}and\stackrel{ˇ}{c}$ are constants.

Let

where ${f}_{1}\left(u\right)$ and ${f}_{2}\left(u\right)$ are positive functions of $u$. In the application procedure of the proposed fractional scheme, a difference equation is constructed first in the form,

The next step is to use the fractional Taylor series for the term ${u}_{i}^{n+1}$ in the form,

Substituting the fractional Taylor series (23) into Eq (22) and rewriting the resulting equation in the form,

Comparing coefficients of ${u}_{i}^{n}$, ${\left(\frac{{\partial }~~^{\alpha }u}{\partial {t}~~^{\alpha }}\right)}_{i}^{n}$ and ${\left(\frac{{\partial }^{2\alpha }u}{\partial {t}^{2\alpha }}\right)}_{i}^{n}$ on both sides of Eq (24) yields,

Solving Eqs (25) and (26) yields the values of $\stackrel{ˇ}{a}$ and $\stackrel{ˇ}{b}$ and the difference equation using the values of $\stackrel{ˇ}{a}$ and $\stackrel{ˇ}{b}$ will be the discretized equation for the second-order fractional implicit scheme for solving Eq (20).

4.   Construction of numerical scheme for epidemic disease model

To construct the fractional numerical scheme for the model (15)–(19), one of the equations in the model will be chosen. For doing so, consider the Eq (15) with discretization is given as,

In Eq (27), $‘‘a’’$ and $‘‘b’’$ are unknown to be determined. Their values will be found by applying Taylor series expansion in the fractional derivative. So, for doing so, consider the following Taylor series in fractional derivatives,

Substituting Taylor series expansion (28) into (27), it is obtained

Comparing coefficients of ${S}_{i}^{n}, {D}_{t}~~^{\alpha }{S}_{i}^{n}$ and ${D}_{t}^{2\alpha }{S}_{i}^{n}$ on both sides of Eq (29) yields the following two equations in two unknowns,

Equation (30) can be simplified as

Upon solving Eqs (31) and (32), values for $"a"$ and $"b"$ can be found. After finding $"a"$ and $"b"$, the discretized Eq (27) can be used to find the unknown ${S}_{i}^{n+1}$ and this value is found by the use of an additional iterative method.

The positivity of the numerical solution obtained by applying the proposed scheme for the first equation in a model (15)–(19) depends on two conditions. First, it depends on found values of parameters $"a"$ and $"b"$ by solving Eqs (31) and (32), and second, it depends on those values obtained by the remaining difference equations discretized by the proposed scheme to the rest of the equations in the model (15)–(19). If all the parameters in the difference equation are positive, the positivity would be guaranteed. Since there are conditions for a positive solution, the scheme is a conditionally positivity-preserving scheme. If the parameters $a$ and $b$ are positive, then by looking at the difference Eq (27) and collecting the terms of ${S}_{i}^{n+1}$ on the left-hand side of the equation, the resulting equation will produce a positive solution because each term in the resulting equation will be positive with the assumption of choosing positive initial conditions, i.e.,

where $\theta = 1+\frac{a{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left\{\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{E}_{i}^{n+1}+\beta {I}_{i}^{n+1}\right\}+\frac{b{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left\{\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{E}_{i}^{n}+\beta {I}_{i}^{n}\right\}$,

provided that $a > 0~~and~~b > 0$, exposed, infected, quarantined, and recovered individuals at time level $n+1$ and grid point $i$ are non-negative. The positivity of other variables depends on the positivity of included parameters in their difference equations.

5.   Stability analysis of proposed numerical scheme

To find the stability condition consider the first general type of linear fractional partial differential equation of the form,

for $0 < \alpha \le 1$. The stability analysis can be found only for parabolic equations without diffusion, and in that case, one may consider the case when,

To discretize parabolic fractional equation with souse term, Eq (33) using the proposed fractional scheme, the difference equation is expressed as:

Using fractional Taylor series expansion (23) for ${u}_{i}^{n+1}$ and substituting it into Eq (35) and rewriting the resulting in the following manner,

Comparing coefficients of ${u}_{i}^{n}, {D}_{t}~~^{\alpha }{u}_{i}^{n}$ and ${D}_{t}^{2\alpha }{u}_{i}^{n}$ on both sides of Eq (36), it is obtained

The expressions for $"a"$ and $"b"$ can be found by solving Eqs (37) and (38).

Theorem. Proposed fractional numerical scheme is unconditionally stable.

Proof. Consider Eq (33) and the following transformations

Substituting transformations (39) into Eq (35), it is obtained

Dividing both sides of Eq (40) by ${e}^{iI\psi }$ gives,

Using trigonometric identities in Eq (41) and collecting the coefficients of ${\widehat{E}}^{n+1}$ of the left side, it is obtained

Let $\stackrel{~}{d} = \frac{{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\frac{\bar{d}}{{\left(\mathrm{\Delta }x\right)}^{2}}$, then Eq (42) can be expressed as,

The amplification factor can be expressed as,

The stability condition can be expressed as,

Here special cases will be discussed with the specific value of $\alpha$.

Let $\alpha = 0.9$, then the stability condition (45) can be expressed as,

where $a\approx \frac{0.889627}{1.61239-1.50299\bar{d}\frac{{\left(\mathrm{\Delta }t\right)}~~~^{0.9}}{{\left(\mathrm{\Delta }x\right)}~~^{2}}-0.751497\bar{\beta }{\left(\mathrm{\Delta }t\right)}~~~^{0.9}}$ and $b\approx \frac{0.722764}{1.61239-1.50299\bar{d}\frac{{\left(\mathrm{\Delta }t\right)}~~~^{0.9}}{{\left(\mathrm{\Delta }x\right)}~~^{2}}-0.751497\bar{\beta }{\left(\mathrm{\Delta }t\right)}~~~^{0.9}}$.

Condition (46) can be expressed as

Clearly, for any value of $cos\psi$ inequality (47) holds. Similarly, the stability can be checked for any fixed value of $\alpha,$ and it can be concluded that the proposed fractional scheme is unconditionally stable for any value of $\alpha$ in the interval $\left(\mathrm{0, 1}\right].$

4.1.   Convergence of proposed numerical scheme

Theorem. The proposed fractional numerical scheme converges for non-linear epidemic fractional diffusive equation (15) if the following condition holds,

where $\delta = 1+a\frac{{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left(\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{H}_{1}^{1}+\beta {H}_{2}^{1}\right)+b\frac{{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left(\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{H}_{3}^{1}+\beta {H}_{4}^{1}\right)$,

where ${F}_{1} = {F}_{1}\left({S}_{i}^{n+1}, {E}_{i}^{n+1}\right)$ , ${G}_{1} = {G}_{1}\left({S}_{i}^{n+1}, {I}_{i}^{n+1}\right), {F}_{2} = {F}_{2}\left({S}_{i}^{n+1}, {E}_{i}^{n}\right)\&{G}_{2} = {G}_{2}\left({S}_{i}^{n+1}, {I}_{i}^{n}\right)$ .

Proof. The discretized form of Eq (15) can be expressed as,

The expression for $"a"$ and $"b"$ can be found by comparing the coefficients of fractional Taylor series expansions, and their values can be found by solving Eqs (31) and (32).

Let the exact difference equation for Eq (15) be expressed as,

Consider,

By using the Mean value Theorem, the following relationship can be established

where ${C}_{1} = {\bar{S}}_{i}^{n+1}+{\varepsilon }_{1}{e}_{1, i}^{n+1}$, ${C}_{2} = {\bar{E}}_{i}^{n+1}+{\varepsilon }_{2}{e}_{2, i}^{n+1}$, ${C}_{3} = {\bar{S}}_{i}^{n+1}+{\varepsilon }_{3}{e}_{1, i}^{n+1}$, ${C}_{4} = {\bar{I}}_{i}^{n+1}+{\varepsilon }_{4}{e}_{3, i}^{n+1}$, ${C}_{5} = {\bar{S}}_{i}^{n+1}+{\varepsilon }_{5}{e}_{1, i}^{n+1}$, ${C}_{6} = {\bar{E}}_{i}^{n}+{\varepsilon }_{6}{e}_{2, i}^{n}$, ${C}_{7} = {\bar{S}}_{i}^{n+1}+{\varepsilon }_{7}{e}_{1, i}^{n+1}$, ${C}_{8} = {\bar{I}}_{i}^{n}+{\varepsilon }_{8}{e}_{3, i}^{n}$

Let ${\boldsymbol{e}}_{1, i}^{n+1} = {\left[{e}_{1, i}^{n+1}, {e}_{2, i}^{n+1}\right]}^{t}, {\boldsymbol{e}}_{2, i}^{n+1} = {\left[{e}_{1, i}^{n+1}, {e}_{3, i}^{n+1}\right]}^{t}, {\boldsymbol{e}}_{3, i}^{n+\frac{1}{2}} = {\left[{e}_{1, i}^{n+1}, {e}_{2, i}^{n}\right]}^{t}, {\boldsymbol{e}}_{4, i}^{n+\frac{1}{2}} = {\left[{e}_{1, i}^{n+1}, {e}_{3, i}^{n}\right]}^{t}$.

Subtracting Eq (49) from Eq (48) yields an equation of the form

where $\nabla {F}_{1} = {\left[{H}_{1}^{1}, {H}_{1}^{2}\right]}^{t}, \nabla {G}_{1} = {\left[{H}_{2}^{1}, {H}_{2}^{2}\right]}^{t}, \nabla {F}_{2} = {\left[{H}_{3}^{1}, {H}_{3}^{2}\right]}^{t}, \nabla {G}_{2} = {\left[{H}_{4}^{1}, {H}_{4}^{2}\right]}^{t}$.

Collecting the coefficients of ${e}_{1, i}^{n+1}$ on the left-hand side of Eq (54), it is obtained

Let ${e}^{n+1} = \mathrm{max}\left\{\underset{1\le i\le N}{\mathrm{max}}\left|{e}_{1, i}^{n+1}\right|, \underset{1\le i\le N}{\mathrm{max}}\left|{e}_{2, i}^{n+1}\right|, \underset{1\le i\le N}{\mathrm{max}}\left|{e}_{3, i}^{n+1}\right|\right\}$ and applying absolute on both sides of Eq (55), it is obtained

where $\delta = 1+a\frac{{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left(\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{H}_{1}^{1}+\beta {H}_{2}^{1}\right)+b\frac{{\left(\mathrm{\Delta }t\right)}~~^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}\left(\frac{2{d}_{1}}{{\left(\mathrm{\Delta }x\right)}^{2}}+\bar{\alpha }{H}_{3}^{1}+\beta {H}_{4}^{1}\right)$.

Collecting coefficients of ${e}^{n+1}$, it is obtained

Let,

Then, inequality (57) can be expressed as,

Let $n = 0$ in (58), it is obtained that

Since initial condition is exact, so ${e}^{0} = 0$, so Eq (59) can be expressed as,

Let $n = 1$ in (58), which yields

If this continued, then

Since the series $1+{\delta }_{3}+{\delta }_{3}^{2}+\dots +{\delta }_{3}^{n-1}+\dots$ is an infinite geometric series that will converge if $\left|{\delta }_{3}\right| < 1$. Therefore the convergence condition is $\left|\frac{{\delta }_{2}}{{\delta }_{1}}\right| < 1$. Similarly, the convergence condition can be found for the remaining equations in the system (16)–(18).

4.2.   Comparison

The proposed scheme results can be compared with results obtained by an existing scheme. One of the existing schemes considered for the epidemic model is the non-standard finite difference method. The fractional non-standard finite difference method is given for fractional diffusive epidemic Eq (15).

The fractional non-standard finite difference method for Eq (15) is given as,

The explicit form of the non-standard fractional scheme is expressed as,

The main advantage of applying the fractional non-standard finite difference method (NSFD) is the stability and positivity of the solution. But the main drawback is its accuracy. It's not even first-order accurate in both cases when $\alpha = 1$ and $\alpha \ne 1$, its lack of accuracy can be proved by applying standard and fractional Taylor series expansions. So, for the diffusive epidemic model, it can be proved that NSFD is neither first-order accurate nor consistent. So, due to these drawbacks, the obtained solution is also not accurate compared with the standard first-order method. This contribution also compared the proposed scheme with the fractional NSFD scheme and the first-order fractional Backward Euler method. The fractional Backward Euler method can be used for those cases when the solution is positive, and the method stays stable. The fractional Backward Euler method for Eq (15) is given as:

Now the drawback of using fractional NSFD is discussed. To become the method first-order accurate, it should satisfy first-order fractional Taylor series expansion. For this reason, expand ${S}_{i}^{n+1}$ in fractional Taylor series form for Eq (62) as,

Comparing coefficients of ${S}_{i}^{n}, {D}_{t}~~^{\alpha }{S}_{i}^{n}$ and ${D}_{t}^{2\alpha }{S}_{i}^{n}$ on both sides of Eq. (65), it is obtained

From Eq (66), it is observed that the left-hand side is not equal to the right-hand side, so it means that coefficients of ${D}_{t}~~^{\alpha }{S}_{i}^{n}$ on both sides are not canceled with each other, or reminder term has a fractional derivative of the form ${D}_{t}~~^{\alpha }{S}_{i}^{n}$ which means that fractional NSFD is not first-order accurate.

4.3.   Order of convergence

A parabolic equation with a classical derivative is tested to find the order of convergence of the proposed scheme. For doing so, consider the following parabolic equation

subject to the boundary conditions,

and using an initial condition,

Equation (67) using boundary conditions (68) and initial condition (69) are solved with three different numerical schemes. The order of convergence of the three schemes is shown in Table 1. Theoretically, the employed schemes for Eqs (67)–(69) are second-order accurate as proved by applying the Taylor series, but according to Table 1, their convergence order is nearly one.

One of the advantages of using the Non-standard finite difference method is its consumption of timeless than the proposed scheme. But it has the drawback of lacking accuracy. Since the considered NSFD is an explicit method, it takes less time than the proposed implicit scheme if the iterations are carried out by Gauss-Seidel iterative method for linear PDEs. So, the time consumption of any implicit scheme depends on the considered way of solving the difference equations obtained by applying the implicit scheme to some particular type of differential equation. The computation time of the implicit scheme can be reduced by using second-order Newton's method for non-linear problems. For linear problems, computation time can be reduced by directly solving the matrix-vector form equation instead of using another iterative solver. Matlab has the facility to solve the system of equations, so the difference equation in matrix-vector form can be solved directly by employing that Matlab facility.

For the accuracy and computation time of the proposed and non-standard finite difference schemes, their comparison is made in solving Eqs (67)–(69). Figure 2 compares two schemes using 50 grid points and 21-time levels. The time consumed by NSFD is 0.596988 seconds, and the proposed scheme consumed 1.9425 seconds to perfume the same task using the Gauss-Seidel iterative method. So due to the considered iterative method, the proposed scheme consumed the mentioned time, which may have been reduced by adopting a Matlab solver for solving matrix-vector equations for the linear PDE.

6.   Results and discussions

The fractional proposed scheme with an accuracy of two in both time and space has been implemented, and the results are compared with existing classical and fractional schemes. As mentioned in the previous section, NSFD is not even first-order accurate. Its comparison is made with the second-order proposed and first-order Euler methods. Figure 3 compares the proposed scheme with the Backward Euler method for faster convergence. The maximum norm considered for drawing Figure 3 is the difference of solutions computed over two consecutive iterations. The utilization of the proposed scheme and an iterative method on the difference equation (27) yields

So, the solution component ${S}_{i}^{n+1, k+1}$ on each grid point and at each time level is found iteratively. The iterative method will be stopped if the required criterion is met. The norm between consecutive two iterations is computed as

where $tol$ is to a very small number nearly equal to zero. In inequality (71), the stopping criterion of susceptible individuals is given. Similarly, stopping criteria for each category of people can be expressed, and their maximum values are computed for each iteration. Figure 3 shows the maximum of all norms on each iteration for both the proposed and Backward Euler methods. Figure 3 shows that the proposed scheme converges faster than the Backward Euler method with the particular value of $\alpha$. Figures 4 and 5 offer the comparison of three schemes over time for the classical case $\alpha = 1$. The difference between results obtained by three schemes can be seen in these Figures 4 and 5. Since it was mentioned that NSFD is not first-order accurate, its drawback in finding an accurate solution can be seen in these two figures. The results were computed for 50 grid points and 50-time levels. The solution obtained by NSFD will be worse if results are computed using a smaller number of grid points and the smaller number of time levels. Figures 6 and 7 also compare three schemes, but this comparison is made over spatial variable $x$. The solution for both NSFD and the Backward Euler method is near because NSFD is second-order accurate in space. Still, the proposed scheme produced better results verified by applying another second-order scheme, the Crank-Nicolson method. But, the remaining categories of the individual solution obtained by the proposed scheme are near to the solution obtained by the Backward Euler method. Figures 8 and 9 show a comparison of three fractional schemes. Similar to the classical case, the fractional results are also not accurately obtained by fractional NSFD because, again, it is not first-order accurate using the fractional Taylor series approach.

7.   Conclusions

A mathematical model of the COVID-19 epidemic disease has been constructed in the form of a system of linear and non-linear differential equations. Its linear stability conditions are found by considering Ruth-Hurwitz criteria, and later on, a time-fractional diffusive epidemic model was constructed. The fractional numerical scheme has been developed by considering the fractional Taylor series. Its stability and convergence conditions have been established, and it has been compared to the present numerical scheme. It was proved that the existing non-standard scheme is not first-order accurate from Taylor series analysis and drawn obtained solution for both integer and fractional derivative cases.

On the other hand, the numerical results give compelling evidence for the advantages of a proposed approach over non-standard approaches. The proposed numerical scheme can be applied to epidemic models to get a conditionally positivity solution with unconditional stability. The proposed numerical scheme had advantages over some existing fractional numerical schemes that provided the solution for the classical case. Also, the proposed numerical scheme can be further applied in epidemiological disease models and other problems ^[52,53] in fractional calculus.

Acknowledgments

The authors wish to express their gratitude to Prince Sultan University for facilitating the publication of this article through the Theoretical and Applied Sciences Lab. This work was supported by the research grants Seed project; SEED-2022-CHS-100.

Conflict of interest

The authors declare no conflict of interest.