A numerical approach for an epidemic SIR model via Morgan-Voyce series

As the ordinary differential system of equations is considerably used to explain mathematical properties, they are used to represent the real-world problems arising in the field of science such as biology, physics, engineering, finance, or sociology. It is well-known that many models in these areas can be described by using a system of ordinary differential equations. Analytically, it is generally very hard or impossible to solve these models. Moreover, as the exact solutions of these models cannot be found, we need to use the numerical methods to find approximate solutions. Many researchers used methods for solving such equations. Hojjati et al. [1] investigated the multistep method. Mastorakis [2] presented the collocation method. Shawagfed et al. [3] proposed the Adomian decomposition method. While Yüzbaşı ve Karaçayır [4] handled the Galerkin method, Yüzbaşı [5] proposed the exponential collocation method. Al-Omari et al. [6] used the Galerkin finite element method to observe numerical properties of models.

In modern times, many experts studied the epidemic models. Iwa et al. studied malaria and cholera models by using the Caputo-Fabrizio derivative of order a ∈ (0,1) varied with some notable parameters in the fractional system [7]. Atede et al. considered a fractional-order vaccination model for the novel Coronavirus 2019 (COVID-19) incorporating environmental transmission and analyzed using tools of fractional calculus [8]. Joshi et al. proposed a non-singular SIR model with the Mittag-Leffler law [9]. Odionyenma et al. investigated the fractional derivative model of Chlamydia-Gonorrhea co-infection via Caputo derivative [10]. Sabbar proposed an epidemic model with three intervention measures: media coverage, isolation, and medical therapy [11]. Joshi and Yavuz presented the fractional-order coinfection model for the transmission dynamics of COVID-19 and tuberculosis [12].

This paper considers the SIR model being an epidemic model. This model is used to study the spread of infectious diseases in a population. Unlike the SIR model, there is an additional factor for vaccination. In a vaccination-inclusive model, controlling the epidemic can be achieved by forecasting its future direction. Because biologists can use it to deal with natural initial conditions to design and conduct experiments in observing the spreading of diseases, with the help of these experiments and modeling, the way of controlling the spreading of the epidemic can be learned. Getting exact solutions for the problems representing such natural events is extremely difficult. Searching for convenient methods is an excellent mission for the science community.

The SIR model is composed of three members: those who are susceptible to disease, those who might be infected with the disease, and those who recovered from or become immune to it. Kermack and McKendrick earlier [13] discovered the SIR model.

Rachah and Torres improved the SIR model [14]. They expressed the susceptible-infectious-recovery (SIR) model with vaccination and the initial conditions as follows, (1) dSdt=−βStIt−νS(t)dIdt=βStIt−γI(t)dRdt=γIt+νS(t) \[\begin{align} & \frac{dS}{dt}=-\beta S\left( t \right)I\left( t \right)-\nu S(t) \\ & \frac{dI}{dt}=\beta S\left( t \right)I\left( t \right)-\gamma I(t) \\ & \frac{dR}{dt}=\gamma I\left( t \right)+\nu S(t) \\ \end{align}\] with initial conditions given by S(0)=NS,I(0)=NI,R0=NR \[S(0)={{N}_{S}},\ I(0)={{N}_{I}},\ R\left( 0 \right)={{N}_{R}}\] where ν is the percentage of individuals vaccinated every day [15].

A person can only be in one of the three groups in any time.

Recently, some researchers have analyzed the SIR model using various methods to obtain an approximate solution. Argub and El-Ajou [16] applied the Homotopy analysis method for different parameter values for the SIR model. Biazar [17] presented the Adomian decomposition method. While Rafai et al. [18] studied the homotopy perturbation method, İbrahim and Ismail [19] used a differential transformation approach. This model was also solved using the Laplace-Adomian decomposition method [20]. Harman ve Johnston [21] solved by using the stochastic Galerkin method. This epidemic model was solved by Kousar et al. [22] by using the 4th-order Runge Kutta method. It was also solved by the Hussain et al. [23] Runge Kutta-2 ve Runge Kutta-4 methods. Kousar solved this problem via the Euler method [22].

In this paper, the solutions of the SIR model have been found using the Morgan-Voyce collocation method (M-VCM). Many researchers have recently utilized the collocation method to solve this differential equation. Taylor [24,25,26,27,28], Pell-Lucas [ 29], Morgan-Voyce [ 30,31,32,33,34], Dickson [35], ortho-exponential [36], and Bernoulli [37] matrix methods have been used to solve the linear differential, integro-differential equations, fractional differential equations [38], nonlinear differential equations, multi-pantograph [39], delay differential equations and difference and Fredholm integro-difference equation [40]. This method can significantly aid in reducing the complexity of solutions for large-scale epidemic models and can contribute to obtaining results that closely approximate the exact solutions of ordinary differential equation systems.

This study is organized as follows. Section 2 describes the method. We give the error estimation of the method in Section 3. Several applications are provided in Section 4. Finally, we introduce the results of this paper in Conclusion 5.

This section presents the M-VCM in clear. This method is used to obtain the approximate solution of the equation (1). Equation (1) is expressed in the truncated Morgan-Voyce series form as follows, (2) S(t)=∑l=0LA1,lBl(t),I(t)=∑l=0LA2,lBl(t),R(t)=∑l=0LA3,lBl(t). \[S(t)=\overset{L}{\mathop{\underset{l=0}{\mathop{\sum{}}}\,}}\,{{A}_{1,l}}{{B}_{l}}(t),\ I(t)=\overset{L}{\mathop{\underset{l=0}{\mathop{\sum{}}}\,}}\,{{A}_{2,l}}{{B}_{l}}(t),\ R(t)=\overset{L}{\mathop{\underset{l=0}{\mathop{\sum{}}}\,}}\,{{A}_{3,l}}{{B}_{l}}(t).\]

Here, A_1,l, A_2,l and A_3,l (l = 0, 1, 2, · , L) are the unknown Morgan-Voyce coefficients, L is an arbitrary positive integer such as (L ≥ 2). Swamy earlier defined the Morgan-Voyce polynomials (B_l(t), (l = 0, 1, 2, · , L)) [41]. B_l(t), (l = 0, 1, 2, · , L) are presented by, BL(t)=∑Lk=0L+k+1L−ktk. \[{{B}_{L}}(t)=\underset{k=0}{\mathop{\overset{L}{\mathop{\sum{}}}\,}}\,\left( \begin{matrix} L+k+1 \\ L-k \\\end{matrix} \right){{t}^{k}}.\ \ \ \]

To write equation (1) in matrix form, we can first write the approximate solutions as follows, (3) S(t)=B(t)A1,I(t)=B(t)A2,R(t)=B(t)A3, \[S(t)=B(t){{A}_{1}},\ I(t)=B(t){{A}_{2}},\ R(t)=B(t){{A}_{3}},\] where, B(t)=[B1(t)B2(t)⋯BL−1(t)BL(t)],A1=[a1,0a1,1⋯a1,L]T,A2=[a2,0a2,1⋯a2,L]T,A3=[a3,0a3,1⋯a3,L]T,BT(t)=[ B0(t)B1(t)B2(t)⋮BL(t) ]=[(10)00⋯0(21)(30)0⋯0(32)(41)(50 )⋯0⋮⋮⋮⋯⋮(L+1L )(L+2L−1 )(L+3L−2 )⋯(2L+10 ) ][1tt2 ⋮tn ], \[\begin{gathered} B(t) = \left[ {\begin{array}{*{20}c} {B_1 (t)} \hfill & {B_2 (t)} \hfill & \cdots \hfill & {B_{L - 1} (t)} \hfill & {B_L (t)} \hfill \\ \end{array} } \right],\;A_1 = \left[ {\begin{array}{*{20}c} {a_{1,0} } \hfill & {a_{1,1} } \hfill & \cdots \hfill & {a_{1,L} } \hfill \\ \end{array} } \right]^T , \hfill \\ A_2 = \left[ {\begin{array}{*{20}c} {a_{2,0} } \hfill & {a_{2,1} } \hfill & \cdots \hfill & {a_{2,L} } \hfill \\ \end{array} } \right]^T ,\;A_3 = \left[ {\begin{array}{*{20}c} {a_{3,0} } \hfill & {a_{3,1} } \hfill & \cdots \hfill & {a_{3,L} } \hfill \\ \end{array} } \right]^T , \hfill \\ B^T (t) = \left[ {\begin{array}{*{20}c} {B_0 (t)} \hfill \\ {B_1 (t)} \hfill \\ {B_2 (t)} \hfill \\ \vdots \hfill \\ {B_L (t)} \hfill \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\left( {\begin{array}{*{20}c} 1 \\ 0 \\ \end{array} } \right)} & 0 & 0 & \cdots & 0 \\ {\left( {\begin{array}{*{20}c} 2 \\ 1 \\ \end{array} } \right)} & {\left( {\begin{array}{*{20}c} 3 \\ 0 \\ \end{array} } \right)} & 0 & \cdots & 0 \\ {\left( {\begin{array}{*{20}c} 3 \\ 2 \\ \end{array} } \right)} & {\left( {\begin{array}{*{20}c} 4 \\ 1 \\ \end{array} } \right)} & {\left( {\begin{array}{*{20}c} 5 \\ 0 \\ \end{array} } \right)} & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ {\left( {\begin{array}{*{20}c} {L + 1} \\ L \\ \end{array} } \right)} & {\left( {\begin{array}{*{20}c} {L + 2} \\ {L - 1} \\ \end{array} } \right)} & {\left( {\begin{array}{*{20}c} {L + 3} \\ {L - 2} \\ \end{array} } \right)} & \cdots & {\left( {\begin{array}{*{20}c} {2L + 1} \\ 0 \\ \end{array} } \right)} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} 1 \\ t \\ {t^2 } \\ \vdots \\ {t^n } \\ \end{array} } \right], \hfill \\ \end{gathered} \] here we can write, F(t)=1tt2⋯tL,R=1000⋯021300⋯0324150⋯0⋮⋮⋮⋯⋮L+1LL+2L−1L+3L−2⋯2L+10. \[F(t)=\left[ \begin{array}{*{35}{l}} 1 & t & {{t}^{2}} & \cdots & {{t}^{L}} \\\end{array} \right],\,R\ =\left[ \begin{matrix} \left( \begin{matrix} 1 \\ 0 \\\end{matrix} \right) & 0 & 0 & \cdots & 0 \\ \left( \begin{matrix} 2 \\ 1 \\\end{matrix} \right) & \left( \begin{matrix} 3 \\ 0 \\\end{matrix} \right) & 0 & \cdots & 0 \\ \left( \begin{matrix} 3 \\ 2 \\\end{matrix} \right) & \left( \begin{matrix} 4 \\ 1 \\\end{matrix} \right) & \left( \begin{matrix} 5 \\ 0 \\\end{matrix} \right) & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ \left( \begin{matrix} L+1 \\ L \\\end{matrix} \right) & \left( \begin{matrix} L+2 \\ L-1 \\\end{matrix} \right) & \left( \begin{matrix} L+3 \\ L-2 \\\end{matrix} \right) & \cdots & \left( \begin{matrix} 2L+1 \\ 0 \\\end{matrix} \right) \\\end{matrix} \right].\]

Therefore, we can write the following equations, (4) S(t)=F(t)RTA1,I(t)=F(t)RTA2,R(t)=F(t)RTA3. \[S(t)=F(t){{R}^{T}}{{A}_{1}},\ I(t)=F(t){{R}^{T}}{{A}_{2}},\ R(t)=F(t){{R}^{T}}{{A}_{3}}.\]

The relationship between matrix F(t) and its derivative matrix F⁽¹⁾(t) is as follows, (5) F(1)(t)=F(t)TT \[{{F}^{(1)}}(t)=F(t){{T}^{T}}\] and here, TT=010⋯0002⋯0⋮⋮⋮⋱⋮000⋯L000⋯0. \[{{T}^{T}}=\left[ \begin{array}{*{35}{l}} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & L \\ 0 & 0 & 0 & \cdots & 0 \\\end{array} \right].\]

From equations (4) and (5), we can get the following equations, S(1)(t)=F(t)TTRTA1,I(1)(t)=F(t)TTRTA2,R(1)(t)=F(t)TTRTA3. \[{{S}^{(1)}}(t)=F(t){{T}^{T}}{{R}^{T}}{{A}_{1}},\ {{I}^{(1)}}(t)=F(t){{T}^{T}}{{R}^{T}}{{A}_{2}},\ {{R}^{(1)}}(t)=F(t){{T}^{T}}{{R}^{T}}{{A}_{3}}.\]

Thus, we can obtain the matrices y(t) and y⁽¹⁾(t) as follows, y(t)=F¯R¯A,y1(t)=F¯T¯R¯A, \[y(t)=\overline{F}\overline{R}A,\ {{y}^{1}}(t)=\overline{F}\overline{T}\overline{R}A,\] where, y(t)=S(t)I(t)R(t),y(1)(t)=S(1)(t)I(1)(t)R(1)(t),F¯(t)=F(t)000F(t)000F(t),T¯=TT000TT000TT, \[y(t)=\left[ \begin{matrix} S(t) \\ I(t) \\ R(t) \\\end{matrix} \right],\ {{y}^{(1)}}(t)=\left[ \begin{matrix} {{S}^{(1)}}(t) \\ {{I}^{(1)}}(t) \\ {{R}^{(1)}}(t) \\\end{matrix} \right],\ \overline{F}(t)=\left[ \begin{matrix} F(t) & 0 & 0 \\ 0 & F(t) & 0 \\ 0 & 0 & F(t) \\\end{matrix} \right],\ \overline{T}=\left[ \begin{matrix} {{T}^{T}} & 0 & 0 \\ 0 & {{T}^{T}} & 0 \\ 0 & 0 & {{T}^{T}} \\\end{matrix} \right],\] and, R¯=RT000RT000RT,A=A1A2A3. \[\overline{R}=\left[ \begin{matrix} {{R}^{T}} & 0 & 0 \\ 0 & {{R}^{T}} & 0 \\ 0 & 0 & {{R}^{T}} \\\end{matrix} \right],\ A=\left[ \begin{matrix} {{A}_{1}} \\ {{A}_{2}} \\ {{A}_{3}} \\\end{matrix} \right].\]

Equation (1) can be written in matrix form as follows, (6) y1(t)−Dy(t)−Ny(t)−Ey1,2(t)=g, \[{{y}^{1}}(t)-Dy(t)-Ny(t)-E{{y}_{1,2}}(t)=g,\] where, g=000,D=0000−γ00γ0,E=−ββ0,N=−ν00000ν00,y1,2=S(t)I(t). \[g=\left[ \begin{matrix} 0 \\ 0 \\ 0 \\\end{matrix} \right],\ D=\left[ \begin{matrix} 0 & 0 & 0 \\ 0 & -\gamma & 0 \\ 0 & \gamma & 0 \\\end{matrix} \right],\ E=\left[ \begin{matrix} -\beta \\ \beta \\ 0 \\\end{matrix} \right],\ N=\left[ \begin{matrix} -\nu & 0 & 0 \\ 0 & 0 & 0 \\ \nu & 0 & 0 \\\end{matrix} \right],\ {{y}_{1,2}}=\left[ S(t)I(t) \right].\]

Now, to determine the unknown coefficients of A_1,I, A_2,I and A_3,I, we can define the collocation points for the interval a ≤ t ≤ b as follows, ti=a+b−aLi,i=0,1,⋯,L. \[{{t}_{i}}=a+\frac{b-a}{L}i,\ \ \ i=0,\ 1,\ \cdots ,\ L.\]

Substituting the collocation points in equation (6), we get the following system, y(1)(ti)−Dy(ti)−Ny(ti)−Ey1,2(ti)=g, \[{{y}^{(1)}}({{t}_{i}})-Dy({{t}_{i}})-Ny({{t}_{i}})-E{{y}_{1,2}}({{t}_{i}})=g,\] where, (7) y(ti)=F¯(ti)R¯A, \[y({{t}_{i}})=\overline{F}({{t}_{i}})\overline{R}A,\] (8) y(1)(ti)=F¯(ti)T¯R¯A, \[{{y}^{(1)}}({{t}_{i}})=\overline{F}({{t}_{i}})\overline{T}\overline{R}A,\] here, F=F¯(t0)F¯(t1)⋯F¯(tL)T,F¯(ti)=F(ti)000F(ti)000F(ti). \[F={{\left[ \begin{matrix} \overline{F}({{t}_{0}}) & \overline{F}({{t}_{1}}) & \cdots & \overline{F}({{t}_{L}}) \\\end{matrix} \right]}^{T}},\ \overline{F}({{t}_{i}})=\left[ \begin{matrix} F({{t}_{i}}) & 0 & 0 \\ 0 & F({{t}_{i}}) & 0 \\ 0 & 0 & F({{t}_{i}}) \\\end{matrix} \right].\]

So, equation (1) can be written in the following matrix form utilizing the upper matrices, (9) Y(1)−D¯Y−N¯Y−E¯Y¯˜=G, \[{{Y}^{(1)}}-\overline{D}Y-\overline{N}Y-\overline{E}\widetilde{\overline{Y}}=G,\] where, Y(1)=y(1)(t0)y(1)(t1)⋮y(1)(tL),D¯=D0⋯00D⋯0⋮⋮⋱⋮00⋯D(L+1)(L+1),N¯=N0⋯00N⋯0⋮⋮⋱⋮00⋯N(L+1)(L+1),Y=y(t0)y(t1)⋮y(tL),G=gg⋮g(L+1)1,Y¯˜=y1,2(t0)y1,2(t1)⋮y1,2(tL),E¯=E0⋯00E⋯0⋮⋮⋱⋮00⋯E(L+1)(L+1). \[\begin{align} & {{Y}^{(1)}}=\left[ \begin{matrix} {{y}^{(1)}}({{t}_{0}}) \\ {{y}^{(1)}}({{t}_{1}}) \\ \vdots \\ {{y}^{(1)}}({{t}_{L}}) \\\end{matrix} \right],\ \overline{D}={{\left[ \begin{matrix} D & 0 & \cdots & 0 \\ 0 & D & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & D \\\end{matrix} \right]}_{(L+1)(L+1)}},\ \overline{N}={{\left[ \begin{matrix} N & 0 & \cdots & 0 \\ 0 & N & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & N \\\end{matrix} \right]}_{(L+1)(L+1)}}, \\ & Y=\left[ \begin{matrix} y({{t}_{0}}) \\ y({{t}_{1}}) \\ \vdots \\ y({{t}_{L}}) \\\end{matrix} \right],\ G={{\left[ \begin{matrix} g \\ g \\ \vdots \\ g \\\end{matrix} \right]}_{(L+1)1}},\ \widetilde{\overline{Y}}=\left[ \begin{matrix} {{y}_{1,2}}({{t}_{0}}) \\ {{y}_{1,2}}({{t}_{1}}) \\ \vdots \\ {{y}_{1,2}}({{t}_{L}}) \\\end{matrix} \right],\ \overline{E}={{\left[ \begin{matrix} E & 0 & \cdots & 0 \\ 0 & E & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & E \\\end{matrix} \right]}_{(L+1)(L+1)}}. \\ \end{align}\]

Similarly, the following matrix form is obtained by substituting the collocation points in y_1,2(t), Y¯˜=y1,2(t0)y1,2(t1)⋮y1,2(tL)=I(t0)0⋯00I(t1)⋯0⋮⋮⋱⋮00⋯I(tL)S(t0)S(t1)⋮S(tL)=I¯S¯, \[\widetilde{\overline{Y}}=\left[ \begin{matrix} {{y}_{1,2}}({{t}_{0}}) \\ {{y}_{1,2}}({{t}_{1}}) \\ \vdots \\ {{y}_{1,2}}({{t}_{L}}) \\\end{matrix} \right]=\left[ \begin{matrix} I({{t}_{0}}) & 0 & \cdots & 0 \\ 0 & I({{t}_{1}}) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & I({{t}_{L}}) \\\end{matrix} \right]\left[ \begin{matrix} S({{t}_{0}}) \\ S({{t}_{1}}) \\ \vdots \\ S({{t}_{L}}) \\\end{matrix} \right]=\overline{I}\overline{S},\] where, (10) I¯=F˜R˜A2¯ and S¯=T˜˜R˜˜A, \[\overline{I}=\widetilde{F}\widetilde{R}\overline{{{A}_{2}}}\ \ \ \text{and}\ \ \ \overline{S}=\widetilde{\widetilde{T}}\widetilde{\widetilde{R}}A,\] and, F˜=F(t0)0⋯00F(t1)⋯0⋮⋮⋱⋮00⋯F(tL),A2¯=A20⋯00A2⋯0⋮⋮⋱⋮00⋯A2(L+1)(L+1),R˜=RT0⋯00RT⋯0⋮⋮⋱⋮00⋯RT(L+1)(L+1),F˜˜=F(t0)F(t1)⋮F(tL),R˜˜=RTKK,K=00⋯000⋯0⋮⋮⋱⋮00⋯0(L+1)(L+1). \[\begin{align} & \widetilde{F}=\left[ \begin{matrix} F({{t}_{0}}) & 0 & \cdots & 0 \\ 0 & F({{t}_{1}}) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & F({{t}_{L}}) \\\end{matrix} \right],\ \overline{{{A}_{2}}}={{\left[ \begin{matrix} {{A}_{2}} & 0 & \cdots & 0 \\ 0 & {{A}_{2}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & {{A}_{2}} \\\end{matrix} \right]}_{(L+1)(L+1)}},\ \widetilde{R}={{\left[ \begin{matrix} {{R}^{T}} & 0 & \cdots & 0 \\ 0 & {{R}^{T}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & {{R}^{T}} \\\end{matrix} \right]}_{(L+1)(L+1)}}, \\ & \widetilde{\widetilde{F}}=\left[ \begin{matrix} F({{t}_{0}}) \\ F({{t}_{1}}) \\ \vdots \\ F({{t}_{L}}) \\\end{matrix} \right],\ \widetilde{\widetilde{R}}=\left[ \begin{matrix} {{R}^{T}} & K & K \\\end{matrix} \right],\ K={{\left[ \begin{matrix} 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \\\end{matrix} \right]}_{(L+1)(L+1)}}. \\ \end{align}\]

The fundamental matrix equation from equations (7)–(10) is as follows, (11) FB¯R¯−D¯FR¯−N¯FR¯−E¯F˜R˜A2¯F˜˜R˜˜A=G. \[\left\{ F\overline{B}\overline{R}-\overline{D}F\overline{R}-\overline{N}F\overline{R}-\overline{E}\widetilde{F}\widetilde{R}\overline{{{A}_{2}}}\widetilde{\widetilde{F}}\widetilde{\widetilde{R}} \right\}A=G.\]

Briefly, equation (11) can be written as, WA=G or [W;G], \[WA=G\ \ {\rm or}\ \ [W;G],\] where, (12) W=FT¯R¯−D¯FR¯−N¯FR¯−E¯F˜R˜A2¯F˜˜R˜. \[W=F\overline{T}\overline{R}-\overline{D}F\overline{R}-\overline{N}F\overline{R}-\overline{E}\widetilde{F}\widetilde{R}\overline{{{A}_{2}}}\widetilde{\widetilde{F}}\widetilde{R}.\]

Equation (12) refers to a nonlinear system of 3(L + 1) algebraic equations with the unknown Morgan-Voyce coefficients A_1,I, A_2,I and A_3,I. The matrix form of the initial conditions for t → 0 in equation (3) can be written as follows, S(t)=B(0)A1=NSI(t)=B(0)A2=NIR(t)=B(0)A3=NR. \[\begin{align} & S(t)=B(0){{A}_{1}}=\left[ {{N}_{S}} \right] \\ & I(t)=B(0){{A}_{2}}=\left[ {{N}_{I}} \right] \\ & R(t)=B(0){{A}_{3}}=\left[ {{N}_{R}} \right]. \\ \end{align}\]

Thus, these matrix forms can be written as, (13) U1=S(0)=a1,0a1,1⋯a1,LU2=I(0)=a2,0a2,1⋯a2,LU3=R(0)=a3,0a3,1⋯a3,L. \[\begin{align} & {{U}_{1}}=S(0)=\left[ \begin{matrix} {{a}_{1,0}} & {{a}_{1,1}} & \cdots & {{a}_{1,L}} \\\end{matrix} \right] \\ & {{U}_{2}}=I(0)=\left[ \begin{matrix} {{a}_{2,0}} & {{a}_{2,1}} & \cdots & {{a}_{2,L}} \\\end{matrix} \right] \\ & {{U}_{3}}=R(0)=\left[ \begin{matrix} {{a}_{3,0}} & {{a}_{3,1}} & \cdots & {{a}_{3,L}} \\\end{matrix} \right]. \\ \end{align}\]

Finally, we replace the row matrices in equation (13) with any three rows of the matrix in equation (12), and we get a solution for equation (1) under initial conditions. Finally, we get the augmented matrix as, W˜A=G. \[\widetilde{W}A=G.\]

First of all, this system is solved and the unknown coefficients a_i,0, a_i,1, · , a_i,L,(i = 1, 2, 3) are found. Then, these coefficients are substituted in equation (2) and, approximate solutions are found.

In this section, we investigate the accuracy of the presented method. S_L(t), I_L(t), R_L(t) are the approximate solutions of equation (1). Therefore, these functions approximately satisfy the equation (1). Once we substitute these functions and the first derivatives of these functions in the equation (1), we get the error estimations for t = t_r ∈ [0,R], r = 0, 1, ... as following, E1,L(tr)=|S′(tr)+βS(tr)I(tr)+νS(tr)|≅0E2,L(tr)=|I′(tr)−βS(tr)I(tr)+γI(tr)|≅0E3,L(tr)=|R′(tr)−γI(tr)−νS(tr)|≅0, \[\begin{align} & {{E}_{1,L}}({{t}_{r}})=|{S}'({{t}_{r}})+\beta S({{t}_{r}})I({{t}_{r}})+\nu S({{t}_{r}})|\cong 0 \\ & {{E}_{2,L}}({{t}_{r}})=|{I}'({{t}_{r}})-\beta S({{t}_{r}})I({{t}_{r}})+\gamma I({{t}_{r}})|\cong 0 \\ & {{E}_{3,L}}({{t}_{r}})=|{R}'({{t}_{r}})-\gamma I({{t}_{r}})-\nu S({{t}_{r}})|\cong 0, \\ \end{align}\] and E_i,L ≤ 10^−k_r, i = 1, 2, 3 (k_r is a positive number). If we get max10^−k_r = 10^−k, as L increases, error estimation E_i,L(t_r), (t = 1, 2, 3) at each of the points decreases and becomes smaller than 10^−k.

This section supports the accuracy and efficiency of the presented method. We offer two numerical examples, one comparing the errors and solutions with the other methods. In the second sample, the results for different vaccination rates are compared. The tables and figures show the comparisons at specific points within the given interval. Numerical computations are performed using MATLABR2021a software.

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Fig. 8

In this paper, we used the M-VCM to find the approximate solution of the governing model. Several applications were presented to show the efficiency and accuracy of the proposed method. To display the accuracy of the approximate solutions, these solutions are once again substituted into equation (1) using MATLAB software. Thus, this work offers a different measure for the reliability of the approximate solutions. The obtained approximate solutions and errors have been compared by the HPM [18] and the LADM [20]. These comparisons have shown that the presented method to find approximate solutions to the SIR model from epidemic models is more efficient and valuable. Tables 2, 4, and 6 illustrated that the numerical solutions of HPM [18] and LADM [20] are nearly the same. As can be seen from the figures and tables, minimal errors were obtained in this study, indicating that the results are very close to the exact solutions. In particular, in the second example, we can accurately determine the vaccination rate needed to control the epidemic. In conclusion, all these findings show that our method is efficient and feasible for solving the nonlinear ODE systems. The most significant advantage of the proposed method is that all calculations can be easily and quickly applied with the help of MATLAB code. In addition, the fact that the SIR model is a real-life problem increases the importance of our study. One of the most significant advantages of the method is changing the system equation (1) into an easily solved nonlinear algebraic equation system. As a further study, this method can be applied to other epidemic models for longer time periods.