Stability and bifurcation analysis of $SIQR$ for the COVID-19 epidemic model with time delay

1.   Introduction

Mathematical modelling is the best way to formulating problems from an application area and it is well known that several mathematical characterization of numerous growth in chemical and physical sciences is described by differential equations (DEs). In chemistry, chemical kinetics problem and $CO_2$ with PGE problems are described by system of nonlinear (DEs) with different kind of Neumann boundary and Drichlet type conditions in different published work such as chemistry problem by Jawary and Raham ^[1], chemistry problem by Abbasbandy and Shirzadi ^[2], $CO_2$ absorbed into PGE problem by Jawary et al. ^[3], Choe et al. ^[4], Singha et al. ^[5], $CO_2$ absorbed into PGE problem by Robertson ^[6], chemistry problem by Matinfar et al. ^[7], chemistry problem by Ganji et al. ^[8], Dokoumetzidis et al. ^[9]. In the past few years, fractional calculus (FC) has found many diverse and robust applications in various research areas such as fluid dynamics, image processing, viscoelasticity and other physical phenomena. Many definitions of fractional derivatives are discovered by several mathematicians but two most famous definitions of fractional derivative are Riemann $-$ Liouville and Caputo. Some interesting and fundamental works on various direction of the FC is given in several famous books such as by Mainardi ^[10], fractional differential equations by Podlubny ^[11], Diethelm ^[12], Kilbas et al. ^[13] and Das ^[14].

In the past few years, wavelets have become an increasingly newly developed famous mechanism in the several research areas of physical, chemical, computational sciences, Image manipulation, signal analysis, data compression, numerical analysis and several others research areas such as a primer on wavelets and their scinentific applications by Walker ^[15], wavelet: mathematics and applications by Benedetto ^[16], a mathematical tool for signal analysis by Chui ^[17], wavelet methods for dynamical problems by Gopalakrishnan and Mitra ^[18], Wang ^[19] and a wavelet operational matrix method by Wu ^[20]. Due to this reason, wavelets have been applied for the solution of differential equations (DEs) since the 1980s. The interesting features in this method are possibility to find-out singularities, irregular structure and transient phenomena exhibited by the analysed equations such as by Heydari et al. ^[21], Wang and Fan ^[19], Balaji ^[22], Rehman and Khan ^[23], Hosseininia ^[24], Pirmohabbati et al. ^[25], Hosseininia ^[26], Heydari ^[27] and Kumar et al. ^[28].

Among the several wavelet families most simple are the Haar wavelets and it has been successfully applied to several linear and nonlinear problems of physical science and other research areas such as fractional order stationary neutron transport equation, neutron point kinetics equation, fractional order nonlinear oscillatory van der pol system and fractional bagley torvik equation by Ray and Patra ^{[29,30,31,32]}, a comparative study on haar wavelet and hybrid functions, nonlinear integral and integro $-$differential equation of first and higher order and parabolic differential equatons by Aziz et al. ^[33,34,35], burgers equation by Jiwari ^[36], fractional integral equations by Lepik ^[37], Poisson and biharmonic equations by Shi and Cao ^[38], delamination detection in composite beams by Hein and Feklistova ^[39], fractional order integral equations by Gao and Liao ^[40], lumped and distributed parameters systems by Chen and Hsiao ^[41], FDEs by Chen et al. ^[42], free vibration analysis by Xie et al. ^[43], fractional nonlinear differential equations by Saeed and Rehman ^[44], magnetohydrodynamic flow equations by Celik and Brahin ^[45,46], fishers equations by Hariharan et al. ^[47], FPDEs by Wang et al. ^[48], nonlinear oscillators equations Kaur et al. ^[49], poisson and biharmonic equations by Shi et al. ^[50] and free vibration analysis of functionally graded cylindrical shells by Jin et al. ^[51].

It is compulsory to note that the fractional chemical kinetics and condensations of $CO_2$ and PGE problems is the first one to be solved by the Haar wavelet and generalization of Adams– Bashforth$-$Moulton method by us. It is also noted that there are no similar works with these methods for fractional chemical kinetics and condensations of $CO_2$ and PGE problems available in any present published literature. It is well known by the several published research papers that the Caputo and Riemann-Liouville is most popular definition of fractional calculus.

The complete work is systematized in the following sections: Overview of basic FC are provided in section 2. Fractional Model of both Chemical Kinetics and $CO_2$ absorbed into PHE problems are provided in section 3. In section 4, a haar wavelet and Adam Bashforth's-Moulton methods are discussed and presented for both chemistry problems. The proposed methods for solutions of both chemistry problem are provided in section 5. Numerical result and discussions are provided in section 6. Conclusion and future scope are given in sections 7.

2.   Overview of FC

There are numerous definition of derivative and integration are available in literature ^{[52,53,54,55,56,57,58,59,60,61]}.

Definition 1. The (left sided) Riemann$-$Liouville fractional integral of order $\alpha > 0$ of a function $\Theta(t)\in C_\alpha, \alpha \geq-1$ is defined as,

where $\Gamma(.)$ is well known Gamma function.

Definition 2. The next two equations define Riemann $–$ Liouville and Caputo fractional derivatives of order a, respectively,

${}^{RL}D_{t}^\alpha \Theta(t) = \frac{{{d^m}}}{{d{t^m}}}\left({I^{m-\alpha} _t}\Theta(t)\right) = \left\{ \begin{array}{l} \frac{{{d^m}\Theta(t)}}{{d{t^m}}}, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \alpha = m \in N, \\ \frac{1}{{\Gamma (m - \alpha)}}\frac{{{d^m}}}{{d{t^m}}}\int\limits_0^t {\frac{{\Theta(\xi)}}{{{{(t - \xi)}^{\alpha - m + 1}}}}} \;\;\;\;\; d\xi, \, \, 0 \le m - 1 < \alpha < m, \end{array} \right.\, \,$

and,

${}^{C}D_{t}^\alpha \Theta(t) = {I^{m-\alpha} _t}\left(\frac{{{d^m}}}{{d{t^m}}} \Theta(t)\right) = \left\{ \begin{array}{l} \frac{{{d^m}\Theta(t)}}{{d{t^m}}}, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \alpha = m \in N, \\ \frac{1}{{\Gamma (m - \alpha)}}\int\limits_0^t {\frac{{\Theta^m(\xi)}}{{{{(t - \xi)}^{\alpha - m + 1}}}}}\;\;\;\;\; d\xi, \, \, 0 \le m - 1 < \alpha < m, \end{array} \right.\, \,$

where $t > 0$ and $m$ is an integer. Two basic properties for $m-1 < \alpha\leq m$ and $\Theta\in L_{1}[a, b]$ are given as

3.   Fractional model of Chemical Kinetics and $CO_2$ absorbed into PGE problems

3.1.   Fractional model of Chemical Kinetics (CK) problem

Let $\mathcal{D, E}$ and $\mathcal{H}$ are different location of a model of chemical process then the reactions are presented as

The concetrations of all three spaces of $\mathcal{D, E}$ and $\mathcal{H}$ are denoted by $\Theta_1, \Theta_2$ and $\Theta_3$ respectively. Let $r_1, r_2$ and $r_3$ denotes the reaction rate of Eqs (3.1), (3.2) and (3.3) respectively. We consider an integer order model of chemical kinetics problem as ^[1,2,6,7,8]

with the initial conditions, $\Theta_1(0) = 1$, $\Theta_2(0) = 0,$ $\Theta_3(0) = 0.$ The main target of this section is converted above inter order CK problem into fractional order CK problem. The fractional model of CK problem is presented as

with the initial conditions, $\Theta_1(0) = 1$, $\Theta_2(0) = 0,$ $\Theta_3(0) = 0$ where, $D_t^\alpha = \frac{d^\alpha}{dt^\alpha}, D_t^\beta = \frac{d^\beta}{dt^\beta}, D_t^\gamma = \frac{d^\gamma}{dt^\gamma}$ are fractional derivative with $0 < \alpha, \beta, \gamma\leq 1.$ If $r_1 = 1, r_2 = 0,$ and $r_3 = 1$ then

with the initial conditions, $\Theta_1(0) = 1$, $\Theta_2(0) = 0,$ $\Theta_3(0) = 0.$ The above system is representing a nonlnear reaction which was taken from litrature ^[2,7,8,62].

3.2.   Fractional model of condensations of $CO_{2}$ and PGE problem

The $CO_{2}$ causes in ocean acidification because it dissolves in water to form carbonic acid ^[63].The mathematical formulation of the concentration of $CO_{2}$ and PGE is shown in Muthukaruppan et al. ^[64]. Now, the two nonlinear reactions equations in normalized form is presented as

with boundary conditions, $\Upsilon_1(0) = 0$, $\Upsilon_1(1) = \frac{1}{m}$, ${{\Upsilon'_2}}(0) = \frac{1}{m}$, $\Upsilon_2(1) = \frac{1}{m}$. The whole chemistry of the above problem is given in several litratures ^[1,3,4]. The fractional model of the condensation of $CO_{2}$ and PGE in operator form is given as,

with the same boundary conditions $\Upsilon_1(0) = 0$, $\Upsilon_1(1) = \frac{1}{m}$ and ${{\Upsilon'}_2}(0) = \frac{1}{m}, \, \, {\Upsilon_2}(1) = \frac{1}{m}$; where $m \ge 3$ and fractional operator is taken in Caputo sence.

4.   Overview of Haar wavelet and Adam Bashforth's$-$Moulton (ABM) methods

4.1.   Haar wavelet

The Haar functions have been discovered by Alfred Haar in 1910 and Haar wavelets are the simplest wavelet among all wavelet. The Haar sequence was also introduced by itself Alfred Haar in 1909 which is recognised as wavelet basis. The Haar wavelets are the mathematical operations which are known as Haar transform. These wavelets are build up by piecewise constant function on the real line. We used Haar wavelet operational matrix method because of its flexibility, simplicity and require very less effort of computation. Usually Haar wavelet is defined for [0, 1) but in general case we extend it up to certain interval. Haar functions are very useful in many applications as image coding, extraction of edge, binary logic design etc ^{[20,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51]}. The Haar scaling function is defined as

The Haar wavelet mother function is defined as

The orthogonal set of Haar wavelet functions for $t \in [0, 1]$ are defined as

where $i = 0, 1, 2, ..., m - 1$, $m = {2^{r + 1}}$ and $r$ is positive integer known as resolution of Harr wavelet. Also $j$ and $k$ represent integer decomposition of $i = {2^j} + k - 1$.

4.2.   Function approximation

Any function $\Theta(t) \in {L^2}([0, 1))$ can be expanded in terms of Haar wavelet by

If we approximated as piecewise constant during each interval, Eq. 4.4 will terminated at finite terms as ^[65]:

where $C_m = {[{c_0}, {c_1}, {c_2}, ..., {c_{m-1}}]^T}$ and ${H_m}(t) = {[{h_0}(t), {h_1}(t), {h_2}(t), ..., {h_{m-1}}(t)]^T}$,

Using collocation points ${t_l} = \frac{{(l - 0.5)}}{m}$, where $l = 0, 1, ..., m - 1,$ we obtained the discrete form as

4.3.   Haar wavelet operational matrix (HWOM) of fractional order integration

The HWOM of fractional order integration without using block pulse functions we integrate ${H_m}(t)$ using Reimann-Liouville integration operator ^[41,66]. Then the HWOM of fractional order integration ${Q^\alpha }$ is given by

where

$Q{h_0}(t) = \, \, \frac{1}{{\sqrt m }}\frac{{{t^\alpha }}}{{\Gamma (1 + \alpha)}}, 0 \le t \le 1,$

$Q{h_i}(t) = \, \, \frac{1}{{\sqrt m }}\left\{ \begin{array}{l} 0, \, \, \, \, \, \, \, \, 0 \le t < \frac{{k - 1}}{{{2^j}}}, \\ {2^{j/2}}{\zeta _1}(t)\, \, \, \, \, \, \frac{{k - 1}}{{{2^j}}}\, \, \le t < \frac{{k - 0.5}}{{{2^j}}}, \\ {2^{j/2}}{\zeta _2}(t)\, \, \, \, \, \, \frac{{k - 0.5}}{{{2^j}}} \le t < \frac{k}{{{2^j}}}, \\ {2^{j/2}}{\zeta _3}(t)\, \, \, \, \, \, \frac{k}{{{2^j}}} \le t < 1, \, \end{array} \right.$

where

${\zeta _1}(t)\, = \, \frac{1}{{\Gamma (1 + \alpha)}}{\left({t - \frac{{k - 1}}{{{2^j}}}} \right)^\alpha },$

${\zeta _2}(t)\, = \, \frac{1}{{\Gamma (1 + \alpha)}}{\left({t - \frac{{k - 1}}{{{2^j}}}} \right)^\alpha } - \frac{2}{{\Gamma (1 + \alpha)}}{\left({t - \frac{{k - 0.5}}{{{2^j}}}} \right)^\alpha },$

${\zeta _3}(t) = \, \frac{1}{{\Gamma (1 + \alpha)}}{\left({t - \frac{{k - 1}}{{{2^j}}}} \right)^\alpha } - \frac{2}{{\Gamma (1 + \alpha)}}{\left({t - \frac{{k - 0.5}}{{{2^j}}}} \right)^\alpha } + \, \frac{1}{{\Gamma (1 + \alpha)}}{\left({t - \frac{k}{{{2^j}}}} \right)^\alpha }.$

If we take, $\alpha = 1/2, \, m = 8,$ then we have the operational matrix as given below:

$Q^{1/2}H_8$ = $\begin{bmatrix} 0.0997 & 0.1727 & 0.2230 & 0.2639 & 0.2992 & 0.3308 & \, \, \, \, 0.3596 & \, \, \, \, 0.3863 \\ 0.0997 & 0.1727 & 0.2230 & 0.2639 & 0.0997 & -0.0147 & -0.0864 & -0.1415 \\ 0.2443 & 0.0333 & -0.1154 & -0.0666 & -0.0343 & -0.0223 & -0.0193 & -0.0132 \\ 0 & 0 & 0 & 0 & 0.1410 & 0.2443 & 0.033 & -0.1154 \\ 0.1995 & -0.0534 & -0.0455 & -0.0188 & -0.0111 & -0.0075 & -0.0055 & -0.0043 \\ 0 & 0 & 0.1995 & -0.0534 & -0.0455 & -0.0188 & -0.0111 & -0.0075 \\ 0 & 0 & 0 & 0 & 0.1995 & -0.0534 & -0.0455 & -0.0188 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0.1995 & -0.0534 \\ \end{bmatrix}$.

The above matrix is the operational matrix of Haar wavelets.

4.4.   Adam Bashforth's-Moulton method

In this section we discuss about Predictor-Corrector scheme (PECE), which is the genralization of (ABM) mehod ^[67,68]. We obtain the numerical solution of nonlinear FDES as

where derivative in Caputo's sense. which is equivalent to the Volterra integral equation

Assume $h = T/N$, $t_n = nh$, $n = 0, 1, 2, ..., \mathrm{N}$ $\in{\mathrm{Z}^+}$ then the discrete form for the above equation will be

The corrector values for chemistry problem is

The corresponding predictor values are,

From Eqs (4.12) and (4.14) we can calculate {$\alpha_{j, n+1}$}, {$\beta_{j, n+1}$}, {$\gamma_{j, n+1}$}, and $B_{j, n+1}$, $C_{j, n+1}$, $D_{j, n+1}$.

5.   Proposed methods for solution of both chemistry problems

Example: 1 We assume a fractional model of chemical kinetics problem is given as

with the initial conditions, $\Theta_1(0) = 1$, $\Theta_2(0) = 0$ $\Theta_3(0) = 0$, where $r_1$, $r_2$ and $r_3$ are reaction rates. Let us assume higher derivatives in the terms of haar wavelet series.

where ${C} = {[{c_0}, {c_1}, {c_2}, ..., {c_{m-1}}]^T}$, ${G} = {[{g_0}, {g_1}, {g_2}, ..., {g_{m-1}}]^T}$ and ${K} = {[{k_0}, {k_1}, {k_2}, ..., {k_{m-1}}]^T}$ are unknown vectors. Applying Riemann-Liouville fractional integral in Eq. $(5.2)$ and using initial conditions, we obtained

Now substituting the values of $\Theta_1$, $\Theta_2$ and $\Theta_3$ into the Eq. $(5.1)$, we obtained.

Let $r_1 = 0.1$, $r_2 = 0.02$ and $r_3 = 0.009$ as given in Aminikhah ^[69]. Now disperse the Eq. $(5.4)$ at the collocation points ${t_l} = \frac{{(l - 0.5)}}{m}$, where $l = 1, 2, 3, ..., m$. We obtained $3m$ nonlinear algebraic equations which can be solved by Newton iteration method, after solving we obtained the coefficients $c_i$, $g_i$ and $k_i$. Substitute these coefficients into the Eq. $(5.3)$ we get desired solutions $\Theta_1$, $\Theta_2$ and $\Theta_3$.

Example 2: Consider the system of condensations of $CO_2$ and PGE problem of arbitrary order.

with boundary conditions $\Upsilon_1(0) = 0$, $\Upsilon_1(1) = \frac{1}{m}$, and ${{\Upsilon'}_2}(0) = \frac{1}{m}$, $\Upsilon_2(1) = \frac{1}{m}$ where ${\Upsilon}_1^\alpha (t) =$ ${}^CD_{t}^\alpha {\Upsilon_1}(t)$. Here for simplicity we have taken $m = 3$ and we will take the value of ${\alpha _1} = 1$, ${\alpha _2} = 2$, ${\beta _1} = 1$ and ${\beta _2} = 3$ as given in Duan et al. ^[70], AL-jawary ad Radhi ^[71]. Further, we assume the higher derivative in terms of Haar wavelet series.

applying Riemann-Liouville integral operator on the above equation and using boundary conditions, we obtained

substituting $t = 1$ into Eq. $(5.7)$ we obtained

${\Upsilon_1}(1) - {{\Upsilon'}_1}(0) = {C^T}{Q^\alpha }{H_m}(1)$

and

Therefore,

similarly

Substituting the values of $\Upsilon_1$, $\Upsilon_2$ into the Eq. $(5.5)$ and using Eq. $(5.6)$ we obtained

Now disperse the Eqs $(5.12)$ and $(5.13)$ at the collocation points ${t_l} = \frac{{(l - 0.5)}}{m}$, where $l = 0, 1, ..., m - 1$. We obtained a system of nonlinear algebraic equations which can be easily solved by Newton-Iteration method using mathematical softwares, after solving we obtained the unknowns coefficients $c_i$ and $k_i$. Substituting these coefficients into the Eqs $(5.10)$ and $(5.11)$ we get desired solutions $\Upsilon_1$ and $\Upsilon_2$.

6.   Numerical results and discussion

All numerical simulation and graphical results of both examples are depicted through the Figures 1–14 where Figures 1–6 and Figures 7–14 are depicted for examples 1 and 2 respectively. We have depicted a comparison between numerical obtained solutions using by Haar wavelet and Adam's-Bashforth-Moulton predictor-corrector schemes through the Figures 1–3 and these figures are depicted for the values of $m = 64$. It is clear from all figures that both obtained solutions by HWM and ABM are identical. The obtained solutions $\Theta_1,$ $\Theta_2$ and $\Theta_3$ are plotted through the Figures 3–6 where the nature of solution $\Theta_1$ is of decreasing nature while other solutions $\Theta_2$ and $\Theta_3$ is of increasing nature. We plotted the resolutions Figures 7–14 for better understanding the nature of obtained solution of example 2. We plotted resolutions figures due to non-availability of its exact solution.

7.   Conclusion

In this work, Haar wavelet operational matrix and Adam Bashforth's Moulton scheme are proposed to solve fractional chemical kinetics and another problem that relates the condensations of carbon dioxide $CO_2$ numerically. A comparative study between fractional chemical kinetics and another problem that relates the conden sations of carbon dioxide $CO_2$ has been done for $m = 64$ in this work. Our tabulated and graphical results indicate that the solution will ameliorate if we will take more collocation points, i.e greater values of $m$. The essential advantage of HWM is that it converts problems into the system of linear or nonlinear algebraic equations so that the computation is facile and computer-oriented. Furthermore, wavelet method is much easier than other numerical methods for system of FDEs. Again, we have solved the chemistry problems at different resolutions, which produced the same results at each resolution. The precision of the solution will ameliorate if we increase the resolution. This new comparative study between the Haar wavelet operational matrix and Adam Bashforth's Moulton scheme for fractional chemical kinetics and another problem that relates the condensations of carbon dioxide $CO_2$ indicates that both approaches can be applied successfully to the chemistry problems of chemistry science.

Acknowledgement

The first author Dr. Sunil Kumar would like to acknowledge the financial support received from the National Board for Higher Mathematics, Department of Atomic Energy, Government of India (Approval No. 2/48(20)/2016/ NBHM(R.P.)/R and D II/1014). The authors are also grateful to the editor and anonymous reviewers for their constructive comments and valuable suggestions to improve the quality of article.

Conflict of interest

The authors declare no conflict of interest in this manuscript.