A higher order Galerkin time discretization scheme for the novel mathematical model of COVID-19

1.   Introduction

Throughout human history, viral infections have always played a vital role, and there have been a bunch of new epidemics that resulted in the devastation of thousands of individuals. For example, consider the outbreak known as the Spanish Flu, which killed millions of people all over the world. Many of the diseases that kill thousands of people each year, such as HIV/AIDS, become endemic. Likewise, coronavirus pandemics have been reported in the last couple of decades ^[1,2,3,4,5]. An infectious disease due to a new variant called COVID-19 viral infection, popularly known as Respiratory Syndrome COVID-19-2, formerly called 2019-nCoV, was reportedly noticed as a respiratory disorder in the city of Hubei province of China in 2019 ^[6,7,8]. More than 800 people died when SARS, a type of coronavirus outbreak, struck, resulting in more than 8000 positive cases. MERS, or Middle East Respiratory Syndrome, has reportedly migrated from its origin in the Kingdom of Saudi Arabia to many surrounding and distant countries, especially countries around the Persian Gulf. MERS is already a factor in a few cases ^[9]. On December 31, 2019, the World Health Organization (WHO) received the report for the first time. It was officially announced by WHO as COVID-19 on January 30th, 2020, as a pandemic that had drastically been affecting the world ^[10,11]. According to reports, the disease first spread from animals to people. According to some studies, bats carry the coronavirus, which has been spread to mankind, as demonstrated by many of the patients who were originally detected working in Wuhan's city ^[12]. The pandemic gained exponential momentum in the months of January and February 2020, when the entire origin was sealed in Wuhan. Later on, positive cases were also reported in the US and across the Pacific and Atlantic. Then it was observed that the disease is communicable and can be transferred from one person to the other by physical contact. The WHO declared the virus to be an epidemic in mid-March 2020. According to reports on March 7th, 2021, the pandemic has expanded to many countries around the world, resulting in more than 116.17 million confirmed cases and 2.58 million deaths so far ^[13].

Depending on how it is affected by sunlight, climate change, and the surface material, COVID-19 can exist on a surface for hours or days. COVID-19 is spread by interacting with a virus-infected material or item, then touching one's own lips, nose, or eyes. Unfortunately, that is not the only way that the virus can spread. Outside the home, social separation reduces the risks of coming into contact with infected surfaces or contagious humans ^[14]. In order to effectively minimize COVID-19 transmission, governments developed a variety of control measures, including strict and mandatory lockdowns, as well as social distancing, avoiding crowded meetings, imposing a minimum number of participants in meetings, as well as using face masks. In order to efficiently contain the outbreak of COVID-19, the majority of member nations have fortified interaction verification of known infected individuals, and incidences are immediately isolated for medical care ^[15].

The World Health Organization denotes March 11, 2020, as the official start of the COVID-19 pandemic. Though millions of people worldwide had already been living under severe restrictions by that point, this date marked the start of a global lock down period that was often described as unprecedented. This period was marked by pain-the loss of life and suffering of millions, but also fear of unabated transmission, stressed healthcare systems, and strained and scarce resources. Further, the pandemic intersected with pre-existing vulnerabilities and axes of inequality, including by race, class, geography, and (dis)ability status; the impacts of this tragedy are highly uneven ^[16]. Sunthrayuth et al. ^[17] studied to predict the COVID-19 infection fifth wave in South Africa using the Gaussian mixture model for the available data of the early four waves. The quantification data is considered, and the time unit is used in days. They presented the modeling of COVID-19 in South Africa and predict the future fifth wave in the country. Initially, they used the Gaussian mixture model to characterize the coronavirus infection to fit the early reported cases of four waves and then to predict the future wave. Actual data and the statistical analysis using the Gaussian mixture model are performed which give close agreement with each other. The differential equations model is simulated for various values of the model parameters in order to determine the disease's possible eliminations. Vijayalakshmi et al. ^[18] studied the effectiveness of vaccination in COVID-19 pandemic disease by modelling three compartments susceptible, vaccinated and infected (SVI) of Atangana Baleanu of Caputo (ABC) type derivatives in non-integer order. The disease dynamics is analysed and its stability is performed. Numerical approximation is derived using Adam's Moulton method and simulated to forecast the results for controllability of pandemic spread. Cerqueti et al. ^[19] discussed the cluster analysis of selected countries based on COVID-19 new deaths per million data. They implemented a statistical procedure that combines a rank-size exploration and a k-means approach for clustering. Specifically, they first carried out a best-fit exercise on a suitable polynomial rank-size law at an individual country level; then, they cluster the considered countries by adopting a k-means clustering procedure based on the calibrated best-fit parameters. On January 30, 2020, China had 7734 tested positive incidents of CIVID-19, in addition to 90 international incidents exported around the world to approximately 13 countries, like India, France, Germany, Canada, the United Arab Emirates, and the United States. By October 31st, 2020, over 4,667,780 positive cases were found around the world (Asia: 13,461,293 cases, Africa: 1,776,595 cases, Europe: 9,840,736 cases, America: 20,546,580 cases, Oceania: 41,880 cases, and others: 696 cases), with 1,189,499 deaths (Asia: 239,675 deaths, Africa: 42,688 deaths, Europe: 265,565 deaths, America: 640,513 deaths, and others: 7 deaths) ^[19,20].

Mathematical models are beneficial for understanding and comprehending the transmission of contagious and infectious diseases. It is a valuable tool for evaluating real-world phenomena and processes ^{[21,22,23,24,25]}. Bernoulli was the first mathematician to propose mathematical modelling of contagious and infectious disease transmission in 1760. Later on, the issue piqued the attention of a large number of scientists and researchers. These models simplify the understanding of a wide variety of physical and biological phenomena and their processes. This field has advanced significantly, with models ranging from simple to elaborate and convoluted. Numerous infectious and non-contagious disorders have been investigated using mathematical models (see for example ^[22,23]). Researchers utilize computational equations to understand and analyze the dynamics of a disease (see ^[26,27] for more information). The researchers used nonlinear numerical analytic techniques to determine the global and local stability of the endemic and negative pool equilibria (disease-free-equilibria). In a similar fashion, researchers lately conducted an extensive investigation into the novel COVID-19 using mathematical models from different perspectives. Stability theory, numerical simulation, and global and local dynamics are all part of the research. In this regard, we have alluded to some excellent work, such as ^{[28,29,30,31,32,33,34]}.

Researchers have been applying and constructing mathematical models to develop a clear understanding of the mechanism of pandemic spread, transmission, effect, prevention, and handling of the pandemic since the arrival of COVID-19. Okuonghae and Omame ^[35] investigated COVID-19's transmission and social behavior in Lagos, Nigeria. Roseline et al. ^[36] analyzed the pandemic's mortality rate in Nigeria using the linear regression approach. Abdo et al. ^[37] investigated the solution for the model of nonlinear fractional differential equations (FDEs) describing the deadly and widespread coronavirus (COVID-19). The mathematical model based on fourteen nonlinear FDEs is presented, and the numerical results are investigated using the fractional Adams Bashforth (AB) method. In order to realize more effectively, a recently introduced fractional nonlocal operator known as Atangana-Baleanu (AB) is used. The fixed point theorems of Krasnoselskii and Banach are used to demonstrate the model's existence, uniqueness, and stability. The behavior of the approximate solution is presented in terms of graphs through various fractional orders for numerical simulations. The earlier stage of the emergence of COVID-19 in Nigeria was studied by Adegboye et al. ^[38]. Ajisegiri et al. ^[39] investigated the COVID-19 pandemic situation in Nigeria. In Hubei, China, Anastassopoulou et al. ^[40] conducted research on the diagnosis of serious epidemiological restrictions along with demonstrating and predicting the spread of the COVID-19 pandemic. Fanelli and Piazza ^[41] investigated the dynamic behavior of the COVID-19 outbreak in China, Italy, and France. In Wuhan, China, Roda et al. ^[42] connected the traditional SIR and SEIR structures to analyze COVID-19. Al-Ganess et al. ^[43] improved the Adaptive-Neuro-Fuzzy-Inference-System (ANFIS) by using the Salp-Swarm Algorithm to establish an Enhanced Flower Pollination Algorithm to estimate the number of verified COVID-19 situations in China. Using publically available datasets, Wang et al. ^[44] developed the Patient Information Based Algorithm for estimating the COVID-19 death rate in real-time. Nave et al. ^[45] considered the COVID-19 model and implemented the Laplace Adomian decomposition method (LADM) for finding the approximate solutions. Schiøler et al. ^[46] described the Markov model for infection spread and an approximation of a two-stage sampling scheme to infer the probability of extinction. The potential of the method is illustrated via a simulation study. The model is used to assess the Danish containment strategy when SARS-CoV-2 spread from mink to man during the summer of 2020, including the Cluster-5 sub-type. The previously stated COVID-19 research piqued our curiosity in adopting a novel approach to the COVID-19. The following are the key objectives of the present study:

1.1.   The main objectives

The primary contributions of the present article are as follows:

1) To develop a novel mathematical model that is composed of a system of differential equations that represents the population dynamics of susceptible, healthy, infected, quarantined, and recovered individuals.

2) To use the next generation technique and examined the boundedness, local and global behavior of equilibria, and the threshold quantity.

3) To determine the conditions for stability analysis, the linearization and Lyapunov function theories are used.

4) The partial rank correlation coefficient (PRCC) technique is used to perform a sensitivity analysis on the basic reproduction number in order to identify the most significant parameter for regulating the threshold values of the model.

5) Additionally, numerically solved the model using the well-known continuous Galerkin Petrov time discretization approach. In comparison, solved the model using the conventional fourth-order Runge-Kutta approach to determine the validity and reliability of the aforementioned scheme. Finally, showed the computational cost of the suggested scheme in term of time.

1.2.   The outline of the paper

The remaining of the article is summarized as follow: In Section 2, we present the COVID-19 mathematical model and investigate its appropriate mathematical characteristics. In Section 3, we performed a stability analysis, computed the basic reproductive number, and discussed the local and global characteristics of the various equilibria of the proposed model. Implementation of the novel numerical scheme for COVID-19 model is also described in this section. Section 4 discusses comparison and validation of the results obtained through suggested scheme. Section 5 includes the manuscript conclusions and future recommendations.

2.   Basic mathematical model of COVID-19

The model proposed by Ahmad et al. ^[47] consists of four population, that is, susceptible S(t), healthy individuals H(t), infected I(t) and quarantined Q(t). All the model's parameters are non-negative. Susceptible change to the infectious category when the flow rate remains constant. Suspected or infectious people are transfer to quarantine, whereas verified cases are returned to the infectious section for additional hospitalization. The following autonomous ordinary differential equations represent the model under investigation.

Where state-variables' initial conditions are defined as ^[47].

We propose a novel mathematical model based on susceptible individuals S(t), healthy or resistant individuals H(t), infected-I(t), quarantine individuals Q(t), and recovered population R(t). According to new model, when the flow rate remains constant, the susceptibility to infection increases. People who are suspected of being contagious are quarantined, while confirmed patients are sent back to the infectious area for further treatment. Researchers also considered that post-COVID-19 normalcy does not provide lifelong protection, consequently recovered people could be reclassified as susceptible. The new model under consideration is represented by the autonomous ordinary differential equations listed below:

The Models (1)–(5) will be investigated in the biologically feasible initial condition listed below.

If $P\left(t\right)$ represents the total population at time t, then

We get the following when we take the derivative of $P\left(t\right)$ with regard to t.

after simplification, we have

In the absence of infected patients in the quarantine class are treated at a high rate, infected individual treatment rate and Immunity loss rate (${\rm{ \mathsf{ ϱ} }} = 0, {\rm{ \mathsf{ ξ} }} = 0, \; {\rm{ \mathsf{ ϑ} }} = 0$). The Eq (7) become

The Eq (8) is exact differential equation with solution

As a result, it can be concluded that Eq (9) has nonnegative solutions for every $\mathrm{t}\in \left[0, \;\mathrm{\infty }\right).$

Theorem 1. The stochastic Models (1)–(5) demonstrates boundedness with all non-negative initial conditions which were not all exactly zero in the whole region given by

Proof. We assume that ${\Phi } = \left\{\left(\mathrm{S}, \mathrm{H}, \mathrm{I}, \mathrm{Q}, \mathrm{R}\right){\rm{\epsilon }} {R}_{+}^{5}\right\}$ be any solution set of Models (1)–(5) with some non-negative initial condition such as

For $S, H, I, Q\;and\;R,$ this corresponds to any other non-negative initial condition. Because υ is a positive parameter, one may deduce Eq (8)

To solve the above equation, we get

where ${P}_{0}$ is the whole population of the dynamical system at the start. As a result, for $t\to \mathrm{\infty }$, we have

As a result, $P\left(t\right)$ is positive and bounded, and ${\Phi }$ is the greatest collection of positive and bounded solutions. This completes the proof of the theorem.

3.   Stability analysis

We will concentrate on estimating stability and finding the system's probable stationary states in the Models (1)–(5). We analyses a situation in which no one in the community has been infected with the disease. This state is called disease-free-equilibria (DFE). We will refer to such a state as ${E}^{0}$ in the following.

3.1.   DFE state

Substituting $I = Q = R = 0$ in the specific system and solving the autonomous differential Eqs (1) and (2) for $S\left(t\right), H\left(t\right)$ then we get

To compute the endemic equilibrium, we should first calculate the ${R}_{0}$ threshold quantity, which is important to figure out the global behavior of a dynamical system.

3.2.   Basic reproductive number

The basic reproductive number is denoted by ${R}_{0}$ and define as the average amount of infections that occur as a result of the primary infection which can be brought on by an infected person. To derive ${R}_{0}$ we can use the next generation method ^[50]. In the Models (1)–(5), let $X = (I, Q)$ represent the infectious class. We have the ability to write

The Jacobians of the Matrices (12) are

The inverse of the matrix V multiplicatively is

The next generation matrix of the DFE of the suggested is follows:

The eigenvalue of this matrix $F{V}^{-1}$ gives us the threshold parameter ${R}_{0}$. Thus ${R}_{0}$ is given by

3.3.   Local stability at the DFE state

We can derive the stability criteria for Models (1)–(5) before finding the threshold quantity.

Theorem 2. With threshold quantity ${R}_{0} < 1$, the DFE of the Systems (1)–(5) is locally monotonically stable and

For ${R}_{0} > 1$ and in contrast to ${L}_{1}$, the DFE is unstable.

Proof: We taking the Jacobian matrix ($J$) for our model to determine the local stability of ${E}^{0}$, we have

Convert the Matrix (14) into echelon form the last matrix is

where

The Matrix (15) has eigenvalues ${\omega }_{1} = {\omega }_{2} = –\left(d+\upsilon \right) < 0, {\omega }_{5} = 0$ the other two eigenvalues that is ${\omega }_{3} = \frac{\rho \partial +\varphi \omega }{\left(d+\upsilon \right)}-\left(d+\upsilon +\xi +\varrho \right)$ and ${\omega }_{4} = \frac{\propto \xi }{{c}_{0}}-\left(d+\upsilon +\propto \right)$ are negative if

The inequality (16) implies that ${R}_{0} < 1$ and the condition ${L}_{1}$ holds along ${R}_{0} < 1$. Which shows that the Systems (1)–(5) is asymptotically stabile. However, we could see it if ${R}_{0} > 1$ because ${{\rm{ \mathsf{ ω} }}}_{4}$ contains a non-negative real portion, the system is instable. The proof is now complete.

3.4.   Global stability of the DFE state

This section discusses the disease-free state's global stability. Using Castillo-Chaves ^[50], we construct requirements for the global stability of the disease-free-equilibrium. The following two conditions define the global stability of the DFE.

where X represent the susceptible and healthy classes while Y represent the infected, quarantine and recovered class at the DFE point ${E}^{0} = ({X}^{0}, 0)$ i.e.,

The validity of ${E}^{0}$ global asymptotic stability is dependent on the following. ${\chi }_{1}:$ If $\frac{dX}{dt} = F\left(X, 0\right)$ then ${X}^{0}$ is globally stable.${\chi }_{2}:U\left(X, Y\right) = {B}_{Y}-{U}^{*}\left(X, Y\right),$ where ${U}^{*}\left(X, Y\right)\ge 0,$ for $\left(X, Y\right)\in {\Phi }$.

In the case of ${\chi }_{2},$ ${D}_{Y}-U\left({X}^{0}, 0\right)$ is an M-matrix with non-negative off diagonal entries, at which physiologically viable zone is denoted by. As a result, we can prove the following lemma.

Lemma 1. When conditions $\left({\chi }_{1}\right)$ and $\left({\chi }_{2}\right)$ are satisfied, and if ${R}_{0} <$ the dynamical system's equilibrium point ${E}^{0} = ({X}^{0}, 0)$ is globally asymptotically stable.

Theorem 3. The Systems (1)–(5) shows globally asymptotically stable at the DFE point if ${R}_{0} < 1.$

Proof: The $\left({\chi }_{1}\right)and\left({\chi }_{2}\right)$ requirements can be verified to prove shown above theorem. Let the symbol $X = \left(S, H\right)\;and \;\;Y = \left(I, Q, R\right)$ and identify ${E}^{0} = ({X}^{0}, 0)$ where ${X}^{0} = \left(\frac{{\rm{ \mathsf{ ω} }}}{\mathrm{d}+{\rm{ \mathsf{ υ} }}}, \frac{{\rm{ \mathsf{ ρ} }}}{\mathrm{d}+{\rm{ \mathsf{ υ} }}}, \mathrm{0, 0}, 0\right)$ with the use of the Systems (1)–(5), we have

Now if $S = {S}^{0}, H = {H}^{0}, I = Q = R = 0$ and $U\left(X, 0\right) = \left({\rm{ \mathsf{ ω} }}+{\rm{ \mathsf{ ρ} }}\right)–\left(\mathrm{d}+{\rm{ \mathsf{ υ} }}\right)\left(\frac{{\rm{ \mathsf{ ρ} }}+{\rm{ \mathsf{ ω} }}}{\mathrm{d}+{\rm{ \mathsf{ υ} }}}\right) = 0.$ When $\mathrm{t}\to \mathrm{\infty }$ then $X\to {X}^{0}$ this means that $X = {X}^{0}$ is globally asymptotically stable. The second condition is as follows:

where ${B}_{Y}$ is the Jacobian matrix of quarantine and infected groups at $Y = {\left(IQ\right)}^{T}$ and ${U}^{*}\left(X, Y\right)\ge 0$, therefore we can write matrix B in the form of

By taking

we have

Hence the entire population of stochastic Models (1)–(5) is restricted by ${S}^{0}$ and ${H}^{0}$ that is $S, H, I, Q, R\le {S}^{0},$ ${H}^{0}$ from which we have $(\varphi I+\partial I)(S+H)\le \left(\varphi I+\partial I\right)\left({S}^{0}+{H}^{0}\right)$. Therefore ${U}^{*}$ is positive which follows that M-matrix are non-negative. As a result condition $\left({\chi }_{1}\right)$ and $\left({\chi }_{2}\right)$ are satisfied. According to Lemma (1) the ${E}^{0}$ is globally asymptotically stable.

3.5.   Backward bifurcation and endemic equilibrium point

This portion determines the endemic equilibrium point of Models (1)–(5), whenever the number of infected people in the system is greater than zero. The Systems (1)–(5)'s arbitrary endemic equilibrium point is supposed to be ${E}^{*} = \left({S}^{*}, {H}^{*}, {I}^{*}, {Q}^{*}, {R}^{*}\right).$ Then, at the steady-state condition, we can solve the dynamical Systems (1)–(5) simultaneously.

We obviously have ${I}^{*}\ne 0$ for the endemic state. The following quadratic equation in term of ${I}^{*}$ is obtained by substituting ${S}^{*}, {H}^{*}, {Q}^{*}$ and ${R}^{*}$ in the third Eq (3) of the Models (1)–(5) at stationary state.

where ${z}_{1} = \partial {\rm{ \mathsf{ φ} }}, {z}_{2} = \left(\partial +{\rm{ \mathsf{ φ} }}\right)\left(d+{\rm{ \mathsf{ υ} }}\right)-\partial \left({\rm{ \mathsf{ φ} }}{\rm{ \mathsf{ ω} }}+{\rm{ \mathsf{ ϱ} }}\right)$ and ${z}_{3} = {\left(d+{\rm{ \mathsf{ υ} }}\right)}^{2}(1-{R}_{0})$. As the parameter ${\rm{ \mathsf{ φ} }}$ and $\partial$ are positive then surly ${z}_{1}$ is always positive. And ­${z}_{3}$ is negative if ${R}_{0} > 1$ and positive if ${R}_{0} < 1$. There ${z}_{2}$ and ${z}_{3}$ are most important because Eq (23) depends upon the sings of these to get the positive solution. If ${R}_{0} > 1$ then the roots of Eq (23) are positive and real which means that the endemic state is unique. Now if ${R}_{0} = 1$ then ${z}_{3} = 0$ and hence the endemic is either trivial or no endemic state. From the above discussion we conclude that the endemic equilibrium of the model depends on ${R}_{0}$. We have the following interval

If ${z}_{3}\ge 0$ and either ${z}_{2}^{2} < 4{z}_{1}{z}_{3}$ or ${z}_{2}\ge 0$ then there is no positive solution of Eq (23) which mean that there is no endemic equilibrium state. For various range of the parameter, we determine the following.

Theorem 4. The proposed Models (1)–(5) has the following facts:

1) If ${z}_{3} < 0\iff {R}_{0} > 1$ then the Models (1)–(5) have unique state of endemic equilibrium.

2) If ${z}_{2} < 0$ and ${z}_{3} = 0$ or ${z}_{2}^{2}-4{z}_{1}{z}_{3} = 0$ then the Models (1)–(5) have unique state of endemic equilibrium.

3) If ${z}_{2} < 0$, ${z}_{3} > 0$ and ${z}_{2}^{2}-4{z}_{1}{z}_{3} > 0$ then the model have two unique endemic equilibria.

4) Otherwise the Models (1)–(5) has no endemic.

In Theorem (4) case 3 is illustrates the existence of phenomenon of backward bifurcation, when ${R}_{0} < 1$ then the local asymptotic stabilities of both DFE and endemic equilibrium co-exist ^[51,52]. If we set the discriminant ${z}_{2}^{2}-4{z}_{1}{z}_{3} = 0$ to obtain the critical value ${R}_{c}$ of ${R}_{0}$ to probe the backward bifurcation. We have

When ${z}_{2}^{2}-4{z}_{1}{z}_{3} > 0$ or ${R}_{c} < {R}_{0} < 1$ then the backward bifurcation occurs. The backward bifurcation and its epidemiological significance needs the requirement of ${R}_{0} < 1.$ Thus this condition is necessary but not sufficient to eliminate the disease. In similar way the elimination of the disease depends on the initial state variable or the initial size of the sub population. When ${R}_{0} > 1$ the backward bifurcation suggests the feasibility of controlling the disease and depends on the initial stature of the sub-population of our system.

Lemma 2. If case 3 of the Theorem (4) holds along with ${R}_{c} < {R}_{0} < 1$ then the Systems (1)–(5) endures backward bifurcation.

Lemma 3. At ${R}_{0} = 1$ the dynamical Systems (1)–(5) undergoes backward bifurcation $\iff {z}_{2} < 0$.

Proof. We assume that the graph of $f\left(I\right) = {z}_{1}{I}^{2}+{z}_{2}I+{z}_{3}$ for sufficient condition. Now if ${R}_{0} = 1$ then ${z}_{3} = 0$ and hence $f\left(0\right) = 0$ which means that the graph of function passes through the origin. Further-more if ${z}_{2} < 0$ then $f\left(I\right) = 0$ has a positive root i.e $I = \frac{-{z}_{2}}{{z}_{1}}$. If ${z}_{3} > 0$ this means that there exist an open interval $(0, \epsilon$) which contain ${z}_{3}$ where $f\left(I\right) = 0$ has two real and positive root. In short when ${R}_{0} < 1$ then there are two endemic equilibrium states. If ${z}_{2}\ge 0$ the sufficient condition is conspicuous and (23) has no real and positive solution when ${R}_{0} < 1.$

3.6.   Local stability of endemic equilibrium state

Now we have to show that the dynamical Systems (1)–(5) has local stability at the endemic equilibrium point ${E}^{*}$ for this have to prove the following theorem.

Theorem 5. If ${R}_{0} > 1$ and ${\rm{ \mathsf{ ϱ} }} > (\mathrm{d}+{\rm{ \mathsf{ υ} }}+{\rm{ \mathsf{ ξ} }})$ then the unique endemic equilibrium state ${E}^{*}$ of the dynamical Systems (1)–(5) is locally asymptotically stable.

Proof. The Jacobian matrix of the Systems (1)–(5) for unique endemic equilibrium state ${E}^{*} = ({S}^{*}, {H}^{*}, {I}^{*}, {Q}^{*}, {R}^{*})$ is

Now we have to determine the nature of the eigenvalue of the above matrix we convert (26) in echelon form we obtain

where ${C}_{1} = \left(\mathrm{d}+{\rm{ \mathsf{ υ} }}+{\rm{ \mathsf{ ξ} }}+{{\rm P}}\right)-\left({\rm{ \mathsf{ φ} }}{S}^{*}+\partial {H}^{*}+\frac{{\rm{ \mathsf{ φ} }}{S}^{*}{I}^{*}}{{\rm{ \mathsf{ φ} }}{I}^{*}+d+{\rm{ \mathsf{ υ} }}}+\frac{{\partial I}^{*}}{\partial {I}^{*}+d+{\rm{ \mathsf{ υ} }}}\right).$ The eigenvalue of J are ${{\rm{ \mathsf{ ω} }}}_{1} = -\left({\rm{ \mathsf{ φ} }}{I}^{*}+d+{\rm{ \mathsf{ υ} }}\right) < 0, {{\rm{ \mathsf{ ω} }}}_{2} = -\left(\partial {I}^{*}+d+{\rm{ \mathsf{ υ} }}\right) < 0, {{\rm{ \mathsf{ ω} }}}_{3} = \left(\mathrm{d}+{\rm{ \mathsf{ υ} }}+{\rm{ \mathsf{ ξ} }}+{\rm{ \mathsf{ ϱ} }}\right)-\left({\rm{ \mathsf{ φ} }}{S}^{*}+\partial {H}^{*}+\frac{{\rm{ \mathsf{ φ} }}{S}^{*}{I}^{*}}{{\rm{ \mathsf{ φ} }}{I}^{*}+d+{\rm{ \mathsf{ υ} }}}+\frac{{\partial I}^{*}}{\partial {I}^{*}+d+{\rm{ \mathsf{ υ} }}}\right), {{\rm{ \mathsf{ ω} }}}_{4} = -\frac{\propto \;{\rm{ \mathsf{ ξ} }}}{{C}_{1}}-\left(\mathrm{d}+{\rm{ \mathsf{ υ} }}+\propto \right) < 0$ and ${{\rm{ \mathsf{ ω} }}}_{5} = 0$. As ${{\rm{ \mathsf{ ω} }}}_{3} < 0$ if and only if $\left(\mathrm{d}+{\rm{ \mathsf{ υ} }}+{\rm{ \mathsf{ ξ} }}+{\rm{ \mathsf{ ϱ} }}\right) < \left({\rm{ \mathsf{ φ} }}{S}^{*}+\partial {H}^{*}+\frac{{\rm{ \mathsf{ φ} }}{S}^{*}{I}^{*}}{{\rm{ \mathsf{ φ} }}{I}^{*}+d+{\rm{ \mathsf{ υ} }}}+\frac{{\partial I}^{*}}{\partial {I}^{*}+d+{\rm{ \mathsf{ υ} }}}\right).$ Now we use algebra to simplify, we have

Thus if ${\rm{ \mathsf{ ϱ} }} > \left(\mathrm{d}+{\rm{ \mathsf{ υ} }}+{\rm{ \mathsf{ ξ} }}\right)$ and ${\mathrm{R}}_{0} > 1$ then clearly all the coefficient of (28) are positive. Hence, the condition of the theorem is satisfied which means that the system is locally asymptotically stable.

3.7.   Global stability of endemic equilibrium state

In this section we determined the global stability of the endemic equilibrium point of the dynamical Systems (1)–(5) in term of the basic reproduction number ${\mathrm{R}}_{0}$.

Theorem 6. If ${\mathrm{R}}_{0} > 1$ then the endemic equilibrium state ${\mathrm{E}}^{0}$ of Models (1)–(5) is globally asymptotically stable, otherwise unstable.

Proof. For studying of the global stability of the Systems (1)–(5) we construct the Lyapunov function which is the following:

Obviously at the endemic state E, the function $\mathit{\Psi } \left(t\right) = 0$ for the values of $S = {S}^{*},$ $H = {H}^{*},$ $I = {I}^{*},$ $Q = {Q}^{*},$ $R = {R}^{*}$. Clearly when $S > {S}^{*},$ $H > {H}^{*},$ $I > {I}^{*},$ $Q > {Q}^{*},$ $R > {R}^{*}$ the given function is strictly positive which means that the function $\mathit{\Psi } \left(t\right)$ for the given E is positive semi-definite. Now we take the time derivative of the (29) we have

Substituting values from Models (1)–(5) into (30) and then simplifying, we have

Obviously $, \frac{d\mathit{\Psi } \left(t\right)}{dt} < 0$, for $S < 1 = P, H < 1 = P$ and $\left(\frac{{\rm{ \mathsf{ ϱ} }}}{H}+\frac{{\rm{ \mathsf{ ξ} }}}{{\rm{ \mathsf{ ϱ} }}}+d+{\rm{ \mathsf{ υ} }}\right) > \left(\partial +{\rm{ \mathsf{ φ} }}\right)$. Which means that $\frac{d\mathit{\Psi } }{dt}$ of the definite function is negative semi definite. Thus E is the largest invariant sub set and hence the dynamical Systems (1)–(5) is globally asymptotically stable under the condition ${R}_{0} > 1$.

Through a number of assumptions, sensitivity analysis is a technique for evaluating how the values of independent variables change a specific dependent variable. One or even more input variables within the given parameters, as well as the effect of interest rate rises on the price of a bond, will decide its use. The sensitivity analysis outlined in the preceding part aids in the development of an effective control strategy to combat the pandemic. Local and global sensitivity analysis are used to calculate the impact of input values on the output of a system. The relevance of model inputs and their interactions in relation to model output is quantified via global sensitivity analysis. In contrast to local sensitivity analysis, which offers a local perspective of partial derivatives, it provides an overall picture of the effect of inputs on outcomes. Here, we will utilize the most well-known partial rank correlation coefficient (PRCC) method for sensitivity-analysis of the basic reproduction number of the model.

We perform PRCC sensitivity test in order to calculate the influence of input factors on the output of threshold parameter ${\mathrm{R}}_{0}$. The results of our significant test is illustrated in Table 2 and Figure 2, which shows that φ and ω are the most sensitive factors with PRCC values 0.8037 and 0.7295 followed by ${\rm{ \mathsf{ ρ} }}$ and d with PRCC values 0.7602 and 0.6045. Figure 2 and Table 2 illustrate the results of PRCC significant test and show the importance of each parameter in the structure of ${\mathrm{R}}_{0}$.

3.8.   The numerical schemes

This section discusses the numerical schemes used to solve the mentioned above model. Nowadays, the Galerkin approach has been used effectively to handle a wide range of complicated problems in engineering and science, seen for example ^{[52,53,54,55,56,57,58]}. In addition, we employ the conventional Runge Kutta approach to assess the precision of the Galerkin scheme findings.

3.9.   The continuous Galerkin-Petrov (cGP) technique

The system of ODE's for the considered model can be written follows: Find $\widetilde {\mathit{\boldsymbol{u}}}:[0, {t}_{max}]\to \mathit{\boldsymbol{V}} = {\mathbb{R}}^{d}$ like as follows

where ${d}_\mathit{\boldsymbol{t}}$ shows the time derivative of $\widetilde {\mathit{\boldsymbol{u}}}\left(t\right), I = \left[0, T\right]$ is the total interval, $\widetilde {\mathit{\boldsymbol{u}}}\left(t\right) = \left({\widetilde {u}}_{1}\left(0\right), {\widetilde {u}}_{2}\left(0\right), {\widetilde {u}}_{3}\left(0\right)\right)\in V\Rightarrow \widetilde {\mathit{\boldsymbol{u}}}\left(0\right) = \left({\widetilde {u}}_{1}\left(0\right), {\widetilde {u}}_{2}\left(0\right), {\widetilde {u}}_{3}\left(0\right)\right)\in V$ are the initial values of $\widetilde {\mathit{\boldsymbol{u}}}\left(t\right)\;\mathrm{a}\mathrm{t}\;t = 0.$ We also assume that $({\widetilde {u}}_{1}\left(0\right), {\widetilde {u}}_{2}\left(0\right), {\widetilde {u}}_{3}\left(0\right)) = \left(T\right(t), I(t), V(t\left)\right)$ which implies that $({\widetilde {u}}_{1}\left(0\right), {\widetilde {u}}_{2}\left(0\right), {\widetilde {u}}_{3}\left(0\right)) = \left(T\right(0), I(0), V(0\left)\right)$. The function $F = ({f}_{1}, {f}_{2}, {f}_{3})$ is nonlinear and is described as $F:I\times K\to K.$

The weak formulation (see ^{[52,53,54,55,56,57,59]} for explanations) of the problem (32) is follows: find $\widetilde {u}\in X$ such that $\widetilde {\mathit{\boldsymbol{u}}}\left(0\right) = {\widetilde {\mathit{\boldsymbol{u}}}}_{0}$ and

where X and Y represent the solution and test space respectively. To explain the time discretization of a variational type problem (32).

Characterize the function $t\to \widetilde {\mathit{\boldsymbol{u}}}\left(t\right),$ we describe the space $C(I, K) = {C}^{0}(I, K)$ that is the space of continuous functions $\widetilde {u}:I\to K$ equipped with the norm as follows:

We will use the space ${L}^{2}(I, K)$ as the space of discontinuous functions which is given by

In time discretization, we split the intervals I into N subintervals ${I}_{\tau } = [{t}_{\tau -1}, {t}_{\tau }],$ where $\tau = 1, ···, N$ and $0 = {t}_{0} < {t}_{1} < {t}_{2} < ··· < {t}_{N-1} < {t}_{N} = T$. The parameter j indicates the time discretization parameter, as well as the maximum time step size $j = {}_{1\le \tau \le N}{}^{max}{j}_{\tau }$, where ${j}_{\tau } = {t}_{\tau }{-t}_{\tau -1},$ the length of nth time interval ${I}_{\tau }.$ The following set of time intervals ${M}_{j} = \{{I}_{1}, ···, {I}_{N}\}$ will be called the time-mesh. We find out the solution $\widetilde {\mathit{\boldsymbol{u}}}:I\to K$ on each time interval ${I}_{\tau }$ by a function ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}:I\to K$ which is a piecewise polynomial of some order l w.r.t time. The time-discrete solution space for ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}$ is ${X}_{j}^{l}\subset X$ and is defined by

where

The discrete test space for ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}$ is ${Y}_{j}^{l}\subset Y$ and is defined by

which is composed of $l–1$ piecewise polynomials (see ^{[52,53,54,55,56,57,59]} for details) and is discontinuous at the time step end nodes. We multiply the Eq (32) by the test function ${{v}_{j}\in Y}_{j}^{k}$ and integrate over the interval I. We get the discrete-time problem: Find $\widetilde {\mathit{\boldsymbol{u}}}\in {X}_{j}^{k}$ such that ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}\left(0\right) = 0$ and

where $\left\langle{\cdot , \cdot }\right\rangle$ represents the usual inner product in ${L}^{2}(I, K)$. This discretization is known as the exact continuous Galerkin-Petrov method or simply the "exact cGP(l)-scheme" of order k. The Galerkin-Petrov name is due to the fact that the solution space ${X}_{j}^{l}$ is different from the test space ${Y}_{j}^{l}$. The term "exact" denotes that the time integral on the right hand side of the Eq (35) is determined exactly. Because the discrete test space ${Y}_{j}^{l}$ is discontinuous, Eq (35) could be computed by a time marching technique in which local problems on the interval are handled successively. Therefore, we select the test function ${v}_{j}\left(t\right) = \nu \psi \left(t\right)$ with arbitrary time independent ν∈K and a scalar function $\psi :I\to \mathrm{R}$ which is zero on $I|{I}_{\tau }$ and a polynomial of order less than or equal to $l–1$ on the time interval ${I}_{\tau } = [{t}_{\tau -1}, {t}_{\tau }]$ Then, we get from (35) the ${I}_{\tau }$-problem of the exact cGP-scheme of order l: find ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}{ǀ}_{{I}_{\tau }}\in {\mathbb{p}}_{k-1}\left({I}_{\tau }, K\right)$ such that

with the initial condition ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}{ǀ}_{{I}_{\tau }}\left({t}_{\tau -1}\right) = {\widetilde {\mathit{\boldsymbol{u}}}}_{j}{ǀ}_{{I}_{\tau -1}}\left({t}_{\tau -1}\right)$ for $\tau \ge 2$ and ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}{ǀ}_{{I}_{\tau }}\left({t}_{\tau -1}\right) = {\widetilde {u}}_{0}$ for $\tau = 1$. In case of a nonlinear function $F\left\langle{\cdot , \cdot }\right\rangle$, we need to calculate the integrals numerically on the right hand side 6 of the Eq (36). The $(l+1)-\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}$ Gauß-Lobatto formula is exact if the function to be integrated has a polynomial of degree less than or equal to $2l-1$. As a result, this formula is applied to the integral on the left hand side of (36) will give the exact value. Then, the "${I}_{\tau }$–problem of the numerically integrated cGP(l) method" is: Find ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}{ǀ}_{{I}_{\tau -1}}\in {\mathbb{p}}_{l}\left({I}_{\tau }, K\right)$ such that ${\widetilde {\mathit{\boldsymbol{u}}}}_{j}\left({t}_{\tau -1}\right) = {\widetilde {\mathit{\boldsymbol{u}}}}_{\tau -1},$

where ${\widehat{w}}_{s}$ are the weights and $t\mathrm{ˆ}\in [-1, 1], s = 0, 1, 2, 3, . . . , l$ represent the nodes on the reference interval. To find ${\widetilde {u}}_{j}$ on each time interval ${I}_{\tau },$ we use a polynomial ansatz to illustrate it as follows:

where the coefficients ${U}_{\tau }^{s}$ are the components of K and the functions ${\phi }_{\tau , s}\in {\mathbb{p}}_{l}\left({I}_{\tau }\right)$ are the Lagrange basis functions (see ^{[52,53,54,55,56,57,59]} for more details) with respect to $l+1$ suitable nodal points ${t}_{\tau , s}\in {I}_{\tau }$ satisfying conditions mentioned below:

where ${\delta }_{\tau , s}$ is the Kronecker delta that is defined as:

Like in ^[57], the ${t}_{\tau , s}$ have been defined as the quadrature points of $(l+1)-\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}$ Gauß-Lobatto formula (for details information (see ^{[53,54,55,56,57,58,60]}) on the interval ${I}_{\tau }$. For the selection of initial conditions, we can set ${t}_{\tau , 0} = {t}_{\tau -1}$ which denotes the initial condition for Eq (36),

We define the basis functions ${\phi }_{\tau , s}\in {\mathbb{p}}_{l}\left({I}_{\tau }\right)$ via the affine reference transformation (see ^{[52,53,54,55,56,57,59]} for details explanations) $\overline {T}:\widehat{I}\to {I}_{\tau }$ where $\widehat{I} = [-1, \;1]$ and

Let ${\widehat{\phi }}_{s}\in {\mathbb{p}}_{l}\left(\widehat{I}\right)$, $\mathrm{s} = 0, \;\;\;\;1, \;.\;.\;.\;, l,$ demonstrate the basis functions that meet the requirements

where ${\widehat{t}}_{0} = -1$ and ${\widehat{t}}_{r}, r = 1, 2, . . . , l,$ are the quadrature points for the reference interval $\widehat{I}$. Then, we define the basis functions on the given interval ${I}_{\tau }$ by the mapping

Likewise, the test basis functions ${\varphi }_{\tau , r}$ are described by the suitable reference basis functions $\widehat{\varphi }\in {\mathbb{p}}_{l-1}\left(\widehat{I}\right),$ i.e.,

From the representation (37), we get for ${d}_{t}{\widetilde {\mathit{\boldsymbol{u}}}}_{j}$

By putting (40) in (36), we get

The integral is now transformed into the reference interval $\widehat{I}$ and computed using the $(l+1)-\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}$ Gauß-Lobatto quadrature formula which leads, for each test basis function $\phi \in {\mathbb{p}}_{l-1}$ and for all ν∈K,

where ${\widehat{w}}_{\mu }$ are the weights and ${\widehat{t}}_{\mu }\in [-1, 1]$ (see ^[60] for more details) are the points of integration with ${\widehat{t}}_{0} = -1$ and ${\widehat{t}}_{l} = 1.$ If we choose the test functions ${\varphi }_{\tau , i}\in {\mathbb{p}}_{l-1}\left({I}_{\tau }\right)$ such that

Now find the coefficients that are unknown ${U}_{\tau }^{s}\in K$ for s = 1, ..., l,

where ${U}_{\tau }^{s} = {U}_{\tau -1}^{s}\;for\;\tau > 1\;and\;{U}_{1}^{0} = {\widetilde {\mathit{\boldsymbol{u}}}}_{0}\;for\;\tau = 1 \;and$

We will discuss the cGP(k) method for the cases $l = 1\;and\;l = 2$ in the following.

3.9.1.   The cGP(1) method

We used the two point Gauß-Lobatto formula with ${t}_{\tau , 0} = {t}_{\tau -1}, {t}_{\tau , 1} = {t}_{\tau }$, and weights ${\widehat{w}}_{0} = {\widehat{w}}_{1} = 1$ which gives the well-known Trapezoidal rule. We obtain ${\alpha }_{\mathrm{1, 0}} = -1, {\alpha }_{\mathrm{1, 1}} = 1and{\beta }_{1} = 1.$ For the single coefficient ${U}_{\tau }^{1} = {\widetilde {u}}_{j}\left({t}_{\tau }\right)\in K$ the problem leads to the following block equation:

3.9.2.   The cGP(2) method

Three-point Gauß-Lobatto formula (Simpson rule) is used to define the quadratic basis functions with weights ${\widehat{w}}_{0} = {\widehat{w}}_{2} = 1/3, {\widehat{w}}_{1} = 4/3$ and ${\widehat{t}}_{0} = -1, {\widehat{t}}_{1} = 0, {\widehat{t}}_{2} = 1$ then, we get

Thus, the system to be solved for $\mathit{\boldsymbol{U}}_{\tau }^{1}$, $\mathit{\boldsymbol{U}}_{\tau }^{2}\in K$ from the known $\mathit{\boldsymbol{U}}_{\tau }^{0} = \mathit{\boldsymbol{U}}_{\tau -1}^{2}$ becomes:

where ${U}_{\tau }^{0}$ indicate the initial conditions at the current time interval.

4.   Comparison and computational analysis

Here, we apply the Galerkin and RK4 approaches and find out the numerical simulations based on the data shown in Table 1. One thousand people were identified as susceptible in the first hundred days of a community, with 170 getting infected and 790 being certified healthy. A total of 450 people died or recovered from the virus. Individuals who had been infected were placed in quarantine. We simulated the new developed model using such data and discovered that the virus will expand exponentially in the next few months, resulting in more people being quarantined if residents and government officials fail to follow SOPs. Therefore, as a consequence, the resistant individual falls, causing classes $I\left(t\right)$ and $Q\left(t\right)$ to expand. Susceptibility is reducing, as shown in Figure 3, and the healthy population is reducing, as shown in Figure 4. As a consequence, the infected and quarantined groups are rising (see Figures 5 and 6). Whereas the recovered class is strictly increasing, see Figure 7. Compared the outcomes of both methods and observed that the proposed schemes yield more accurate findings of the model. Figures 3–8 clearly demonstrate that the new techniques provide reasonably good agreement with the RK4-method results used for the COVID-19 model. This validate that the Galerkin and RK4 frameworks are capable of determining the behavior of variables in the region under investigation. The mesh graphs for the Galerkin-method and RK4-method solutions are displayed in Figures 8 and 9. Additionally, the computational cost of both method are carried out at various time intervals presented in Table 3 and displayed graphically in Figure 10. The simulations are carried with several time step sizes, and the time consuming is verified at various time points. The comparison of computational cost in terms of time for different step sizes indicate that the Galerkin technique is more time-consuming than the RK4-method for $t\in \left[\mathrm{0, 1}\right]$. Ultimately, we may deduce that when applied to similar problems, the numerical approach described in this article can be depended on to give relatively versatile and accurate solutions. The model, as mentioned earlier, is vital in the field of mathematical modeling of COVID-19 infection. This will be used to analyze the population dynamics of COVID-19 infected and susceptible population in the presence and absence of virus, helpful to observe the symptoms of COVID-19 seen clinically and valuable to hold back the disease. This study will be a valuable addition to the current literature on biomathematics.

5.   Conclusions and future recommendations

In the present article, we established a novel mathematical model that is composed of a system of differential equations that represents the population dynamics of susceptible, healthy, infected, quarantined, and recovered individuals. The next generation technique is used and examined the boundedness, local and global behaviour of equilibria, and the threshold quantity. The disease fee equilibria (DFE) are locally asymptotically stable when ${R}_{0} < 1$ and unstable when ${R}_{0} > 1$. The partial rank correlation coefficient (PRCC) technique is used to perform a sensitivity analysis on the basic reproduction number in order to identify the most significant parameter for regulating the threshold values of the model. To determine the conditions for stability analysis, the linearization and Lyapunov function theories are used. Additionally, we solved the model numerically using the well-known continuous Galerkin-Petrov time discretization scheme. In comparison, we solved the model using the conventional fourth-order Runge-Kutta approach to determine the validity and reliability of the aforementioned schemes. All the outcomes are displayed graphically. Furthermore, showed the computational cost in terms of time and concluded that the RK4 method is more efficient than the cGP (2) method.

Fractional analysis is a well-known fact that has recently become a popular research topic. Fractional calculus has been proven to be particularly useful in simulating a variety of real-world phenomena. Researchers have studied infectious diseases such as HIV, AIDS, and others by using fractional order derivatives and integrals. In the future, we have plane to implement the aforementioned numerical scheme to various fractional order differential equations.

Acknowledgments

This research received funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, (grant number B05F650018)

Conflict of interest

The authors declare no conflicts of interest.