Numerical Analysis of the Susceptible Exposed Infected Quarantined and Vaccinated (SEIQV) Reaction-Diffusion Epidemic Model

In this paper, two structure-preserving nonstandard finite difference (NSFD) operator splitting schemes are designed for the solution of reaction diffusion epidemic models. The proposed schemes preserve all the essential properties possessed by the continuous systems. These schemes are applied on a diffusive SEIQV epidemic model with a saturated incidence rate to validate the results. Furthermore, the stability of the continuous system is proved, and the bifurcation value is evaluated. A comparison is also made with the existing operator splitting numerical scheme. Simulations are also performed for numerical experiments.

1. Introduction

Mathematical modeling has a prominent role in describing physical phenomena in various disciplines of mathematics, physical sciences, social sciences, engineering, life sciences, and many more [1–6]. The transmission of infectious diseases and the control of their spread can be studied effectively by constructing mathematical models for various strategies like vaccination and quarantine. The word quarantine denotes forced isolation or being cut off from interactions with others. Quarantine is an effective intervention process for restraining the spread of infection by isolating individuals who are affected by the disease. Such isolation has been adopted to decrease the communication of infectious diseases like dengue, measles, smallpox, cholera, leprosy, tuberculosis, and many more.

Epidemic models, that is, mathematical models of infectious diseases, are a simplified way to illustrate the transmission dynamics of the complicated nonlinear processes and complex behavior of an infectious disease in individuals within a population. These are deterministic models that are used to allocate the population to different subclasses or compartments, describing a particular stage of the epidemic. The incidence rate, which is proportional to the number of susceptible and infected persons, is an important parameter of compartment-based epidemic models. The mathematical models of infectious diseases are often based on bilinear incidence rate βSI, but a more concise approach to use the saturated incidence rate rather than the bilinear incidence rate. In the saturated incidence rate $\frac{\beta SI}{1+\alpha I}$, if number of infected individuals I becomes very large, $\frac{\beta SI}{1+\alpha I}$ approaches the saturation level. The infection force is measured by βI, which describes the penetration of the disease into a fully susceptible population. $\frac{1}{1+\alpha I}$ is used to measure the inhibition effect of behavioral change of susceptible persons. Liu and Yang [7] proposed the SEIQV epidemic model, which uses the saturated incidence rate. The model is expressed as:

$\begin{array}{ll}\frac{dS(t)}{dt}=b-\frac{\beta S(t)I(t)}{1+\alpha I(t)}-(\omega +\mu +q_3)S(t)& (1)\end{array}$

$\begin{array}{ll}\frac{dE(t)}{dt}=\frac{\beta S(t)I(t)}{1+\alpha I(t)}-(\mu +\sigma +q_2)E(t)& (2)\end{array}$

$\begin{array}{ll}\frac{dI(t)}{dt}=\sigma E(t)-(\mu +ϵ+\gamma +q_1)I(t)& (3)\end{array}$

$\begin{array}{ll}\frac{dQ(t)}{dt}=q_{3}S(t)+q_{2}E(t)+q_{1}I(t)-(\mu +\varphi )Q(t)& (4)\end{array}$

$\begin{array}{ll}\frac{dV(t)}{dt}=\omega S(t)+\varphi Q(t)+\gamma I(t)-\mu V(t)& (5)\end{array}$

The variables and parameters of the model are defined as:

S(t) = Susceptible persons at time t,

E(t) = Exposed persons at time t,

I(t) = Infected persons at time t,

Q(t) = Quarantined persons at time t,

V(t) = Vaccinated persons at time t,

b = Rate of recruitment,

β = Rate of transmission,

μ = Rate of natural death,

ϵ = Rate of death due to disease in infected compartment,

α = Parameter that measures psychological or inhibitory effects,

γ = Rate at which infected individuals are being vaccinated infected persons,

σ = Rate at which exposed persons become infected,

ω = Rate at which infected individuals are being vaccinated susceptible persons,

ϕ = Rate at which infected individuals are being vaccinated quarantined persons,

q₁, q₂, q₃ = Effective quarantine probabilities.

The above model (1)–(5) assumes a homogeneous population, where the population mixes in such a way that there is no difference between person in one place and person in another place. However, in actual scenarios, the disease may spread faster in one place than in another because of different circumstances like different weather conditions, etc. Hence, it is essential for the variables to depend on space also. Therefore, we extend system (1)–(5) to make it a reaction-diffusion system by adding a diffusion term.

$\begin{array}{l}\frac{\partial S(x,t)}{\partial t}=b-\frac{\beta S(x,t)I(x,t)}{1+\alpha I(x,t)}-(\omega +\mu +q_3)S(x,t)\\ +d_1\frac{{\partial }^{2}S(x,t)}{\partial x^2}& (6)\end{array}$

$\begin{array}{ll}\frac{\partial E(x,t)}{\partial t}=\frac{\beta S(x,t)I(x,t)}{1+\alpha I(x,t)}-(\mu +\sigma +q_2)E(x,t)+d_2\frac{{\partial }^{2}E(x,t)}{\partial x^2}& (7)\end{array}$

$\begin{array}{ll}\frac{\partial I(x,t)}{\partial t}=\sigma E(x,t)-(\mu +ϵ+\gamma +q_1)I(x,t)+d_3\frac{{\partial }^{2}I(x,t)}{\partial x^2}& (8)\end{array}$

$\begin{array}{l}\frac{\partial Q(x,t)}{\partial t}=q_{3}S(x,t)+q_{2}E(x,t)+q_{1}I(x,t)-(\mu +\varphi )Q(x,t)\\ +d_4\frac{{\partial }^{2}Q(x,t)}{\partial x^2}& (9)\end{array}$

$\begin{array}{l}\frac{\partial V(x,t)}{\partial t}=\omega S(x,t)+\varphi Q(x,t)+\gamma I(x,t)-\mu V(x,t)\\ +d_5\frac{{\partial }^{2}V(x,t)}{\partial x^2}& (10)\end{array}$

with the initial conditions:

$\begin{array}{ll}S(x,0)=g_1(x)0\le x\le L& (11)\end{array}$

$\begin{array}{ll}E(x,0)=g_2(x)0\le x\le L& (12)\end{array}$

$\begin{array}{ll}I(x,0)=g_3(x)0\le x\le L& (13)\end{array}$

$\begin{array}{ll}Q(x,0)=g_4(x)0\le x\le L& (14)\end{array}$

$\begin{array}{ll}V(x,0)=g_5(x)0\le x\le L& (15)\end{array}$

The boundary conditions are no flux,

$\begin{array}{ll}{S}_x(0,t)=S_x(L,t)=0& (16)\end{array}$

$\begin{array}{ll}{E}_x(0,t)=E_x(L,t)=0& (17)\end{array}$

$\begin{array}{ll}{I}_x(0,t)=I_x(L,t)=0& (18)\end{array}$

$\begin{array}{ll}{Q}_x(0,t)=Q_x(L,t)=0& (19)\end{array}$

$\begin{array}{ll}{V}_x(0,t)=V_x(L,t)=0& (20)\end{array}$

Epidemic models always demonstrate two equilibrium points: the disease-free equilibrium (DFE) point and the endemic equilibrium (EE) point. The DFE point exists if R₀ < 1, where R₀ is the reproductive number, which basically measures the occurrence of disease. The EE point exists if R₀ > 1. This implies that the SEIQV reaction-diffusion system (6–10) always converges to the DFE point or EE point if R₀ < 1 or R₀ > 1, respectively. Analytical solution of the SEIQV epidemic system is not possible, so we have to use numerical techniques to find its solution. Note that the numerical technique must show the same behavior as is possessed by the continuous SEIQV reaction-diffusion epidemic system.

In this work, we propose two operator-splitting NSFD methods, one explicit and one implicit. These methods are used to solve the SEIQV epidemic model with diffusion. As S, E, I, Q, and V are population sizes and evaluated in absolute scale, we propose NSFD methods because they give a positive solution. Also, the convergence of the proposed NSFD operator splitting methods toward the equilibrium points is the same as the convergence of continuous an SEIQV reaction-diffusion epidemic system. The proposed splitting methods are designed with the aid of rules given by Mickens [8]. In the recent era, positivity preserving FD methods have gained importance, as many physical dynamical systems possess the positivity property [9–11]. The NSFD method presented by Mickens [8, 12, 13] has becomes an effective and important structure-preserving FD method for solving differential equations. In epidemic models, population dynamics and population size cannot be negative, so the numerical technique must be a positivity-preserving technique. Various authors have used different positivity-preserving numerical techniques for the approximate solution of epidemic models: see, for example [14–22].

In this work, we also show the numerical stability of the SEIQV epidemic model with diffusion and evaluate the bifurcation value of the vaccination parameter ω with the aid of the Routh-Hurwitz method.

2. Equilibrium Points

The model (6–10) has two equilibrium points [7], the DFE point and EE point. The DFE point is:

$\begin{array}{l}DFE=(S_0,E_0,I_0,Q_0,V_0)\\ =(\frac{b}{\omega +\mu +q_3},0,0,\frac{q_3{S}_0}{\mu +\varphi },\frac{\omega S_0+\varphi Q_0}{\mu })& (21)\end{array}$

and the EE point is:

$\begin{array}{ll}EE=(S_{*},E_{*},I_{*},Q_{*},V_{*})& (22)\end{array}$

where,

$\begin{array}{l}{S}_{*}=\frac{(\mu +\sigma +q_2)(\mu +ϵ+\gamma +q_1)(1+\alpha I)}{\gamma \beta }\\ E_{*}=\frac{(\mu +ϵ+\gamma +q_1)I}{\sigma }\\ Q_{*}=\frac{q_{3}S+q_{2}E+q_{1}I}{\mu +\varphi }\\ V_{*}=\frac{\omega S+\varphi Q+\gamma I}{\mu }\\ I_{*}=\frac{\sigma \beta b-(\omega +\mu +q_3)(\mu +\sigma +q_2)(\mu +ϵ+\gamma +q_1)}{(\mu +\sigma +q_2)(\mu +ϵ+\gamma +q_1)(\beta +\alpha (\omega +\mu +q_3))}\end{array}$

Reproductive number R₀ is given as:

$\begin{array}{l}{R}_0=\frac{\sigma \beta b}{(\omega +\mu +q_3)(\mu +\sigma +q_2)(\mu +ϵ+\gamma +q_1)},\\ \text{ when},d_1=d_2=d_3=d_4=d_5=0\end{array}$

R₀ is the reproductive value. Now, if R₀ < 1, the model acquires a DFE point, and if R₀ > 1, the model acquires an EE point.

3. Numerical Stability of the SEIQV Model at Equilibrium Point

We evaluated the small perturbation S₁(x, t), E₁(x, t), I₁(x, t), Q₁(x, t), and V₁(x, t) so that (6)–(10) is linearized at the EE point (S_*, E_*, I_*, Q_*, V_*), as discussed in Chakrabrty et al. [23].

$\begin{array}{ll}\frac{\partial S_1}{\partial t}=a_{11}{S}_1+a_{12}{E}_1+a_{13}{I}_1+a_{14}{Q}_1+a_{15}{V}_1+d_1\frac{{\partial }^2{S}_1}{\partial x^2}& (23)\end{array}$

$\begin{array}{ll}\frac{\partial E_1}{\partial t}=a_{21}{S}_1+a_{22}{E}_1+a_{23}{I}_1+a_{24}{Q}_1+a_{25}{V}_1+d_2\frac{{\partial }^2{E}_1}{\partial x^2}& (24)\end{array}$

$\begin{array}{ll}\frac{\partial I_1}{\partial t}=a_{31}{S}_1+a_{32}{E}_1+a_{33}{I}_1+a_{34}{Q}_1+a_{35}{V}_1+d_3\frac{{\partial }^2{I}_1}{\partial x^2}& (25)\end{array}$

$\begin{array}{ll}\frac{\partial Q_1}{\partial t}=a_{41}{S}_1+a_{42}{E}_1+a_{43}{I}_1+a_{44}{Q}_1+a_{45}{V}_1+d_4\frac{{\partial }^2{Q}_1}{\partial x^2}& (26)\end{array}$

$\begin{array}{ll}\frac{\partial V_1}{\partial t}=a_{51}{S}_1+a_{52}{E}_1+a_{53}{I}_1+a_{54}{Q}_1+a_{55}{V}_1+d_5\frac{{\partial }^2{V}_1}{\partial x^2}& (27)\end{array}$

Suppose a Fourier series solution is demonstrated for Equations (23)–(27) of the form:

$\begin{array}{ll}{S}_1(x,t)={\displaystyle {\displaystyle \sum }_k}{S}_k{e}^{\lambda t}cos(kx)& (28)\end{array}$

$\begin{array}{ll}{E}_1(x,t)={\displaystyle {\displaystyle \sum }_k}{E}_k{e}^{\lambda t}cos(kx)& (29)\end{array}$

$\begin{array}{ll}{I}_1(x,t)={\displaystyle {\displaystyle \sum }_k}{I}_k{e}^{\lambda t}cos(kx)& (30)\end{array}$

$\begin{array}{ll}{Q}_1(x,t)={\displaystyle {\displaystyle \sum }_k}{Q}_k{e}^{\lambda t}cos(kx)& (31)\end{array}$

$\begin{array}{ll}{V}_1(x,t)={\displaystyle {\displaystyle \sum }_k}{V}_k{e}^{\lambda t}cos(kx)& (32)\end{array}$

Here, k = nπ/2, (n = 1, 2, 3,…) exhibits the value of the wave number for the node n. Substituting Equations (28)–(32) in the system (23)–(27), the system is converted into:

$\begin{array}{l}{\displaystyle {\displaystyle \sum }_k}(a_{11}-d_1{k}^2-\lambda )S_k+{\displaystyle {\displaystyle \sum }_k}{a}_{12}{E}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{13}{I}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{14}{Q}_k\\ +{\displaystyle {\displaystyle \sum }_k}{a}_{15}{V}_k=0& (33)\end{array}$

$\begin{array}{l}{\displaystyle {\displaystyle \sum }_k}{a}_{21}{S}_k+{\displaystyle {\displaystyle \sum }_k}(a_{22}-d_2{k}^2-\lambda )E_k+{\displaystyle {\displaystyle \sum }_k}{a}_{23}{I}_k\\ +{\displaystyle {\displaystyle \sum }_k}{a}_{24}{Q}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{25}{V}_k=0& (34)\end{array}$

$\begin{array}{l}{\displaystyle {\displaystyle \sum }_k}{a}_{31}{S}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{32}{E}_k+{\displaystyle {\displaystyle \sum }_k}(a_{33}-d_3{k}^2-\lambda )I_k\\ +{\displaystyle {\displaystyle \sum }_k}{a}_{34}{Q}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{35}{V}_k=0& (35)\end{array}$

$\begin{array}{l}{\displaystyle {\displaystyle \sum }_k}{a}_{41}{S}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{42}{E}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{43}{I}_k\\ +{\displaystyle {\displaystyle \sum }_k}(a_{44}-d_4{k}^2-\lambda )Q_k+{\displaystyle {\displaystyle \sum }_k}{a}_{45}{V}_k=0& (36)\end{array}$

$\begin{array}{l}{\displaystyle {\displaystyle \sum }_k}{a}_{51}{S}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{52}{E}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{53}{I}_k+{\displaystyle {\displaystyle \sum }_k}{a}_{54}{Q}_k\\ +{\displaystyle {\displaystyle \sum }_k}(a_{55}-d_5{k}^2-\lambda )V_k=0& (37)\end{array}$

The variational matrix ${V}$ for the system (33)–(37) is:

$\begin{array}{ll}{V}=\left(\begin{array}{ccccc}{a}_{11}-d_1{k}^2& a_{12}& a_{13}& a_{14}& a_{15}\\ a_{21}& a_{22}-d_2{k}^2& a_{23}& a_{24}& a_{25}\\ a_{31}& a_{32}& a_{33}-d_3{k}^2& a_{34}& a_{35}\\ a_{41}& a_{42}& a_{43}& a_{44}-d_4{k}^2& a_{45}\\ a_{51}& a_{52}& a_{53}& a_{54}& a_{55}-d_5{k}^2\end{array}\right)& (38)\end{array}$

where,

$\begin{array}{l}{a}_{11}=-\frac{\beta I_{*}}{(1+\alpha I_{*})}-(\omega +\mu +q_3),a_{12}=0,a_{13}=-\frac{\beta S_{*}}{{(1+\alpha I_{*})}^2}\\ a_{14}=0,a_{15}=0,a_{21}=\frac{\beta I_{*}}{(1+\alpha I_{*})},a_{22}=-(\sigma +\mu +q_2)\\ a_{23}=\frac{\beta S_{*}}{{(1+\alpha I_{*})}^2},a_{24}=0,a_{25}=0,a_{31}=0,a_{32}=\sigma \\ a_{33}=-(\mu +ϵ+\gamma +q_1),a_{34}=0,a_{35}=0,a_{41}=q_3,a_{42}=q_2,\\ a_{43}=q_1\\ a_{44}=-(\mu +\varphi ),a_{45}=0,a_{51}=\omega ,a_{52}=0,a_{53}=\gamma ,a_{54}=\varphi ,\\ a_{55}=-\mu \end{array}$

The characteristics equation for matrix ${V}$ is:

$\begin{array}{l}{\lambda }^5+{\xi }_1{\lambda }^4+{\xi }_2{\lambda }^3+{\xi }_3{\lambda }^2+{\xi }_4\lambda +{\xi }_5=0\end{array}$

The expressions for ξ₁, ξ₂, ξ₃, ξ₄, and ξ₅ with diffusion and without diffusion are mentioned in Islam and Haider [24].

The Routh-Hurwitz stability criterion gives:

$\begin{array}{l}{\xi }_1>0,{\xi }_2>0,{\xi }_3>0,{\xi }_4>0,{\xi }_5>0\end{array}$

$\begin{array}{l}{t}_1={\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4>0\end{array}$

and

$\begin{array}{l}{t}_2=({\xi }_1{\xi }_4-{\xi }_5)({\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4)-{\xi }_5{({\xi }_1{\xi }_2-{\xi }_3)}^2\\ -{\xi }_1{\xi }_5^2>0.\end{array}$

The Table 2 reflects the numerical stability of the equilibrium point against the various cases, as discussed in Table 1.

TABLE 1

Table 1. Values of parameters.

TABLE 2

Table 2. Stability of equilibrium point.

4. Bifurcation Value of Vaccination Parameter ω Independent of Diffusion

Considering the vaccination parameter ω, to find its bifurcation value, a₁₁, a₁₂, are used instead of S_*, E_*, I_*, Q_*, and V_*.

$\begin{array}{l}{a}_{11}=-0.6893640968\omega -0.2309369724,\\ a_{13}=-1.1198107335\omega -0.3751365957\\ a_{21}=0.0709369724-0.3106359032\omega ,a_{22}=-0.96\\ a_{23}=1.1198107335\omega +0.3751365957,a_{32}=0.7,\\ a_{33}=-0.46,a_{41}=0.1,a_{42}=0.2\\ a_{43}=0.2,a_{44}=-0.46,a_{51}=\omega ,a_{53}=0.15,a_{54}=0.4,\\ a_{55}=-0.06\\ a_{12}=a_{14}=a_{15}=a_{24}=a_{25}=a_{31}=a_{34}=a_{35}=a_{45}=\\ a_{52}=0\end{array}$

The Routh-Hurwitz criterion for stability gives:

$\begin{array}{l}{\xi }_1=0.6893640968\omega +2.1709369724=f_1(\omega )\\ {\xi }_2=0.5534988343\omega +1.3930221095=f_2(\omega )\\ {\xi }_3=-0.7838675135{\omega }^2+0.0368505571\omega \\ +0.3691384683=f_3(\omega )\\ {\xi }_4=-0.4076111070{\omega }^2-0.0380846274\omega \\ +0.0451739663=f_4(\omega )\\ {\xi }_5=-0.0216347434{\omega }^2-0.0023071181\omega \\ +0.00165507=f_5(\omega )\\ {\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4=-0.7198363952{\omega }^4-0.3846861900{\omega }^3\\ +0.4409094061{\omega }^2+0.9265604777\omega \\ +0.7671683354=f_6(\omega )\end{array}$

$\begin{array}{l}({\xi }_1{\xi }_4-{\xi }_5)({\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4)-{\xi }_5{({\xi }_1{\xi }_2-{\xi }_3)}^2-{\xi }_1{\xi }_5^2\hfill \\ =0.2022686014{\omega }^7+0.7777858386{\omega }^6+0.3637035355{\omega }^5\hfill \\ -0.4633365077{\omega }^4-0.8380985505{\omega }^3-0.5245415951{\omega }^2\hfill \\ +0.0491674986\omega +0.0622935244=f_7(\omega )\hfill \end{array}$

where the values of ξ₁, ξ₂, ξ₃, ξ₄, and ξ₅ are obtained from the expression of the characteristic equation (without diffusion) given in paper [24].

f₅(ω) = 0 gives the value of bifurcation for ω. This value transfers the stability of a continuous model from stable to unstable. f₅(ω) = 0 provides the bifurcation value ω = 0.228360507. The EE point is stable for ω less than ω = 0.228360507.

5. Bifurcation Value of Vaccination Parameter ω With Diffusion

For the bifurcation value of vaccination parameter ω, the values of S_*, E_*, I_*, Q_*, and V_* are replaced into a₁₁, a₁₂, to

$\begin{array}{l}{a}_{11}=-0.6893640968\omega -0.2309369724,\\ a_{13}=-1.1198107335\omega -0.3751365957\\ a_{21}=0.0709369724-0.3106359032\omega ,a_{22}=-0.96\\ a_{23}=1.1198107335\omega +0.3751365957,a_{32}=0.7,a_{33}=-0.46,\\ a_{41}=0.1\\ a_{42}=0.2,a_{43}=0.2,a_{44}=-0.46,a_{51}=\omega ,a_{53}=0.15,\\ a_{54}=0.4,a_{55}=-0.06\\ a_{14}=a_{12}=a_{15}=a_{24}=a_{25}=a_{31}=a_{34}=a_{35}\\ =a_{45}=a_{52}=0\end{array}$

The Routh-Hurwitz criterion for stability gives:

$\begin{array}{l}{\xi }_1=0.6893640968\omega +2.3707964615=f_1(\omega )\\ {\xi }_2=0.60622790397\omega +1.7722067841=f_2(\omega )\\ {\xi }_3=-0.7838675135{\omega }^2-0.0208144715\omega \\ +0.5625469972=f_3(\omega )\\ {\xi }_4=-0.4462934183{\omega }^2-0.0884716755\omega \\ +0.0784487226=f_4(\omega )\\ {\xi }_5=-0.0321693682{\omega }^2-0.0070094221\omega \\ +0.0035676470=f_5(\omega )\\ {\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4=-0.7299468904{\omega }^4-0.6247499708{\omega }^3\\ +0.5281660018{\omega }^2+1.6725898825\omega \\ +1.6061706299=f_6(\omega )\end{array}$

$\begin{array}{l}({\xi }_1{\xi }_4-{\xi }_5)({\xi }_1{\xi }_2{\xi }_3-{\xi }_3^2-{\xi }_1^2{\xi }_4)-{\xi }_5{({\xi }_1{\xi }_2-{\xi }_3)}^2-{\xi }_1{\xi }_5^2\hfill \\ =0.2245744816{\omega }^7+1.0320435772{\omega }^6+0.8416645931{\omega }^5\hfill \\ -0.5793147772{\omega }^4-1.7894315920{\omega }^3-1.3916904163{\omega }^2\hfill \\ +0.0896889602\omega +0.2457209648=f_7(\omega )\hfill \end{array}$

f₅(ω) = 0 gives ω = 0.24144152. The EE point is stable therefore for any value less than ω = 0.24144152.

It can be seen that the value of bifurcation of ω for the system with diffusion is greater than value of bifurcation of ω for the system without diffusion.

6. Numerical Methods

In this section, we apply two proposed and classical splitting methods to the SEIQV reaction-diffusion epidemic model with diffusion. Operator-splitting techniques very effectively handle the nonlinearity and complexity of reaction-diffusion equations. Therefore, these techniques are used frequently by several researchers for the solution of nonlinear differential equations [23, 25–33]. The SEIQV epidemic model with diffusion is split in two ways. The nonlinear reaction equations are split in the first step as,

$\begin{array}{ll}\frac{1}{2}\frac{\partial S}{\partial t}=b-\frac{\beta SI}{1+\alpha I}-(\omega +\mu +q_3)S& (39)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial E}{\partial t}=\frac{\beta SI}{1+\alpha I}-(\mu +\sigma +q_2)E& (40)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial I}{\partial t}=\sigma E-(\mu +ϵ+\gamma +q_1)I& (41)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial Q}{\partial t}=q_{3}S+q_{2}E+q_{1}I-(\mu +\varphi )Q& (42)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial V}{\partial t}=\omega S+\varphi Q+\gamma I-\mu V& (43)\end{array}$

and the diffusion equations are split in the second step as:

$\begin{array}{ll}\frac{1}{2}\frac{\partial S}{\partial t}=d_1\frac{{\partial }^{2}S}{\partial x^2}& (44)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial E}{\partial t}=d_2\frac{{\partial }^{2}E}{\partial x^2}& (45)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial I}{\partial t}=d_3\frac{{\partial }^{2}I}{\partial x^2}& (46)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial Q}{\partial t}=d_4\frac{{\partial }^{2}Q}{\partial x^2}& (47)\end{array}$

$\begin{array}{ll}\frac{1}{2}\frac{\partial V}{\partial t}=d_5\frac{{\partial }^{2}V}{\partial x^2}& (48)\end{array}$

Now, we apply forward and backward Euler methods with operator splitting on the system (6)–(7).

$\begin{array}{ll}{\overline{S}}_i^{m+\frac{1}{2}}=S_i^m+\tau \left(b-\frac{\beta S_i^m{I}_i^m}{1+\alpha I_i^m}-(\omega +\mu +q_3)S_i^m\right)& (49)\end{array}$

$\begin{array}{ll}{\overline{E}}_i^{m+\frac{1}{2}}=E_i^m+\tau \left(\frac{\beta S_i^m{I}_i^m}{1+\alpha I_i^m}-(\mu +\sigma +q_2)E_i^m\right)& (50)\end{array}$

$\begin{array}{ll}{\overline{I}}_i^{m+\frac{1}{2}}=I_i^m+\tau \left(\sigma E_i^m-(\mu +ϵ+\gamma +q_1)I_i^m\right)& (51)\end{array}$

$\begin{array}{ll}{\overline{Q}}_i^{m+\frac{1}{2}}=Q_i^m+\tau (q_3{S}_i^m+q_2{E}_i^m+q_1{I}_i^m-(\mu +\varphi )Q_i^m)& (52)\end{array}$

$\begin{array}{ll}{\overline{V}}_i^{m+\frac{1}{2}}=V_i^m+\tau (\omega S_i^m+\varphi Q_i^m+\gamma I_i^m-\mu V_i^m)& (53)\end{array}$

where $S_i^m,E_i^m,I_i^m,Q_i^m$ and $V_i^m$ at mτ, m = 0, 1, … and 0 + ih, i = 0, 1, … reflects difference approximations of S, E, I, Q, and V. The values of ${\overline{S}}_i^{m+\frac{1}{2}}$, ${\overline{E}}_i^{m+\frac{1}{2}}$, ${\overline{I}}_i^{m+\frac{1}{2}}$, ${\overline{Q}}_i^{m+\frac{1}{2}}$, and ${\overline{V}}_i^{m+\frac{1}{2}}$ are the values at the half time step. Both forward and backward Euler methods have same process at first step, but, at the second half step of time, they have different procedures. Since the forward Euler operator splitting method is explicit, we use:

$\begin{array}{ll}{S}_i^{m+1}={\overline{S}}_i^{m+\frac{1}{2}}+{\lambda }_1\left({\overline{S}}_{i-1}^{m+\frac{1}{2}}-2{\overline{S}}_i^{m+\frac{1}{2}}+{\overline{S}}_{i+1}^{m+\frac{1}{2}}\right)& (54)\end{array}$

$\begin{array}{ll}{E}_i^{m+1}={\overline{E}}_i^{m+\frac{1}{2}}+{\lambda }_2\left({\overline{E}}_{i-1}^{m+\frac{1}{2}}-2{\overline{E}}_i^{m+\frac{1}{2}}+{\overline{E}}_{i+1}^{m+\frac{1}{2}}\right)& (55)\end{array}$

$\begin{array}{ll}{I}_i^{m+1}={\overline{I}}_i^{m+\frac{1}{2}}+{\lambda }_3\left({\overline{I}}_{i-1}^{m+\frac{1}{2}}-2{\overline{I}}_i^{m+\frac{1}{2}}+{\overline{I}}_{i+1}^{m+\frac{1}{2}}\right)& (56)\end{array}$

$\begin{array}{ll}{Q}_i^{m+1}={\overline{Q}}_i^{m+\frac{1}{2}}+{\lambda }_4\left({\overline{Q}}_{i-1}^{m+\frac{1}{2}}-2{\overline{Q}}_i^{m+\frac{1}{2}}+{\overline{Q}}_{i+1}^{m+\frac{1}{2}}\right)& (57)\end{array}$

$\begin{array}{ll}{V}_i^{m+1}={\overline{V}}_i^{m+\frac{1}{2}}+{\lambda }_5\left({\overline{V}}_{i-1}^{m+\frac{1}{2}}-2{\overline{V}}_i^{m+\frac{1}{2}}+{\overline{V}}_{i+1}^{m+\frac{1}{2}}\right)& (58)\end{array}$

For the backward Euler method, we use:

$\begin{array}{ll}-{\lambda }_1{S}_{i-1}^{m+1}+(1+2{\lambda }_1)S_i^{m+1}-{\lambda }_1{S}_{i-1}^{m+1}={\overline{S}}_i^{m+\frac{1}{2}}& (59)\end{array}$

$\begin{array}{ll}-{\lambda }_2{E}_{i-1}^{m+1}+(1+2{\lambda }_2)E_i^{m+1}-{\lambda }_2{E}_{i-1}^{m+1}={\overline{E}}_i^{m+\frac{1}{2}}& (60)\end{array}$

$\begin{array}{ll}-{\lambda }_3{I}_{i-1}^{m+1}+(1+2{\lambda }_3)I_i^{m+1}-{\lambda }_3{I}_{i-1}^{m+1}={\overline{I}}_i^{m+\frac{1}{2}}& (61)\end{array}$

$\begin{array}{ll}-{\lambda }_4{Q}_{i-1}^{m+1}+(1+2{\lambda }_4)Q_i^{m+1}-{\lambda }_4{Q}_{i-1}^{m+1}={\overline{Q}}_i^{m+\frac{1}{2}}& (62)\end{array}$

$\begin{array}{ll}-{\lambda }_5{V}_{i-1}^{m+1}+(1+2{\lambda }_5)V_i^{m+1}-{\lambda }_5{V}_{i-1}^{m+1}={\overline{V}}_i^{m+\frac{1}{2}}& (63)\end{array}$

For the proposed NSFD operator splitting methods, we implement the rules constructed by Mickens [8]. The technique for both explicit and implicit schemes is similar at the first half time step:

$\begin{array}{ll}{\overline{S}}_i^{m+\frac{1}{2}}=\frac{S_i^m+\tau b}{1+\frac{\tau \beta I_i^m}{1+\alpha I_i^m}+\tau (\omega +\mu +q_3)}& (64)\end{array}$

$\begin{array}{ll}{\overline{E}}_i^{m+\frac{1}{2}}=\frac{E_i^m+\frac{\tau \beta S_i^m{I}_i^m}{1+\alpha I_i^m}}{1+\tau (\mu +\sigma +q_2)}& (65)\end{array}$

$\begin{array}{ll}{\overline{I}}_i^{m+\frac{1}{2}}=\frac{I_i^m+\tau \sigma E_i^m}{1+\tau (\mu +ϵ+\gamma +q_1)}& (66)\end{array}$

$\begin{array}{ll}{\overline{Q}}_i^{m+\frac{1}{2}}=\frac{Q_i^m+\tau (q_3{S}_i^m+q_2{E}_i^m+q_1{I}_i^m)}{1+\tau (\mu +\varphi )}& (67)\end{array}$

$\begin{array}{ll}{\overline{V}}_i^{m+\frac{1}{2}}=\frac{V_i^m+\tau (\omega S_i^m+\varphi Q_i^m+\gamma I_i^m)}{1+\tau \mu }& (68)\end{array}$

A positive solution desires that if:

$\begin{array}{l}{S}_i^m\ge 0,E_i^m\ge 0,I_i^m\ge 0,Q_i^m\ge 0,V_i^m\ge 0\\ \Rightarrow {\overline{S}}_i^{m+\frac{1}{2}}\ge 0,{\overline{E}}_i^{m+\frac{1}{2}}\ge 0,{\overline{I}}_i^{m+\frac{1}{2}}\ge 0,{\overline{Q}}_i^{m+\frac{1}{2}}\ge 0,{\overline{V}}_i^{m+\frac{1}{2}}\ge 0& (69)\end{array}$

The techniques for the implicit and explicit NSFD schemes are not similar for the second half of the time step. The procedure for the explicit NSFD scheme is as follows:

$\begin{array}{ll}{S}_i^{m+1}=(1-2{\lambda }_1){\overline{S}}_i^{m+\frac{1}{2}}+{\lambda }_1\left({\overline{S}}_{i-1}^{m+\frac{1}{2}}+{\overline{S}}_{i+1}^{m+\frac{1}{2}}\right)& (70)\end{array}$

$\begin{array}{ll}{E}_i^{m+1}=(1-2{\lambda }_2){\overline{E}}_i^{m+\frac{1}{2}}+{\lambda }_2\left({\overline{E}}_{i-1}^{m+\frac{1}{2}}+{\overline{E}}_{i+1}^{m+\frac{1}{2}}\right)& (71)\end{array}$

$\begin{array}{ll}{I}_i^{m+1}=(1-2{\lambda }_3){\overline{I}}_i^{m+\frac{1}{2}}+{\lambda }_3\left({\overline{I}}_{i-1}^{m+\frac{1}{2}}+{\overline{I}}_{i+1}^{m+\frac{1}{2}}\right)& (72)\end{array}$

$\begin{array}{ll}{Q}_i^{m+1}=(1-2{\lambda }_4){\overline{Q}}_i^{m+\frac{1}{2}}+{\lambda }_4\left({\overline{Q}}_{i-1}^{m+\frac{1}{2}}+{\overline{Q}}_{i+1}^{m+\frac{1}{2}}\right)& (73)\end{array}$

$\begin{array}{ll}{V}_i^{m+1}=(1-2{\lambda }_5){\overline{V}}_i^{m+\frac{1}{2}}+{\lambda }_5\left({\overline{V}}_{i-1}^{m+\frac{1}{2}}+{\overline{V}}_{i+1}^{m+\frac{1}{2}}\right)& (74)\end{array}$

We use an implicit procedure for the second NSFD method at the second half of the time step:

$\begin{array}{ll}-{\lambda }_1{S}_{i-1}^{m+1}+(1+2{\lambda }_1)S_i^{m+1}-{\lambda }_1{S}_{i-1}^{m+1}={\overline{S}}_i^{m+\frac{1}{2}}& (75)\end{array}$

$\begin{array}{ll}-{\lambda }_2{E}_{i-1}^{m+1}+(1+2{\lambda }_2)E_i^{m+1}-{\lambda }_2{E}_{i-1}^{m+1}={\overline{E}}_i^{m+\frac{1}{2}}& (76)\end{array}$

$\begin{array}{ll}-{\lambda }_3{I}_{i-1}^{m+1}+(1+2{\lambda }_3)I_i^{m+1}-{\lambda }_3{I}_{i-1}^{m+1}={\overline{I}}_i^{m+\frac{1}{2}}& (77)\end{array}$

$\begin{array}{ll}-{\lambda }_4{Q}_{i-1}^{m+1}+(1+2{\lambda }_4)Q_i^{m+1}-{\lambda }_4{I}_{i-1}^{m+1}={\overline{Q}}_i^{m+\frac{1}{2}}& (78)\end{array}$

$\begin{array}{ll}-{\lambda }_5{V}_{i-1}^{m+1}+(1+2{\lambda }_5)V_i^{m+1}-{\lambda }_5{V}_{i-1}^{m+1}={\overline{V}}_i^{m+\frac{1}{2}}& (79)\end{array}$

where,

$\begin{array}{l}{\lambda }_1=d_1\frac{\tau }{h^2},{\lambda }_2=d_2\frac{\tau }{h^2},{\lambda }_3=d_3\frac{\tau }{h^2},{\lambda }_4=d_4\frac{\tau }{h^2},{\lambda }_5=d_5\frac{\tau }{h^2}\end{array}$

6.1. Stability and Accuracy of Splitting Schemes

In finite difference operator splitting techniques, the step involving the reaction term is unconditionally stable because it is solved exactly [25, 26]. On the other hand, the step involving the diffusion term has different stability in different techniques. The explicit procedure has conditional stability in the region:

$\begin{array}{ll}{\lambda }_i\le \frac{1}{2},(i=1,2,3,4,5).& (80)\end{array}$

while the implicit procedure has unconditionally stability [25, 26]. The accuracy of both schemes is O(τ) and O(h²) for all the methods under study.

6.2. Positivity of Proposed Schemes

Equations (64)–(65) in the reaction step of both proposed methods preserve the property of positivity depicted by the continuous SEIQV model, as there is no negative term involved in (64)–(65).

As far as the diffusion step is concerned, the proposed explicit scheme (70)–(74) demonstrates the positivity of the solution if:

$\begin{array}{l}1-2{\lambda }_i\ge 0,i=1,2,3,4,5\end{array}$

so,

$\begin{array}{l}{\lambda }_i\le \frac{1}{2},(i=1,2,3,4,5)\end{array}$

which is the condition of stability for the explicit operator-splitting NSFD scheme (70)–(74). This verifies that the explicit NSFD scheme retains the positive solution in the region of stability. M matrix theory has been used for the verification of the positivity of the implicit NSFD method (75)–(79). For more details [34] is referred.

6.2.1. Theorem [21, 22]

For any positive τ and h, the system described by (75)–(79) is also positive, i.e., S^m > 0, E^m > 0, Q^m > 0 and V^m > 0, ∀m = 0, 1, 2…

Proof

Equations (75)–(79) can be written as:

$\begin{array}{ll}A{S}^{m+1}=S^m& (81)\end{array}$

$\begin{array}{ll}B{E}^{m+1}=E^m& (82)\end{array}$

$\begin{array}{ll}C{I}^{m+1}=I^m& (83)\end{array}$

$\begin{array}{ll}D{Q}^{m+1}=Q^m& (84)\end{array}$

$\begin{array}{ll}G{V}^{m+1}=V^m& (85)\end{array}$

In Equations (81)–(85), the letters A, B, C, D, and G represent the square matrices. Where,

$\begin{array}{ll}A=\left(\begin{array}{cccccccc}{a}_3& a_1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ a_2& a_3& a_2& \ddots & & & & ⋮\\ 0& a_2& a_3& a_2& \ddots & & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & a_2& a_3& a_2& 0\\ ⋮& & & & \ddots & a_2& a_3& a_2\\ 0& \cdots & \cdots & \cdots & \cdots & 0& a_1& a_3\end{array}\right)& (86)\end{array}$

$\begin{array}{ll}B=\left(\begin{array}{cccccccc}{b}_3& b_1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ b_2& b_3& b_2& \ddots & & & & ⋮\\ 0& b_2& b_3& b_2& \ddots & & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & b_2& b_3& b_2& 0\\ ⋮& & & & \ddots & b_2& b_3& b_2\\ 0& \cdots & \cdots & \cdots & \cdots & 0& b_1& b_3\end{array}\right)& (87)\end{array}$

$\begin{array}{ll}C=\left(\begin{array}{cccccccc}{c}_3& c_1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ c_2& c_3& c_2& \ddots & & & & ⋮\\ 0& c_2& c_3& c_2& \ddots & & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & c_2& c_3& c_2& 0\\ ⋮& & & & \ddots & c_2& c_3& c_2\\ 0& \cdots & \cdots & \cdots & \cdots & 0& c_1& c_3\end{array}\right)& (88)\end{array}$

$\begin{array}{ll}D=\left(\begin{array}{cccccccc}{d}_3& d_1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ d_2& d_3& d_2& \ddots & & & & ⋮\\ 0& d_2& d_3& d_2& \ddots & & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & d_2& d_3& d_2& 0\\ ⋮& & & & \ddots & d_2& d_3& d_2\\ 0& \cdots & \cdots & \cdots & \cdots & 0& d_1& d_3\end{array}\right)& (89)\end{array}$

$\begin{array}{ll}G=\left(\begin{array}{cccccccc}{g}_3& g_1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ g_2& g_3& g_2& \ddots & & & & ⋮\\ 0& g_2& g_3& g_2& \ddots & & & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& & & \ddots & g_2& g_3& g_2& 0\\ ⋮& & & & \ddots & g_2& g_3& g_2\\ 0& \cdots & \cdots & \cdots & \cdots & 0& g_1& g_3\end{array}\right)& (90)\end{array}$

The off-diagonal entries of A are a₁ = −2λ₁, a₂ = −λ₁, and the diagonal entries are a₃ = 1 + 2λ₁. The entries of B in the off-diagonal are b₁ = −2λ₁, b₂ = −λ₁, and the diagonal entries are b₃ = 1 + 2λ₂. The entries of C in the off-diagonal are c₁ = −2λ₃, c₂ = −λ₃, and the diagonal entries are c₃ = 1 + 2λ₃. The off-diagonal entries of D are d₁ = −2λ₄, d₂ = −λ₄, and the diagonal entries are d₃ = 1 + 2λ₄. The off-diagonal entries of G are g₁ = −2λ₅, g₂ = −λ₅, and the diagonal entries are g₃ = 1 + 2λ₅. Thus, A, B, C, D, and G are M-matrices, and Equations (81), (82), (83), (84), and (85) are:

$\begin{array}{ll}{S}^{m+1}=A^{-1}{S}^m& (91)\end{array}$

$\begin{array}{ll}{E}^{m+1}=B^{-1}{E}^m& (92)\end{array}$

$\begin{array}{ll}{I}^{m+1}=C^{-1}{I}^m& (93)\end{array}$

$\begin{array}{ll}{Q}^{m+1}=D^{-1}{I}^m& (94)\end{array}$

$\begin{array}{ll}{V}^{m+1}=G^{-1}{I}^m& (95)\end{array}$

If we consider that S^m > 0, E^m > 0, I^m > 0, Q^m > 0, and V^m > 0, then the M-matrix along with (60) implies that the values of all of the state variables, i.e., S^m+1, E^m+1, I^m+1, Q^m+10, and V^m+1 are positive. Hence, the theorem is done by using the principle of mathematical induction.

This theorem is applied for drawing the conclusion that the proposed scheme, which is implicit in nature, guarantees the positive solution unconditionally.

7. Numerical Experiment and Simulations

A numerical test is performed on both the points of equilibrium for all the schemes under consideration. The set of parametric values considered for the test problem at disease-free equilibrium point [7] is given as:

b = 0.7, μ = 0.06, σ = 0.7, β = 0.35, γ = 0.15, ω = 0.3, q₁ = 0.2, q₂ = 0.2, q₃ = 0.1, ϕ = 0.4, ϵ = 0.05, α = 2, d₁ = 0.05, d₂ = 0.01, d₃ = 0.001, d₄ = 0.01, d₅ = 0.01.

For the endemic equilibrium point, the following parametric values are used:

b = 0.7, μ = 0.06, σ = 0.7, β = 0.35, γ = 0.15, ω = 0.06, q₁ = 0.2, q₂ = 0.2,

q₃ = 0.1, ϕ = 0.4, ϵ = 0.05, α = 2, d₁ = 0.05, d₂ = 0.01, d₃ = 0.001, d₄ = 0.01, d₅ = 0.01.

The initial condition for the model (6)–(10) is:

$\begin{array}{ll}S(x,0)=\{\begin{array}{cc}0.7x\hfill & \text{if }x\in [0,0.5)\\ 0.7(10-x)& \text{if }x\in \left[0.5,1\right]\end{array}& (96)\end{array}$

$\begin{array}{ll}E(x,0)=\{\begin{array}{cc}0.5x\hfill & \text{if }x\in [0,0.5)\\ 0.5(10-x)& \text{if }x\in \left[0.5,1\right]\end{array}& (97)\end{array}$

$\begin{array}{ll}I(x,0)=\{\begin{array}{cc}0.3x\hfill & \text{if }x\in [0,0.5)\\ 0.3(10-x)& \text{if }x\in \left[0.5,1\right]\end{array}& (98)\end{array}$

$\begin{array}{ll}Q(x,0)=\{\begin{array}{cc}0.1x\hfill & \text{if }x\in [0,0.5)\\ 0.1(10-x)& \text{if }x\in \left[0.5,1\right]\end{array}& (99)\end{array}$

$\begin{array}{ll}V(x,0)=\{\begin{array}{cc}0.1x\hfill & \text{if }x\in [0,0.5)\\ 0.1(10-x)& \text{if }x\in \left[0.5,1\right]\end{array}& (100)\end{array}$

7.1. Disease-Free Equilibrium Point

In this section, graphs of all the state variables against time are presented (for DFE) to illustrate the behavior of the schemes. In Figures 1, 2, we consider h = 0.5, λ₁ = 0.3, λ₂ = 0.06, λ₃ = 0.006, λ₄ = 0.06, and λ₅ = 0.06.

FIGURE 1

Figure 1. The explicit operator splitting NSFD scheme is used to simulate the graphs (A–E). (A) Mesh graph of S; (B) Mesh graph of E; (C) Mesh graph of I; (D) Mesh graph of Q; (E) Mesh graph of V.

FIGURE 2

Figure 2. The implicit operator splitting NSFD scheme is used to simulate graphs (A–E). (A) Mesh graph of S; (B) Mesh graph of E; (C) Mesh graph of I; (D) Mesh graph of Q; (E) Mesh graph of V.

Figures 1, 2 validate the preservation of the positivity property in both of the proposed operator splitting NSFD schemes. Also, all the graphs in Figures 1, 2 show that both proposed NSFD schemes achieve convergence to the DFE point. Next, we examine the behavior of forward Euler and backward Euler splitting schemes at different values of h and τ.

In parts (a) and (b) of Figure 3, we take the same values of h and τ as in Figures 1, 2 for the graphs of forward Euler operator splitting scheme and backward Euler operator splitting scheme. The graphs clearly show the failure of the positivity property of both classic schemes. In parts (c) and (d) of Figure 3, both existing splitting schemes converge to the false DFE equilibrium point for susceptible individuals.

FIGURE 3

Figure 3. (A) The forward Euler FD operator splitting method is used to simulate the graph of exposed persons for the DFE point at h = 0.5, λ₁ = 0.3. (B) The backward Euler FD operator splitting method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.3. (C) The forward Euler FD operator splitting method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4. (D) The backward Euler FD operator splitting method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4.

In parts (a) and (b) of Figure 4, we take the same values of h and τ as given in parts (c) and (d) of Figure 3 for the explicit and implicit NSFD operator splitting schemes, respectively. The graphs clearly show that both of the proposed NSFD schemes not only preserve the positivity property but also achieve convergence to the true equilibrium point.

FIGURE 4

Figure 4. (A) The explicit operator splitting NSFD method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.3. (B) The implicit operator splitting NSFD method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4.

7.2. Endemic Equilibrium Point

In this section, we present simulations of the SEIQV epidemic model at the EE point using all of the operator splitting FD schemes. In Figures 5, 6, we consider h = 0.5, λ₁ = 0.3, λ₂ = 0.06, λ₃ = 0.006, λ₄ = 0.06, and λ₅ = 0.06.

FIGURE 5

Figure 5. (A) The explicit operator splitting NSFD scheme is used to simulate graphs (A–E). (A) Mesh graph of S; (B) Mesh graph of E; (C) Mesh graph of I; (D) Mesh graph of Q; (E) Mesh graph of V.

FIGURE 6

Figure 6. The implicit operator splitting NSFD scheme is used to simulate graphs (A–E). (A) Mesh graph of S; (B) Mesh graph of E; (C) Mesh graph of I; (D) Mesh graph of Q; (E) Mesh graph of V.

Figures 5, 6 depict the graphs of susceptible, exposed, infected, quarantined, and vaccinated individuals for the EE point using the explicit operator splitting NSFD scheme and implicit operator splitting NSFD scheme, respectively. All the graphs in Figures 5, 6 demonstrate that both of the proposed operator splitting NSFD schemes preserve the property of positivity. These graphs also show that both proposed schemes converge to the EE point.

Again, both the forward and backward Euler FD schemes fail to preserve the positivity property and converge to the false EE point, as shown in Figure 7.

FIGURE 7

Figure 7. (A) The forward Euler FD operator splitting method is used to simulate the graph of exposed persons for the DFE point at h = 0.5, λ₁ = 0.3. (B) The backward Euler FD operator splitting method is used to simulate the graph of exposed persons for the DFE point at h = 0.5, λ₁ = 0.3. (C) The forward Euler FD operator splitting method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4. (D) The backward Euler FD operator splitting method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4.

Figure 8 shows that the proposed NSFD operator splitting methods are consistent with the continuous reaction-diffusion system as they not only preserve the positivity property but converge to the EE point.

FIGURE 8

Figure 8. (A) The explicit operator splitting NSFD method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4. (B) The implicit operator splitting NSFD method is used to simulate the graph of susceptible persons for the DFE point at h = 0.5, λ₁ = 0.4.

8. Conclusion

In this work, we consider the SEIQV reaction-diffusion epidemic model. The stability of the SEIQV model is guaranteed numerically by using criteria defined by Routh-Hurwitz. We also find the bifurcation value of the important vaccination parameter ω of SEIQV epidemic systems with diffusion and without diffusion. We design two novel and efficient operator splitting NSFD schemes for the SEIQV reaction-diffusion system. The NSFD schemes put forth, which are technically operator splitting schemes, possess the same behavior as is possessed by the SEIQV epidemic system. To conclude regarding the designed methods, we present two novel numerical schemes, one of which is explicit and the other of which is implicit in nature. The explicit scheme is more computationally efficient than the implicit scheme, but it has conditional stability while the implicit scheme is stable unconditionally. Both schemes employ structural splitting, due to which they deal adroitly with the nonlinearity of the reaction-diffusion system. These schemes avoid the false chaos that is a part of many existing methods. The positive solution of the SEIQV model is sustained by both schemes. Also, the nature of the stability of equilibria is preserved effectively by the proposed NSFD schemes. It is also shown that classical schemes, in parallel to our proposed schemes, produce chaos, leading to inconsistencies and instabilities numerically. The currently designed schemes are a valuable contribution for finding the solutions of nonlinear dynamical systems comprising differential equations. These NSFD schemes will become very efficient for the solution of one- and multi-dimensional reaction-diffusion population models, auto-catalytic chemical reaction models, and many more.

Author Contributions

NA, MR, and MAR designed the study. NA, MF, and MR developed the methodology. NA and MF collected the data. NA, DB, KN, and IK performed the analysis. NA, MF, DB, MA, and IK wrote the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors are grateful to the referees for their suggestions and useful comments on this paper.

References

1. Ahmed N, Wei Z, Baleanu D, Rafiq M, Rehman MA. Spatio-temporal numerical modeling of reaction-diffusion measles epidemic system. Chaos. (2019) 29:103101. doi: 10.1063/1.5116807

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Abro KA, Khan I, Nisar KS. Novel technique of Atangana and Baleanu for heat dissipation in transmission line of electrical circuit. Chaos Solitons Fractals. (2019) 129:40–5 doi: 10.1016/j.chaos.2019.08.001

CrossRef Full Text | Google Scholar

3. Ahmed N, Rafiq M, Baleanu D, Rehman MA. Spatio-temporal numerical modeling of auto-catalytic Brusselator model. Roman J Phys. (2019) 64:110.

Google Scholar

4. Gill V, Singh J and Singh Y. Analytical solution of generalized space-time fractional advection-dispersion equation via coupling of sumudu and fourier transforms. Front Phys. (2019) 6:151. doi: 10.3389/fphy.2018.00151

CrossRef Full Text | Google Scholar

5. Kumar D, Singh J, Al Qurashi M, Baleanu D. A new fractional SIRS-SI malaria disease model with application of vaccines, antimalarial drugs, and spraying. Adv Differ Equ. (2019) 278. doi: 10.1186/s13662-019-2199-9

CrossRef Full Text | Google Scholar

6. Saqib M, Khan I, Shafie S. Application of fractional differential equations to heat transfer in hybrid nanofluid: modeling and solution via integral transforms. Adv Differ Equ. (2019) 52: doi: 10.1186/s13662-019-1988-5

CrossRef Full Text | Google Scholar

7. Liu X, Yang L. Stability analysis of an SEIQV epidemic model with saturated incidence rate. Nonlinear Anal Real World Appl. (2012) 13:2671–9. doi: 10.1016/j.nonrwa.2012.03.010

CrossRef Full Text | Google Scholar

8. Mickens RE. Nonstandard finite difference models of differential equations world scientific (1994).

Google Scholar

9. Ahmed N, Rafiq M, Rehman MA, Iqbal MS, Ali M. Numerical modelling of three dimensional Brusselator reaction diffusion system. AIP Adv. (2019) 9:015205. doi: 10.1063/1.5070093

CrossRef Full Text | Google Scholar

10. Mickens RE. Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations. J Differ Equ Appl. (2005) 11:645–53. doi: 10.1080/10236190412331334527.

CrossRef Full Text | Google Scholar

11. Macias-Diaz JE, Puri A. An explicit positivity-preserving finite-difference scheme for the classical Fisher-Kolmogorov-Petrovsky-Piscounov equation. Appl Math Comput. (2012) 218:5829–37. doi: 10.1016/j.amc.2011.11.064

CrossRef Full Text | Google Scholar

12. Mickens RE. A nonstandard finite difference scheme for a Fisher PDE having nonlinear diffusion. Comput Math Appl. (2003) 45:429–36. doi: 10.1016/S0898-1221(03)80028-7

CrossRef Full Text | Google Scholar

13. Mickens RE. A nonstandard finite difference scheme for an advection-reaction equation J Diff Eq Appl. (2004) 10:307–1312. doi: 10.1080/10236190410001652838

CrossRef Full Text | Google Scholar

14. Ahmed N, Rafiq M, Rehman MA, Ali M, Ahmad MO. Numerical modeling of SEIR measles dynamics with diffusion. Commun Math Appl. (2018) 9:315–26.

Google Scholar

15. Ahmed N, Shahid N, Iqbal Z, Jawaz M, Rafiq M, Tahira SS, et al. Numerical modeling of SEIQV epidemic model with saturated incidence rate. J Appl Environ Biol Sci. (2018) 8:67–82.

Google Scholar

16. Ahmed N, Jawaz M, Rafiq M, Rehman MA, Ali M, Ahmad MO. Numerical treatment of an epidemic model with spatial diffusion. J Appl Environ Biol Sci. (2018) 8:17–29.

Google Scholar

17. Chinviriyasit S, Chinviriyasit W. Numerical modeling of SIR epidemic model with diffusion. Appl Math Comp. (2010) 216:395–409. doi: 10.1016/j.amc.2010.01.028

CrossRef Full Text | Google Scholar

18. Fatima U, Ali M, Ahmed N, Rafiq M. Numerical modeling of susceptible latent breaking-out quarantine computer virus epidemic dynamics. Heliyon. (2018) 4:e00631. doi: 10.1016/j.heliyon.2018.e00631

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Jansen H, Twizell EH. An unconditionally convergent discretization of the SEIR model. Math Comput Simulat. (2002) 58:147–58. doi: 10.1016/S0378-4754(01)00356-1

CrossRef Full Text | Google Scholar

20. Rafiq M. Numerical modeling of infectious diseases dynamics. Ph.D. thesis, Lahore, Pakistan: University of Engineering and Technology (2016).

Google Scholar

21. Qin W, Wang L, Ding X. A non-standard finite difference method for a hepatitis B virus infection model with spatial diffusion. J Differ Equ Appl. (2014) 20:1641–51. doi: 10.1080/10236198.2014.968565

CrossRef Full Text | Google Scholar

22. Manna K, Chakrabarty SP. Global stability and a non-standard finite difference scheme for a diffusion driven HBV model with capsids. J Differ Equ Appl. (2015) 21:918–33. doi: 10.1080/10236198.2015.1056524

CrossRef Full Text | Google Scholar

23. Chakrabrty A, Singh M, Lucy B, Ridland P. Predator-prey model with prey- taxis and diffusion. Math Comput Model. (2007) 46:482–98. doi: 10.1016/j.mcm.2006.10.010

CrossRef Full Text | Google Scholar

24. Islam S, Haider N. Numerical solution of compartmental model by meshless finite difference methods. Appl Math Comput. (2014) 238:408–35. doi: 10.1016/j.amc.2014.04.014

CrossRef Full Text | Google Scholar

25. Harwood RC. (2011) Operator splitting method and applications for semilinear parabolic partial differential equations. Ph.D. dissertation. Pullman, WA: Dept. Math., Washington State Univ.

Google Scholar

26. Harwood RC, Manoranjan VS, Edwards DB. Lead-acid battery model under discharge with a fast splitting method. IEEE Trans Energy Convers. (2011) 26:1109–17. doi: 10.1109/TEC.2011.2162093

CrossRef Full Text | Google Scholar

27. Yanenko NN. The method of fractional steps. Berlin; Heidelberg: Springer-Verlag (1971).

28. Zharnitsky V. Averaging for split-step scheme. Nonlinearity. (2003) 16:1359–66. doi: 10.1088/0951-7715/16/4/310

CrossRef Full Text | Google Scholar

29. Ansarizadeh F, Singh M, Richards D. Modelling of tumor cells regression in response to chemotherapeutic treatment. Appl Math Model. (2017) 48:96–112. doi: 10.1016/j.apm.2017.03.045

CrossRef Full Text | Google Scholar

30. Naheed A, Singh M, Lucy D. Numerical study of SARS epidemic model with the inclusion of diffusion in the system. Appl Math Comput. (2014) 229:480–98. doi: 10.1016/j.amc.2013.12.062

CrossRef Full Text | Google Scholar

31. Naheed A, Singh M, Lucy D. Effect of treatment on transmission dynamics of SARS epidemic. J Infect Non Infect Dis. (2016) 2:16. doi: 10.24966/INID-8654/100016

CrossRef Full Text | Google Scholar

32. Samsuzzoha MD. A Study on numerical solutions of epidemic models. Ph.D. thesis, Australia: Swinburne University of Technology (2012).

Google Scholar

33. Wang H. Q., Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations. Appl Math Comput. (2005) 170:17–35. doi: 10.1016/j.amc.2004.10.066

CrossRef Full Text | Google Scholar

34. Fujimoto T, Ranade R. Two characterizations of inverse-positive matrices: the Hawkins-Simon condition and the Le Chatelier-Braun principle. Electron J Linear Algebra. (2004) 11:59–65. doi: 10.13001/1081-3810.1122

CrossRef Full Text | Google Scholar