An Efficient Numerical Simulation for the Fractional COVID-19 Model Using the GRK4M Together with the Fractional FDM

Abstract

In this research, we studied a mathematical model formulated with six fractional differential equations to characterize a COVID-19 outbreak. For the past two years, the disease transmission has been increasing all over the world. We included the considerations of people with infections who were both asymptomatic and symptomatic as well as the fact that an individual who has been exposed is either quarantined or moved to one of the diseased classes, with the chance that a susceptible individual could also migrate to the quarantined class. The suggested model is solved numerically by implementing the generalized Runge–Kutta method of the fourth order (GRK4M). We discuss the stability analysis of the GRK4M as a general study. The acquired findings are compared with those obtained using the fractional finite difference method (FDM), where we used the Grünwald–Letnikov approach to discretize the fractional differentiation operator. The FDM is mostly reliant on correctly converting the suggested model into a system of algebraic equations. By applying the proposed methods, the numerical results reveal that these methods are straightforward to apply and computationally very effective at presenting a numerical simulation of the behavior of all components of the model under study.

1. Introduction

Many models have been used to explain a wide range of biological phenomena [1,2] with the first appearance of COVID-19 [3,4]. Mathematical models are constructed as a strategy for gaining insight into the pandemic’s mode of effect, transmission, spread, prevention, and control, and the impact of preventive measures such as hand-washing with a disinfectant such as hand sanitizer, increasing social distancing, and wearing face masks [5,6]. In [7,8], the authors looked at the early stages of analyzing COVID-19’s distribution in Nigeria. However, a more thorough examination of the most recent models as well as a justifiable and satisfactory assessment are required. The proposed model is considered by many researchers to have different forms, i.e., it is different in the number of equations and in the nature of the given derivatives (ordinary/fractional). Additionally, it is solved numerically by using some numerical methods, and the existence, the convergence, and the stability were given through some articles. For more details about the proposed model, see [9,10,11,12].

Most of the models described in the abovementioned studies were based on regular derivatives, with some limitations in the order of the involved DEs. To circumvent these constraints, numerous academics have turned to fractional calculus, which is a rapidly growing branch of mathematics. Non-integer or fractional order differential operators are employed in fractional calculus. They have many properties and are important in demonstrating many scientific phenomena and facts with nonlocal dynamical behavior. The fact that these operators’ kernels are nonlocal and nonsingular is their primary advantages. Many studies containing fractional and variable order operators were published to describe many daily problems [13,14].

Many numerical and approximate methods are implemented to solve numerically the COVID-19 in different forms; for example, in [15] published in 2021, the authors used the Atangana–Toufik numerical scheme to numerically solve the proposed COVID-19 model and used the least square curve fitting method to obtain the best values for some of the unknown biological parameters involved in the proposed model; in [16], the authors applied the FDM to study the COVID-19 model; and in [17], the authors applied the nonstandard FDM to study the COVID-19 model.

The main aim and the novelty of this work are to extend the COVID-19 system, which is given in [15], to a fractional-order model. To achieve this aim, the numerical solution for the proposed model is obtained by using the GRK4M and the fractional FDM with the approach of Grünwald-Letnikov. The FDM reduces the presented problem to a system of algebraic equations, substantially simplifying it and subsequently optimizing it to obtain the unknown coefficients and, as a result, the solution. We examine the disease-free equilibrium’s locally asymptotically-stability, the unique endemic equilibrium point, the locally asymptotically-stability of its unique endemic equilibrium point, and the globally asymptotically-stability of its endemic equilibrium to provide a qualitative analysis of the proposed model. Finally, we give comparative studies. We highlight the main limitation in this study, which is due to the lack of real data, so we study the model from the theoretical side with the help of its corresponding system of ODE.

The following is how the paper is structured: Section 2 contains the model formulation as well as an explanation of the model’s parameters. Section 3 contains some preliminary information and notations on the proposed methods. We provide a qualitative examination of the suggested model in Section 4. In Section 5, we numerically solve the proposed problem by using the GRK4M and the fractional FDM and discuss the stability analysis of the GRK4M. In addition, in Section 6, the numerical simulation and the validation of the GRK4M and the fractional FDM of the model are presented. Finally, in Section 7, we present the conclusion and the planned future work.

2. Model Construction

This study presents a COVID-19 model that uses fundamental transmission rates. Let $N\left(t\right)$ be the whole human population. Susceptible individuals $S\left(t\right)$, exposed individuals $E\left(t\right)$, asymptotically infected individuals $I_A\left(t\right)$, symptomatic infected individuals $I_S\left(t\right)$, quarantined persons $Q\left(t\right)$, and individuals who have recovered/removed from COVID-19 $R\left(t\right)$ are the seven classifications in this population. Taking this into account, the total population is [15]

$$N\left(t\right)=S\left(t\right)+E\left(t\right)+Q\left(t\right)+I_A\left(t\right)+I_S\left(t\right)+R\left(t\right).$$

Assume that $\Lambda$ and $\mu$ are the natural human natality and death rates, respectively. $\left(S\right)$ becomes infected after contact with $\left(E\right)$ or migration to $\left(Q\right)$ at a rate of $\beta$. $\left(E\right)$ may be the first to be moved to $\left(Q\right)$ or they may get sick without displaying any symptoms $\left(I_A\right)$ or $\left(I_S\right)$ with rates of $\gamma ,\sigma$, and $\eta$, respectively. Individuals who have been quarantined (Q) may also be confirmed infected through a test with symptoms $\left(I_S\right)$ or without $\left(I_A\right)$ with rates v and $\theta$, respectively. The persons in $\left(I_A\right)$ may recover at the rate $r_1$ and the persons in $\left(I_S\right)$ may recover at the rate of $r_2$ [15].

In the following, some notations will be useful to formulate the studied model:

• $\tau$: The number of people who are transferred to quarantine from those who are susceptible;
• $\beta$: Contact rate between people who are susceptible and those who are exposed;
• $\delta$: Coronavirus-related death rate in the symptomatically infected individual class;
• $\gamma$: Rate of transfer of exposed people to quarantine;
• $\eta$: Rate of transfer of persons from the exposed to the symptomatically infected classes;
• $\theta$: The proportion of quarantined people who are asymptomatically infected;
• $\mu$: Natural death rate;
• $\nu$: Rate of transfer quarantined people to the symptomatically infected people class;
• $\sigma$: Rate of exposed people to asymptomatically infected people class;
• $\Lambda$: Recruitment (natality) rate;
• $r_1$: Asymptomatically infected people’s recovery rate;
• $r_2$: Symptomatically infected people’s recovery rate.

Each of these groups may drop as a consequence of natural death $\mu$, while $\left(I_S\right)$ may decrease as a result of disease death at the rate of $\delta$. The mortality rate as a result of the disease is not regarded in $I_A$. This model does not account for the risk of reinfection following recovery.

The schematic diagram in Figure 1 is used to build a system, which is discussed below:

$$\begin{array}{l}\frac{dS\left(t\right)}{dt}=\Lambda -(\tau +\mu )S\left(t\right)-\beta S\left(t\right)E\left(t\right),\\ \frac{dE\left(t\right)}{dt}=\beta S\left(t\right)E\left(t\right)-(\gamma +\mu +\eta +\sigma )E\left(t\right),\\ \frac{dQ\left(t\right)}{dt}=\tau S\left(t\right)+\gamma E\left(t\right)-(\mu +v+\theta )Q\left(t\right),\\ \frac{d{I}_A\left(t\right)}{dt}=\sigma E\left(t\right)+\theta Q\left(t\right)-\left(\mu +r_1\right)I_A\left(t\right),\\ \frac{d{I}_S\left(t\right)}{dt}=\eta E\left(t\right)+vQ\left(t\right)-\left(\delta +\mu +r_2\right)I_S\left(t\right),\\ \frac{dR\left(t\right)}{dt}=r_1{I}_A\left(t\right)+r_2{I}_S\left(t\right)-\mu R\left(t\right).\end{array}$$

(1)

A system of fractional differential equations (FDEs) is offered as a mathematical model to characterize the dynamics of propagation for COVID-19 ($\forall t\ge 0$):

$$\begin{array}{l}{D}^{\alpha }S\left(t\right)=\Lambda -(\tau +\mu )S\left(t\right)-\beta S\left(t\right)E\left(t\right),\\ D^{\alpha }E\left(t\right)=\beta S\left(t\right)E\left(t\right)-(\gamma +\mu +\eta +\sigma )E\left(t\right),\\ D^{\alpha }Q\left(t\right)=\tau S\left(t\right)+\gamma E\left(t\right)-(\mu +v+\theta )Q\left(t\right),\\ D^{\alpha }{I}_A\left(t\right)=\sigma E\left(t\right)+\theta Q\left(t\right)-\left(\mu +r_1\right)I_A\left(t\right),\\ D^{\alpha }{I}_S\left(t\right)=\eta E\left(t\right)+vQ\left(t\right)-\left(\delta +\mu +r_2\right)I_S\left(t\right),\\ D^{\alpha }R\left(t\right)=r_1{I}_A\left(t\right)+r_2{I}_S\left(t\right)-\mu R\left(t\right),\end{array}$$

(2)

where $\alpha$ refers to the order of the Liouville–Caputo fractional derivative.

Figure 1 shows the employed coefficients (parameters) in the model together. Additionally, the initial conditions are as follows:

$$S\left(0\right)=S_0,E\left(0\right)=E_0,Q\left(0\right)=Q_0,I_A\left(0\right)=I{A}_0,I_S\left(0\right)=I{S}_0,R\left(0\right)=R_0.$$

(3)

With this model (2) in its fractional form, we can more precisely characterize the effect of the spread of the pandemic of COVID-19 in the future and in history based on the memory effect of fractional derivatives. As it is also known, although the mathematical models with integer derivatives play an important role and have their significance to understand the dynamics of epidemiological systems, they have some limitations, as these systems do not have memory or non-local effects; therefore, these models in the usual form are sometimes inappropriate, and thus, it is necessary to convert several epidemiological models into FDEs to have an effective ability to explore many natural phenomena, facts that have non-local dynamic behavior. As a general case, FDEs are frequently used in the study of anomalous phenomena in nature and in the theory of complex systems as well as when taking into account the properties of the curve over a large extent. Finally, it can explain the time delay, fractal properties, etc. For more details about this model, see [18,19].

3. Preliminaries and Notations

In this part, we go over some of the fundamental concepts and useful properties that are used in this work.

Definition 1.

The Liouville–Caputo and the Riemann–Liouville derivatives of order α are, respectively, defined by [20,21]:

$$D^{\alpha }f\left(t\right)=\frac{1}{\Gamma (m-\alpha )}{\int }_0^t\frac{f^{\left(m\right)}\left(\xi \right)}{{(t-\xi )}^{\alpha -m+1}}d\xi ,m-1\le \alpha \le m,m\in \mathbb{N}.$$

$$D_R^{\alpha }f\left(t\right)=\frac{1}{\Gamma (m-\alpha )}\frac{d^m}{d{t}^m}{\int }_0^t\frac{f\left(\xi \right)}{{(t-\xi )}^{\alpha -m+1}}d\xi ,m-1\le \alpha \le m,m\in \mathbb{N}.$$

Definition 2.

The shifted Grünwald–Letnikov fractional derivative is defined as follows [22]:

$$\frac{d^{\alpha }y\left(t\right)}{d{t}^{\alpha }}=\underset{h⟶0}{lim}\frac{1}{h^{\alpha }}\sum _{i=0}^{\left[\frac{t}{h}\right]+1}{(-1)}^i\left(\begin{array}{c}\alpha \\ i\end{array}\right)y(t-(i-1)h),$$

(4)

for some constant h.

Lemma 1.

For $\mid t\mid <\rho$, with arbitrary constant, assume that ρ$y\left(t\right)$ can be written as a power series. The Grünwald–Letnikov formula holds for each $0<\alpha <\rho$ and a series of step size h with $\frac{\tau }{h}\in \mathbb{N}$,

$$D_R^{\alpha }y\left(\tau \right)=\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }y\left(nh\right)+O\left(h\right),(h⟶0),$$

(5)

where $D_R^{\alpha }$ refers to the Riemann–Liouville definition of the fractional derivative and

$${\Delta }_h^{\alpha }y\left(nh\right)=\sum _{i=0}^n{(-1)}^i\left(\begin{array}{c} \alpha \\ i \end{array}\right)y\left(t_{n-i}\right).$$

Therefore, by using the relation between both definitions (Liouville–Caputo and the Riemann–Liouville) [23] for $0<\alpha \le 1$, we can find the following in the case of the Liouville–Caputo sense:

$$D^{\alpha }y\left(\tau \right)=\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }y\left(nh\right)-\frac{{\tau }^{-\alpha }}{\Gamma (1-\alpha )}y\left(0\right)+O\left(h\right),(h⟶0).$$

(6)

4. Qualitative Analysis of the Proposed Model

As we can see in [24], one of the most important parameters in studying the infectious disease model is the basic reproduction number ${\mathfrak{R}}_0=1$, which is defined and derived by the following formula [15,24]:

$${\mathfrak{R}}_0=\frac{\beta \Lambda }{(\gamma +\mu +\eta +\sigma )(\tau +\mu )}>0.$$

Additionally, in this section, we present the invariant region of (1).

4.1. Invariant Region

It is very important to note that $S\left(t\right),E\left(t\right),Q\left(t\right),I_A\left(t\right),I_S\left(t\right),R\left(t\right)$ are nonnegative for all $t\ge 0$ in the model (1). For all $t\ge 0$, solution via initial positive data remains bounded and positive. From system (1), we can see that $\frac{dN\left(t\right)}{dt}=\Lambda -\mu N\left(t\right)-\delta I_S\left(t\right)$ and ${sup}_{t\to +\infty }N\left(t\right)⩽\frac{\Lambda }{\mu }.$ The system (1) will be studied in the feasible region:

$$\Omega =\left\{\left(S\left(t\right),E\left(t\right),Q\left(t\right),I_A\left(t\right),I_S\left(t\right),R\left(t\right)\right)\in {\mathbb{R}}_{+}^6:0⩽N\left(t\right)⩽\frac{\Lambda }{\mu }\right\}.$$

(7)

Equation (7) is now positive invariant in relation to (1), which means that (1) is epidemiologically well-posed.

The disease-free equilibrium (DFE) is obtained by letting $E=Q=I_A=I_S=R=0$. Therefore, the equations in (1) indicate that this point is given by the following:

$${\overline{\mathbf{E}}}_0=\left(S_0,0,0,0,0,0\right)=\left(\frac{\Lambda }{\tau +\mu },0,0,0,0,0\right).$$

4.2. Existence of Endemic Equilibrium Point

Let us denote the endemic equilibrium by ${\overline{\mathbf{E}}}_1=$$\left(S^{*},E^{*},Q^{*},I_A^{*},I_S^{*},R^{*}\right),$ where the values of all the components $S^{*},E^{*},Q^{*},I_A^{*},I_S^{*},R^{*}$ can be found in [15].

Theorem 1.

${\overline{\mathit{E}}}_0$ is locally asymptotically stable if ${\mathfrak{R}}_0<1$.

Theorem 2.

The system (1) has a unique ${\overline{\mathit{E}}}_1$.

Theorem 3.

${\overline{\mathit{E}}}_1$ is locally asymptotically stable if ${\mathfrak{R}}_0>1$.

Theorem 4.

The system (1) has no periodic orbits.

Theorem 5.

${\overline{\mathit{E}}}_1$ for (1) is globally asymptotically stable whenever ${\mathfrak{R}}_0>1$.

The definition of ${\overline{\mathbf{E}}}_1$ and the proof of the five Theorems 1–5 can be given in [15].

5. Numerical Studies

In this section, we use GRK4M and fractional FDM to solve the suggested model numerically. In addition, as part of general research, we cover the stability analysis of the GRK4M.

5.1. Procedure of Solution by the GRK4M

To show how to implement GRK4M in solving the proposed model (2), we consider the following general form of the FDE [25]:

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)=\psi (t;u\left(t\right)),0<\alpha \le 1,0<t⩽T,\\ u\left(0\right)=u_0.\end{array}$$

(8)

Now to construct the numerical scheme of the GRK4M [26], we divide the interval $[0,T]$ to n of equal subintervals with the nodes $\left\{t_0,t_1,\dots ,t_n\right\}$, where $t_0=0,t_j=jh,$$j=1,2,\dots n$, and $t_n=T$, and $h=\frac{T}{n}$ is the step size.

The numerical scheme of the GRK4M takes the following form:

$$\begin{array}{ll}{u}_{n+1}=u_n+\frac{1}{6}\left(K_1+2K_2+2K_3+K_4\right),& \\ K_1=\hslash \psi \left(t_n,u_n\right),& K_2=\hslash \psi \left(t_n+\frac{1}{2}\hslash ,u_n+\frac{1}{2}{K}_1\right),\\ K_3=\hslash \psi \left(t_n+\frac{1}{2}\hslash ,u_n+\frac{1}{2}{K}_2\right),& K_4=\hslash \psi \left(t_n+\hslash ,u_n+K_3\right),\end{array}$$

(9)

where $\hslash =\frac{h^{\alpha }}{\Gamma (\alpha +1)}$.

To study the stability of the GRK4M, we take for simplicity the function $\psi (t;u(t\left)\right)=\zeta u\left(t\right)$, for some constant $\zeta <0.$

Theorem 6.

The numerical scheme (9) of the FDE (8) with the assumption $\psi (t;u(t\left)\right)=\zeta u\left(t\right),\zeta <0$ is conditionally stable under the following criterion:

$$0<h^{\alpha }<12\left|\frac{\Gamma (\alpha +1)}{\zeta }\right|.$$

(10)

Proof. 

We can rewrite Equation (8) by using GRK4M as follows:

$$u\left(t_{j+1}\right)=u\left(t_j\right)+\frac{1}{6}\frac{h^{\alpha }\zeta }{\Gamma (\alpha +1)}u\left(t_j\right),j=0,1,\dots ,n-1.$$

When the terms of this equation are regrouped, we can obtain the following recurrence formula [27]:

$$u\left(t_{j+1}\right)={\left(1+\frac{1}{6}\frac{h^{\alpha }\zeta }{\Gamma (\alpha +1)}\right)}^j{u}_0,j=0,1,\dots ,n-1.$$

Then, the stability condition of this recurrence formula is given as follows:

$$\left|\left(1+\frac{1}{6}\frac{h^{\alpha }\zeta }{\Gamma (\alpha +1)}\right)\right|<1,$$

under some simplifications, we can obtain the condition (10), and this ends the proof. □

5.2. Procedure of Solution by the Fractional FDM

First, for system (2), we consider the following:

$$\begin{array}{l}{D}^{\alpha }{S}_n\left(t\right)=\Lambda -(\tau +\mu )S_n\left(t\right)-\beta S_n\left(t\right)E{(}_{n}t),\\ D^{\alpha }{E}_n\left(t\right)=\beta S_n\left(t\right)E_n\left(t\right)-(\gamma +\mu +\eta +\sigma )E_n\left(t\right),\\ D^{\alpha }{Q}_n\left(t\right)=\tau S_n\left(t\right)+\gamma E_n\left(t\right)-(\mu +v+\theta )Q_n\left(t\right),\\ D^{\alpha }I{A}_n\left(t\right)=\sigma E_n\left(t\right)+\theta Q_n\left(t\right)-\left(\mu +r_1\right)I{A}_n\left(t\right),\\ D^{\alpha }I{S}_n\left(t\right)=\eta E_n\left(t\right)+v{Q}_n\left(t\right)-\left(\delta +\mu +r_2\right)I{S}_n\left(t\right),\\ D^{\alpha }{R}_n\left(t\right)=r_{1}I{A}_n\left(t\right)+r_{2}I{S}_n\left(t\right)-\mu R_n\left(t\right).\end{array}$$

(11)

Second, we use $t_n=nh$, where $n=0,1,...,M,Mh=T$ and the abbreviations $S_n,E_n,Q_n,I{A}_n,I{S}_n,$ and $R_n$ for approximation of the true solutions $S\left(t_n\right),E\left(t_n\right),Q\left(t_n\right)$, $IA\left(t_n\right),IS\left(t_n\right),$ and $R\left(t_n\right),$ in the grid point $t_n$.

By applying the definition of the approximate fractional derivative (6)–(11), we obtain the following:

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{S}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{S}_0=\Lambda -(\tau +\mu )S_n-\beta S_n{E}_n,$$

(12)

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{E}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{E}_0=\beta S_n{E}_n-(\gamma +\mu +\eta +\sigma )E_n,$$

(13)

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{Q}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{Q}_0=\tau S_n+\gamma E_n-(\mu +v+\theta )Q_n,$$

(14)

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }I{A}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}I{A}_0=\sigma E_n+\theta Q_n-\left(\mu +r_1\right)I{A}_n,$$

(15)

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }I{S}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}I{S}_0=\eta E_n+v{Q}_n-\left(\delta +\mu +r_2\right)I{S}_n,$$

(16)

$$\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\theta }_i^{\alpha }{R}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{R}_0=r_{1}I{A}_n+r_{2}I{S}_n-\mu R_n,$$

(17)

where ${\psi }_i^{\alpha }=\left(\begin{array}{c} \alpha \\ i \end{array}\right)$.

We express the problem defined by (12)–(17) as a constrained optimization problem with the following cost functions (CFs):

$$CF1=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{S}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{S}_0-\Lambda +(\tau +\mu )S_n+\beta S_n{E}_n,$$

(18)

$$CF2=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{E}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{E}_0-\beta S_n{E}_n+(\gamma +\mu +\eta +\sigma )E_n,$$

(19)

$$CF3=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{Q}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{Q}_0-\tau S_n-\gamma E_n+(\mu +v+\theta )Q_n,$$

(20)

$$CF4=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }I{A}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}I{A}_0-\sigma E_n-\theta Q_n+\left(\mu +r_1\right)I{A}_n,$$

(21)

$$CF5=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }I{S}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}I{S}_0-\eta E_n-\upsilon Q_n+\left(\delta +\mu +r_2\right)I{S}_n,$$

(22)

$$CF6=\frac{1}{h^{\alpha }}\sum _{i=0}^{n+1}{(-1)}^i{\psi }_i^{\alpha }{R}_{n+1-i}-\frac{{\left(nh\right)}^{-\alpha }}{\Gamma (1-\alpha )}{R}_0-r_{1}I{A}_n-r_{2}I{S}_n+\mu R_n.$$

(23)

6. Numerical Simulation

6.1. Validation of the GRK4M

The goal of this subsection is to demonstrate the accuracy and reliability of the numerical values acquired through the GRK4M. Clearly, the comparison is carried out for the values of the numerical solution in the interval $(0,140)$ by using the GRK4M with the standard RK4M (SRK4M) at the order of the fractional derivative $\alpha =1$ and $h=0.05$, with $S_0=0.2,E_0=Q_0=0.1,I{A}_0=I{S}_0=0.2,R_0=0.0$, where in this case, ${\mathfrak{R}}_0=0.359916<1$. This verification is suggested in Figure 2, which shows that the current results are in excellent accord with those obtained by the SRK4M. This also confirms that the proposed method (GRK4M) is good implemented. In addition, the numerical scheme by the introduced method, and the given models (1) and (2) are stable.

We consider the systems (1) and (2) with the following parameters [15]:

$$\tau =0.0002,\beta =0.0805,\delta =0.000016728,\gamma =0.00020138,$$

$$\eta =0.4478,\theta =0.0101,mu=0.0106,\upsilon =0.0003208,$$

$$\sigma =0.0668,\Lambda =0.02537,r_1=0.00005734,r_2=0.00001673.$$

6.2. Numerical Results Using GRK4M and the Fractional FDM

We consider (2) with the following parameters [15]:

$$\tau =0.0002,\beta =0.08,\delta =0.000016,\gamma =0.0002,\eta =0.44,\theta =0.01,$$

$$\mu =0.01,\upsilon =0.0003,\sigma =0.06,\Lambda =0.025,r_1=0.000057,r_2=0.000016,$$

and the different values of $S_0,E_0,Q_0,I{A}_0,I{S}_0,R_0,$ and distinct values of the order of the fractional derivative $\alpha$.

In Figure 3, the solution for all components of the model by using GRK4M via distinct values of $\alpha =1.0,0.92,0.82,0.72,$ is presented in the interval ($0,60$) and $h=0.01$, with $S_0=0.5,E_0=0.2,Q_0=0.1,I{A}_0=0.2,I{S}_0=0.1,R_0=0.0$, where in this case, ${\mathfrak{R}}_0=0.359916<1$, and in the view of Theorem 1, we can note that ${\xi }_0$ is locally asymptotically stable.

In Figure 4, the solution for all components of the model is presented using the fractional FDM via the same values in Figure 3. From the Figure 3a,b and Figure 4a,b, we can observe that if $\alpha$ decreases, the amounts $S\left(t\right)$ and $R\left(t\right)$, respectively, will be reduced. Hence, a small memory of the infection effect maximizes the number of COVID-19 healthy individuals. In Figure 5, the solution via distinct values of the initial conditions with $h=0.05,\alpha =0.95$ is presented in the interval ($0,120$), where the components of solution $S\left(t\right),E\left(t\right),Q\left(t\right),I_A\left(t\right),I_S\left(t\right),R\left(t\right)$ are presented in Figure 5a–f, respectively. There are three cases:

• $S_0=0.5,E_0=0.1,Q_0=0.1,I{A}_0=0.3,I{S}_0=0.3,R_0=0.0$;
• $S_0=0.3,E_0=0.2,Q_0=0.2,I{A}_0=0.1,I{S}_0=0.2,R_0=0.0$;
• $S_0=0.4,E_0=0.3,Q_0=0.3,I{A}_0=0.2,I{S}_0=0.1,R_0=0.0$.

We note that, in all these cases, the basic reproduction number ${\mathfrak{R}}_0<1$.

Finally, in Figure 6, we presented a comparison between the numerical solution obtained by the GRK4M and the fractional FDM at $\alpha =0.94$ and $h=0.05$. With the same parameters included in Figure 4, and the following initial conditions:

$$S_0=0.2,E_0=0.1,Q_0=0.1,I{A}_0=0.2,I{S}_0=0.2,R_0=0.0.$$

From the COVID-19 dynamics in this figure, we can observe and confirm that the majority of the population recover because the recovered population grows significantly; see Figure 6f. In addition, we can observe and confirm that the populations of those infected and exposed decrease significantly; see Figure 6b,d,e. This means that the majority of the population will be recovered, which results in a decrease in the deaths caused by COVID-19.

We can confirm that the expected behavior of the disease has occurred, and this presents a clear simulation of the model. Additionally, the good physical interpretation of these numerical results is devoted to knowing the behavior of all of the components (stats) of the model with different values of the derivative, not only with $\alpha =1$.

7. Conclusions and Remarks

The studied problem regarding the COVID-19 model is solved numerically by employing the GRK4M and the fractional FDM with the help of the Grünwald–Letnikov’s formula of the fractional derivative. The Caputo-derivative was used here since it simply requires the initial circumstances to be expressed in terms of integer-order derivatives. The proposed model’s qualitative analysis is demonstrated by studying the locally asymptotically stable equilibrium of its disease-free equilibrium, the globally asymptotically stable equilibrium of its endemic equilibrium, and the locally asymptotically stable equilibrium of its disease-free equilibrium. Two equilibrium points, the disease-free equilibrium point ${\overline{\mathbf{E}}}_0$ and the endemic equilibrium point ${\overline{\mathbf{E}}}_1$, were found for the studied model. The stability analysis shows that ${\overline{\mathbf{E}}}_0$ is locally asymptotically stable whenever ${\mathfrak{R}}_0<1$ and ${\overline{\mathbf{E}}}_1$ is globally asymptotically stable whenever ${\mathfrak{R}}_0>1.$ Using the fitted values of the parameters showing the spread of the disease, a sensitivity analysis was performed on the parameters in the ${\mathfrak{R}}_0$ as well as the profile of each state variable. The contact rates between susceptible individuals and transmission rates from the exposure to the symptomatically infected classes are key parameters of the ${\mathfrak{R}}_0$.

Besides computing numerical solutions with different values of $\alpha$ and m, the initial conditions, and the REF, we perform a verification of the proposed approach. We can conclude that the Caputo-operator is a very efficient tool for investigating the numerical solution for the considered model in our study. The method presents an efficient way of investigating the numerical solution for such a model. Additionally, from the presented comparison, we can see the excellent agreement between our numerical solution for the fractional model obtained by using the GRK4M and the numerical solution obtained using the fractional FDM at various values of $\alpha$. In contrast, the numerical analysis of our work provides field meanings. We performed our numerical analysis using Mathematica. Finally, in the future, the authors are interested extending the model with studies of its optimal control and considering different forms to cover all of the affected parameters of the disease as well as the concerned fractional COVID-19 model can also be analyzed.

In addition, we can recommend that, by increasing the awareness of the population and by adhering to the government’s actions, the amount of susceptible individuals is reduced and the number of infected population is increased, and this in turn plays a major role in reducing the severity of infection. As a final conclusion, as we know that most numerical methods have some limitations, the two proposed methods have some shortness since we cannot evaluate the solution at any point but only at the nodes of the given domain, and they are not unconditionally stable. We plan to achieve some future work such as the following:

• Solving the same model but with different techniques such as the finite element method, the finite volume method, and others;
• Optimal control of the resulting solutions;
• Some deepness studies in the theoretical part to describe the COVID-19 model;
• Change the sense of the fractional derivative to be, for example, Atangana–Baleanu–Caputo or to be variable-order.