The risk-return relationship in South Africa: tail optimization of the GARCH-M approach

Nitesha Dwarika; Nitesha Dwarika

doi:10.3934/DSFE.2022020

Data Science in Finance and Economics

2022, Volume 2, Issue 4: 391-415. doi: 10.3934/DSFE.2022020

Previous Article Next Article

Research article

The risk-return relationship in South Africa: tail optimization of the GARCH-M approach

Nitesha Dwarika ^,

Independent Researcher, 15 Sisal Place, Crossmoor, Chatsworth, Durban, 4092, South Africa

Received: 16 September 2022 Revised: 27 October 2022 Accepted: 07 November 2022 Published: 14 November 2022
JEL Codes: C32, C58, G11, G32

The risk-return relationship is of fundamental significance in the field of economics and finance. It is used to structure investment strategies, allocate resources, as well as assist in the construction of policy and regulatory frameworks. The accurate forecast of the risk-return relationship ensures sound financial decisions, whereas an inaccurate one can underestimate risk and thus lead to losses. The GARCH-M approach is one of the foremost models used in South African literature to investigate the risk-return relationship. This study made a novel and significant contribution, on a local and international level, as it was the first study to investigate GARCH-M type models with different innovation distributions. This study analyzed the JSE ALSI returns of the South African market for the sample period from 05 October 2004 to 05 October 2021. Results revealed that the EGARCH (1, 1)-M with the Skewed Student-t distribution (Skew-t) is optimal relative to the standard GARCH, APARCH and GJR. However, the EGARCH-M Skew-t failed to capture the financial data's asymmetric, volatile and random nature. To improve forecast accuracy, this study applied different nonnormal innovation distributions: the Pearson Type Ⅳ distribution (PIVD), Generalized Extreme Value distribution (GEVD), Generalized Pareto distribution (GPD) and Stable. Model diagnostics revealed that the nonnormal innovation distributions adequately captured asymmetry. The Value at Risk and backtesting procedure found that the PIVD, followed by Stable, outperformed the Extreme Value Theory distributions (GEVD and GPD). Thus, investors, risk managers and policymakers would opt to use the EGARCH-M in combination with the PIVD when modelling the risk-return relationship. The main contribution of this study was to confirm that applying GARCH type models with the conventional and normal type distributions, to a volatile emerging market, is considered ineffective. Therefore, this study recommended the exploration of other innovation distributions, that were not included in the scope of this study, for future research purposes.

Keywords:

Citation: Nitesha Dwarika. The risk-return relationship in South Africa: tail optimization of the GARCH-M approach[J]. Data Science in Finance and Economics, 2022, 2(4): 391-415. doi: 10.3934/DSFE.2022020

Related Papers:

[1]	Mustafa Gürbüz, Yakup Taşdan, Erhan Set . Ostrowski type inequalities via the Katugampola fractional integrals. AIMS Mathematics, 2020, 5(1): 42-53. doi: 10.3934/math.2020004
[2]	Shuang-Shuang Zhou, Saima Rashid, Muhammad Aslam Noor, Khalida Inayat Noor, Farhat Safdar, Yu-Ming Chu . New Hermite-Hadamard type inequalities for exponentially convex functions and applications. AIMS Mathematics, 2020, 5(6): 6874-6901. doi: 10.3934/math.2020441
[3]	Hong Yang, Shahid Qaisar, Arslan Munir, Muhammad Naeem . New inequalities via Caputo-Fabrizio integral operator with applications. AIMS Mathematics, 2023, 8(8): 19391-19412. doi: 10.3934/math.2023989
[4]	Dazhao Chen . Weighted boundedness for Toeplitz type operator related to singular integral transform with variable Calderón-Zygmund kernel. AIMS Mathematics, 2021, 6(1): 688-697. doi: 10.3934/math.2021041
[5]	Haoliang Fu, Muhammad Shoaib Saleem, Waqas Nazeer, Mamoona Ghafoor, Peigen Li . On Hermite-Hadamard type inequalities for $n$ -polynomial convex stochastic processes. AIMS Mathematics, 2021, 6(6): 6322-6339. doi: 10.3934/math.2021371
[6]	Tekin Toplu, Mahir Kadakal, İmdat İşcan . On n-Polynomial convexity and some related inequalities. AIMS Mathematics, 2020, 5(2): 1304-1318. doi: 10.3934/math.2020089
[7]	Muhammad Imran Asjad, Waqas Ali Faridi, Mohammed M. Al-Shomrani, Abdullahi Yusuf . The generalization of Hermite-Hadamard type Inequality with exp-convexity involving non-singular fractional operator. AIMS Mathematics, 2022, 7(4): 7040-7055. doi: 10.3934/math.2022392
[8]	Gültekin Tınaztepe, Sevda Sezer, Zeynep Eken, Sinem Sezer Evcan . The Ostrowski inequality for $s$ -convex functions in the third sense. AIMS Mathematics, 2022, 7(4): 5605-5615. doi: 10.3934/math.2022310
[9]	Suriyakamol Thongjob, Kamsing Nonlaopon, Sortiris K. Ntouyas . Some (p, q)-Hardy type inequalities for (p, q)-integrable functions. AIMS Mathematics, 2021, 6(1): 77-89. doi: 10.3934/math.2021006
[10]	Chunyan Luo, Yuping Yu, Tingsong Du . Estimates of bounds on the weighted Simpson type inequality and their applications. AIMS Mathematics, 2020, 5(5): 4644-4661. doi: 10.3934/math.2020298

Abstract

1. Introduction

The term "infectious diseases" refers to illnesses caused by living organisms, usually parasites and microorganisms, that spread from one host to another either directly or indirectly. Infectious diseases have been a significant concern in our society. These illnesses can sometimes pose severe threats and lead to epidemics. Medical research has rapidly advanced in identifying and controlling three primary characteristics: Infectivity, epidemic potential, and uncertainty, which have helped reduce the burden of such diseases. Infectious diseases are modelled mathematically as a dynamic transmission cycle involving interactions between susceptible and diseased hosts. These interactions are usually expressed as a coupled ordinary differential equation. The SIR model takes into account several variables while examining infectious disease systems. Several models that examine significant aspects in the research of infectious diseases include SIS, SIRV, SEIR, SVEIRS, SIQR, SEIRV, SEIQS, and others ^{[1,2,3,4,5,6,7]}. Understanding disease dynamics is complex due to factors like environmental fluctuations, making it challenging to collect precise data. For example, temperature variations can significantly impact the reproduction rate of disease-spreading bacteria, influencing the disease's dynamics. As a result, it's essential to consider how uncertainty in biological parameters affects disease behavior. If incorrect parameter values are used, model outputs may be biased, so sensitivity analysis is necessary for understanding how model outputs change in response to parameter variations. When studying the transmission of infectious diseases, sensitivity analysis is often used to develop strategies for slowing infection spread by targeting key model parameters, such as the primary reproduction number $R_0$ . The literature provides information on dynamic parameter sensitivities in infectious spread for traditional SIR and SEIR models to maximize the usefulness of observed data, with the primary source of uncertainty being the random variability of the included parameters.

Interactions between species can impact the quality of life for all involved. These interactions are crucial for an organism's survival since its primary objective is to stay alive. Without these relationships, life as we know it would be impossible. It's remarkable how fascinating organism interactions can be and can be studied using mathematical models, an exciting field within ecology research ^[5,7]. Different factors, such as environmental conditions, impact the growth and balance of living populations, and modeling can be an invaluable tool for understanding their dynamics. Organisms, typically proto-microorganisms and parasites can cause infectious illnesses through direct or indirect transmission from one host to another. Due to their reason to periodically spread like epidemics and cause serious problems, infectious diseases have recently become a significant worry in our culture. The three major components that lessen the severity of these disorders are infectivity, epidemic potential, and uncertainty. The medical community has made great strides in recognising and managing these components. The dynamic transmission cycle is a mathematical model of an infectious illness that involves interactions between susceptible and infected hosts. Typically, these interactions are expressed as a linked set of ordinary differential equations.

Fractional calculus is a mathematical analysis that extends beyond differentiation and integration conventions. Fractional calculus concerns derivatives and integrals of non-integer orders, such as fractions or decimals, while classical calculus concentrates on the powers of integers. The Caputo operator is a valuable tool for problem-solving and is commonly used to handle real-world challenges. The main difference between the Riemann-Liouville operator and other fractional operators is the singularity quality of the kernels. They are applied in roughly related ways to observe different model dynamics. It follows that using these operations would inevitably result in better outcomes. Researchers attempting to understand better a model's dynamics advise using fractional operators with non-singular kernels. Establishing the C-F derivative has addressed certain standard limits on fractional derivatives ^[8]. The C-F derivative outperforms the most commonly used techniques for assessing human knowledge ^[9]. The operator's performance and its benefits were validated using various approaches. Studies on numerous scientific, engineering, and mathematical models have proved their usefulness. However, experts from real-world cases are a more meaningful example ^[10,11]. Atangana and Baleanu introduced a derivative operator with a non-singular kernel. It utilizes the Mittag-Leffer function and is in a convenient area for modellers of real-world situations because of its non-local and non-singular kernel. Over the past few decades, the Atangana-Baleanu variant has obtained strong research value and applicability in diverse fields, particularly in biological models ^[12]. This fractional derivative introduced non-local and complex dynamic behavior and numerous natural results. In addition, since biological models are expected to have a memory effect and hereditary features, which can be defined more precisely using fractional calculus, the solution of the fractional order system is expected to obtain how to control epidemic diseases ^[13]. The authors of this study investigated the sensitivity of the model solution in delay differential systems using variational and direct methods ^[14]. Corruption and terrorism have become significant issues in many nations worldwide. However, more needs to be written about the relationship between the two. To address this, the author has developed a novel fractional-order mathematical model to explore the coexistence of terrorism and corruption ^[15]. This article examines the authors' mathematical model of COVID-19, which incorporates a fractional-order system and considers the effectiveness of vaccination ^[16].

Fractional-order models have become increasingly popular in control engineering, physics, neural networks, and medicine. They can predict the condition of a system at any future moment by considering its current state and all of its past states. Due to its practicality, researchers are interested in fractional-order calculus. The concept of fractional-order differentiation and integration was first introduced by two prominent figures in mathematics history, Riemann and Liouville. The R-L and Caputo fractional derivatives are examples of fractional differential operators that use a power law kernel. These operators have precise definitions and have been widely researched and applied across various disciplines ^[17]. However, their limited applicability and inability to represent other natural and artificial systems have hindered their development. To overcome these limitations, novel C-F fractional-order derivatives have been developed based on the exponential kernel, successfully describing various real-world events ^[18]. Several mathematical models with C-F fractional order have been discussed in the literature.

Numerous numerical techniques have been developed to solve biological models and fractional-order differential equations ^[19]. The Newton interpolation formula, the Toufik-Atangana approach, the Adam-Bashforth method, the predictor-corrector method, and the Runge-Kutta technique are a few of the numerous numerical algorithms that have been the subject of extensive literature ^{[17,18,19,20,21,22,23]}. To solve non-linear systems numerically, the Toufik-Atangana approach is widely used ^[24,25]. Recently, a novel numerical system called the Toufik-Atangana numerical system has emerged as a potential solution to the problems with the Adams-Bashforth approach ^[26,27]. Combining the fundamental theorem of fractional calculus with the two-step Lagrange polynomial, this new methodology develops a novel and very efficient numerical method ^[26]. Using this approach, problems can be resolved quickly and accurately. A biological model has been subjected to many fractional operators using the Toufik-Atangana numerical framework to get the required findings. Stability analysis is an essential method in numerical analysis that yields effective results. When expressed mathematically, stability is necessary for spectrum analysis of real-world problems ^[28]. Hyers expounded on Ulam's 1940 introduction of the U-H stability hypothesis in 1941. As previously stated, the optimal approximation or exact solution to the issue is the basis for computing the stability. Furthermore, it is simple to implement and evaluate the recommended stability ideas ^[28,29].

The following subjects in this work: The definitions and basic ideas of fractional calculus are covered in Section 2. We examine a non-integer C-F model for infectious diseases in Section 3. The positivity and sensitivity analysis of the model under discussion is covered in detail in Section 4. In Section 5, we focuses on the proposed model solution's existence, uniqueness, and stability in the U-H. In Section 6, we present the numerical method for the C-F operator Adams-Bashforth scheme. Section 7 features numerical findings and a graphic analysis. Finally, in Section 8, we discuss the findings of the investigation.

2. Definitions and basic concepts

The essential definitions and theorems relevant to the C-F have been provided in this section.

Definition 2.1. ^[18,30] Suppose that $\vartheta$ $\in$ $H^1(q_1, q_2),$ $q_1 < q_2,$ and $\varsigma \in (0, 1),$ so the C-F fractional differential operator. Then, we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{B(\varsigma)}{1-{\varsigma}}\int_{q_1}^{\mathfrak{k}} \vartheta^{'}(\mathfrak{z}){exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}, \end{align}$

(2.1)

where $B(\varsigma)$ is the normalization function with $B(0) = B(1) = 1.$

However, if $\vartheta$ $\notin$ $H^1(q_1, q_2)$ , then we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{B(\varsigma)}{1-{\varsigma}}\int_{q_1}^{\mathfrak{k}} \left( \vartheta({{\mathfrak{k}}}) - \vartheta(\mathfrak{z})\right) {exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}. \end{align}$

(2.2)

Definition 2.2. ^[18,31] Suppose that $\vartheta$ $\in$ $H^1(q_1, q_2),$ $q_1 < q_2,$ and $\varsigma \in (0, 1),$ so the C-F fractional differential operator. Then, we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{(2-\varsigma)B(\varsigma)}{2(1-{\varsigma})}\int_{q_1}^{\mathfrak{k}} \vartheta^{'}(\mathfrak{z}){exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}. \end{align}$

(2.3)

Definition 2.3. ^[18,32] Suppose that $\varsigma$ is the order of integral in the C-F integral operator. Then, we have

$\begin{align} {}^{CF}{}_{\; }I^{\varsigma}_{{\mathfrak{k}}}[\vartheta({\mathfrak{k}})] = \frac{2(1-\varsigma)}{(2-\varsigma)B({\varsigma})}\vartheta ({\mathfrak{k}}) + \frac{2\varsigma}{(2-\varsigma)B({\varsigma})}\int_{q_1}^{{\mathfrak{k}}}\vartheta (\mathfrak{z})d\mathfrak{z}. \end{align}$

(2.4)

3. Fractional order of infectious diseases model

This section presents the generalized form of the infectious spread transmission dynamics with specific parameter values for a given population. In this paper, we build upon the fundamental SEIR contagious disease model by expanding it to the SEQIRDV model ^[33], which considers different disease stages and traits. There are seven subpopulations within the overall population $N$ : susceptible $(S)$ , infected $(I)$ , exposed $(E)$ , recovered $(R)$ , quarantined $(Q)$ , dead $(D)$ , and vaccinated $(V)$ . Realistically, the natural death rate caused for each subpopulation is represented by a value $\mu$ . A second parameter, $\Pi$ , indicates the recruitment of susceptible individuals in any infected population at any time $\mathfrak{k}$ . At the disease transmission rate $\beta$ , susceptible people can be exposed and become members of the diseased class. We assume that specific individuals receive vaccinations for a particular infectious disease at a rate of $\nu$ . The rate at which people contract the disease is $\beta$ , and they move on to the affected group after a latent period of $\gamma$ . Based on the vaccination's efficiency measure $\sigma$ , the interaction of those who have had vaccinations falls into this category. With a predetermined period until death $\rho$ , the infected population is confined for a proposed average length of $\delta$ , either entering the dead population with a disease mortality rate of $\tau$ or the recovered class with a recovery rate of $\omega$ . Depending on the availability of vaccines for the particular condition, the class "vaccinated" is added. The flowchart of the infectious dynamical disease system (3.1) is represented by Figure 1. After accounting for these variables, the dynamical system that results is as follows ^[33]:

$\begin{align} \frac{dS}{d{\mathfrak{k}}} & = \Pi - \beta SI - \upsilon S - \mu S, \\ \frac{dE}{d{\mathfrak{k}}} & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ \frac{dI}{d{\mathfrak{k}}} & = \gamma E - \delta I - \mu I, \\ \frac{dQ}{d{\mathfrak{k}}} & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ \frac{dR}{d{\mathfrak{k}}} & = (1-\tau)\omega Q -\mu R, \\ \frac{dD}{d{\mathfrak{k}}} & = \tau \rho Q, \\ \frac{dV}{d{\mathfrak{k}}} & = \upsilon S - \sigma \beta VI - \mu V, \end{align}$

(3.1)

Figure 1. A diagrammatic graph that illustrates the transmission of the disease.

DownLoad: Full-Size Img PowerPoint

with initial conditions $S(0) \geq 0,$ $E(0) \geq 0,$ $I(0) \geq 0,$ $Q(0) \geq 0,$ $R(0) \geq 0,$ $D(0) \geq 0$ and $V(0) \geq 0.$

Fractional derivatives in mathematical biology offer a versatile framework for modeling various biological phenomena, including vaccination processes and recovery dynamics. By extending beyond traditional epidemic and noninfectious disease models, researchers gain a deeper understanding of the complex behaviors exhibited by biological systems. This enables them to develop innovative solutions to critical vaccination research and healthcare issues by incorporating fractional calculus into their mathematical models. Initially, we transformed the model with integer order into one with fractional order. We apply the C-F derivative in the system (3.1), and then we get

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}S & = \Pi - \beta SI - \upsilon S - \mu S, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}E & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}I & = \gamma E - \delta I - \mu I, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}Q & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}R & = (1-\tau)\omega Q -\mu R, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}D & = \tau \rho Q, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}V & = \upsilon S - \sigma \beta VI - \mu V. \end{align}$

(3.2)

4. Preliminary analysis of the model

4.1. Positivity and boundedness

In fractional order modeling of biological systems, the population remains positive and bounded over time to represent biological constraints accurately. Positivity and boundedness are essential in fractional order models to ensure stability, realistic behavior, and relevance to real-world situations. These features are necessary due to the complex dynamics introduced by fractional calculus. Positive states play a crucial role in maintaining a system's stability and meaningful behavior. In contrast, positive outputs are vital for applying and understanding the model's results in practical applications. Boundedness ensures that the system does not exhibit unlimited growth, which is essential for stability and control. Furthermore, boundedness guarantees that the behavior of a system remains predictable and controllable, thus facilitating the analysis of system performance and reaction. The proposed model has a built-in feature where its solutions are always positive and bounded. We ensure that all state variables have non-negative values for any time ${\mathfrak{k}} > 0$ . This ^[34] means a trajectory starting with a positive initial condition will stay favorable for ${\mathfrak{k}} > 0$ . Therefore, system (3.2) provides that

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}S|_{S = 0} & = \Pi \geq 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}E|_{E = 0} = \sigma \beta VI \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}I|_{I = 0} & = \gamma E \geq 0, \; \; \; \; \; \; \; \; \; \; \; \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}Q|_{Q = 0} = \delta I \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}R|_{R = 0} & = (1-\tau)\omega Q \geq 0, \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}D|_{D = 0} = \tau \rho Q \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}V|_{V = 0} & = \upsilon S \geq 0. \end{align}$

(4.1)

Since $\mathbb{N}(t) = S({\mathfrak{k}}) + I({\mathfrak{k}}) + E({\mathfrak{k}}) + Q({\mathfrak{k}}) + D({\mathfrak{k}}) + R({\mathfrak{k}}) + V({\mathfrak{k}})$ is total population, we have

$\begin{align*} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}\mathbb{N}({\mathfrak{k}}) & = \Pi -\mu S - \mu E - \mu I - \mu Q - \mu R - \mu V \nonumber \\ & \leq \Pi - \mu S, \end{align*}$

then one has

$\begin{align*} \mathbb{N}({\mathfrak{k}}) \leq \left( \mathbb{N}(0) - \frac{\Pi}{\mu}\right)E_{\varsigma}(-\mu {\mathfrak{k}}) + \frac{\Pi}{\mu}. \end{align*}$

Therefore, we have

$\begin{align*} \Omega = \left\lbrace \left( S({\mathfrak{k}}), I({\mathfrak{k}}), E({\mathfrak{k}}), Q({\mathfrak{k}}), D({\mathfrak{k}}), R({\mathfrak{k}}), V({\mathfrak{k}})\right) \in \Re^7_+ : 0 \leq \mathbb{N}({\mathfrak{k}}) \leq \frac{\Pi}{\mu} \right\rbrace, \end{align*}$

which provides the feasible region for the infectious diseases model, and $\Omega$ is positively invariant. Thus, the proposed model (3.2) is well-posed mathematically.

4.2. Sensitivity analysis

This study gives us insights into each parameter's importance in the disease's spread. This information is crucial not only for planning experiments but also for integrating data and simplifying complex models. Sensitivity analysis is often used to assess the resilience of model predictions to changes in parameter values, as there are sometimes errors in data collection and assumed parameter values. This tool helps identify the characteristics that significantly impact the threshold $R_0$ and should be the focus of intervention strategies. Sensitivity indices provide a way to measure the proportional change in a variable resulting from a change in the parameter. We use a variable's normalized forward sensitivity index concerning a specific parameter to achieve this objective. This index is defined as the ratio of the relative change in the variable to the relative change in the parameter. The sensitivity index can be determined using partial derivatives if a variable is differentiable concerning the parameter. We conducted a sensitivity study in this section to determine the impact of each parameter on the $R_0$ . A sensitivity analysis can also determine the critical parameters vital to illness management. This method ascertains the relative contribution of each parameter value to the $R_0$ . Therefore, we have

$\begin{align*} R_0 = \frac{\gamma \beta \Pi (\mu + \sigma \nu )}{\mu (\nu + \mu)(\delta + \mu)(\gamma + \mu)} . \end{align*}$

The system's qualitative characteristics are being considered while organising the relevant parameters based on their impact on the value of $R_0$ , which is quite beneficial. The results of this investigation will help identify the most crucial disease control factors. We can utilize sensitivity indices to determine the effect of a parameter change on the state variable. These indices are calculated using the definition provided in the paper ^[34,35]. To estimate the sensitivity index, we use partial derivatives, as shown below:

$\begin{align} \Lambda^{R_0}_{\amalg} = \frac{\partial R_0}{\partial \amalg}\times \frac{\amalg}{R_0}. \end{align}$

(4.2)

The sensitivity indices for different parameters are present in . These indices have been calculate using the starting values, except for $\sigma$ . shows the sensitivity indices of $R_0$ for the considered parameters of interest. The results indicate that $\Pi$ , $\beta$ , $\delta$ , and $\mu$ are highly significant characteristics, as seen in and . Based on the sensitivity analysis results, it was found that the $R_0$ values increase or decrease in direct proportion to changes in the values of $\Pi,$ $\beta,$ $\sigma,$ $\delta,$ $\mu,$ $\nu$ and $\gamma$ . The reproduction number within the defined limits has been analyzed, considering the following factors. After reviewing the analyses and graphs, it is evident that precautionary measures should be taken to prevent the spread of the disease by improving adverse conditions and reducing factors that contribute to increased reproduction numbers. is the plot of $R_0$ versus the $\beta$ and $\Pi$ . is the plot of $R_0$ versus the $\gamma$ and $\delta$ . is the plot of $R_0$ versus the $\Pi$ and $\mu$ . is the plot of $R_0$ versus the $\sigma$ and $\nu$ . Notably, we have observed that $\mu$ is quite sensitive, and raising this parameter could result in a notable drop in the $R_0$ value. Consequently, limiting these factors can aid in the prevention of the spread of illness. The variables are considered when examining the relevant reproduction number within the specified limits. Based on the evaluations and visual representations, it is concluded that appropriate measures should be implemented to prevent the spread of the disease. This can be achieved by decreasing the factors that contribute to the positive growth rate of reproduction numbers and raising the factors that have a detrimental effect on it. Hence, it can be concluded that increasing awareness about quarantine and vaccination among affected individuals can significantly reduce the spread of infection within the population.

Table 1. Sensitivity indices.

Parameters	Sensitivity index
$\Pi$	$1$
$\beta$	$1$
$\sigma$	$0.1429$
$\delta$	$-0.9444$
$\mu$	$-1.2506$
$\gamma$	$0.0556$
$\nu$	$0.1395$

| Show Table

DownLoad: CSV

Figure 2. Sensitivity plot for the basic reproduction number

$R_0$ .

DownLoad: Full-Size Img PowerPoint

Figure 3. Sensitivity analysis of

$R_0$ according to the model parameters for (a)

$\beta$ and

$\Pi$ , (b)

$\gamma$ and

$\delta$ (c)

$\Pi$ and

$\mu$ , (d)

$\sigma$ and

$\nu$ .

DownLoad: Full-Size Img PowerPoint

5. Essential existence and uniqueness of the solutions

Here, we examine the uniqueness and existence of solutions using the fixed point theorem, which is essential for the proposed model ^[31,36]. The C-F fractional integral operator, when applied to system (3.2), yields

$\begin{align} S({\mathfrak{k}}) - S(0) = &\frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \Pi - \beta SI - \upsilon S - \mu S \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \Pi - \beta S(\mathfrak{Y})I(\mathfrak{Y}) - \upsilon S(\mathfrak{Y}) - \mu S(\mathfrak{Y}) \right] d(\mathfrak{Y}), \\ E({\mathfrak{k}}) - E(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \beta SI - \gamma E + \sigma \beta VI - \mu E \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \beta S(\mathfrak{Y})I(\mathfrak{Y}) - \gamma E(\mathfrak{Y}) + \sigma \beta V(\mathfrak{Y})I(\mathfrak{Y}) - \mu E(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ I({\mathfrak{k}}) - I(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \gamma E - \delta I - \mu I \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \gamma E(\mathfrak{Y}) - \delta I(\mathfrak{Y}) - \mu I(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ Q({\mathfrak{k}}) - Q(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \delta I(\mathfrak{Y}) - (1-\tau)\omega Q(\mathfrak{Y}) - \tau \rho Q(\mathfrak{Y}) - \mu Q(\mathfrak{Y}) \right] d(\mathfrak{Y}), \\ R({\mathfrak{k}}) - R(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \gamma E - \delta I - \mu I \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \gamma E(\mathfrak{Y}) - \delta I(\mathfrak{Y}) - \mu I(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ D({\mathfrak{k}}) - D(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \tau \rho Q \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \tau \rho Q(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ V({\mathfrak{k}}) - V(0) = &\frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \upsilon S - \sigma \beta VI - \mu V \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \upsilon S(\mathfrak{Y}) - \sigma \beta V(\mathfrak{Y})I(\mathfrak{Y}) - \mu V(\mathfrak{Y}) \right] d(\mathfrak{Y}). \end{align}$

(5.1)

Define the following kernels:

$\begin{align} M_1({\mathfrak{k}}, S) & = \Pi - \beta SI - \upsilon S - \mu S, \\ M_2({\mathfrak{k}}, E) & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ M_3({\mathfrak{k}}, I) & = \gamma E - \delta I - \mu I, \\ M_4({\mathfrak{k}}, Q) & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ M_5({\mathfrak{k}}, R) & = (1-\tau)\omega Q -\mu R, \\ M_6({\mathfrak{k}}, D) & = \tau \rho Q, \\ M_7({\mathfrak{k}}, V) & = \upsilon S - \sigma \beta VI - \mu V, \end{align}$

(5.2)

then we have

$\begin{align} S({\mathfrak{k}}) - S(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S) \right] d\mathfrak{Y}, \\ E({\mathfrak{k}}) - E(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_2({\mathfrak{k}}, E) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_2(\mathfrak{Y}, E)\right] d\mathfrak{Y}, \\ I({\mathfrak{k}}) - I(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_3({\mathfrak{k}}, I) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_3(\mathfrak{Y}, I) \right] d\mathfrak{Y}, \\ Q({\mathfrak{k}}) - Q(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_4({\mathfrak{k}}, Q) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_4(\mathfrak{Y}, Q) \right] d\mathfrak{Y}, \\ R({\mathfrak{k}}) - R(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_5({\mathfrak{k}}, R) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_5(\mathfrak{Y}, R)\right] d\mathfrak{Y}, \\ D({\mathfrak{k}}) - D(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_6({\mathfrak{k}}, D) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_6(\mathfrak{Y}, D)\right] d\mathfrak{Y}, \\ V({\mathfrak{k}}) - V(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_7({\mathfrak{k}}, V) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_7(\mathfrak{Y}, V) \right] d\mathfrak{Y}. \end{align}$

(5.3)

Theorem 5.1. If the dissimilarity listed below holds:

$\begin{equation*} 0 \leq \beta \varpi_2 + \upsilon + \mu < 1. \end{equation*}$

Afterwards, the contraction mapping and Lipschitz condition are convinced by the kernel $M_1$ .

Proof. Assume that $S$ and $S_1$ are any two functions, then we get

$\begin{align*} \|M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_1)\| & = \left\| \left( \Pi - \beta S({\mathfrak{k}})I({\mathfrak{k}}) - \upsilon S({\mathfrak{k}}) - \mu S({\mathfrak{k}}) \right) \right. - \left. \left( \Pi - \beta S_1({\mathfrak{k}})I({\mathfrak{k}}) - \upsilon S_1({\mathfrak{k}}) - \mu S_1({\mathfrak{k}})\right) \right\| \\ & \leq \left( \beta \varpi_2 + \upsilon + \mu \right)\left\|S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| \\ & \leq {\Upsilon}_1\left\|S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| . \end{align*}$

Let ${\Upsilon}_1 = \beta \varpi_2 + \upsilon + \mu$ , we do suppose that $S, E, I, Q, R, D$ and $V$ is bounded functions, i.e., $\|S({\mathfrak{k}})\| \leq \varpi_1,$ $\|E({\mathfrak{k}})\|\leq \varpi_2,$ $\|I({\mathfrak{k}})\|\leq \varpi_3,$ $\|Q({\mathfrak{k}})\|\leq \varpi_4,$ $\|R({\mathfrak{k}})\|\leq \varpi_5,$ $\|D({\mathfrak{k}})\|\leq \varpi_6$ and $\|V({\mathfrak{k}})\|\leq \varpi_7.$

$\begin{align} \|M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_1)\| & \leq \Upsilon_1 \left\| S({\mathfrak{k}}) - S_1({\mathfrak{k}})\right\| . \end{align}$

(5.4)

Hence, kernel $M_1$ satisfies the Lipschitz condition and since

$0 \leq \beta \varpi_2 + \upsilon + \mu < 1,$

then it is also a contraction for $M_1$ .

Similarly, kernels $M_2,$ $M_3,$ $M_4,$ $M_5,$ $M_6$ , and $M_7$ satisfy the Lipschitz condition, respectively, as follows:

$\begin{align} \|M_2({\mathfrak{k}}, E) - M_2({\mathfrak{k}}, E_{1})\| & \leq \Upsilon_2 \left\| E({\mathfrak{k}}) - E_{1}({\mathfrak{k}})\right\|, \\ \|M_3({\mathfrak{k}}, I) - M_3({\mathfrak{k}}, I_{1})\| & \leq \Upsilon_3 \left\| I({\mathfrak{k}}) - I_{1}({\mathfrak{k}})\right\|, \\ \|M_4({\mathfrak{k}}, Q) - M_4({\mathfrak{k}}, Q_{1})\| & \leq \Upsilon_4 \left\| Q({\mathfrak{k}}) - Q_{1}({\mathfrak{k}})\right\|, \\ \|M_5({\mathfrak{k}}, R) - M_5({\mathfrak{k}}, R_{1})\| & \leq \Upsilon_5 \left\| R({\mathfrak{k}}) - R_{1}({\mathfrak{k}})\right\|, \\ \|M_6({\mathfrak{k}}, D) - M_6({\mathfrak{k}}, D_{1})\| & \leq \Upsilon_6 \left\| D({\mathfrak{k}}) - D_{1}({\mathfrak{k}})\right\|, \\ \|M_7({\mathfrak{k}}, V) - M_7({\mathfrak{k}}, V_{1})\| & \leq \Upsilon_7 \left\| V({\mathfrak{k}}) - V_{1}({\mathfrak{k}})\right\|. \end{align}$

(5.5)

Equations (5.4) and (5.5) are used to apply the previously described kernels, which transform system (5.1) into

$\begin{align} S({\mathfrak{k}}) & = S(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_1({\mathfrak{k}}, S) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_1(\mathfrak{Y}, S)d\mathfrak{Y}, \\ E({\mathfrak{k}}) & = E(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_2({\mathfrak{k}}, E) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_2(\mathfrak{Y}, E)d\mathfrak{Y}, \\ I({\mathfrak{k}}) & = I(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_3({\mathfrak{k}}, I) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_3(\mathfrak{Y}, I)d\mathfrak{Y}, \\ Q({\mathfrak{k}}) & = Q(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_4({\mathfrak{k}}, S) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_4(\mathfrak{Y}, Q)d\mathfrak{Y}, \\ R({\mathfrak{k}}) & = R(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_5({\mathfrak{k}}, R) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_5(\mathfrak{Y}, R)d\mathfrak{Y}, \\ D({\mathfrak{k}}) & = D(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_6({\mathfrak{k}}, D) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_6(\mathfrak{Y}, D)d\mathfrak{Y}, \\ V({\mathfrak{k}}) & = V(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_7({\mathfrak{k}}, V) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_7(\mathfrak{Y}, V)d\mathfrak{Y}. \end{align}$

(5.6)

We now introduce the following recursive formulas:

$\begin{align} S_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_1({\mathfrak{k}}, S_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_1(\mathfrak{Y}, S_{r-1})d\mathfrak{Y}, \\ E_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_2({\mathfrak{k}}, E_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_2(\mathfrak{Y}, E_{r-1})d\mathfrak{Y}, \\ I_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_3({\mathfrak{k}}, I_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_3(\mathfrak{Y}, I_{r-1})d\mathfrak{Y}, \\ Q_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_4({\mathfrak{k}}, Q_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_4(\mathfrak{Y}, Q_{r-1})d\mathfrak{Y}, \\ R_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_5({\mathfrak{k}}, R_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_5(\mathfrak{Y}, R_{r-1})d\mathfrak{Y}, \\ D_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_6({\mathfrak{k}}, D_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_6(\mathfrak{Y}, D_{r-1})d\mathfrak{Y}, \\ V_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_7({\mathfrak{k}}, V_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_7(\mathfrak{Y}, V_{r-1})d\mathfrak{Y}, \end{align}$

(5.7)

where the initial condition are

$\begin{align*} S_0({\mathfrak{k}}) & = S(0), \; E_0({\mathfrak{k}}) = E(0), \; I_0({\mathfrak{k}}) = I(0), \; Q_0({\mathfrak{k}}) = Q(0), \\ R_0({\mathfrak{k}}) & = R(0), \; D_0({\mathfrak{k}}) = D(0), \; V_0({\mathfrak{k}}) = V(0). \nonumber \end{align*}$

With regard to the recursive formulas, the differences between successive terms can be represented as

$\begin{align*} \chi_{1, r}({\mathfrak{k}}) & = S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{2, r}({\mathfrak{k}}) & = E_r({\mathfrak{k}}) - E_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_2({\mathfrak{k}}, E_{r-1}) - M_2({\mathfrak{k}}, E_{r-2})\right] \nonumber \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_2(\mathfrak{Y}, E_{r-1})- M_2(\mathfrak{Y}, E_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{3, r}({\mathfrak{k}}) & = I_r({\mathfrak{k}}) - I_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_3({\mathfrak{k}}, I_{r-1}) - M_3({\mathfrak{k}}, I_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_3(\mathfrak{Y}, I_{r-1})- M_3(\mathfrak{Y}, I_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align} \chi_{4, r}({\mathfrak{k}}) & = Q_r({\mathfrak{k}}) - Q_{r-1}({\mathfrak{k}}) \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_4({\mathfrak{k}}, Q_{r-1}) - M_4({\mathfrak{k}}, Q_{r-2})\right] \\ & \quad+ \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_4(\mathfrak{Y}, Q_{r-1})- M_4(\mathfrak{Y}, Q_{r-2})\right] d\mathfrak{Y}, \end{align}$

(5.8)

$\begin{align*} \chi_{5, r}({\mathfrak{k}}) & = R_r({\mathfrak{k}}) - R_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_5({\mathfrak{k}}, R_{r-1}) - M_5({\mathfrak{k}}, R_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_5(\mathfrak{Y}, R_{r-1})- M_5(\mathfrak{Y}, R_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{6, r}({\mathfrak{k}}) & = D_r({\mathfrak{k}}) - D_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_6({\mathfrak{k}}, D_{r-1}) - M_6({\mathfrak{k}}, D_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_6(\mathfrak{Y}, D_{r-1})- M_6(\mathfrak{Y}, D_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{7, r}({\mathfrak{k}}) & = V_r({\mathfrak{k}}) - V_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_7({\mathfrak{k}}, V_{r-1}) - M_7({\mathfrak{k}}, V_{r-2})\right] \nonumber \\ & \quad+ \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_7(\mathfrak{Y}, V_{r-1})- M_7(\mathfrak{Y}, V_{r-2})\right] d\mathfrak{Y}. \nonumber \end{align*}$

Note that

$\begin{align} S_r({\mathfrak{k}}) & = \sum\limits_{j = 1}^{r} \chi_{1, j}, \; \; E_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{2, j}, \; \; I_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{3, j}, \; \; Q_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{4, j}, \\ R_r({\mathfrak{k}}) & = \sum\limits_{j = 1}^{r} \chi_{5, j}, \; \; D_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{6, j}, \; \; V_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{7, j}. \end{align}$

(5.9)

After solving system (5.8), we apply the usual supremum norm to both sides of the first equation of system (5.8), then we obtain

$\begin{align} \left\| \chi_{1, r}({\mathfrak{k}})\right\| & = \left\| S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\right\| \\ & = \left\| \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right]\right. \\ &\quad \left. + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\mathfrak{Y}\right\| . \end{align}$

(5.10)

Using Eq (5.10) and the triangle inequality, we get

$\begin{align} \left\| S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\right\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| \left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right]\right\| \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\left\| \int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\omega\right\| . \end{align}$

(5.11)

Thus, using the Lipschitz constant $\Upsilon_1$ to propitiate the Lipschitz condition, the kernel $M_1$ enables us to determine

$\begin{align} \|S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\| \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S_{r-1} - S_{r-2}\right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \int_{0}^{{\mathfrak{k}}} \left\|S_{r-1} - S_{r-2}\right\| d\mathfrak{Y} . \end{align}$

(5.12)

Thus, we obtain

$\begin{align} \|\chi_{1, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1\|\chi_{1, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \int_{0}^{{\mathfrak{k}}}\|\chi_{1, r-1}(\mathfrak{Y})\|d\mathfrak{Y}. \end{align}$

(5.13)

In a similar manner, we get

$\begin{align} \|\chi_{2, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_2\|\chi_{2, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_2 \int_{0}^{{\mathfrak{k}}}\|\chi_{2, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{3, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_3\|\chi_{3, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_3 \int_{0}^{{\mathfrak{k}}}\|\chi_{3, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{4, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_4\|\chi_{4, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_4 \int_{0}^{{\mathfrak{k}}}\|\chi_{4, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{5, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_5\|\chi_{5, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_5 \int_{0}^{{\mathfrak{k}}}\|\chi_{5, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{6, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_6\|\chi_{6, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_6 \int_{0}^{{\mathfrak{k}}}\|\chi_{6, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{7, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_7\|\chi_{7, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_7 \int_{0}^{{\mathfrak{k}}}\|\chi_{7, r-1}(\mathfrak{Y})\|d\mathfrak{Y}. \end{align}$

(5.14)

Theorem 5.2. If there $\exists$ a time ${\mathfrak{k}}_0 > 0$ , then the ensuing disparities are valid:

$\begin{equation} \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}}_0 < 1, \end{equation}$

(5.15)

solutions exist for the infectious disease system.

Proof. Let $S({\mathfrak{k}}),$ $E({\mathfrak{k}}),$ $I({\mathfrak{k}}),$ $Q({\mathfrak{k}}),$ $R(t),$ $D({\mathfrak{k}})$ , and $V({\mathfrak{k}})$ be bounded functions and use the Lipschitz condition. Now, in Eqs (5.13) and (5.14), using the recursive method, we obtain

$\begin{align} \left\| \chi_{1, r}({\mathfrak{k}})\right\| & \leq \left\| S_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_1 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_1{\mathfrak{k}}\right]^r , \\ \left\| \chi_{2, r}({\mathfrak{k}})\right\| & \leq \left\| E_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_2 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_2 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{3, r}({\mathfrak{k}})\right\| & \leq \left\| I_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_3 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_3 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{4, r}({\mathfrak{k}})\right\| & \leq \left\| Q_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_4 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_4 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{5, r}({\mathfrak{k}})\right\| & \leq \left\| R_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_5 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_5 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{6, r}({\mathfrak{k}})\right\| & \leq \left\| D_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_6 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_6 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{7, r}({\mathfrak{k}})\right\| & \leq \left\| V_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_7 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_7 {\mathfrak{k}}\right]^r . \end{align}$

(5.16)

Thus, it is proven that the aforementioned solutions exist and continue to exist. To demonstrate the function is a solution of system (3.2), we make the following assumptions:

$\begin{align} S({\mathfrak{k}}) - S(0) & = S_r({\mathfrak{k}}) - \mathfrak{A}_{1, r}({\mathfrak{k}}), \\ E({\mathfrak{k}}) - E(0) & = E_r({\mathfrak{k}}) - \mathfrak{A}_{2, r}({\mathfrak{k}}), \\ I({\mathfrak{k}}) - I(0) & = I_r({\mathfrak{k}}) - \mathfrak{A}_{3, r}({\mathfrak{k}}), \\ Q({\mathfrak{k}}) - Q(0) & = Q_r({\mathfrak{k}}) - \mathfrak{A}_{4, r}({\mathfrak{k}}), \\ R({\mathfrak{k}}) - R(0) & = R_r({\mathfrak{k}}) - \mathfrak{A}_{5, r}({\mathfrak{k}}), \\ D({\mathfrak{k}}) - D(0) & = D_r({\mathfrak{k}}) - \mathfrak{A}_{6, r}({\mathfrak{k}}), \\ V({\mathfrak{k}}) - V(0) & = V_r({\mathfrak{k}}) - \mathfrak{A}_{7, r}({\mathfrak{k}}). \end{align}$

(5.17)

Then, we have

$\begin{align} \left\| \mathfrak{A}_{1, r}({\mathfrak{k}})\right\| & = \left\| \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{r-1})\right]\right. \\ &\quad \left. + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S)- M_1(\mathfrak{Y}, S_{r-1})\right] d\mathfrak{Y}\right\| \\ & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{r-1})\right\| \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| M_1(\mathfrak{Y}, S)- M_1(\mathfrak{Y}, S_{r-1})\right\| d\mathfrak{Y} \\ & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S - S_{r-1}\right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S - S_{r-1}\right\| {\mathfrak{k}}. \end{align}$

(5.18)

Repeating this process recursively, it becomes

$\begin{align} \left\| \mathfrak{A}_{1, r}({\mathfrak{k}})\right\| \leq \left[ \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)} {\mathfrak{k}} \right]^{r+1} {\Upsilon}^{r+1}_1 a . \end{align}$

(5.19)

At the point ${\mathfrak{k}}_0$ , we get

(5.20)

As $r$ $\rightarrow \infty$ in Eq (5.20), then we get

$\begin{equation*} \|\mathfrak{A}_{1, r}({\mathfrak{k}})\| \rightarrow 0. \end{equation*}$

Similarly, we get

$\begin{equation*} \|\mathfrak{A}_{2, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{3, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{4, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{5, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{6, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{7, r}({\mathfrak{k}})\| \rightarrow 0 . \end{equation*}$

We establish that a system of system (3.2) solutions is unique. Consider the possibility that there is another set of model (3.2) to determine the uniqueness of the solution. We have $S_1({\mathfrak{k}}),$ $E_1({\mathfrak{k}}),$ $I_1({\mathfrak{k}}),$ $Q_1({\mathfrak{k}}),$ $R_1({\mathfrak{k}}),$ $D_1({\mathfrak{k}}),$ and $V_1({\mathfrak{k}});$ then

$\begin{align} S({\mathfrak{k}}) - S_1({\mathfrak{k}}) = &\frac{2(1-{\varsigma})}{(2-{\varsigma})B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{1})\right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S) - M_1(\mathfrak{Y}, S_{1})\right] d\mathfrak{Y}. \end{align}$

(5.21)

When we apply the norm to Eq (5.21), we obtain

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\| & \leq \frac{2(1-{\varsigma})}{(2-\varsigma) B(\varsigma)}\left\| \left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{1})\right] \right\| \\ & \quad+ \frac{2\varsigma}{(2- \varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| M_1(\mathfrak{Y}, S) - M_1(\mathfrak{Y}, S_{1})\right\| d\mathfrak{Y}. \end{align}$

(5.22)

Using the kernel's Lipschitz condition, one can achieve

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S({\mathfrak{k}}) - S_{1}({\mathfrak{k}}) \right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \left\| S({\mathfrak{k}}) - S_{1}({\mathfrak{k}}) \right\|. \end{align}$

(5.23)

By reducing the complexity of Eq (5.23), we get

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\|\left[ 1 - \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 - \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \right] \leq 0. \end{align}$

(5.24)

Theorem 5.3. The proposed model has a unique solution if the following conditions are hold:

$\begin{equation} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\|\left[ 1 - \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 - \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \right] > 0. \end{equation}$

(5.25)

Proof. Using Eq (5.24), we have

(5.26)

Implying that

$\begin{align} \left\| S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| = 0, \end{align}$

(5.27)

we obtain

$\begin{align} S({\mathfrak{k}}) = S_1({\mathfrak{k}}). \end{align}$

(5.28)

Continuing in the same manner, we have

$\begin{equation} \begin{array}{ll} E({\mathfrak{k}}) & = E_{1}({\mathfrak{k}}), \; \; I({\mathfrak{k}}) = I_{1}({\mathfrak{k}}), \; \; Q({\mathfrak{k}}) = Q_{1}({\mathfrak{k}}), \\ R({\mathfrak{k}}) & = R_{1}({\mathfrak{k}}), \; \; D({\mathfrak{k}}) = D_{1}({\mathfrak{k}}), \; \; V({\mathfrak{k}}) = V_{1}({\mathfrak{k}}). \end{array} \end{equation}$

(5.29)

As a result, we proved that system (3.2)'s system of solutions is unique.

Stability analysis

We employ nonlinear functional analysis to investigate the Ulam-Hyers (U-H) stability ^[28,37] of the proposed fractional model (3.2). Eighty-four years have passed since Professor Ulam presented the stability problem to the University of Wisconsin Mathematics Club. Ulam asked whether a proposition remains valid or approximately valid when the hypothesis is slightly modified ^[28,38]. This subject is intriguing and significant in multiple scientific disciplines, motivating numerous individuals to pursue its investigation. The current name for this topic is the U-H stability problem. Initially, the focus was on the stability of group homomorphisms. One year later, Hyers resolved the inquiry of the additive mappings across Banach spaces using the "contraction mapping theorem". After this significant advancement, numerous investigations were conducted on Ulam's challenge using diverse approaches and variations. Rassias extended the findings of Hyers ^[39]. Instead of using a positive constant, he employed a dominating function to regulate the estimate. This problem is commonly referred to as the Ulam-Hyers-Rassias stability problem or the generalized U-H stability problem.

In recent decades, numerous research articles have been published on U-H stability, focusing on ordinary differential equations (ODEs) ^[40,41]. ODEs result in productive outcomes, encompassing both linear and nonlinear equations ^[42,43]. Many recent studies have focused on the U-H stability concerning fractional systems ^[44]. The concepts of existence, uniqueness, and U-H stability have been established in various functional systems. U-H stability measures the difference in solutions of differential equations with fractional order. This stability ensures that controlled systems remain stable and perform predictably, even with slight variations in the system dynamics. When modeling physical processes with fractional-order dynamics, U-H stability guarantees that the model stays accurate even when minor changes or approximations are made to the governing equations. U-H stability is used to analyze the sensitivity of solutions to small changes in boundary data for fractional differential equations with boundary conditions. This analysis ensures the reliability of the solutions in practical applications. The concept of U-H stability is beneficial for analyzing and implementing fractional-order systems. It provides a framework to ensure that solutions remain stable even when subjected to minor disturbances, which is essential in numerous practical situations. For simplicity, we consider the proposed model (3.2) to be as follows:

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{U}({{\mathfrak{k}}}) & = \mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})), \\ \mathfrak{U}({0}) & = \mathfrak{U}_0 \geq 0, \end{align}$

(5.30)

where

$\begin{align*} &\mathfrak{U}({{\mathfrak{k}}}) = \left( S({\mathfrak{k}}), E({\mathfrak{k}}), I({\mathfrak{k}}), Q({\mathfrak{k}}), R({\mathfrak{k}}), D({\mathfrak{k}}), V({\mathfrak{k}})\right)^T, \\ &\mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})) = \left( M_1, M_2, M_3, M_4, M_5, M_6, M_7\right)^T. \end{align*}$

Applying the fractional integral on (5.30), we obtain

$\begin{align} \mathfrak{U}({\mathfrak{k}}) - \mathfrak{U}(0) = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \mathfrak{U}(\mathfrak{Y}))d\mathfrak{Y} . \end{align}$

(5.31)

Definition 5.1. The proposed model (3.2) is UH stable if there exists $\varsigma > 0$ and along with the condition for any $\epsilon > 0$ and $\bar{\mathfrak{U}} \in \mathbb{B},$ if

$\begin{align} \left| {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{U}({{\mathfrak{k}}}) - \mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}}))\right| \leq \epsilon, \end{align}$

(5.32)

then $\exists$ $\mathfrak{U} \in \mathbb{B}$ , and $\mathbb{B}$ in the model (3.2),

$\begin{align*} \mathfrak{U}(0) = \bar{\mathfrak{U}}(0) = \bar{\mathfrak{U}_0}, \end{align*}$

such that

$\begin{align*} \left\| \bar{\mathfrak{U}} - \mathfrak{U} \right\| \leq \varsigma \epsilon, \end{align*}$

where

$\begin{align*} \begin{cases} & \bar{\mathfrak{U}}({{\mathfrak{k}}}) = \left( \bar{S}({\mathfrak{k}}), \bar{E}({\mathfrak{k}}), \bar{I}({\mathfrak{k}}), \bar{Q}({\mathfrak{k}}), \bar{R}({\mathfrak{k}}), \bar{D}({\mathfrak{k}}), \bar{V}({\mathfrak{k}})\right)^T, \\ & \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) = \left( \bar{M_1}, \bar{M_2}, \bar{M_3}, \bar{M_4}, \bar{M_5}, \bar{M_6}, \bar{M_7}\right)^T, \\ & \epsilon = \max(\epsilon_i)^T; \; i = 1, 2, ..., 7, \\ & \varsigma = \max(\varsigma_i)^T; \; i = 1, 2, ..., 7. \end{cases} \end{align*}$

Remark 5.1. Let a small perturbation $\kappa \in \mathbb{C}[0, \mathbb{T}],$ such that $\kappa(0) = 0$ with the following condition: $\| \kappa({\mathfrak{k}}) \| \leq \bar{\epsilon},$ for ${\mathfrak{k}} \in [0, \mathbb{T}]$ and $\bar{\epsilon} > 0.$

Lemma 5.1. Let the solution $\bar{\mathfrak{U}}_{\kappa}({{\mathfrak{k}}})$ of the perturbed system

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \bar{\mathfrak{U}}({{\mathfrak{k}}}) & = \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{{\mathfrak{k}}}})) + \kappa({\mathfrak{k}}), \\ \bar{\mathfrak{U}}(0) & = \bar{\mathfrak{U}}_0 \end{align}$

(5.33)

hold the condition

$\begin{align*} \left\| \bar{\mathfrak{U}}_\kappa({{\mathfrak{k}}}) - \bar{\mathfrak{U}}({{\mathfrak{k}}})\right\| \leq \Phi \bar{\epsilon}, \end{align*}$

where

$\begin{align*} \Phi & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2 \varsigma}{(2-\varsigma)B(\varsigma)} \mathbb{T}, \\ \kappa({\mathfrak{k}}) & = \left( \kappa_1({\mathfrak{k}}), \kappa_2({\mathfrak{k}}), \kappa_3({\mathfrak{k}}), \kappa_4({\mathfrak{k}}), \kappa_5({\mathfrak{k}}), \kappa_6({\mathfrak{k}}), \kappa_7({\mathfrak{k}})\right)^T. \end{align*}$

Proof. Applying the fractional integral in Eq (5.33), we obtain

$\begin{align} \bar{\mathfrak{U}}_{\kappa}({\mathfrak{k}}) - \bar{\mathfrak{U}}(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \bar{\mathfrak{U}}(\mathfrak{Y}))d\mathfrak{Y} \\ & \quad+ \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\kappa({\mathfrak{k}}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\kappa(\mathfrak{Y})d\mathfrak{Y}. \end{align}$

(5.34)

Also,

$\begin{align} \bar{\mathfrak{U}}({\mathfrak{k}}) & = \bar{\mathfrak{U}}(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \bar{\mathfrak{U}}(\mathfrak{Y}))d\mathfrak{Y}. \end{align}$

(5.35)

Using Remark 5.1, we obtain

$\begin{align*} \left\| \bar{\mathfrak{U}}_{\kappa}({\mathfrak{k}}) - \bar{\mathfrak{U}}({\mathfrak{k}})\right\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left| \kappa({\mathfrak{k}})\right| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left| \kappa(\mathfrak{Y})\right| d\mathfrak{Y} \nonumber \\ & \leq \left( \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2 \varsigma}{(2-\varsigma)B(\varsigma)} \mathbb{T}\right) \bar{\epsilon} \\ & \leq \Phi \bar{\epsilon}. \end{align*}$

This completes the proof.

Theorem 5.4. The proposed model (3.2) is UM stable if

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| \leq \bar{\varsigma \epsilon}. \end{align*}$

Proof. Let $\bar{\mathfrak{U}}$ be the solution of Eq (5.32) and with help of uniqueness, $\mathfrak{U}$ be a unique solution of the system (5.30), then we have

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| & \leq \Phi \bar{\epsilon} + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) - \mathfrak{V}({{{\mathfrak{k}}}}, {\mathfrak{U}}({{\mathfrak{k}}})) \right\| \\ & \; \; \; + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| \mathfrak{V}({\mathfrak{Y}}, \bar{\mathfrak{U}}({\mathfrak{Y}})) - \mathfrak{V}({\mathfrak{Y}}, {\mathfrak{U}}({\mathfrak{Y}})) \right\|d\mathfrak{Y} + \Phi \bar{\epsilon} \\ & \leq 2 \Phi \bar{\epsilon} + \Phi \bar{\delta} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\|. \end{align*}$

Simplifying the above equation, we obtain

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| & \frac{2 \Phi \bar{\epsilon}}{1 - \Phi \bar{\delta}} = \bar{\varsigma \epsilon}, \end{align*}$

where

$\begin{align*} \varsigma = \frac{2 \Phi}{1 - \Phi \bar{\delta}}. \end{align*}$

Hence, the proposed model (3.2) is UM stable.

6. Numerical scheme

We present a numerical solution for the model (3.2) using the Adams-Bashforth (A-B) technique. Owolabi and Atangana et al. ^[45] introduced a three-step A-B technique with the C-F fractional derivative, which we used to determine the numerical scheme for the fractional-order system (3.2). Take into consideration the C-F derivative of the fractional differential equation

$\begin{equation} \begin{array}{ll} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{E}({{\mathfrak{k}}}) = \mathfrak{F}({{{\mathfrak{k}}}}, \mathfrak{E}({{{\mathfrak{k}}}})). \end{array} \end{equation}$

(6.1)

By utilising the C-F fractional integral on both sides of Eq (6.1), we obtain

$\begin{align} \mathfrak{E}({\mathfrak{k}}) - \mathfrak{E}(0) = \frac{(1-\varsigma)}{B(\varsigma)}{\mathfrak{F}}({{{\mathfrak{k}}}}, \mathfrak{E}({{{\mathfrak{k}}}})) + \frac{\varsigma}{B(\varsigma)}\int_{0}^{{\mathfrak{k}}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y}. \end{align}$

(6.2)

We have divided the time interval into smaller intervals with steps of h, then we get ${\mathfrak{k}}_{e+1} = {\mathfrak{k}}_e + h$ , ${\mathfrak{k}}_0 = 0$ , $e = 0, 1, 2, ..., e-1.$ We put ${\mathfrak{k}} = {\mathfrak{k}}_{e+1}$ and ${\mathfrak{k}} = {\mathfrak{k}}_{e}$ in Eq (6.2). We have found the difference in the resulting equations by performing a computation, then we get

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) - \mathfrak{E}({\mathfrak{k}}_{e}) & = \frac{(1-\varsigma)}{B(\varsigma)}\left[ {\mathfrak{F}}({{{\mathfrak{k}}}_{e+1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e+1}})) - {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}}))\right] + \frac{\varsigma}{B(\varsigma)} \int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y} . \end{align}$

(6.3)

We are estimated the integral $\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}({\mathfrak{Y}}, \mathfrak{E}({\mathfrak{Y}}))d\mathfrak{Y}$ with the approximation of $\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}\mathfrak{Q}_2({\mathfrak{Y}})d\mathfrak{Y},$ where $\mathfrak{Q}_2({\mathfrak{Y}})$ is the Lagrange polynomial of points $({\mathfrak{k}}_{e-2}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-2}}))),$ $({\mathfrak{k}}_{e-1}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-1}})))$ and $({\mathfrak{k}}_{e}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}})))$ . Thus,

$\begin{align} \mathfrak{Q}_2(\mathfrak{Y}) = \sum\limits_{j = 0}^{j = 2}{\mathfrak{F}}({{{\mathfrak{k}}}_{e-j}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-j}}))\mathfrak{L}_{j}(\mathfrak{Y}), \end{align}$

(6.4)

where $\mathfrak{L}_{j}(\mathfrak{Y})$ is the Lagrange basis polynomials on the point ${\mathfrak{k}}_e, {\mathfrak{k}}_{e-1}, {\mathfrak{k}}_{e-2}.$

Let $\mathfrak{E}_{e} = \mathfrak{E}({{{\mathfrak{k}}}_{e}}),$ and $v = \frac{{\mathfrak{k}}_{e+1}-\mathfrak{Y}}{h},$ after replacing the Lagrange basis polynomials and performing integration, we obtain the following result:

$\begin{align} \int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y} & = h \int_{0}^{1}\left( {\mathfrak{F}}({\mathfrak{k}}_e, \mathfrak{E}_{e})\frac{(v-2)(v-3)}{(1-2)(1-3)}\right. \\ & \quad\left. + {\mathfrak{F}}({\mathfrak{k}}_{e-1}, \mathfrak{E}_{e-1})\frac{(v-1)(v-3)}{(2-1)(2-3)} + {\mathfrak{F}}({\mathfrak{k}}_{e-2}, \mathfrak{E}_{e-2})\frac{(v-2)(v-1)}{(3-2)(3-1)}\right) dv \\ & = \frac{23h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e}, \mathfrak{E}_{e}) - \frac{16h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e-1}, \mathfrak{E}_{e-1}) + \frac{5h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e-2}, \mathfrak{E}_{e-2}). \end{align}$

(6.5)

Using Eq (6.5) in Eq (6.3), we have

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) - \mathfrak{E}({\mathfrak{k}}_{e}) = & \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}})) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-1}})) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-2}})). \end{align}$

(6.6)

The error with this technique is

$\begin{align} \mathfrak{R}^{\varsigma}_{e}({\mathfrak{k}}) = \frac{\varsigma}{B(\varsigma)}\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}\frac{3h^3}{8}\mathfrak{F}^4(s)ds = \frac{3\varsigma h^3}{8B(\varsigma)}{\mathfrak{F}^3}(x_e, {\mathfrak{E}(x_e)}), x_e \in ({\mathfrak{k}}_e, {\mathfrak{k}}_{e+1}). \end{align}$

(6.7)

The proposed model's numerical simulations are generated using the three-step A-B approach for the C-F fractional derivative in Eq (6.6).

Let the vectors

$\mathfrak{E}({\mathfrak{k}}) = \left[ S({\mathfrak{k}}), E({\mathfrak{k}}), I({\mathfrak{k}}) Q({\mathfrak{k}}), R({\mathfrak{k}}), D({\mathfrak{k}}), V({\mathfrak{k}})\right],$

and

$\mathfrak{F}({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) = \left[ \mathfrak{F}_1({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_2({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_3({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_4({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_5({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_6({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_7({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}}))\right],$

where the following are the functions specified by system (3.2), we have

$\begin{align} \begin{cases} \mathfrak{F}_1({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \Pi - \beta SI - \upsilon S - \mu S, \\ \mathfrak{F}_2({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ \mathfrak{F}_3({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \gamma E - \delta I - \mu I, \\ \mathfrak{F}_4({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ \mathfrak{F}_5({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; (1-\tau)\omega Q -\mu R, \\ \mathfrak{F}_6({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \tau \rho Q, \\ \mathfrak{F}_7({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \upsilon S - \sigma \beta VI - \mu V. \end{cases} \end{align}$

(6.8)

We can express model (3.2) in vector form as shown below:

(6.9)

Equation (6.6) is utilized to obtain the solution of model (3.2), which is expressed through the iterative formula given below:

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) = & \mathfrak{E}({\mathfrak{k}}_{e}) + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.10)

Let $\mathfrak{E}_0 = \mathfrak{E}({\mathfrak{k}}_0) = \left[S({\mathfrak{k}}_0), E({\mathfrak{k}}_0), I({\mathfrak{k}}_0), Q({\mathfrak{k}}_0), R({\mathfrak{k}}_0), D({\mathfrak{k}}_0), V({\mathfrak{k}}_0) \right]^{T},$ $\mathfrak{E}_{e-2} = \mathfrak{E}({\mathfrak{k}}_{e-2}),$ $\mathfrak{E}_{e-1} = \mathfrak{E}({\mathfrak{k}}_{e-1}),$ $\mathfrak{E}_{e} = \mathfrak{E}({\mathfrak{k}}_{e})$ and $\mathfrak{E}_{e+1} = \mathfrak{E}({\mathfrak{k}}_{e+1}),$ then

$\begin{align} \mathfrak{E}_{e+1} = & \mathfrak{E}_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.11)

We have obtained the iterative formula:

$\begin{align*} S_{e+1} = & S_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} E_{e+1} = & E_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} I_{e+1} = & I_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align} Q_{e+1} = & Q_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align}$

(6.12)

$\begin{align*} R_{e+1} = & R_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} D_{e+1} = & D_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} V_{e+1} = & V_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}). \end{align*}$

Stability analysis of numerical scheme

Here, we are analysing the stability of the fractional A-B scheme

(6.13)

Now, rewrite Eq (6.11), then we have

$\begin{align} \mathfrak{E}_{e+1} & = \mathfrak{E}_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ &\quad - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}) \\ & = \mathfrak{E}_{e} + A \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) - B \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + C \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.14)

Using Eq (5.30), we have

$\begin{align} \mathfrak{E}_{e+1} = (1+A) \mathfrak{E}_{e} - B \mathfrak{E}_{e-1} + C \mathfrak{E}_{e-2}. \end{align}$

(6.15)

Afterward, we apply the Von-Neumann stability analysis for the terms in the aforementioned equation.

$\begin{align} \mathfrak{E}_{e+1} & = \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{(e+1)i \Delta \mathfrak{k}}, \; \; \mathfrak{E}_{e} = \bar{\mathfrak{E}_{e}}\mathfrak{e}^{(e)i \Delta \mathfrak{k}}, \\ \mathfrak{E}_{e-1} & = \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{(e-1)i \Delta \mathfrak{k}}, \; \; \mathfrak{E}_{e-2} = \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{(e-2)i \Delta \mathfrak{k}}. \end{align}$

(6.16)

So that

$\begin{align} \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{(e+1)i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{e}}\mathfrak{e}^{e i \Delta \mathfrak{k}} - B \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{(e-1)i \Delta \mathfrak{k}} + C \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{(e-2)i \Delta \mathfrak{k}}, \end{align}$

(6.17)

which reduces to

$\begin{align} \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{e}} - B \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{- i \Delta \mathfrak{k}} + C \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{-2i \Delta \mathfrak{k}}. \end{align}$

(6.18)

A recursive formula is applied, when $e = 0,$ we have

$\begin{align} \bar{\mathfrak{E}_{1}}\mathfrak{e}^{i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{0}}. \end{align}$

(6.19)

Simplifying the previous equation (6.15), we have

$\begin{align*} \left\| \bar{\mathfrak{E}_{1}}\mathfrak{e}^{i \Delta \mathfrak{k}}\right\| & = \left\| (1 + A)\bar{\mathfrak{E}_{0}} \right\|, \\ \implies \left\| \bar{\mathfrak{E}_{1}}\right\| & = (1 + A) \left\| \bar{\mathfrak{E}_{0}} \right\|. \end{align*}$

$\begin{align*} h \leq 0.55 \; \; \text{and} \; \; \varsigma > \frac{1}{1 - 1.9167h}, \end{align*}$

then, we suppose that

$\begin{align*} \forall e \geq \bar{\mathfrak{E}_{e}} < \bar{\mathfrak{E}_{0}}, \end{align*}$

we have

$\begin{align} \left| \bar{\mathfrak{E}_{e+1}}\right| & \leq (1 + A)\left| \bar{\mathfrak{E}_{e}}\right|\left| \mathfrak{e}^{- i \Delta \mathfrak{k}}\right| - B \left| \bar{\mathfrak{E}_{e-1}}\right| \left| \mathfrak{e}^{- 2i \Delta \mathfrak{k}}\right| + C \left| \bar{\mathfrak{E}_{e-2}}\right| \left| \mathfrak{e}^{-3i \Delta \mathfrak{k}}\right| \\ & < (1 + A)\left| \bar{\mathfrak{E}_{e}}\right| - B \left| \bar{\mathfrak{E}_{e-1}}\right| + C \left| \bar{\mathfrak{E}_{e-2}}\right|. \end{align}$

(6.20)

Simplifying the previous equation (6.15), we have

$\begin{align} \left| \bar{\mathfrak{E}_{e+1}}\right| & < (1 + A)\left| \bar{\mathfrak{E}_{0}}\right| - B \left| \bar{\mathfrak{E}_{0}}\right| + C \left| \bar{\mathfrak{E}_{0}}\right| \\ & < (1 + A - B + C)\left| \bar{\mathfrak{E}_{0}}\right| \\ & < (1 + A - B + C) < 1. \end{align}$

(6.21)

Hence, when applied to Eq (5.30), the three-step A-B method with the C-F derivative is conditionally stable.

7. Numerical results and graphical analysis

In this section, we simulate the non-linear fractional proposed model using a novel numerical technique that uses the C-F non-integer operator. We displayed graphical findings for various fractional orders and infection rates during the simulation. To comprehend the system's behavior and operation, tracking and controlling the infection rate is imperative. Utilizing the numerical approach outlined in Section 6, we followed the instructions and finished the simulation results. We examined numerical findings for various fractional orders and time intervals using MATLAB R2016a. As such, we can see how adjusting the parameters and starting conditions affects the model's predictions through the simulations and improves our grasp of the dynamics of the model. The fractional order analysis's graphical results have been more illuminating and broadly applicable than other related efforts. The following parameter values have been applied to the simulation: $S(0) = 100,$ $E(0) = 50,$ $I(0) = 30,$ $Q(0) = 20,$ $R(0) = 18,$ $D(0) = 2,$ $V(0) = 60,$ $\Pi = 20,$ $\beta = 2 \times 10^{-4},$ $\nu = 1 \times 10^{-3},$ $\mu = 0.3,$ $\gamma = 5.1,$ $\sigma = 0.05,$ $\delta = 5.1,$ $\tau = 0.04,$ $\omega = 0.03$ and $\rho = 15.$ For this study, we have chosen specific orders to represent the complex nature of the system under examination, as arbitrary order derivatives offer greater degrees of freedom and provide a wide range of geometries. depicts the state variable plots of time graphically in the suggested model using the ABM method for different orders $(\varsigma)$ : $0.99, 0.96, 0.93$ and $0.90$ . The fractional suggested model (3.2) for various values of $\varsigma$ has been numerically simulated and shown in Figures 5 and 6.

Figure 4. Solution trajectories for the proposed model (3.2) when the derivative order varies.

DownLoad: Full-Size Img PowerPoint

Figure 5. The plots for the proposed model (3.2) when varies values of

$\varsigma$ .

DownLoad: Full-Size Img PowerPoint

Figure 6. The plots for the proposed model (3.2) when varies values of

$\varsigma$ .

DownLoad: Full-Size Img PowerPoint

When we place a classical order, the infected population is higher, but when we place an order with $\varsigma = 0.85,$ the infected population is smaller, as illustrated in . When placing a classical order, the recovered and quarantined population is smaller. However, when placing an order with $\varsigma = 0.85,$ the recovered and quarantined population is more significant.

We consider a small order, leading to a more significant impact on the dead population than a more substantial higher order. depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\mu$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\mu$ varies. When we place the value of $\mu = 0.45$ , the susceptible, quarantined, and vaccinated population is smaller. However, when we set $\mu = 0.30$ , the susceptible, quarantined, and vaccinated population is more significant, as shown in . If we choose the smaller value for $\mu$ , we have large populations of these populations. The dead population becomes more extensive if we choose the smaller value for $\mu$ . depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ , $\delta$ varies and $\nu$ varies. The infected and quarantined populations are sensitive to parameter $\delta$ increasing order. The exposed and infected populations are sensitive to parameter $\nu$ . depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\beta$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\beta$ varies. When we place the value of $\beta$ is small, we get the susceptible, infected, and recovered population is more significant compared to the more considerable value of $\beta$ . The dead population becomes more extensive if we choose the smaller value for $\beta$ . depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\tau$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\tau$ varies. The exposed and recovered populations are not very sensitive to parameter $\tau$ . We consider a small value for $\tau$ , leading to a more significant impact on the quarantined population than a more substantial value of $\tau$ .

Figure 7. Three-dimensional numerical results for proposed model (3.2) when

$\mu$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 8. The plots for the proposed model (3.2) when

$\mu$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 9. The plots for the proposed model (3.2) when

$\delta$ and

$\nu$ are varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 10. Three-dimensional numerical results for proposed model (3.2) when

$\beta$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 11. The plots for the proposed model (3.2) when

$\beta$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 12. Three-dimensional numerical results for proposed model (3.2) when

$\tau$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 13. The plots for the proposed model (3.2) when

$\tau$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

The numerical simulation suggests we gain more reliable insights into behavior within specific ranges with more precise parameter information. We can explore the model's behavior under different input parameter values. Therefore, this study highlights the advantages of using numerical techniques in real-world disease models.

8. Conclusions

In this paper, we discuss the SEIQRDV model, which focuses on the spread of infectious diseases and considers positivity, boundedness, and sensitivity. We have created a fractional order model based on a described integer order contagious disease model. Using fixed-point theory, we have proven the uniqueness and existence of solutions to this model under specific circumstances. We have also used a novel numerical technique supported by the A-B scheme to find the numerical solution. By plotting the model's output at different fractional orders and parameter values, we have observed how the population curves respond to these variations. However, it is important to note that the uncertainty resulting from random fluctuations in the model parameters cannot be eliminated. Additionally, determining which parameters are essential in the model can be challenging, causing sensitivity analysis for such a model to be difficult. The sensitivity analysis of the basic reproduction number has shown that the model's parameters have a significant impact. The results indicate that the recruitment and transmission rates are the most influential factors contributing to a substantial increase in the $R_0$ . These findings provide valuable insights into the mechanisms that affect the dynamics of infectious disease transmission and offer strategies for reducing disease transmission.

We observe that even slight changes in parameters such as $\beta$ , $\delta$ , $\mu$ , $\nu$ , and $\tau$ result in different population dynamics. These efforts will enhance our understanding of this paradigm and enable us to obtain more information successfully. Moreover, the results indicate that as the number of infectious individuals in the community decreases, the number of recovered individuals in the system increases. This suggests a correlation between the decrease in the number of sick individuals and the increase in the number of those who have effectively recuperated from the illness. The proposed model provides valuable insights into disease transmission dynamics using fractional-order derivatives. The results of our investigation are crucial for improving the accuracy of the infectious disease model and formulating successful strategies. Furthermore, our investigation indicated that meaningful parameters can be found and examined even in uncertainty. These research results suggest expanding this work to encompass control theory and other fractional operators.

Author contributions

Parveen Kumar: Formal analysis, Writing–original draft, Methodology, Software, Validation; Sunil Kumar: Supervision, Conceptualization, Investigation, Writing–review and editing; Badr Saad T Alkahtani: Validation, Writing–review and editing, Conceptualization, Visualization; Sara S Alzaid: Investigation, Formal analysis, Visualization, Writing–original draft. All authors have read and agreed to the published version of the manuscript.

Acknowledgement

This research was funded by Researchers Supporting Project number (RSPD2024R526), King Saud University, Riyadh, Saudi Arabia.

Funding

This research was funded by Researchers Supporting Project number (RSPD2024R526), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

All authors declare no conflicts of interest.

References

[1]	Adu G, Alagidede P, Karimu A (2015) Stock return distribution in the BRICKS. Rev Dev Financ 5: 98–109. https://doi.org/10.1016/j.rdf.2015.09.002 doi: 10.1016/j.rdf.2015.09.002
[2]	Alqaralleh H, Abuhommous AA, Alsaraireh A (2020) Modelling and Forecasting the Volatility of Cryptocurrencies: A Comparison of Nonlinear GARCH-Type Models. Intl J Financ Res 11: 346–356. https://doi.org/10.5430/ijfr.v11n4p346 doi: 10.5430/ijfr.v11n4p346
[3]	Ayed WB, Fatnassi I, Maatoug AB (2020) Selection of Value-at-Risk models for MENA Islamic indices. J Islamic Account Bus Res 11: 1689–1708. 10.1108/JIABR-07-2019-0122 doi: 10.1108/JIABR-07-2019-0122
[4]	Bautista RS, Mata LM (2020) A conditional heteroscedastic VaR approach with alternative distributions. EconoQuantum 17: 81–98. https://doi.org/10.18381/eq.v17i2.7125 doi: 10.18381/eq.v17i2.7125
[5]	Bautista RS, Mora JAN (2020) Value at risk in the oil sector: an analysis of the efficiency in the measurement of the risk of the α-Stable distribution versus the generalized asymmetric Student-t and normal distributions. Contaduría y Administración 65: 1–19. https://doi.org/10.1016/j.physa.2020.124876 doi: 10.1016/j.physa.2020.124876
[6]	Brooks C (2014) Introductory Econometrics for Finance. Cambridge University Press: New York.
[7]	Christoffersen P, Pelletier D (2004) Backtesting Value-at-Risk: A Duration-Based Approach. J Financ Econ 2: 84–108. http://dx.doi.org/10.1093/jjfinec/nbh004 doi: 10.1093/jjfinec/nbh004
[8]	Christoffersen PF (1998) Evaluating interval forecasts. Int Econ Rev 39: 841–862. https://doi.org/10.2307/2527341 doi: 10.2307/2527341
[9]	Delis MD, Savva CS, Theodossiou P (2021) The impact of the coronavirus crisis on the market price of risk. J Financ Stab 53(100840): 1–12. https://doi.org/10.1016/j.jfs.2020.100840 doi: 10.1016/j.jfs.2020.100840
[10]	Ding X, Granger CWJ, Engle, RF (1993) A long memory property of stock market returns and a new model. J Empir Financ 1: 83–106. https://doi.org/10.1016/0927-5398(93)90006-D doi: 10.1016/0927-5398(93)90006-D
[11]	Dwarika N, Moores-Pitt P, Chifurira R (2021) Volatility dynamics and the risk-return relationship in South Africa: A GARCH approach. Investment Manage Financ Innovations 18: 106–117. http://dx.doi.org/10.21511/imfi.18(2).2021.09 doi: 10.21511/imfi.18(2).2021.09
[12]	Echaust K, Just M (2020) Value at Risk Estimation Using the GARCH-EVT Approach with Optimal Tail Selection. Math 8: 1–24. http://dx.doi.org/10.3390/math8010114 doi: 10.3390/math8010114
[13]	Emenogu NG, Adenomon MO, Nweze NO (2020) On the volatility of daily stock returns of Total Nigeria Plc: evidence from GARCH models, value-at-risk and backtesting. Financ Innovation 6: 2–25. https://doi.org/10.1186/s40854-020-00178-1 doi: 10.1186/s40854-020-00178-1
[14]	Engle RF (1982) Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation. Econometrica 50: 987–1007. https://doi.org/10.2307/1912773 doi: 10.2307/1912773
[15]	Engle RF, Lilien D, Robins R (1987) Estimating time varying risk premia in the term structure: the ARCH-M model. Econometrica 55: 391–407. https://doi.org/10.2307/1913242 doi: 10.2307/1913242
[16]	Glosten LR, Jagannathan R, Runkle DE (1993) On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks. J Financ 48: 1779–1801. https://doi.org/10.1111/j.1540-6261.1993.tb05128.x doi: 10.1111/j.1540-6261.1993.tb05128.x
[17]	Ilupeju YE (2016) Modelling South Africa's market risk using the APARCH model and heavy-tailed distributions. Master's thesis. University of KwaZulu-Natal. Durban. South Africa.
[18]	IRESS (2022) Available from: https://www.iress.com/.
[19]	Khan M, Aslam F, Ferreira P (2021) Extreme Value Theory and COVID-19 Pandemic: Evidence from India. Econo Res Guardian 11: 2–10.
[20]	Kuang W (2020) Dynamic VaR forecasts using conditional Pearson type Ⅳ distribution. J Forecasting 2020:1–12. https://doi.org/10.1002/for.2726 doi: 10.1002/for.2726
[21]	Kupiec PH (1995) Techniques for verifying the accuracy of risk measurement models. J Deriv 3: 73–84. https://doi.org/10.3905/jod.1995.407942 doi: 10.3905/jod.1995.407942
[22]	Ma C, Mi X, Cai Z (2020) Nonlinear and time-varying risk premia. China Econ Rev 62: 1–29. https://doi.org/10.1016/j.chieco.2020.101467
[23]	Mandimika NZ, Chinzara Z (2012) Risk–return trade-off and behaviour of volatility on the South African stock market: evidence from both aggregate and disaggregate data. South African J Econ 80: 345–365.
[24]	Mangani R (2008) Modelling return volatility on the JSE Securities Exchange of South Africa. African Financ J 10: 55–71.
[25]	Markowitz H (1952) Portfolio selection. J Financ 7: 77–91.
[26]	Naradh K, Chinhamu K, Chifurira R (2021) Estimating the value-at-risk of JSE indices and South African exchange rate with Generalized Pareto and Stable distributions. Investment Manage Financ Innovations 18: 151–165. http://dx.doi.org/10.21511/imfi.18(3).2021.14 doi: 10.21511/imfi.18(3).2021.14
[27]	Nelson DB (1991) Conditional heteroscedasticity in asset returns: A new approach. Econometrica 59: 347–70. https://doi.org/10.2307/2938260 doi: 10.2307/2938260
[28]	Omari C, Mundia S, Ngina I (2020) Forecasting Value-at-Risk of Financial Markets under the Global Pandemic of COVID-19 Using Conditional Extreme Value Theory. J Math Financ 10: 569–597. https://doi.org/10.4236/jmf.2020.10403 doi: 10.4236/jmf.2020.10403
[29]	Saddah S, Sitanggang ML (2020) Value at risk estimation of exchange rate in banking industry. Jurnal Keuangan dan Perbankan 24: 474–484. https://doi.org/10.26905/jkdp.v24i4.4808 doi: 10.26905/jkdp.v24i4.4808
[30]	Tsay RS (2013) An introduction to analysis of financial data with R. Wiley Series in Probability and Statistics. Wiley and Sons. Springer. New Jersey.

Reader Comments

Your name:*

Email:*
© 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)