The standard numerical methodology will be discussed in the proceeding section for approximation of stochastic model solutions. This is agreed by us that \({I}_{n}=\{\text{0,1},\text{2,3},\dots n\}\). Consider \(N\in {\mathbb{N}}\), as well as the consistent division of the temporal interval \([0,T]\) with uniform partition equal to \(\tau =\frac{T}{N}\), and their respective nodes are presented as \(0={t}_{0}.
For each \(n\in {I}_{N}\). Further, this will be agreed by us \({U}^{n}=U({t}_{n})\), however \(n\in {I}_{N}\) and \(U(t)={(S,I,R)}^{t}\). \(\Delta {W}_{n}=W\left({t}_{n}+1\right)-W({t}_{n})\).
Stochastic delayed model
The phrase “time delay” is used in transmission terms to denote the duration of viral infection’s incubation or the interval between infection and the start of symptoms. The addition of time delay in the model could result in many periodic solutions for varying time delay \(\tau\) values. Considered that \(cdt=cdt+\sigma dB(t)\) with parametric perturbation technique for the Eqs. (1–3) as follows28,
$$\left\{\begin{array}{c}dS\left(\text{t}\right)=\left[\text{a}-\text{cS}\left(\text{t}\right)\text{I}\left(\text{t}\right)\left(1+{\gamma {\rm I}}\left(\text{t}\right)\right){\text{e}}^{-{\mu \tau }}-{\mu {\rm S}}\left(\text{t}\right)+{\alpha {\rm R}}\left(\text{t}\right)\right]dt-\sigma c\left[\text{S}\left(\text{t}\right)\text{I}\left(\text{t}\right)\left(1+{\gamma {\rm I}}\left(\text{t}\right)\right)\right]dB(t), \\ dI\left(\text{t}\right)=[cI(t)S(t)\left(1+{\gamma {\rm I}}\left(\text{t}\right)\right){\text{e}}^{-{\mu \tau }}-(\beta +\mu +\delta -b)I(t)]dt+\sigma cS\left(\text{t}\right)I\left(\text{t}\right)\left(1+{\gamma {\rm I}}\left(\text{t}\right)\right)dB(t),\\ dR\left(\text{t}\right)=\left[{\beta {\rm I}}\left(\text{t}\right)-\left({\alpha }+\upmu \right)\text{R}\left(\text{t}\right)\right]dt, \end{array}\right.$$
(16)
In the above equations, \(\tau\) is the delay effect, \(\upsigma\) is the randomness, \(\forall \text{ t}\ge 0\) and \(\text{B}(\text{t})\) is Brownian motion. These equations are non-integrable due to Brownian motion.
Positivity and boundedness
Let’s introduce a probability space which is represented as \((\Omega , F, P)\) and \({\{{F}_{t}\}}_{t\epsilon R}\) represented as a filtration. For satisfying the right continuous conditions and increasing29, \({F}_{0}\) contains P-null sets. Symbolically represented as \(U\left(t\right)=(S\left(t\right), I\left(t\right), R\left(t\right))\) and norm as \(\left|U(t)\right|=\sqrt{{S}^{2}\left(t\right)+{I}^{2}\left(t\right)+{R}^{2}(t)}\).
Theorem 8
For model (16) and any given initial value \(\left(S\left(0\right),I\left(0\right),R\left(0\right)\right)\in {R}_{+}^{3}\), there is a unique solution \(\left(S\left(t\right),I\left(t\right),R\left(t\right)\right)\), which is defined on \(t\in [0,\infty ]\) and the solution of the system will remain in \({\mathbb{R}}_{+}^{3}\) with probability one.
Proof
The initial value \(\left(S\left(0\right),I\left(0\right),R\left(0\right)\right)\epsilon {\mathbb{R}}_{+}^{3}\) of the model (16) satisfies the local Lipschitz condition, and has a unique local solution existing on \(t\in \left[0,{t}_{e}\right],{t}_{e}\) shows the explosion time. Next to show that the solution of the system is global, \({t}_{e}=\infty\).
Suppose \({d}_{0}\ge 1\) for initial values, they all lie within the interval [\(\frac{1}{{d}_{0}},{d}_{0}]\). For each integer, \(d\ge {d}_{0}\), this sequence is known as stopping time.
$${\tau }_{d}=\text{inf}\left\{t\in \left[0,{t}_{e}\right]:\text{min}\left(S\left(t\right),I\left(t\right),R\left(t\right)\right)\le \frac{1}{d} or\text{max}\left(S\left(t\right),I\left(t\right),R\left(t\right)\right)\ge d\right\}.$$
(17)
where inf ∅ = ∞ (∅ is an empty set) \(.\)
By the definition of stopping time, \({\tau }_{d}\) is increasing function as \(d\to \infty\).
$$\text{So}, {\tau }_{\infty }={lim}_{d\to \infty }{\tau }_{d}$$
(18)
$${\tau }_{\infty }\le {\tau }_{e}$$
If \({\tau }_{\infty }=\infty\), then \({\tau }_{e}=\infty\).
If this argument is false, then there exists a constant \(T>0\) and \(\varepsilon \epsilon (\text{0,1})\) such that,
$$P\{{\tau }_{\infty }\le T\}>\varepsilon .$$
(19)
There exists an integer \({d}_{1}\ge {d}_{0}\).
$$P\left\{{\tau }_{d}\le T\right\}>\varepsilon , \forall d\ge {d}_{1}.$$
(20)
For \(t\le {\tau }_{d}\), we get
$$d\left(S+I+R\right)=\left[a-cSI\left(1+\gamma I\right){e}^{-\mu \tau }-\mu S+\alpha R+cSI\left(1+\gamma I\right){e}^{-\mu \tau }-\left(\beta +\mu +\delta -b\right)I+\beta I-\left(\alpha +\mu \right)R\right].$$
$$=\left[a-\mu \left(S+I+R\right)-\delta I+bI\right]dt.$$
$$\le \left[a-\mu \left(S+I+R\right)\right]dt.$$
After calculation, we see that
$$S\left(t\right)+I\left(t\right)+R\left(t\right)\le \left\{\begin{array}{c}\begin{array}{cc}\frac{a}{\mu } & if \,S\left(0\right)+I\left(0\right)+R\left(0\right)\le \frac{a}{\mu }\end{array}\\ \begin{array}{cc}otherwise& S\left(0\right)+I\left(0\right)+R\left(0\right)\ge \frac{a}{\mu }\end{array}\end{array}\right.$$
Define a \({C}^{3}\)-function \(V:{R}_{+}^{3}\to {R}_{+}\) by
$$V\left(S,I,R\right)=\left(S-1-lnS\right)+\left(I-1-lnI\right)+\left(R-1-lnR\right).$$
(21)
By using Ito’s formula, we can obtain
$$dV\left(S,I,R\right)=\left(1-\frac{1}{S}\right)dS+\left(1-\frac{1}{I}\right)dI+\left(1-\frac{1}{R}\right)dR+\frac{{\sigma }^{2}}{2}dt.$$
Putting values of \(dS,dI,dR.\)
$$=\left[a+2\mu +\beta +\delta +b+\frac{{\sigma }^{2}}{2}\right]dt+\sigma c\left(S-I\right)\left(1+\gamma I\right){e}^{-\mu \tau }dB\left(t\right).$$
Let, \(M=a+2\mu +\beta +\delta +b+\frac{{\sigma }^{2}}{2}\); where M is a constant. Then above equation as
$$dV\left(S,I,R\right)\le Mdt+\sigma c\left(S-I\right)\left(1+\gamma I\right){e}^{-\mu \tau }dB\left(t\right).$$
(22)
Integrating above equation from limit \(0 to {\tau }_{d}\wedge \tau\). Now,
$${\int }_{0}^{{\tau }_{d}\wedge \tau }dV\left(S,I,R\right)\le {\int }_{0}^{{\tau }_{d}\wedge \tau }MdS+{\int }_{o}^{{\tau }_{d}\wedge \tau }\sigma c\left(S-I\right)\left(1+\gamma I\right){e}^{-\mu \tau }dB\left(t\right)$$
(23)
where \({\tau }_{d}\wedge \tau =\text{min}({\tau }_{d},T)\). Now, taking expectations.
$$EV\left(S\left({\tau }_{d}\wedge \tau \right),I\left({\tau }_{d}\wedge \tau \right),R\left({\tau }_{d}\wedge \tau \right)\right)\le V\left(S\left(0\right),I\left(0\right),R\left(0\right)\right)+MT.$$
(24)
Let, \({\Omega }_{d}=\{{\tau }_{d}\le \tau \}\)
$$P=\left\{{\Omega }_{d}\right\}>\varepsilon .$$
For every, \(u \epsilon {\Omega }_{d}\). For some \(i, {v}_{i}({\tau }_{d},u)\) equals either \(d or \frac{1}{d}\) for i = 1, 2, 3.
Hence,
\(V(S\left({\tau }_{d},u\right),I\left({\tau }_{d},u\right),R\left({\tau }_{d},u\right))\) is less \(\text{min}\left\{d-1-lnd, \frac{1}{d}-1-ln\frac{1}{d}\right\},\) we obtain
$$V\left(S\left(0\right),I\left(0\right),R\left(0\right)\right)+MT\ge E\left[{I}_{{\Omega }_{d}u}V\left(S\left({\tau }_{d},u\right),I\left({\tau }_{d},u\right),R\left({\tau }_{d},u\right)\right)\right] \ge \varepsilon [(d-1-lnd)\wedge (\frac{1}{d}-1-ln\frac{1}{d} )]$$
(25)
Letting \({\varvec{d}}\to \boldsymbol{\infty },\) leads to the contradiction,
\(\infty =V\left(S\left(0\right),I\left(0\right),R\left(0\right)\right)+MT as desired.
Extinction and persistence
Consider \(B(t)\) represents the Brownian motion and \(I(t)\) be the Ito drift–diffusion process that satisfies the stochastic model as:
$$dI\left(t\right)=\mu \left(I\left(t\right),t\right)dt+\sigma \left(I\left(t\right),t\right)dB\left(t\right).$$
If \(k(I,t)\in {C}^{2}({\mathcal{R}}^{2},\mathcal{R})\) then \(k(I,t)\) is also Ito drift diffusion. So,
$$d\left(k\left(I,t\right)\right)=\frac{dk}{dt}\left(I\left(t\right),t\right)dt+\frac{dk}{dt}\left(I\left(t\right),t\right)dB\left(t\right)+\frac{1}{2}\frac{{d}^{2}k}{{dt}^{2}}\left(I\left(t\right),t\right)dB{\left(t\right)}^{2}.$$
Theorem 9
The model (16) has contained \((S(t), I(t), R(t))\) be the solution with the initial value \((S(0), I(0), R(0))\), then we have two conditions
\(\left(\text{i}\right)\quad If {\sigma }^{2}>\frac{ca{e}^{-\mu \tau }}{2\mu \left(\beta +\mu +\delta -b\right)}-\left(\beta +\mu +\delta -b\right)\) holds, then \({lim}_{t\to \infty }\frac{lnI(t)}{t}>0\).
\((ii) \quad If {\sigma }^{2}holds, then \({lim}_{t\to \infty }sup\frac{lnI(t)}{t}\le 0\).
Proof
Consider \(the I\) component of SDDE’s as:
$$\text{dI}\left(\text{t}\right)=[\text{cI}(\text{t})\text{S}(\text{t})\left(1+{\gamma{\rm I}}\left(\text{t}\right)\right){\text{e}}^{-{\mu \tau }}-(\upbeta +\upmu +\updelta -\text{b})\text{I}(\text{t})]\text{dt}+{\sigma {\rm cS}}\left(\text{t}\right)\text{I}\left(\text{t}\right)\left(1+{\gamma {\rm I}}\left(\text{t}\right)\right)\text{dB}(\text{t})$$
Applying Ito’s formula, we get
$$dln\left(I\right)={k}^{{{\prime}}}\left(I\right)dI+\frac{1}{2}{k}^{{{\prime\prime}}}\left(I\right){I}^{2}{\sigma }^{2}dt.$$
$$dln\left(I\right)=\left(cS\left(1+\gamma I\right){e}^{-\mu \tau }-\left(\beta +\mu +\delta -b\right)-\frac{1}{2}{\sigma }^{2}\right)dt+\sigma cS\left(1+\gamma I\right){e}^{-\mu \tau }dB.$$
Integrate with respect to \(I\) on both sides from \(0 to t,\) we get
$$\text{ln}\left(I\right)=\text{ln}\left(0\right)+{\int }_{0}^{t}\left(cS\left(1+\gamma I\right){e}^{-\mu \tau }-\left(\beta +\mu +\delta -b\right)-\frac{1}{2}{\sigma }^{2}\right)dt+{\int }_{0}^{t}\sigma cS\left(1+\gamma I\right){e}^{-\mu \tau }dB.$$
\(If {\sigma }^{2}>\frac{ca{e}^{-\mu \tau }}{2\mu \left(\beta +\mu +\delta -b\right)}-\left(\beta +\mu +\delta -b\right)\), then
$$\frac{\text{ln}(I)}{t}>\left(\frac{ca{e}^{-\mu \tau }}{2\mu \left(\beta +\mu +\delta -b\right)}-\left(\beta +\mu +\delta -b\right)\right)+\frac{\text{M}\left(t\right)}{t}+\frac{\text{lnI}(0)}{t}$$
$${lim}_{t\to \infty }\frac{\text{ln}(I)}{t}>\left(\frac{ca{e}^{-\mu \tau }}{2\mu \left(\beta +\mu +\delta -b\right)}-\left(\beta +\mu +\delta -b\right)\right)>0$$
\(If {\sigma }^{2} , then
$${lim}_{t\to \infty }sup\frac{\text{ln}(I)}{t}
Therefore, \({R}_{0}^{s} \({R}_{0}^{s}={R}_{0}^{d}-\frac{{\sigma }^{2}}{2\left(\beta +\mu +\delta -b\right)}
$${lim}_{t\to \infty }sup\frac{I(t)}{t}\le 0.$$
\({lim}_{t\to \infty }I\left(t\right)=0\) as desired.
Numerical analysis
In this section, we will analyze some standard and non-standard numerical methodologies of the solution of the stochastic model and present some numerical results. It is important to note that this method has the benefits of offering limited computational cost and extremely precise estimates. Furthermore, it is important to note that applying the non-standard method may ensure the preservation of the number of intriguing aspects of appropriate solutions.
Stochastic non-standard computational method
Our parametric perturbation model of Eq. (16) can be expressed in the non-standard computational method which is named as stochastic non-standard finite difference method.
$$\frac{dS}{dt}=\left[a-cSI\left(1+\gamma I\right){e}^{-\mu \tau }-\mu S+\alpha R\right]-\sigma c[SI\left(1+\gamma I\right]\frac{dB}{dt}.$$
(26)
Equation (16) decomposed in the stochastic non-standard finite difference method as shown:
$$\frac{{S}^{n+1}-{S}^{n}}{h}=\left[a-c{S}^{n+1}{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }-\mu {S}^{n+1}+\alpha {R}^{n}\right]dt-\sigma c\left[{S}^{n+1}{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }\right]\Delta {B}_{n}.$$
(27)
Same as Eq. (27), which is in the stochastic NSFD process, we write Eq. (16) in this way accordingly.
$${S}^{n+1}=\frac{{S}^{n}+ha+h\alpha {R}^{n}}{1+hc{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }+h\mu +\sigma hc{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }\Delta {B}_{n}},$$
(28)
$${I}^{n+1}=\frac{{I}^{n}+hc{{S}^{n}I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }+hb{I}^{n}+\sigma hc{{S}^{n}I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }\Delta {B}_{n}}{1+h(\beta +\mu +\delta -b)},$$
(29)
$${R}^{n+1}=\frac{{R}^{n}+h\beta {I}^{n}}{1+h(\alpha +\mu )},$$
(30)
where, \(n=\text{0,1},2,\dots\) and \(\Delta {\text{B}}_{\text{n}}=\Delta {\text{B}}_{{t}_{n+1}}-\Delta {\text{B}}_{{t}_{n}}\) is standardized distribution. i.e. \(\Delta {\text{B}}_{\text{n}}\sim \text{N}(0, 1\)).
Convergence analysis
In this section, we present some theorems related to convergence analysis.
Theorem 10
The region \(\Gamma =\{{S}^{n}\ge 0,{I}^{n}\ge 0,{R}^{n}\ge 0; {S}^{n}+{I}^{n}+{R}^{n}\le \frac{a}{\mu }\}\) is feasible for Eq. (30–32), \(\forall n\ge 0\), where n is the positive function.
Proof
The Eq. (27–30) decayed as shown below.
$$\frac{{S}^{n+1}-{S}^{n}}{\text{h}}=(\left[a-c{{S}^{n}I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }-\mu {S}^{n}+\alpha {R}^{n}\right]dt-\sigma c[{S}^{n}{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }]\Delta {B}_{n}$$
Same as for \(I \,and \,R\),
$$\frac{{I}^{n+1}-{I}^{n}}{\text{h}}=(\left[c{{I}^{n}S}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }-(\beta +\mu +\delta -b){I}^{n}\right]dt-\sigma c[{S}^{n}{I}^{n}\left(1+\gamma {I}^{n}\right){e}^{-\mu \tau }]\Delta {B}_{n}$$
$$\frac{{R}^{n+1}-{R}^{n}}{\text{h}}=(\beta {I}^{n}-\left(\alpha +\mu \right){R}^{n}).$$
These equations are added up then we get,
$$\frac{{(S}^{n+1}+{I}^{n+1}+{R}^{n+1})-({S}^{n}+ {I}^{n}+{R}^{n})}{h}=a-\mu \left({S}^{n}+ {I}^{n}+{R}^{n}\right)-(\delta -b){I}^{n}$$
$$\frac{{(S}^{n+1}+{I}^{n+1}+{R}^{n+1})-({S}^{n}+ {I}^{n}+{R}^{n})}{h}\le a-\mu \left({S}^{n}+ {I}^{n}+{R}^{n}\right)$$
$${(S}^{n+1}+{I}^{n+1}+{R}^{n+1})\le \frac{a}{\mu }.$$
The result is bounded \(\forall n\ge 0\).
Theorem 11
(for stability) Our computational analysis method is stable if eigenvalues lie in the unit circle for any \(n\ge 0\) and \(\Delta {B}_{n}=0.\)30
Proof
Let, \(A=\frac{S+ha+h\alpha R}{1+hcSI\left(1+\gamma I\right){e}^{-\mu \tau }+k\mu }, B=\frac{1+hcSI(1+\gamma I){e}^{-\mu \tau }}{1+h\left(\beta +\delta +\mu -b\right)I}, C=\frac{R+h\beta I}{1+h(\alpha +\mu )}\).
The system (27–30) converges to equilibrium points of the model, it satisfies the condition \(\rho (J), where \(\rho (J)\) is spectral radius of the Jacobean. If \(\rho (J)>1\), it is unstable then \(\rho \left(J\right)=1\) is neutrally stable. The Jacobian matrix as
$$J=\left[\begin{array}{ccc}\frac{\partial A}{\partial S}& \frac{\partial A}{\partial I}& \frac{\partial A}{\partial R}\\ \frac{\partial B}{\partial S}& \frac{\partial B}{\partial I}& \frac{\partial B}{\partial R}\\ \frac{\partial C}{\partial S}& \frac{\partial C}{\partial I}& \frac{\partial C}{\partial R}\end{array}\right]$$
$$\frac{\partial A}{\partial S}=\frac{\partial }{\partial S}\left(\frac{S+ha+h\alpha R}{1+hcSI\left(1+\gamma I\right){e}^{-\mu \tau }+h\mu }\right), \frac{\partial A}{\partial S}=\frac{1}{1+hcI\left(1+\gamma I\right){e}^{-\mu \tau }+k\mu }$$
$$\frac{\partial A}{\partial I}=\frac{\partial }{\partial I}\left(\frac{S+ha+h\alpha R}{1+hcSI\left(1+\gamma I\right){e}^{-\mu \tau }+h\mu }\right), \frac{\partial A}{\partial I}=\frac{(S+ha+h\alpha R)(hc+2hc\gamma I{e}^{-\mu \tau })}{{(1+hcI\left(1+\gamma I\right){e}^{-\mu \tau })}^{2}}$$
$$\frac{\partial A}{\partial R}=\frac{\partial }{\partial R}\left(\frac{S+ha+h\alpha R}{1+hcSI\left(1+\gamma I\right){e}^{-\mu \tau }+h\mu }\right), \frac{\partial A}{\partial R}=\frac{h\alpha }{1+hcI(1+\gamma I)}$$
$$\frac{\partial B}{\partial S}=\frac{hcI(1+\gamma I){e}^{-\mu \tau }}{1+h(\beta +\mu +\delta -b)},\frac{\partial B}{\partial I}=\frac{1+hcS+2hcS\gamma I+hb}{1+h(\beta +\mu +\delta -b)}, \frac{\partial B}{\partial R}=0$$
$$\frac{\partial C}{\partial S}=0,\frac{\partial C}{\partial I}=\frac{h\beta }{1+h(\alpha +\beta )},\frac{\partial C}{\partial R}=\frac{1}{1+h(\alpha +\beta )}.$$
At corona free equilibrium point \({E}_{1}=(\frac{a}{\mu },\text{0,0})\). The Jacobian matrix is
$$J({E}_{1})=\left[\begin{array}{ccc}\frac{1}{1+h\mu }& -\left(\frac{a}{\mu }+ha\right)hc& h\alpha \\ 0& \frac{1+\frac{hac}{\mu }+hb}{1+h(\beta +\mu +\delta -b)}& 0\\ 0& \frac{h\beta }{1+h(\alpha +\beta )}& \frac{1}{1+h(\alpha +\beta )}\end{array}\right]$$
The eigenvalues are \({\lambda }_{1}=\frac{1}{1+h\mu }, \({\lambda }_{2}=\frac{\mu +hac\mu +h\mu b}{\mu (1+h\left(\beta +\mu +\delta -b\right))} \({\lambda }_{3}=\frac{1}{1+h(\alpha +\beta )}
At endemic equilibrium point \({E}_{2 }=({S}^{*},{I}^{*},{R}^{*})\). The Jacobian matrix is
$$J({E}_{2})=\left[\begin{array}{ccc}\frac{1}{1+hc{I}^{*}\left(1+\gamma {I}^{*}\right){e}^{-\mu \tau }+k\mu }& \frac{({S}^{*}+ha+h\alpha {R}^{*})(hc+2hc\gamma {I}^{*}{e}^{-\mu \tau })}{{(1+hc{I}^{*}\left(1+\gamma {I}^{*}\right){e}^{-\mu \tau })}^{2}}& \frac{h\alpha }{1+hc{I}^{*}(1+\gamma {I}^{*})}\\ \frac{hc{I}^{*}(1+\gamma {I}^{*}){e}^{-\mu \tau }}{1+h(\beta +\mu +\delta -b)}& \frac{1+hc{S}^{*}+2hc{S}^{*}\gamma {I}^{*}+hb}{1+h(\beta +\mu +\delta -b)}& 0\\ 0& \frac{h\beta }{1+h(\alpha +\beta )}& \frac{1}{1+h(\alpha +\beta )}\end{array}\right]$$
By using MATLAB and the values of parameters given in Table 2, verified that the maximum eigenvalues for both Jacobian \(J({E}_{1})\) and \(J({E}_{2})\) are less than one, as desired.
Discover more from PressNewsAgency
Subscribe to get the latest posts sent to your email.