 Research
 Open access
 Published:
Stability analysis and numerical evaluations of a COVID19 model with vaccination
BMC Medical Research Methodology volume 24, Article number: 97 (2024)
Abstract
A novel (nonlinear) mathematical model for the transmission of Coronavirus 19 (COVID19) with eight compartments and considering the impact of vaccination is examined in this manuscript. The qualitative behavior of the system such as the boundedness of solutions, the basic reproduction number, and the stability of the equilibrium points is investigated in detail. Some domestic real data collected from the Kerman University of Medical Science (KUMC) is used to estimate the parameters of the proposed model. We predict the dynamical behavior of the system through numerical simulations based on a combined spectral matrix collocation methodology. In this respect, we first linearize the nonlinear system of equations by the method of quasilinearization (QLM). Hence, the shifted version of Chebyshev polynomials of the second kind (SCPSK) is utilized along with the domainsplitting strategy to acquire the solutions of the system over a long time interval. The uniform convergence and upper bound estimation of the SCPSK bases are proved in a rigorous manner. Moreover, the technique of residual error functions is used to testify the accuracy of the QLMSCPSK method. The presented numerical results justify the robustness and good accuracy of the QLMSCPSK technique. The achieved numerical orders of convergence indicate that the QLMSCSK algorithm has exponential rate of convergence. Using the linearization technique in one hand and the domainsplitting strategy on the other hand, enable us to predict the behaviour of similar disease problems with high accuracy and maximum efficiency on an arbitrary domain of interest.
Introduction
Infectious diseases have always been a limitation of human life. Many epidemics throughout history have killed millions of people. The last of them was the COVID19 virus, the first case of which was observed in China in December 2019. The disease quickly spread throughout the world within a few months. In the mid of Mach 2020, WHO announced that the infection was an outbreak. A prominent feature of the current pandemic is the high persontoperson transmissibility of the virus, with a basic reproduction number \((R_0)\) estimated at 2.22.5 in Wuhan, China. This caused an increase in hospitalizations and as a result deaths due to the disease [1, 2]. Researchers in the field of epidemiology and other branches of biology tried to find treatments or vaccines for this disease. The first vaccines were clinically tested at the beginning of 2021. In the following months, a number of highly effective vaccines entered the market. Mathematical models can describe the process of transmission of infection and its control using vaccination.
Researchers try to understand the dynamics of a disease in the first stage. Hence, they develop control and curing procedures for the diseases. The first modern mathematical models were introduced by Kermak and Mckendrick in 1927 [3]. Since 2019, many mathematical models have been presented. Different aspects of the novel COVID19 have been investigated through mathematical models [4,5,6]. Before vaccine production, the effect of quarantine and control factors were investigated [7, 8].
Butt et al. developed and analyzed mathematically the extended SEQIHR model. The authors determined the possible control strategies and comprehended the longterm dynamics of disease [9]. A SEQIER model introduced to examine the transmission dynamics of COVID19. The focus is on the effectiveness of hospitalization and quarantine strategies [10, 11]. A deterministic SEIHR fractional model is developed in [12]. The authors performed the mathematical analysis and the design of an optimal control strategy for the proposed CaputoFabrizio fractional model. In [13], the authors analyzed a fractional order initial value problem with the AtanganaBaleanu derivative operator to observe transmission dynamics of the infection in a human population. An SVEIR mathematical model was introduced to predict transmission dynamics of COVID19 infection [14]. In [15], the effect of virus rate in the environment in a deterministic model is investigated. Some of the mathematical models of COVID19 are presented by fractionalorder models [16].
Vaccination is very effective to control infectious diseases. Therefore, after producing vaccines for coronavirus, researchers developed dynamical models regarding the effectiveness of vaccines. The authors in [17] developed vaccination strategies for vaccination and investigated the dynamic of their epidemic model. A dynamical model considering the treatment and vaccination saturating function is introduced in [18]. The study of the vaccination model under vaccine immunization has been considered in [19]. The concerned investigations are devoted to stability theory, numerical simulation, and globallocal dynamics [20]. One of the items that is investigated in the dynamic analysis of models is bifurcation. Researchers have investigated the existence of bifurcation in different models [21,22,23].
In this paper, a (novel) mathematical model for COVID19 with an omicron version is presented. In this model, vaccination and its effect on people is considered. In this study, the importance of vaccination to control the epidemic has been investigated. In this model, vaccinated people and those who have not received the vaccine are divided into two groups. Therefore, the rate of infection and hospitalization rate in each group is investigated separately. In this deterministic model, the population is divided into two main parts, people who have been vaccinated and people who have not received any vaccine. It is also assumed that recovered individuals are prepared to reinfection after some time. As follows, some basic dynamic properties of the proposed model are investigated. Using the statistics of the Ministry of Health of Iran in Kerman province, the parameters are estimated over a period of time.
Furthermore, we develop a combined spectral collocation approach based on the shifted version of Chebyshev polynomials of the second kind (SCPSK) to predict the behaviour of COVID19 disease. To get efficacy, we first transform the nonlinear system of ODEs into a set of linearized system of eight equations to be treated iteratively. Hence, the employed SCPSK basis functions together with the domain decomposition strategy are utilized to find the solutions of the linearized systems. We also prove the convergence of the SCPSK bases and an upper bound estimation for them is performed. It should be noted that spectralbased collocation techniques have been used often to tackle various integral and differential equations due to their simplicity of implementations and highorder accuracy. These spectral methods have been benefited by utilizing numerous basis functions like the VietaLucas [24], Chebyschev [25,26,27], Bessel [28, 29], Chelyshkov [30,31,32], and the Schröder [33] polynomials.
Let us illustrate the main achievements of the present work in a concise form as follow:


A novel mathematical model consists of eight equations for the studying of COVID19 transmission is proposed in which the effect of vaccination is considered.


The parameters of the model are estimated from the real data provided by the KUMC.


The dynamic analysis of the proposed COVID19 model is performed from theoretical points of stability, boundedness, and the existence of bifurcations.


A combined efficient method called the QLMSCPSK algorithm based on the technique of quasilinearization along with the Chebyshev spectral collocation approach is designed to solve the given model numerically and confirm the theoretical findings as well.


The strategy of domainsplitting is further employed to keep the accuracy of the proposed QLMSCPSK algorithm at a desired level. The exponential convergence of the employed SCPSK bases in the infinity norm is established in a rigourous analysis.
The structure of the current work is provided as follows. A mathematical model considering vaccination is introduced and some basic dynamical properties such as the boundedness of solutions, the derivation of the basic reproduction number, the stability of the equilibrium points, and the existence of bifurcation are investigated in the subsequent “Mathematical model” section. In the next “The Chebyshev polynomials of the second kind: shifted version on [t_{a}, t_{b}]” section we first give a review of the SCPSK functions. Also, a rigorous mathematical proof is given for the convergence analysis of the SCPSK bases. “QLMSCPSK collocation strategy based on splitting of time interval” section first describes the fundamental ideas of the QLM. Hence, the methodology of the spectral domain decomposition approach is illustrated and finally, the basic steps of the proposed QLMSCPSK algorithm are given in detail. In “Experimental calculations” section some numerical simulations on a longtime domain are carried out to support the theoretical findings. The last “Conclusions” section includes a summary of the presented research study.
Mathematical model
A dynamical model that is introduced in this paper is based on the Kermak and Mckendrick model [3]. The model is introduced for the omicron version of Coronavirus. During the outbreak of the omicron, a large number of people in the community had received the vaccine. Therefore, the population of susceptible, infected, recovered, and isolated are divided into two groups, vaccinated and unvaccinated. Consider, S(t), \(S_v(t)\) be the populations of unvaccinated and vaccinated susceptible at time t respectively. Denote further by I(t) and \(I_v(t)\) as the number of unvaccinated and vaccinated infected people. Also, R(t) and \(R_v(t)\) are unvaccinated recovered and vaccinated recovered persons. Finally, J(t) and \(J_v(t)\) are the number of unvaccinated and vaccinated isolated people.
We consider The fraction \(\beta\) of susceptible people who get infected. This fraction of unvaccinated is \(\beta '\). Also, a fraction \(\delta\) of vaccinated susceptible goes to \(R_v\) class without infection. The recovery rates are \(\gamma _1\) and \(\gamma _2\) for unvaccinated and vaccinated infected people. The unvaccinated and vaccinated infected people are isolated with rates \(\alpha _1\) and \(\alpha _2\) respectively. The recovery rate in J and \(J_v\) groups are \(\eta _1\) and \(\eta _2\). Two parameters \(\mu _1\) and \(\mu _2\) are mortality rates because of disease in unvaccinated and vaccinated isolated populations. In the end, unvaccinated and vaccinated recovered people are reinfected with the rate of \(\theta _1\) and \(\theta _2\). In this model, \(\Lambda\) is the birth rate of the population. The dynamical model is given by
The flowchart of the model is given in Fig. 1.
The above model is subjected to the following initial conditions
Let us derive the equilibrium points of the given model Eq. (1). The first equilibrium point is the diseasefree one. It is given by
The other equilibrium point is obtained as follows
where
Also, we have defined two parameters \(x_1\) and \(x_2\) as
By doing some calculations we find that \(I^*\) satisfies in the following algebraic equation
where the coefficients \(a_i\), for \(i=0,1,2\) are given by
Also, we have
If \(y_1>\beta \Lambda\), the Eq. 5 has a positive root. Hence, the endemic equilibrium point \(E^*\) exists.
Boundedness of system
Our goal here is to examine the boundedness property of the COVID19 system. To proceed, we have
Theorem 2.1
All solutions of the COVID19 system Eq. (1) are uniformly bounded.
Proof
Consider the new variable
Thus, we have
Since, \(\mu < \mu _i\) for \(i=1, 2\), we arrive at
By solving the above linear differential equation in terms of variable X one gets,
Now, by tending \(t\rightarrow \infty\), we obtain \(0\le X\le \frac{\Lambda }{\mu }\). Therefore, all solutions of the system, i.e., \((S, S_v, I, I_v, R, R_v, J, J_v)\) of the system are confined the region
for any \(\epsilon >0\) as \(t\rightarrow \infty\). \(\square\)
Basic reproduction number
Let us present the basic reproduction number by \(R_0\). It is the most important epidemiological parameter. It is defined as the “number of secondary infected individuals caused by a single infected individual in the whole time interval” [34]. In this paper, using the method presented in the paper [34], we calculate \(R_0\) value. The compartments that include infected people are \(I, I_v, J, J_v\). Therefore, we form the following system
The above system can be rewritten as
where, \(y=(I, I_v, J, J_v)^T\), and
Now, let us obtain the Jacobian matrix related to \(\varphi\) and \(\psi\) at \(E_0\) as the disease free equilibrium point. We have
and
As we know that the basic number, \(R_0\), is the eigenvalue of the matrix \(FV^{1}\). It is given by some simple calculations as
Theorem 2.2
Under the restriction \(R_0<1\), the point \(E_0\) (disease free equilibrium point) is asymptotically stable.
Proof
Denote the Jacobian matrix of COVID19 system Eq. (1) by \(\varvec{J}:=\varvec{J}(S,S_v,I,I_v)\). It has the form
where \(I_1:=I+I_v\) and \(\zeta _j:=\gamma _j+\alpha _j+\mu\) for \(j=1,2\). The former Jacobian matrix at the equilibrium point \(E_0\) can be written as
The characteristic equation is \(p(x):=\det (\varvec{J}(E_0)x\varvec{I})=0\). By computing the determinant and after some manipulations we get
where we have
Owing to the fact that we have \(R_0<1\), one can easily check that all roots are negative. Hence, the equilibrium point \(E_0\) is asymptotically stable. \(\square\)
Bifurcation
In this section, we describe the existence of bifurcation at threshold number \(R_0=1\). We select \(\beta\) to be the bifurcation coefficient that yields bifurcation at \(R_0=1\). The model takes the form
Here, the transmission rate, \(\beta '\), is less than \(\beta\). According to the available information, consider \(\beta '= \frac{\beta }{30}\). Hence, \(R_0=1\) gives
Next, we construct characteristic polynomial at \(E_0\) and bifurcation parameter \(\beta ^*\) as follows
where \(q_1(x)\) and \(q_2(x)\) are defined in Eq. (8). We can easily verify that all roots are negative except one of them is zero. Therefore, there exists bifurcation at \(R_0=1\).
Now, we compute the right eigenvector \(w=(w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8)^T\). Corresponding to the zero eigenvalue that satisfies \(J(E_0)w=0\).
By solving the equation, we obtain
To compute the left eigenvector \(V=(v_1, v_2, v_3, v_4, v_5, v_6, v_7, v_8)\), we solve the equation \(V\cdot J(E_0)=\textbf{0}\). We obtain \(v_1=0, v_2=0, v_4=\frac{\xi _1 \beta S^0}{\beta ' S_v^0}v_3,v_5= v_6= v_7= v_8=0.\)
Using Theorem 4.1 in [35], we compute the bifurcation coefficients a and b as follow. For the first one, we have
where, \(x_1= \frac{\xi _1 \beta S^0}{\beta S_v^0}\). One can be easily seen that \(a>0\).
Similarly, for the second parameter b, we have
Thus, if \(\beta =\beta ^*\), the diseasefree equilibrium point \(E_0\) is unstable. If \(R_0<1\), \(E_0\) is locally asymptotically stable. Otherwise, \(E_0\) is unstable.
The Chebyshev polynomials of the second kind: shifted version on [t _{a}, t _{b}]
In order to devise our multidomain collocation scheme, we need to introduce the set of (orthogonal) polynomials known as the Chebyshev polynomials of the second kind (CPSK), see [36]. These polynomials on \([1,1]\) are defined by the formula
where \(\Theta \in [0,\pi ]\) and \(j\ge 0\). The set of CPSK polynomials constitutes an orthogonal system on \([1,1]\) with regard to positive function \(w_1(p)=(1p^2)^{\frac{1}{2}}\). The shifted version of them are obtained by using \(p=(2z1)\) on [0, 1]. They have the explicit representation as [37]
with \(U_0(z)=1\). The second one is obtained as \(U_1(z)=4z2\). In addition to Eq. (10), we may use the next recursive formula to get the remaining polynomials as
Moreover, the orthogonality condition of these shifted polynomials with regard to weight function \(w_2(z)=(zz^2)^{\frac{1}{2}}\) is stated as
We finally notice that the roots of \(U_j(z)\) are located inside the interval (0, 1). Precisely speaking, these zeros are [37]
We will utilize the shifted CPSK on an arbitrary interval \([t_a,t_b]\). Toward this end, we take \(z=(tt_a)/T\) or \(t=T z+t_a\), where \(T=t_bt_a\). So, we have
Definition 3.1
The shifted CPSK (SCPSK) on \([t_a,t_b]\) is defined via relation
From the above change of variable, we are able to represent the summation Eq. (10) as
where
With the preceding transformation, the orthogonality condition Eq. (12) is changed accordingly. The orthogonality of the set \(\{\mathbb {U}_j(t)\}_{j=0}^{\infty }\) will be gotten with regard to the weight function \(w_a(t)=\sqrt{T(tt_a)(tt_a)^2}\). We have indeed the next formula
Finally, based on relations \(t=T z+t_a\) and Eq. (13) we can locate the zeros of SCPSK \(\mathbb {U}_j(t)\) on \((t_a,t_b)\). The associated zeros are
These finite number of points will be utilized as the collocation nodes in the main proposed algorithm, below.
Function approximation: convergence and error analysis of SCPSK
To continue, let us suppose that a given function x(t) is belonged to the weighted \(L^2\) space on \([t_a,t_b]\). It then can be expressed as a summation of SCPSK in the form
To obtain the unknown coefficients \(\nu _j,j\ge 0\), the aim is to exploit the orthogonality property Eq. (16). It gives us
We first find an upper bound for the coefficients \(\nu _j\) in terms of j. Hence, we show that series form Eq. (18) is convergent uniformly on \([t_a,t_b]\). To continue. let assume that \(M_{ab}:=\max _{t\in [t_a,t_b]}\left x''(t)\right\). Moreover, the space \(L^2_w\) stands for the space of squareintegrable functions on \([t_a,t_b]\) with regard to weight function \(w_a(t)\) defined previously. Thus, we have the next assertion.
Theorem 3.1
For the function \(x(t)\in C^{(2)}([t_a,t_b])\cap L^2_w([t_a,t_b])\), which can be written as Eq. (18) we obtain the following estimate
Proof
We begin by the relation Eq. (19) and make the change of variable \(tt_a=\frac{T}{2}(1+\cos \Theta )=:r(\Theta )\) we arrive at
We then integrate by parts two times on the last formula. By introducing the auxiliary function
we get from Eq. (21) the next expression for \(\nu _j\) as
We next utilize the fact that \(\left x''\left( r(\Theta )\right) \right \le M_{ab}\) and \(\sin (\Theta )\le 1\). The resulting inequality is
A simple calculation can be done to obtain the exact values of the integral on the former line. It follows that
By utilizing an odd \(j>1\), the resulting integral’s value is zero obviously. If we choose a \(j>1\) to be even, then we have
We use the relation \(j1\ge \frac{j}{2}\), which holds for all \(j\ge 2\). Hence, we reach at the inequality
By placing Eq. (24) into Eq. (23) yields the result Eq. (20).
To handle the model problem Eq. (1) in practical situations, we require to truncate the infinite series solution Eq. (18) into a finite one. If we use only \((\mathcal {J}+1)\) (\(\mathcal {J}\in \mathbb {N}\)) basis functions, we can approximate x(t) in the way that
So, our goal now is to examine the difference between x(t) and \(x_{\mathcal {J}}(t)\) as an approximation for it. This global error is defined by the formula
To get an upper bound for the global error \(\mathbb {E}_{\mathcal {J}}(t)\), we employ the following relation is valid for the normal CPSK on \([1,1]\) given as
With the aid of transformation \(p=2(tt_a)/T1\), we arrive at the same conclusion for \(\mathbb {U}_j(t)\) on \([t_a,t_b]\). Namely, we have
The next result establishes that the norm of error \(\Vert \mathbb {E}_{\mathcal {J}}\Vert _{\infty }\) converges to 0 if we let \(\mathcal {J}\) goes to infinity.
Theorem 3.2
Under assumptions of Theorem 3.1, the error \(\mathbb {E}_{\mathcal {J}}(t)\) converges to 0 as \(\mathcal {J}\rightarrow \infty\). In fact, we have
Proof
To prove this result, we consider Eq. (26) together with Eq. (27) to obtain
Hence, we apply the upper bound for \(\left \nu _j\right\) derived in Eq. (20). The resultant inequality is
Utilizing the wellknown Integral Test [38] gives us
Our proof is established by inserting the foregoing result into Eq. (28) followed by taking the supremum over all \(t\in [t_a,t_b]\). By tending \(\mathcal {J}\) to infinity, we have done the proof.
QLMSCPSK collocation strategy based on splitting of time interval
Let us emphasize that the spectral matrix collocation approach based on the SCPSK may not yield convergence on a long time interval \([t_a,t_b]\). One remedy is to use a large number of bases on the long domains accordingly to reach the desired level of accuracy. Another approach is to divide the given interval into a sequence of subintervals and employ the proposed collocation scheme on each subinterval consequently.
Towards this end, we split the time interval \([t_a,t_b]\) into \(N\ge 1\) subdomains in the forms
Here, we have \(t_0:=t_a\) and \(t_N:=t_b\). The uniform time step is taken as \(h=t_{n+1}t_n=(t_bt_a)/N\). Note that by selecting \(N=1\), we turn back to the traditional spectral collocation method on the whole domain \([t_a,t_b]\). Therefore, on each subinterval \(K_n\) we take the approximate solution of the model Eq. (1) to be in the form Eq. (25) as
where we utilized the notations
as the vector of unknown coefficients and the vector of SCPSK bases respectively. Once we get the all local approximate solutions for \(n=0,1,\ldots ,N1\), the global approximate solution on the given (large) interval \([t_a,t_b]\) will be constructed in the form
In order to collocate a set of \((\mathcal {J}+1)\) linear equations to be obtained later at some suitable points, we consider the roots of \(\mathbb {U}_{\mathcal {J}+1}(t)\) on the subinterval \(K_n\). By modifying the points given in Eq. (17), we take the collocation nodes as
At the end, we note that in the proposed splitting approach, the given initial conditions of the underlying model problem are prescribed on the first subinterval \(K_0\). Once the approximate solution on \(K_0=[t_0,t_1]\) is determined, we utilize it to assign the initial conditions on the next time interval \(K_1\). To do so, it is sufficient to evaluate the obtained approximation at \(t_1\). We repeat this idea on the next subintervals in order until we arrive at the last subinterval \(K_{N1}\). Below, we illustrate the main steps of our matrix collocation algorithm on an arbitrary subinterval \(K_n\) for \(n=0,1,\ldots ,N1\).
The QLMSCPSK matrix collocation technique
Our chief aim is to solve the nonlinear COVID19 system Eq. (1) efficiently by using the spectral method based on SCPSK basis. Towards this end, we first need to get rid of the nonlinearity of the model. This can be done by employing the Bellman’s quasilinearization method (QLM) [39]. Thus we will get more advantages in terms of running time, especially for large values of J in comparison to the performance of directly applied collocation methods to nonlinear models, see cf. [40,41,42]. By combining the idea of QLM and the splitting of the domain we will obtain more gains in terms of accuracy for the approximate solutions of nonlinear model Eq. (1). Let us first describe the technique of QLM. For more information, we may refer the readers to the abovementioned works.
By reformulating the original COVID19 model Eq. (1) in a compact form we get
where
To begin the QLM process, we assume \(\varvec{z}_0(t)\) is available as an initial rough approximation for the solution \(\varvec{z}(t)\) of the COVID19 system Eq. (31). Through an iterative manner, the QLM procedure reads as follows
Here, the notation \(\varvec{G}_{\varvec{z}}\) stands for the Jacobian matrix of the COVID19 system Eq. (31), which is of size 8 by 8. By performing some calculations we reach the linearized equivalent model form as
where \(\varvec{M}_{s1}(t):=\varvec{J}(S_{s1}(t), (S_v)_{s1}(t), I_{s1}(t), (I_v)_{s1}(t))\) as the Jacobian matrix \(\varvec{J}\) previously constructed in Eq. (7). Also we have
Along with the system Eq. (32) the initial conditions
are given due to Eq. (2). We now are able to solve the family of linearized initialvalue problems Eqs. (32)(33) numerically by our proposed matrix collocation method on an arbitrary (long) domain \([t_a,t_b]\). For this purpose and for clarity of exposition, we restrict our illustrations to a local subinterval \(K_n\) for \(n=0,1,\ldots ,N1\).
In view of Eq. (29) by utilizing only (\(\mathcal {J}+1\)) SCPSK basis functions, we assume that the eight solutions of system Eq. (32) can be represented in terms of Eq. (29). Thus, we take these solutions at iteration \(s\ge 1\) as
for \(t\in K_n\). Moreover, by \(\varvec{W}^{n,s}_{\mathcal {J},i}= \left[ \begin{array}{cccc} \omega ^{n,s}_{0,i}&\omega ^{n,s}_{1,i}&\dots&\omega ^{n,s}_{\mathcal {J},i} \end{array}\right] ^T\) we denote the vectors of unknowns for \(1\le i\le 8\) at the iteration \(s\ge 1\). Also, the vector of SCPSK basis, i.e., \(\varvec{U}_\mathcal {J}(t)\) is defined in Eq. (29). We next provide a decomposition for \(\varvec{U}_\mathcal {J}(t)\) given by
Here, the vector \(\varvec{Q}_\mathcal {J}(t)\) including the powers of \((tt_n)\) introduced by
The next object is the matrix \(\varvec{F}_\mathcal {J}=(f_{i,j})_{i,j=0}^{\mathcal {J}}\) of size \((\mathcal {J}+1)\times (\mathcal {J}+1)\). The entries of the latter matrix are given in Eq. (15). One can also show that \(\det (\varvec{F}_\mathcal {J})\ne 0\) and it is a triangular matrix. It follows that
We then insert the obtained term \(\varvec{U}_\mathcal {J}(t)\) in Eq. (35) into Eq. (34). The resulting expansions are
We then proceed by nothing that the derivative of the vector \(\varvec{Q}_\mathcal {J}(t)\) can be stated in terms of itself. A vivid calculation reveals that
From this relation, we are able to derive a matrix forms of the derivatives of the unknown solutions in Eq. (36).
The exact solutions of the linearized system Eq. (32) can be written in a vectorized form as
We next introduce the following block diagonal matrices of dimensions \(8(\mathcal {J}+1)\times 8(\mathcal {J}+1)\) as
By the aid of the former definitions, the matrix formats of \(\varvec{z}^n_{\mathcal {J},s}(t)\) and \(\dot{\varvec{z}}^n_{\mathcal {J},s}(t)\) will rewrite concisely as
Here, \(\varvec{W}^n\) is the successive vector of eight previously defined vector of unknowns
We now can collocate the linearized Eq. (32) at the zeros of SCPSK given in Eq. (17) on the subdomain \(K_n\). We get
for \(s=1,2,\ldots\). Denote the coefficient matrix by \(\widehat{\varvec{M}}^n_{s1}\) and the righthandside vector as \(\widehat{\varvec{R}}^n_{s1}\). These are defined by
Let us define further the vectors of unknowns as
Consequently, the system of Eq. (41) can be stated briefly as
and with \(s=1,2,\ldots\). Before we talk about the fundamental matrix equation, we need to state two vectors \(\varvec{Z}^n_s\) and \(\dot{\varvec{Z}}^{n}_{s}\) in Eq. (42) in the matrix representation forms. The proof is easy by just considering the definitions of the involved matrices and vectors in Eq. (40).
Lemma 4.1
If two vectors \(\varvec{z}^n_{\mathcal {J},s}(t)\) and \(\dot{\varvec{z}}^n_{\mathcal {J},s}(t)\) in Eq. (40) computed at the collocation points Eq. (30), we arrive at the next matrix forms
where the matrix \(\bar{\widehat{\varvec{Q}}}\) is given by
Moreover, two matrices \(\widehat{\varvec{Q}}, \widehat{\varvec{F}}\) are defined in Eq. (40). Similarly, the vector \(\varvec{W}^n\) is given in Eq. (40).
By turning to relation Eq. (40) we substitute the derived matrix formats into it. Precisely speaking, after replacing \(\varvec{Z}^n_s\) and \(\dot{\varvec{Z}}^n_s\) we gain the socalled fundamental matrix equation (FME) of the form
where
To complete the process of QLMSCPSK approach, it is necessary to implement the initial conditions in Eq. (2) and add them into Eq. (44). So, the next task is to constitute the matrix representation of Eq. (2). Let us approach \(t\rightarrow 0\) in the first relation of Eq. (40). It gives us
We then replace eight rows of the augmented matrix \([\varvec{B}_n;\widehat{\varvec{R}}^n_{s1}]\) by the already obtained row matrix \([\varvec{B}_{0,n};\widehat{\varvec{R}}^n_{s1,0}]\). Denote the modified FME by
This implies that the solution of the model Eq. (1) is obtainable on each subdomain \(K_n\) by iterating \(n=0,1,\ldots ,N1\). On \(K_0\) as the first subdomain, the given initial conditions in Eq. (2) will be used to find the corresponding approximations for the system Eq. (1). Hence, this approximate solutions on \(K_0\) evaluated at the starting point of \(K_1\) will be utilized for the initial conditions on \(K_1\). By repeating this process we acquire all approximations on all \(K_n\) for \(0\le n\le N1\).
REFs and QLMSCPSK technique
Generally, finding the true solutions of the COVID19 system Eq. (1) is not possible practically. In this case, the residual error functions (REFs) help us to measure the quality of approximations obtained by the QLMSCPSK technique. Once we calculate the eight approximations by the illustrated method, we substitute them into the model system Eq. (1). In fact, the REFs are defined as the difference between the lefthand side and the righthand side of the considered equation. On the subdomain \(K_n\) we set the REFs as
for a fixed iteration number s and we have defined \(L^n_{\mathcal {J},s}:=I^{n}_{\mathcal {J},s}(t)+ (I_v)^{n}_{\mathcal {J},s}(t)\) for brevity.
Analogously, at the fixed iteration s, the numerical order of convergence associated with the obtained REFs can be defined in the infinity norm. These are given by
Therefore, the convergence order (Co) for each solution is defined by
Experimental calculations
We now exploit the proposed QLMSCPSK collocation technique to solve the COVID19 system Eq. (1) numerically. We use Matlab version 2021a installed on a personnel laptop to run our algorithm. Different values of the model parameters will be utilized in the experimental results. Setting \(s=5\) is sufficient in the iterated QLM to reach the desired accuracy. For each run, to begin computations we use \(\varvec{z}_0(t)\equiv 0\) as the first rough approximation for the linearized system Eq. (32).
To solve the underlying COVID19 model, we need to determine the involved parameters in the model. Towards this end, the estimations of these parameters are done in accordance with the statistics given by the Kerman University of Medical Sciences (KUMS) during the period between 22 December 2021 and 5 May 2022. The estimated parameters are listed in Table 1.
The initial conditions are taken as
To start computations, we first set \(J=10\) and the computational domain is taken as \([t_a,t_b]=[0,1000]\) in terms of days spent during the pandemic. The number of subintervals used is \(N=100\). Through Figs. 2, 3, 4 and 5 we show the behaviors of the approximate solutions obtained via QLMSCPSK procedure for susceptible, infected, recovered, and isolated populations in both unvaccinated and vaccinated parts. These curves are presented by dashed red lines. Besides, the outcomes of the wellknown Matlab function ode45 are also depicted to validate our approximations. These plots are shown by solid black lines. As one sees the plots obtained by both methods are matched very well. More accuracy is attainable by just increasing the number of bases, namely J, or using a larger number of subintervals N in the computations. Note that for each solution the value of the corresponding population at the endpoint \(t=1000\) is given for completeness.
As you see, the population of susceptible individuals will decrease in both vaccinated and unvaccinated groups. But with consideration, so a large fraction of the population receives the vaccine, a large number of people directly enter the group of recovered people \((R_v)\). Hence, the number of infected and isolated individuals is not high compared to the population and a significant decrease is observed in both vaccinated and unvaccinated groups. According to the prediction of these graphs, with the passage of time, the population goes through an epidemic, which is consistent with the reality that happened in the province of Kerman.
We next show the effect of utilizing a diverse number of bases, \(J=10,20\), on the computed solutions. In this respect, we compute the REFs formulae in Eq. (46) related to the approximate solutions of the COVID19 model Eq. (1). To save space, we only visualize the REFs corresponding to the susceptible and recovered populations. These results for both vaccinated and unvaccinated counterparts are displayed in Figs. 6 and 7. Not that we used the fixed number of subintervals as \(N=100\). It can clearly seen that the desirable form of accuracy is attainable if one increases the number of bases accordingly.
As we mentioned, we can also get the smaller magnitude of REFs by increasing the number of subintervals N in the proposed domain decomposition spectral QLMSCPSK approach. Towards this end, let use take \(N=50,100\) and \(N=200\). For a fixed \(J=15\), we show the results of REFs associated with the recovered populations (vaccinated and unvaccinated) in Fig. 8. From the plots shown in Fig. 8, it can be obviously concluded that highorder accuracy is attained for smaller lengths of intervals.
Finally, we compute \(L^{\infty }_\ell\) error norms and the related numerical order of accuracy as shown by \(\textrm{Co}_{\mathcal {J}}^\ell\) for \(\ell =5,6,7,8\), namely for recovered and isolated populations. The results by using diverse values of \(\mathcal {J}=2^k,\) for \(k=1,2,\ldots ,5\) are listed in Table 2. For this purpose, we use \(N=10\) and the computational interval is set as [0, 100]. By looking at the results tabulated in this table one can infer that the proposed spectral QLMSCPSK approach has an exponential order of accuracy.
Parameter study
In this part, the relationship between \(\beta\), \(\gamma\) and \(R_0\) is examined. As you can see in the next Fig. 9 (left picture), as the \(\gamma\) value increases, \(R_0\) decreases. The lowest value of \(R_0\) is at \(\gamma =0.3\). Therefore, if the recovery rate reaches 0.3, the epidemic will disappear. The relationship between \(\beta\) and \(R_0\) is shown in the same figure but on the right panel. As you see, there is a direct relationship between them.
Conclusions
A novel mathematical model for studying the COVID19 pandemic disease has been suggested by dividing the population into vaccine and nonvaccine groups. From a dynamical point of view, the boundedness of the system was proved and the basic reproduction number was obtained to prove the stability of the equilibrium points. Moreover, the existence of bifurcation for the COVID19 system was discussed. To get the solutions of this system, an efficient spectral method based on the shifted Chebyshev polynomials of the second kind (SCPSK) combined with the quasilinearization methodology (QLM) was used. The methodology of domainsplitting further was implemented to keep the accuracy of the proposed QLMSCPSK approach at a desired level on a long time interval. The convergence analysis of the SCPSK basis functions in the \(L_{\infty }\) norm was performed and an upper bound estimation for the error was done. The involved parameters were estimated by the real data provided by the Kerman University of Medical Sciences during the period between 22 December 2021 and 5 May 2022. Numerous computational experiments based on the proposed QLMSCPSK technique were conducted to predict the behavior of disease over a long time interval. The presented QLMSCPSK technique with high accuracy and efficacy can be generalized to solve similar epidemiological models with integerorder and fractionalorder derivatives and even with more than eight equations.
Availability of data and materials
Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.
References
Li Q, Guan X, Wu P, Wang X, Zhou L, Tong Y, et al. Early transmission dynamics in Wuhan, China, of novel coronavirusinfected pneumonia. N Engl J Med. 2020;382:1199–207.
Wu JT, Leung K, Leung GM. Nowcasting and forecasting the potential domestic and international spread of the 2019nCoV outbreak originating in Wuhan, China: a modelling study. Lancet. 2020;395(10225):689–97.
Kermack WO, Mckendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc Lond A. 1927;115(722):700–21.
Waezizadeh T, Ebrahimi N. Dynamical model for COVID19 in a population. J Mahani Math Res. 2022;11(1):25–34.
Naik PA, Owolabi KM, Zu J, Niak M. Modeling the transmission dynamics of COVID19 pandemic in Caputotype fractional derivative. J Multiscale Model. 2021;12(03):21500062.
Khan MA, Atangana A. Mathematical modeling and analysis of COVID19: A study of new variant Omicron. Phys A. 2022;599:127452.
Khajanchi S, Sakar K, Banerjee S. Modeling the dynamics of COVID19 pandemic with implementation of intervention strategies. Eur Phys J Plus. 2022;137:1–22.
Ullah MS, Higazy M, Kabir KA. Modeling the epidemic control measures in overcoming COVID19 outbreaks: A fractionalorder derivative approach. Chaos Solit Fract. 2022;155:1–21.
Butt AIK, Ahmad W, Rafiq M, Ahmad N, Imran M. Computationally efficient optimal control analysis for the mathematical model of Coronavirus pandemic. Expert Syst Appl. 2023;234:121094.
Butt AIK, Ahmad W. Optimally analyzed fractional Coronavirus model with AtanganaBaleanu derivative. Results Phys. 2023;53:106929.
Butt AIK, Rafiq M, Ahmad W, Ahmad N. Implementation of computationally efficient numerical approach to analyze a Covid19 pandemic model. Alex Eng J. 2023;69:341–62.
Deressa CT, Duressa GF. Investigation of the dynamics of COVID19 with SEIHR nonsingular and nonlocal kernel fractional model. Int J Model Simul. 2022;42(6):1030–48.
Hanif A, Butt AIK, Ahmad W. Numerical approach to solve CaputoFabriziofractional model of corona pandemic with optimal control design and analysis. Math Methods Appl Sci. 2023;46(8):9751–82.
Ahmad W, Abbas M, Rafiq M, Baleanu D. Numerical analysis of AtanganaBaleanu fractional model to understand the propagation of a novel corona virus pandemic. Alex Eng J. 2022;61(9):7007–27.
Naik PA, Zu J, Ghori MB, Naik M. Modeling the effects of the contaminated environments on COVID19 transmission in India. Results Phys. 2021;29:104774.
Naik PA, Yavuz M, Qureshi S, Zu J, Townley S. Modeling and analysis of COVID19 epidemics with treatment in fractional derivatives using real data from Pakistan. Eur Phys J Plus. 2020;795:1–42.
Zheng Q, Wang X, Bao C, Ji Y, Liu H, Meng Q, Pan Q. A multiregional, hierarchicaltier mathematical model of the spread and control of COVID19 epidemics from epicentre to adjacent regions. Transbound Emerg Dis. 2022;69:549–58.
Deng J, Tang S, Shu H. Joint impacts of media, vaccination and treatment on an epidemic Filippov model with application to COVID19. J Theor Biol. 2021;523:1–15.
Han X, Liu H, Lin X, Wei Y, Ming M. Dynamic analysis of a VSEIR model with vaccination efficacy and immune decline. Adv Math Phys. 2022;2022:1–16.
Shaikh AS, Shaikh IN, Nisar KS. A mathematical model of COVID19 using fractional derivative: outbreak in India with dynamics of transmission and control. Adv Differ Equ. 2020;2020(1):373.
Xu C, Liu Z, Pan Y, Cui X, Ou W, Rahman MU, Farman M, Ahmad S, Zeb A. Extended hybrid controller design of bifurcation in a delayed chemostat model. MATCH Commun Math Comput Chem. 2023;90:609–48.
Xu C, Mu D, Pan Y, Aouiti C, Yao L. Exploring bifurcation in a fractionalorder predatorprey system with mixed delays. J Appl Anal Comput. 2023;13:1119–36.
Li P, Lu Y, Xu C, Ren J. Insight into hopf bifurcation and control methods in fractional order BAM neural networks incorporating symmetric structure and delay. Cogn Comput. 2023;15:1825–67.
Izadi M, Yüzbası, Ansari KJ. Application of VietaLucas series to solve a class of multipantograph delay differential equations with singularity. Symmetry. 2021;13(12):2370.
Razavi M, Hosseini MM, Salemi A. Error analysis and Kronecker implementation of Chebyshev spectral collocation method for solving linear PDEs. Comput Methods Differ Equ. 2022;10(4):914–27.
Izadi M, Afshar M. Solving the Basset equation via Chebyshev collocation and LDG methods. J Math Model. 2021;9(1):61–79.
Ahmed HM. Numerical solutions for singular LaneEmden equations using shifted Chebyshev polynomials of the first kind. Contemp Math. 2023;4(1):132–49.
Arora S, Bala I. Numerical study of sineGordon equations using Bessel collocation method. Iran J Numer Anal Optim. 2023;13(4):728–46.
Izadi M, Yüzbası, Cattani C. Approximating solutions to fractionalorder BagleyTorvik equation via generalized Bessel polynomial on large domains. Ricerche Mat. 2023;72(1):235–61.
Talaei Y, Micula S, Hosseinzadeh H, Noeiaghdam S. A novel algorithm to solve nonlinear fractional quadratic integral equations. AIMS Math. 2022;7(7):13237–57.
Izadi M. A novel matrix technique to solve a new singular nonlinear functional LaneEmden model. J Adv Math Model. 2022;12(2):232–47.
Mohseni Moghadam M, Saeedi H, Razaghzadeh N. A spectral Chelyshkov wavelet method to solve systems of nonlinear weakly singular Volterra integral equations. J Mahani Math Res. 2020;9(1):23–43.
Izadi M, Sene N, Adel W, ElMesady A. The Layla and Majnun mathematical model of fractional order: stability analysis and numerical study. Results Phys. 2023;51:106650.
Van den Driessche P, Watmough J. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002;180:29–48.
CastilloChavez C, Song B. Dynamical models of tuberculosis and their applications. Math Biosci Eng. 2004;1:361–404.
Mason JC, Handscomb D. Chebyshev Polynomials. Boca Raton: CRC Press; 2003.
Izadi M, Roul P. A highly accurate and computationally efficient technique for solving the electrohydrodynamic flow in a circular cylindrical conduit. Appl Numer Math. 2022;181:110–24.
Stewart GW. Afternotes on Numerical Analysis. vol 49. SIAM; 1996.
Bellman RE, Kalaba RE. Quasilinearization and Nonlinear BoundaryValue Problems. New York: Elsevier Publishing Company; 1965.
Izadi M. A combined approximation method for nonlinear foam drainage equation. Sci Iran. 2022;29(1):70–8.
Aznam SM, Ghani NA, Chowdhury MS. A numerical solution for nonlinear heat transfer of fin problems using the Haar wavelet quasilinearization method. Results Phys. 2019;14:102393.
Izadi M, Yüzbası, Adel W. Two novel Bessel matrix techniques to solve the squeezing flow problem between infinite parallel plates. Comput Math Math Phys. 2021;61(12):2034–53.
Funding
Not applicable.
Author information
Authors and Affiliations
Contributions
Conceptualization, M.I. and T.W.; methodology, M.I. and T.W.; software, M.I.; validation, M.I. and T.W.; formal analysis, M.I. and T.W.; investigation, M.I. and T.W.; writingoriginal draft preparation, M.I. and T.W.; writingreview and editing, M.I. and T.W. All authors have read and agreed to the published version of the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Izadi, M., Waezizadeh, T. Stability analysis and numerical evaluations of a COVID19 model with vaccination. BMC Med Res Methodol 24, 97 (2024). https://doi.org/10.1186/s12874024022092
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12874024022092