New Coronavirus (2019-nCov) Mathematical Model Using Piecewise Hybrid Fractional Order Derivatives; Numerical Treatments ^†

Abstract

A new mathematical model of Coronavirus (2019-nCov) using piecewise hybrid fractional order derivatives is given in this paper. Moreover, in order to be consistent with the physical model problem, a new parameter $\mu$ is presented. The boundedness, existence, and positivity of the solutions for the proposed model are discussed. Two improved numerical methods are presented in this paper. The Caputo proportional constant nonstandard modified Euler–Maruyama method is introduced to study the fractional stochastic model, and the Grünwald–Letnikov nonstandard finite difference method is presented to study the hybrid fractional order deterministic model. Comparative studies with real data from Spain and Wuhan are presented.

1. Introduction

Infectious diseases continue to impact a vast number of people across the world, despite substantial technical improvements in treatment and prevention. To reduce the prevalence of these diseases, tight management of infection transmissibility, population size, human contact rates, infection duration, and other critical parameters are required to control the spread of infectious diseases. Mathematical models play an important role in understanding the dynamics of pandemic diseases and their control, especially in the initial stages of diseases or when vaccination is not accessible. In recent years, the number of researchers from numerous domains have studied and investigated various types of biological models. A variety of mathematical models have been created in the literature to represent the dynamics of various infectious diseases [1,2,3,4,5,6,7,8,9,10].

Recently, the literature has proposed piecewise fractional differential equations, i.e., deterministic–stochastic differential equations or vice versa. The goal of this piecewise approach is to effectively investigate the model with real data [11,12,13]. Many processes in real-world problems display crossover behavior. Humanity has found it challenging to model processes based on crossover behavior. Real-world challenges have been found in a variety of cases as a result of the transition from Markovian to random processes, such as in epidemiology with the spread of infectious diseases and even some chaos. Deterministic and stochastic methodologies were created independently to improve the future state of the system and predictability [12].

The concepts of piecewise integral and differential operators were used and proposed to describe several complicated real-world problems, such as chaos problems [14]; the authors introduced three different cases for formulating the deterministic–stochastic chaotic models and studied these models numerically. Furthermore, the concept of piecewise differentiation appears to be efficient for the epidemiological problems, when there are modeling problems with crossover behaviors [12,15].

Recently, there have been very interesting papers on piecewise Coronavirus (2019-nCov) mathematical models such as [11,16,17].

In the current paper, we will present a novel Coronavirus (2019-nCov) mathematical model using piecewise fractional differential derivatives. Moreover, we will construct new numerical methods dependent on the modified Euler–Maruyama method (MEM) and the nonstandard finite difference method with the discretization of the Caputo proportional constant (CPC-NMEMM) to solve the hybrid fractional Brownian motion stochastic model and Grünwald–Letnikov nonstandard finite difference approximation, with the discretization of the Caputo proportional constant, to numerically study the hybrid fractional order deterministic models [18]. The mean square stability analysis of the CPC-NMEMM is studied. Some important analysis including boundedness, positivity, uniqueness and existence of the proposed model solutions are studied. In order to demonstrate the effectiveness and broad applicability of the suggested methods, numerical simulations will be provided.

To our knowledge, the proposed piecewise model and CPC-NMEMM given in this paper are new and have not been conducted before.

The rest of this paper is organized as follows: In Section 2, important definitions of hybrid fractional order derivatives, fractional Brownian motion and fractional Gaussian noise are introduced. An analysis of the fractional stochastic–deterministic model is proposed in Section 3; moreover, the existence and uniqueness of the proposed model are proved in Section 4. In Section 5, numerical schemes for solving the proposed model are constructed. Numerical simulations are presented with discussions in Section 6. Finally, the conclusions are offered in Section 7.

2. Fractional Order Operator Definitions

Definition 1. 

The hybrid fractional order proportional Caputo operator is defined as follows [19]:

$$\begin{array}{c}{ }_0^{CP}{D}_t^{\alpha }y\left(t\right)={(\Gamma (1-\alpha ))}^{-1}\left({\int }_0^t{(t-s)}^{-\alpha }\left(y\left(s\right)K_1(\alpha ,s)+y^{\prime }\left(s\right)K_0(\alpha ,s)\right)ds\right),\\ \hfill =(t^{-\alpha }\left(\Gamma \left(1-\alpha \right){)}^{-1}\right)\left(K_0\left(\alpha ,t\right)y^{\prime }\left(t\right)+y\left(t\right)K_1\left(\alpha ,t\right)\right).\end{array}$$

(1)

Furthermore, the hybrid fractional order Caputo proportional constant (CPC) operator is defined as follows [19]:

$$\begin{array}{c}{ }_0^{CPC}{D}_t^{\alpha }y\left(t\right)={(\Gamma (1-\alpha ))}^{-1}\left({\int }_0^t{(t-s)}^{-\alpha }\left(y^{\prime }\left(s\right)K_0\left(\alpha \right)+y\left(s\right)K_1\left(\alpha \right)\right)ds\right)\\ =K_0\left(\alpha \right){}_0^C{D}_t^{\alpha }y\left(t\right){+}_0^{RL}{I}_t^{1-\alpha }y\left(t\right)K_1\left(\alpha \right),\end{array}$$

(2)

where K₁(α) and K₀(α) are constants. Here, we consider K₀(α) = αQ^{(1 − α)}, K₁(α) = (1 − α)Q^α,, 0 < α < 1 and Q as a constant.

In both of these formulas, the domain of functional space is given by requiring that y is differentiable and both y and y0 are locally L¹ functions on the positive reals.

Definition 2. 

The fractional order CPC derivative’s inverse operators are given by [19]:

$$\begin{array}{c}{ }_0^{CPC}{I}_t^{\alpha }y\left(t\right)=\left({\int }_0^t\mathit{exp}{\left[K_1\left(\alpha \left(t\right)\right){\left(K_0\left(\alpha \right)\right)}^{-1}(t-s)\right]}_0^{RL}{D}_t^{1-\alpha }y\left(s\right)ds\right)\frac{1}{K_0\left(\alpha \right)}.\end{array}$$

(3)

2.1. Fractional Gaussian Noise and Fractional Brownian Motion

The integer order Brownian motion $B_t$ is the classical white noise $W_t$:

$$\begin{array}{c}{B}_t={\int }_0^t{W}_{s}ds.\end{array}$$

(4)

For the fractional Brownian motion (FBM), there are various definitions. The most prevalent definition will be indicated as $B_1^{H^{*}}\left(t\right)$ and is referred to as the order’s fractional integral $H^{*}-0.5$ of classical white noise, i.e., [20]:

$$\begin{array}{c}{B}_1^{H^{*}}\left(t\right)=\frac{1}{\Gamma \left(H^{*}-0.5\right)}{\int }_0^t{(t-s)}^{H^{*}-0.5}{W}_{s}ds.\end{array}$$

(5)

When $H^{*}+0.5$, it reverts to the standard Brownian motion definition (4). Another formulation, $B_2^{H^{*}}\left(t\right),$ extends the classical Brownian motion but with a fractional order increment of order $H^{*}$ [21]:

$$\begin{array}{c}{B}_2^{H^{*}}\left(t\right)={\int }_0^{t}W\left(t\right)d{t}^{H^{*}+0.5}.\end{array}$$

(6)

For simplicity, we denote $B_1^{H^{*}}\left(t\right)$ as $B^{H^{*}}\left(t\right)$.

Now, we can write:

$$d{B}^{H^{*}}=W^{H^{*}}dt.$$

(7)

$$\begin{array}{c}{B}_t={\int }_0^t{W}^{H^{*}}\left(s\right)ds.\end{array}$$

(8)

Then, (5) is given as:

$$\begin{array}{c}{B}^{H^{*}}\left(t\right)={\left(\Gamma \left(H^{*}-0.5\right)\right)}^{-1}{\int }_0^t{(t-s)}^{-0.5+H^{*}}{W}_{s}ds,\end{array}$$

(9)

or

$$\begin{array}{c}{B}^{H^{*}}\left(t\right)=\frac{1}{\Gamma \left(H^{*}+1.5\right)}{\int }_0^t{(t-s)}^{0.5+H^{*}}{W}_{s}ds.\end{array}$$

(10)

2.2. Modified Euler–Maruyama Technique

Consider the following stochastic differential equations (SDEs) driven by FBM:

$$d{y}_{\vartheta }={\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)+{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)d{B}_{\eta }^{H^{*}}\left(t\right), 0<t\le T,$$

(11)

$$y_{\vartheta }\left(t_0\right)=y_{\vartheta ,0}, \eta ,\vartheta =1,\dots ,\kappa ,$$

where ${\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)$ and ${\mathsf{{\rm Y}}}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)$ are real continuous functions and $y_{\vartheta ,0}$ is a deterministic initial value. Furthermore, ${\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)$ is the mean rate of change of the system state $y_{\kappa }$ at time t and ${\mathsf{{\rm Y}}}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)d{B}_{\eta }^{H^{*}}\left(t\right)$ is the random perturbation. The term ${\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)$ is used to represent the average or the deterministic part of the problem and ${\mathsf{{\rm Y}}}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t\right)$ is used to represent the intensity of the random part.

In the case of classical Brownian motion, i.e., when $H^{*}-0.5$, the following Euler–Maruyama method (EMM) is a popular solution for Equation (11):

$$\begin{array}{c}{y}_{\vartheta }^{n+1}=y_{\vartheta }^n+{\mathsf{\Lambda }}_{\vartheta }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right)\mathsf{\Delta }t+{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right)\mathsf{\Delta }{B}_{\eta }^n, 0<t\le T,\\ \eta ,\vartheta =1,\dots ,\kappa .\end{array}$$

(12)

The MEMM for solving Equation (11) in the case of the FBM can be written as $H^{*}>0.5$ [22].

$$\begin{gathered} y_{\vartheta }^{n+1}=y_{\vartheta }^n+{\mathsf{\Lambda }}_{\vartheta }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right)\mathsf{\Delta }t+{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right)\mathsf{\Delta }{B}_{\eta }^n \\ +0.5{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right){\mathsf{{\rm Y}}}_{\vartheta }^{\prime }\left(y_1^n,y_2^n,\dots ,y_{\kappa }^n,t^n\right)\mathsf{\Delta }{t}^{2H^{*}}, 0<t\le T, \eta ,\vartheta =1,\dots ,\kappa , \end{gathered}$$

(13)

where $H^{*}$ is the Hurst index. The EMM has a relatively low convergence rate, according to [21,23]; the strong convergence of the EMM for pure Brownian motion $H^{*}=0.5$ is $O\left(\mathsf{\Delta }{t}^{0.5}\right),$ but for FBM with $H^{*}>0.5$ the convergence of the MEMM is $O\left(\mathsf{\Delta }{t}^{H^{*}-0.5}\right),$ meaning that as the Hurst index $H^{*}$ increases, the rate of convergence of the MEMM increases. The biggest challenge with adopting the EMM in the instance of FBM is the time required to generate random samples that include a matrix inversion in addition to the huge sample size required for reliable results [24].

3. Fractional Stochastic–Deterministic Model

In the following, we will extend the Coronavirus (2019-nCov) mathematical model [25] by applying the piecewise differentiation technique. This model takes into account the super-spreading phenomenon of some individuals. Moreover, we consider a fatality compartment related to death due to the virus infection. In doing so, the constant total population size N is subdivided into eight epidemiological classes: susceptible class (S), exposed class (E), symptomatic and infectious class (I), super-spreaders class (P), infectious but asymptomatic class (A), hospitalized class (H), recovery class (R), and fatality class (F). The new hybrid fractional order operator CPC, which is defined as a linear combination of integral Riemann–Liouville and the fractional order derivative of Caputo, is applied to extend the deterministic model. The fractional Brownian motion and the hybrid fractional order operator are applied to extend SDEs. A new parameter, $\mu$, is presented in order to be consistent with the physical model problem. Moreover, we avoid dimensional mismatches by modifying the fractional model with an auxiliary parameter $\mu .$ As a result, the left side possesses the dimension day⁻¹ [26]. Now, the updated nonlinear piecewise fractional order differential mathematical model is given as follows:

$$\{\begin{array}{c}{\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }S=-\beta \frac{IS}{N}-l\beta \frac{HS}{N}-{\beta }^{\prime }\frac{PS}{N},\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }E=\beta \frac{IS}{N}+l\beta \frac{HS}{N}+{\beta }^{\prime }\frac{PS}{N}-KE,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }I=K{\rho }_{1}E-\left({\gamma }_a+{\gamma }_i\right)I-{\delta }_{i}I,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }P=K{\rho }_{2}E-\left({\gamma }_a+{\gamma }_i\right)P-{\delta }_{p}P,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }A=K\left(1-{\rho }_1-{\rho }_2\right)E,0<t\le T_1,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }H={\gamma }_a(I+P)-{\gamma }_{r}H-{\delta }_{h}H,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }R={\gamma }_i(I+P)+{\gamma }_{r}H,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }F={\delta }_{i}I+{\delta }_{p}P+{\delta }_{h}H,\hfill \end{array}$$

(14)

$$\begin{gathered} P\left(0\right)=P_0, E\left(0\right)=E_0, S\left(0\right)=S_0, F\left(0\right)=F_0, A\left(0\right)=A_0, R\left(0\right)=R_0, \\ H\left(0\right)=H_0, I\left(0\right)=I_0. \end{gathered}$$

(15)

$$\{\begin{array}{c}{\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}S=\left(-\beta \frac{IS}{N}-l\beta \frac{HS}{N}-{\beta }^{\prime }\frac{PS}{N}\right)+{\sigma }_{1}S\frac{d{B}_1^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}E=\left(\beta \frac{IS}{N}+l\beta \frac{HS}{N}+{\beta }^{\prime }\frac{PS}{N}-KE\right)+{\sigma }_{2}E\frac{d{B}_2^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}I=\left(-{\delta }_{i}I+K{\rho }_{1}E-\left({\gamma }_a+{\gamma }_i\right)I\right)+{\sigma }_{3}I\frac{d{B}_3^{H^{*}}\left(t\right)}{dt},T_1<t\le T,\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}P=\left(-{\delta }_{p}P-\left({\gamma }_a+{\gamma }_i\right)P+K{\rho }_{2}E\right)+{\sigma }_{4}P\frac{d{B}_4^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}A=\left(\left(-{\rho }_1-{\rho }_2+1\right)KE\right)+{\sigma }_{5}A\frac{d{B}_5^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}H=\left(-{\gamma }_{r}H-{\delta }_{h}H+{\gamma }_a(P+I)\right)+{\sigma }_{6}H\frac{d{B}_6^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}R=\left({\gamma }_{r}H\right)+{\gamma }_i(P+I)+{\sigma }_{7}R\frac{d{B}_7^{H^{*}}\left(t\right)}{dt},\hfill \\ {\frac{1}{{\mu }^{1-{\alpha }_1}}}_0^{CPC}{D}_t^{{\alpha }_1}F=\left({\delta }_{h}H+{\delta }_{i}I+{\delta }_{p}P\right)+{\sigma }_{8}F\frac{d{B}_8^{H^{*}}\left(t\right)}{dt}.\hfill \end{array}$$

(16)

$$E\left(T_1\right)=E_1, P\left(T_1\right)=P_1, I\left(T_1\right)=I_1, H\left(T_1\right)=H_1, R\left(T_1\right)=R_1, A\left(T_1\right)=A_1,$$

$$F\left(T_1\right)=F_1, S\left(T_1\right)=S_1.$$

4. Analysis of the Model

It is well known that it makes no sense for the system to have negative solutions because the model’s purpose is population evolution. As a result, we first show that the state variables are both nonnegative and bounded. Boundedness of the proposed model solution can be verified by adding all of the equations of system (14) as follows:

$$\begin{gathered} { }_0^{CPC}{D}_t^{\alpha }S\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }I\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }E\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }A\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }P\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }H\left(t\right) \\ +{ }_0^{CPC}{D}_t^{\alpha }R\left(t\right)+{ }_0^{CPC}{D}_t^{\alpha }F\left(t\right)={ }_0^{CPC}{D}_t^{\alpha }N\left(t\right). \end{gathered}$$

Then,

$${ }_0^{CPC}{D}_t^{\alpha }N\left(t\right)=0, 0\le N\left(0\right)=A,$$

(17)

where $A$ is a constant and $N$ represents the population. The solution to (17) is provided as follows [19]:

$$N\left(t\right)\ge A{e}^{\left(-\frac{K_1\left(\alpha \right)}{K_0\left(\alpha \right)}t\right)}.$$

(18)

The inequality (18) implies that $0\le N\left(t\right),$ $t⟶\infty ,$ i.e., the solutions of the system (14) are bounded.

Theorem 1. 

All solutions of system (14) are nonnegative and remain in ${ℝ}^8$ for all $t\ge 0$.

Proof. 

Using (3), we have

$$\begin{array}{c}y\left(t\right)=y\left(t_0\right)+\frac{1}{K_0\left(\alpha \right)}{\int }_{t_0}^t\mathit{exp}{\left(-\frac{K_1\left(\alpha \right)}{K_0\left(\alpha \right)}(t-s)\right)}_0^{RL}{D}_t^{1-\alpha }y\left(s\right)ds>0.\end{array}$$

(19)

Using [27], Equations (15) and (19), if $0\le t<\infty$, we can prove the soluations of the following equations are positive:

$$\begin{array}{c}{\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }S{|}_{S=0}=0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }E{|}_{E=0}=\beta \frac{IS}{N}+l\beta \frac{HS}{N}+{\beta }^{\prime }\frac{PS}{N}\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }I{|}_{I=0}=K{\rho }_{1}E\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }P{|}_{P=0}=K{\rho }_{2}E\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }A{|}_{A=0}=K\left(1-{\rho }_1-{\rho }_2\right)E\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }H{|}_{H=0}={\gamma }_a\left(I+P\right)\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }R{|}_{R=0}=\left({\gamma }_i\left(I+P\right)+{\gamma }_{r}H\right)\ge 0,\hfill \\ {\frac{1}{{\mu }^{1-\alpha }}}_0^{CPC}{D}_t^{\alpha }F{|}_{F=0}=({\delta }_{i}I+{\delta }_{p}P+{\delta }_{h}H)\ge 0.\hfill \end{array}$$

(20)

It follows from (20) that the solution is positive and will remain in ${ℝ}^8$ and, thus, we define the biologically feasible region for the model (15), given

$$Z=\left\{\left(S\left(t\right), I\left(t\right),P\left(t\right), A\left(t\right),R\left(t\right), H\left(t\right),F\left(t\right), E\left(t\right)\right)\in {ℝ}^8|S\left(t\right)+I\left(t\right)+P\left(t\right)+A\left(t\right)+R\left(t\right)+H\left(t\right)+F\left(t\right)+E\left(t\right)=N\ge 0\right\}\square$$

4.1. Basic Reproduction Number

In the following, we briefly derive the basic reproduction number $\left(R_0\right)$ for (14) by using the next generation method [28]. Consider the following matrices $F_1$ and $V$, where $F_1$ represents the new infection terms, and $V$ represents the remaining transfer terms [28,29]:

$$F_1={\mu }^{1-\alpha }\left(\begin{array}{llll}0\hfill & \beta \hfill & \beta \hfill & {\beta }^{\prime }\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right),$$

$$V={\mu }^{1-\alpha }\left(\begin{array}{llll}-K\hfill & 0\hfill & 0\hfill & 0\hfill \\ -K\hfill & -\left({\gamma }_a+{\gamma }_i\right)\hfill & 0\hfill & 0\hfill \\ -K{\rho }_1\hfill & 0\hfill & -\left({\gamma }_a+{\gamma }_i+{\delta }_p\right)\hfill & 0\hfill \\ 0\hfill & {\gamma }_a\hfill & {\gamma }_a\hfill & -\left({\gamma }_r+{\delta }_h\right)\hfill \end{array}\right)$$

Then,

$$R_0=\rho \left(F_1{V}^{-1}\right)={\mu }^{1-\alpha }\left[\frac{\beta {\rho }_1\left({\gamma }_{a}l+\left({\gamma }_r+{\delta }_h\right)\right)}{\left({\gamma }_r+{\delta }_h\right)\left({\gamma }_a+{\gamma }_i+{\delta }_i\right)}+\frac{\left({\beta }^{\alpha }{\gamma }_{a}l+{\beta }^{\prime }\left({\gamma }_r+{\delta }_h\right)\right){\rho }_2}{\left({\gamma }_r+{\delta }_h\right)\left({\gamma }_a+{\gamma }_i+{\delta }_p\right)}\right],$$

(21)

where $\rho$ denotes the spectral radius of the $F_1{V}^{-1}$ [28].

4.2. Uniqueness and Existence

The existence and uniqueness of the solutions of the proposed model will be established using.

Banach fixed point theorem. Let $0\le t<\infty ;$ the system (14) can be written as follows:

$${ }_0^{CPC}{D}_t^{\alpha }y\left(t\right)=q\left(y\left(t\right),t\right), y\left(0\right)=y_0\ge 0.$$

(22)

$y\left(t\right)={(S,H,E,I,P,D,R,A)}^T$ and refers to the state variables, and $q$ is a vector of continuous function such that:

$$\left(\begin{array}{c}{q}_1\hfill \\ q_2\hfill \\ q_3\hfill \\ q_4\hfill \\ q_5\hfill \\ q_6\hfill \\ q_7\hfill \\ q_8\hfill \end{array}\right)=\left(\begin{array}{c}{\mu }^{1-\alpha }\left(-\beta \frac{IS}{N}-l\beta \frac{HS}{N}-{\beta }^{\prime }\frac{PS}{N}\right)\hfill \\ {\mu }^{1-\alpha }\left(\beta \frac{IS}{N}+l\beta \frac{HS}{N}+{\beta }^{\prime }\frac{PS}{N}-KE\right)\hfill \\ {\mu }^{1-\alpha }\left(-\left({\gamma }_a+{\gamma }_i\right)I-{\delta }_{i}I+K{\rho }_{1}E\right)\hfill \\ {\mu }^{1-\alpha }\left(-\left({\gamma }_i+{\gamma }_a\right)P+K{\rho }_{2}E-{\delta }_{p}P\right)\hfill \\ {\mu }^{1-\alpha }\left(\left(-{\rho }_2-{\rho }_1+1\right)KE\right)\hfill \\ {\mu }^{1-\alpha }\left(-{\gamma }_{r}H+{\gamma }_a(P+I)-{\delta }_{h}H\right)\hfill \\ {\mu }^{1-\alpha }\left({\gamma }_{r}H+{\gamma }_i(P+I)\right)\hfill \\ {\mu }^{1-\alpha }\left({\delta }_{h}H+{\delta }_{i}I+{\delta }_{p}P\right)\hfill \end{array}\right)$$

with initial condition $y_0.$ Moreover, the Lipschitz condition is satisfied by $q,$ where $q$ is a quadratic vector function, i.e., there exists $M^0\in ℝ,$ such that [30]:

$$\parallel q\left(y_1\left(t\right),t\right)-q\left(y_2\left(t\right),t\right)\parallel \le M^0\parallel y_1\left(t\right)-y_2\left(t\right)\parallel .$$

(23)

Theorem 2. 

The proposed model (14) has a unique solution if the below condition holds:

$$\frac{M^0{{\rm Y}}_{max}^{\alpha }{X}_{max}^{\alpha }}{\mathsf{\Gamma }\left(\alpha -1\right)K_0\left(\alpha \right)}<1, 0\le t<\infty .$$

(24)

Proof. 

Applying (2) in (22), we have:

$$y\left(t\right)=y\left(t_0\right)+\frac{1}{K_0\left(\alpha \right)}{\int }_0^t\mathit{exp}{\left(-\frac{K_1\left(\alpha \left(t\right)\right)}{K_0\left(\alpha \right)}(t-s)\right)}_0^{RL}{D}_t^{1-\alpha }q(y\left(s\right),s)ds.$$

(25)

Let $K=\left(0,T\right)$ and $B:C\left(K,{ℝ}^8\right)⟶C\left(K,{ℝ}^8\right)$ $\Rightarrow$

$$B\left[y\left(t\right)\right]=y\left(t_0\right)+\frac{1}{K_0\left(\alpha \right)}{\int }_0^t\mathit{exp}{\left(-\frac{K_1\left(\alpha \right)}{K_0\left(\alpha \left(t\right)\right)}(t-s)\right)}_0^{RL}{D}_t^{1-\alpha }q(y\left(s\right),s)ds.$$

(26)

It provides:

$$y\left(t\right)=B\left[y\left(t\right)\right].$$

Let $\parallel .{\parallel }_K$ be the supremum norm on $K$. Thus

$$\parallel y\left(t\right){\parallel }_K=\underset{t\in K}{\sup }\parallel y\left(t\right)\parallel , y\left(t\right)\in C\left(K,{ℝ}^8\right).$$

So, $C\left(K,{ℝ}^8\right)$ with $\parallel .{\parallel }_K$ as a Banach space. The relationship below is valid:

$$\parallel {\int }_0^t\phi (s,t)y\left(s\right)ds\parallel \le {\Lambda \parallel \phi (s,t)\parallel }_K{\parallel y\left(s\right)\parallel }_K,$$

with $y\left(t\right)\in C\left(K,{ℝ}^8\right)$, $\phi \left(s,t\right)\in C\left(K^2,{ℝ}^8\right)$ $\Rightarrow$: $\parallel \phi \left(s,t\right){\parallel }_K={\sup }_{t,s\in K}\left|\phi \left(s,t\right)\right|.$

Then (26) is written as:

$$\begin{array}{c}{∥B\left[y_1\left(t\right)\right]-B\left[y_2\left(t\right)\right]∥}_K\le \parallel \frac{1}{K_0\left(\alpha \right)}{\int }_0^t\mathit{exp}\left(-\frac{K_1\left(\alpha \right)}{K_0\left(\alpha \right)}(t-s)\right)\left({ }_0^{RL}{D}_t^{1-\alpha }q\left(y_1\left(s\right),s\right)\right.\\ {-}_0^{RL}{D}_t^{1-\alpha }q\left(y_2\left(s\right),s\right){)ds\Vert }_K.\\ \le \frac{{{\rm Y}}_{\mathit{max}}^{\alpha }}{\Gamma (\alpha -1)K_0\left(\alpha \right)}{\parallel {\int }_0^t{(t-s)}^{\alpha -2}\left(q\left(y_1\left(s\right),s\right)-q\left(y_2\left(s\right),s\right)\right)ds\parallel }_K,\\ \le \frac{{{\rm Y}}_{max}^{\alpha \left(t\right)}{X}_{max}^{\alpha \left(t\right)}}{\Gamma \left(\alpha -1\right)K_0\left(\alpha \right)}{\Vert q\left(y_1\left(t\right),t\right)-q\left(y_2\left(t\right),t\right)\Vert }_K,\\ \le \frac{M^0{{\rm Y}}_{max}^{\alpha }{X}_{max}^{\alpha }}{\Gamma \left(\alpha -1\right)K_0\left(\alpha \right)}{\Vert y_1\left(t\right)-y_2\left(t\right)\Vert }_K.\end{array}$$

(27)

Then, we have:

$$\parallel B\left[y_1\left(t\right)\left]-B\right[y_2\left(t\right)\right]{\parallel }_K\le L\parallel y_1\left(t\right)-y_2\left(t\right){\parallel }_K,$$

(28)

where

$$L=\frac{M^0{{\rm Y}}_{max}^{\alpha \left(t\right)}{X}_{max}^{\alpha }}{\mathsf{\Gamma }\left(\alpha -1\right)K_0\left(\alpha \right)}.$$

If $L<1,$ the operator $B$ is called a contraction. Then, (14) is a unique solution. □

5. Numerical Procedure for Piecewise Model

The numerical procedure can be used to solve the model (14–16) as follows: assuming that a real-world problem exhibits crossover from memory processes (i.e., fractional order) to fractional stochastic processes, the Cauchy problem associated with this problem can be expressed as follows:

$${ }_0^{CPC}{D}_t^{\alpha }{y}_{\vartheta }={\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa }\right), 0<t\le T_1, 0<\alpha \le 1,$$

(29)

$$y_{\vartheta }\left(t_0\right)=y_{\vartheta ,0},$$

$$\begin{gathered} { }_0^{CPC}{D}_t^{{\alpha }_1}{y}_{\vartheta }={\mathsf{\Lambda }}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t^n\right)+{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1,y_2,\dots ,y_{\kappa },t^n\right)d{B}_{\eta }^{H^{*}}\left(t\right), \\ {T}_1<t\le T, 0<{\alpha }_1\le 1, \\ {y}_{\vartheta }\left(t_1\right)=y_{\vartheta ,1}. \end{gathered}$$

(30)

First, by dividing $\left[0,T\right]$ into two parts:

$$0<t\le t_1\le \dots \le t_{n1}=T1\le t_{n1+1}\le t_{n1+2}\le \dots t_{n2}=T.$$

Relation (2) can be written as follows:

$$\begin{array}{c}{ }_0^{CPC}{D}_t^{\alpha }y\left(t\right)=\frac{1}{\Gamma (1-\alpha ))}{\int }_0^t{(t-s)}^{-\alpha }(K_1\left(\alpha \right)y\left(s\right)+K_0\left(\alpha \right)y^{\prime }\left(s\right))ds,\\ =K_1{\left(\alpha \right)}_0^{RL}{I}_t^{1-\alpha }y\left(t\right)+K_0{\left(\alpha \right)}_0^C{D}_t^{\alpha }y\left(t\right),\\ =K_1{\left(\alpha \right)}_0^{RL}{D}_t^{\alpha -1}y\left(t\right)+K_0{\left(\alpha \right)}_0^C{D}_t^{\alpha }y\left(t\right),\end{array}$$

(31)

where $K_1\left(\alpha \right)$, and $K_0\left(\alpha \right)$ are solely dependent on $\alpha$.

One of the most developed methods was the nonstandard finite difference method (NSFDM) that was introduced by Mickens [31]. This method can be constructed easily to improve the discretization of some terms in the differential equations such that, depending on the denominator function and the specific discretization, this method is more accurate and more stable than the standard method; for more details on the advantages of the NSFDM compared to the standard finite difference method see [3].

We can discretize (31) using GL-approximation and the NFDM as follows:

$$\begin{array}{c}{\left.{ }_0^{CPC}{D}_t^{\alpha }y\left(t\right)\right|}_{t=t^{n1}}=\frac{K_1\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha -1}}\left(y_{n1+1}+\sum _{i=1}^{n1+1}{\omega }_i{y}_{n1+1-i}\right)\\ +\frac{K_0\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha }}\left(y_{n1+1}-\sum _{i=1}^{n1+1}{\mu }_i{y}_{n1+1-i}-q_{n1+1}{y}_0\right),\end{array}$$

(32)

$$\begin{array}{c}\frac{K_1\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha -1}}\left(y_{n_1+1}+\sum _{i=1}^{n_1+1}{\omega }_i{y}_{n_1+1-i}\right)+\frac{K_0\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha }}\left(y_{n_1+1}-\sum _{i=1}^{n_1+1}{\mu }_i{y}_{n_1+1-i}-q_{n_1+1}{y}_0\right)\\ ={\Lambda }_{\vartheta }\left(y_1^{n_1},y_2^{n_1},\dots ,y_k^{n_1}\right),\end{array}$$

(33)

where ${\omega }_0=1,$ ${\omega }_i=\left(1-\frac{\alpha }{i}\right){\omega }_{i-1},$ $t^{n_1}=n_1\left(\mathsf{\Theta }\left(\mathsf{\Delta }t\right)\right), \mathsf{\Delta }t=\frac{T_1}{N_n},$ and $N_n$ is a natural number. ${\mu }_i={(-1)}^{i-1}\left({ }_i^{\alpha }\right),$ ${\mu }_1=\alpha ,$ $q_i=\frac{i^{\alpha }}{\mathsf{\Gamma }\left(1-\alpha \right)}$ and $i=1,2,\dots ,n_1+1$.

Also, let us assume that [25]:

$$\begin{gathered} 0<{\mu }_{i+1}<{\mu }_i<\dots <{\mu }_1=\alpha <1, \\ 0<q_{i+1}<q_i<\dots <q_1=\frac{1}{\mathsf{\Gamma }\left(-\alpha +1\right)}. \end{gathered}$$

Remark 1. 

If $K_1\left(\alpha \right)=0$ and $K_0\left(\alpha \right)=1$ in (33), we can obtain the discretization of the nonstandard finite difference technique with the Caputo operator (C-NFDM).

If $H^{*}>0.5,$ the NMEMM for solving Equation (30) in the case of the FBM in $\left[T_1,T\right]$ can be written as [22]:

$$\begin{gathered} y_{\vartheta }^{n2+1}=y_{\vartheta }^{n2}+{\mathsf{\Lambda }}_{\vartheta }\left(y_1^{n2},y_2^{n2},\dots ,y_{\kappa }^{n2},t^{n2}\right)\mathsf{\Delta }t+{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1^{n2},y_2^{n2},\dots ,y_{\kappa }^{n2},t^{n2}\right)\mathsf{\Delta }{B}_{\eta }^{n2} \\ +0.5{\mathsf{{\rm Y}}}_{\vartheta }\left(y_1^{n2},y_2^{n2},\dots ,y_{\kappa }^{n2},t^{n2}\right){\mathsf{{\rm Y}}}_{\vartheta }^{\prime }\left(y_1^{n2},y_2^{n2},\dots ,y_{\kappa }^{n2},t^{n2}\right)\mathsf{\Delta }{t}^{2H^{*}}, \\ 0<t\le T,\eta ,\vartheta =1,\dots ,\kappa , \end{gathered}$$

(34)

Now, to solve FSDEs driven by FBM and CPC fractional derivatives (30) in $T_1<t\le T,$ we constructed a new method to solve $\left(30\right)$; this method is developed from a modified EM technique and Grünwald–Letnikov nonstandard finite difference method with discretization of the CPC operator. The new method, called Caputo proportional constant nonstandard modified Euler–Maruyama method (CPC-NMEMM), is given as follows:

$$\begin{array}{c}\frac{K_1\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha -2}}\left(y_{n_2+1}+{\displaystyle \sum _{i=1}^{n_2+1}}{\omega }_i{y}_{n_2+1-i}\right)+\frac{K_0\left(\alpha \right)}{{(\Theta (\Delta t))}^{\alpha -1}}\left(y_{n_2+1}-{\displaystyle \sum _{i=1}^{n_2+1}}{\mu }_i{y}_{n_2+1-i}-q_{n_2+1}{y}_0\right)\\ ={\Lambda }_{\vartheta }\left(y_1^{n_2},y_2^{n_2},...,y_{\kappa }^{n_2}\right)+{{\rm Y}}_{\vartheta }\left(y_1^{n_2},y_2^{n_2},...,y_{\kappa }^{n_2},t^{n_2}\right)\Delta B_{\eta }^{n_2}\\ +0.5{{\rm Y}}_{\vartheta }\left(y_1^{n_2},y_2^{n_2},...,y_{\kappa }^{n_2},t^{n_2}\right){{\rm Y}}_{\vartheta }^{\prime }\left(y_1^{n_2},y_2^{n_2},...,y_{\kappa }^{n_2},t^{n_2}\right)\Delta t^{2H^{*}},\\ T_1<t\le T,\vartheta =1,...,\kappa ,0<{\alpha }_1\le 1,\end{array}$$

(35)

where ${\omega }_0=1,$ ${\omega }_i=\left(1-\frac{{\alpha }_1}{i}\right){\omega }_{i-1},$ $t^{n_2}=n_2\left(\mathsf{\Theta }\left(\mathsf{\Delta }t\right)\right), \mathsf{\Delta }t=\frac{T_1}{N_n},$ and $N_n$ is a natural number. ${\mu }_i={(-1)}^{i-1}\left({ }_i^{{\alpha }_1}\right),$ ${\mu }_1={\alpha }_1,$ $q_i=\frac{i^{{\alpha }_1}}{\mathsf{\Gamma }\left(1-{\alpha }_1\right)}$ and $i=1,2,\dots ,n_2+1$. Furthermore, let us assume that [32]:

$$0<{\mu }_{i+1}<{\mu }_i<\dots <{\mu }_1={\alpha }_1<1,$$

$$0<q_{i+1}<q_i<\dots <q_1=\frac{1}{\mathsf{\Gamma }\left(-{\alpha }_1+1\right) }$$

Mean Square Stability of the CPC-NMEMM

The stability of the CPC-NMEMM (35) will be proved here. Consider:

$${ }_0^{CPC}{D}_t^{\alpha }y\left(t\right)=\sigma y\left(t\right)\frac{d{B}^{H^{*}}}{dt}+ay\left(t\right), 0<t\le T.{\displaystyle }y\left(0\right)=y_0,$$

(36)

as a test problem.

Theorem 3. 

The CPC-NMEMM (35) is a mean square stable method.

Proof. 

Using the CPC-NMEMM (35) and (36), we have:

$$\begin{array}{c}A{y}_{n_2+1}+A\sum _{i=1}^{n_2+1}{y}_{n_2-i+1}BB+G\left(y_{n_2+1}-\sum _{i=1}^{n_2+1}{\mu }_i{y}_{n_2+1-i}-q_{n_2+1}{y}_0\right)\\ =a\Theta \left(△t\right)y_{n_2}+\sigma y_{n_2}△B_{n_2}+\frac{1}{2}{\sigma }^2{y}_{n_2}{(\Theta \left(△t\right))}^{2H^{*}},\end{array}$$

(37)

where $A=\frac{Q^{\alpha }\left(1-\alpha \right)}{\mathsf{\Theta }{(△t)}^{\alpha -2}\mathsf{\Gamma }\left(2-\alpha \right)},$ $BB=\frac{Q^{\alpha }\left(1-\alpha \right)}{\mathsf{\Theta }{(△t)}^{\alpha -2}\mathsf{\Gamma }\left(2-\alpha \right)},$ $G=\frac{\alpha C^{2\alpha }{Q}^{\left(1-\alpha \right)}}{\mathsf{\Theta }{(△t)}^{\alpha -1}}$.

Then,

$$\begin{array}{c}{y}_{n_2+1}=\left(\frac{\Theta (\Delta t)a+\sigma \Delta W_{n_2}+\frac{1}{2}{\sigma }^2{(\Theta (\Delta t))}^{2H^{*}}}{A+G}\right)y_{n_2}+G\left(\sum _{i=1}^{n_2+1}{\mu }_i{y}_{n_2+1-i}-q_{n_2+1}{y}_0\right)\\ -\left(A\sum _{i=1}^{n_2+1}{y}_{n_2-i+1}BB\right),\end{array}$$

(38)

also,

$$\begin{array}{c}{y}_{n_2+1}=\left(\frac{\Theta (▵t)a+\sigma \Delta B_{n_2}+\frac{1}{2}{\sigma }^2{\left(\Theta \left(▵t_{n_2}\right)\right)}^{2H^{*}}}{A+G}\right)y_{n_2}+G\left({\displaystyle \sum _{i=1}^{n_2+1}}{\mu }_i{y}_{n_2+1-i}-q_{n_2+1}{y}_0\right)\\ -\left(A{\displaystyle \sum _{i=1}^{n_2+1}}{y}_{n_2-i+1}B\right).\end{array}$$

(39)

Then, we can claim:

$$y_{n_2+1}\le \left(\frac{\mathsf{\Theta }\left(△t\right)a+\sigma △B_{n_2}+\frac{1}{2}{\sigma }^2{(△B_{n_2})}^{2H^{*}}}{A+G}\right)y_{n_2},$$

(40)

From [33], the proposed scheme (35) for $\sigma ,\mathsf{\Theta }\left(△t\right),a$ is a mean square stable if

$$1\ge \mathbb{E}(|\frac{a\mathsf{\Theta }\left(△t\right)+\sigma △B_n+\frac{1}{2}{(\mathsf{\Theta }\left(△t\right))}^{2H^{*}}{\sigma }^2}{G+A}{|}^2).$$

□

6. Numerical Simulations

This section studies the piecewise fractional differential equations model (14–16) and its numerical simulations. We consider the reported cases of COVID-19 infection, so we consider the time unit per day. We have chosen certain parameter values from [25] for the real data fitting of the model (14–16), and the remainder values are fitted for the data gathered for Spain from $25$ February to 16 May 2020. Furthermore, since in some parts of Spain there is a more concentrated population and intensive use of public transportation, we consider $N=\frac{\mathrm{47,000,000}}{425}.$ The data collected for Wuhan were from 4 January to 9 March 2020. For all of these cases, we have considered the official data published by the WHO. The total population of Wuhan is about $\mathrm{11,000,000}$. During the COVID-19 outbreak, there was a restriction of movement of individuals due to quarantine in the city. As a consequence, there was a limitation on the spread of the disease. In agreement, in our model, we consider the total population, $N=\frac{\mathrm{11,000,000}}{250} \left[8\right].$ To clarify the importance of the proposed model, we compared the results we obtained with the real data for the states of Spain and Wuhan, as we explained previously. The methods (33 and 35) are constructed to study the proposed model. A comparative study using data from WHO for Spain versus the simulation of the proposed model (14–16) at different values of $H^{*},\alpha ,{\sigma }_i, i=1,2,\dots ,8$ and the real data from Spain are given in Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5. We noted that the obtained results of Figure 1, Figure 2 and Figure 3 are in excellent agreement with real data from Spain. The results in these figures are better than the results obtained in Figure 4 and Figure 5. By comparing our results with the results in [34] and using the same data where the real data in this case come from Spain, we noticed that our results match with the real data better than the results in [34]. Figure 6 shows the comparison between real data from WHO for Wuhan versus the simulation of the proposed model (14–16). Furthermore, by comparing our results with the results in [7,8] and the same data where the real data in this case come from Wuhan, we noticed that our results match with the real data better than the results in [7,8]. Moreover, Figure 7 shows how the behavior of solutions changes according to the different values $\alpha ,$ ${\sigma }_i$ and a new behavior appears. Furthermore, Figure 8 shows how the behavior of solutions changes according to the different values of ${\alpha }_1,$ ${\sigma }_i$ and a new behavior appears. Figure 9, Figure 10 and Figure 11 show how the behavior of solutions changes according to the different values of $H^{*},{\sigma }_i$ and different values of $\alpha ,{\alpha }_1.$

7. Conclusions

A new COVID-19 hybrid fractional order using piecewise derivatives has been formulated in this article. The new hybrid fractional order operator is written as a linear combination of the fractional order integral of Riemann–Liouville, and the fractional order Caputo derivation is applied to extend the deterministic model. Furthermore, the fractional Brownian motion and the hybrid fractional order operator are applied to extend the stochastic differential equations. The basic reproduction number $R_0$ is proved. The existence, uniqueness, boundedness, and positivity of the solutions of the proposed model are established. In order to be compatible with the physical model problem, a new parameter $\mu$ is added. The Caputo proportional constant nonstandard modified Euler–Maruyama method is presented to study the fractional stochastic model, and the Grünwald–Letnikov nonstandard finite difference method is presented to study the hybrid fractional order deterministic model. The real data from Spain and Wuhan for the COVID-19 infection are considered to obtain the parameter’s numerical values. We utilized the Grünwald–Letnikov nonstandard finite difference method for the approximation solutions of the proposed model. The Grünwald–Letnikov nonstandard finite difference method has good properties in stability for solving fractional systems such as these, provides accurate approximation solutions, and saves computational time when the final time is very long. We studied the proposed model numerically and presented the graphical results. Some phase portraits in a stochastic environment are presented. The obtained results lead us to believe that this new approach of model and numerical methods can study complicated real-world problems in the future. In the future, the present study can be extended to optimal control and to examine the impact of multiple vaccination strategies on the dynamics of COVID-19 in a population.