Accurate Approximations for a Nonlinear SIR System via an Efficient Analytical Approach: Comparative Analysis

Abstract

The homotopy perturbation method (HPM) is one of the recent fundamental methods for solving differential equations. However, checking the accuracy of this method has been ignored by some authors in the literature. This paper reanalyzes the nonlinear system of ordinary differential equations (ODEs) describing the SIR epidemic model, which has been solved in the literature utilizing the HPM. The main objective of this work is to obtain a highly accurate analytical solution for this model via a direct technique. The proposed technique is mainly based on reducing the given system to a single nonlinear ODE that can be easily solved. Numerical results are conducted to compare our approach with the previous HPM, where the Runge–Kutta numerical method is chosen as a reference solution. The obtained results reveal that the current technique exhibits better accuracy over HPM in the literature. Moreover, some physical properties are introduced and discussed in detail regarding the influence of the transmission rate on the behavior of the SIR model.

1. Introduction

Modeling the outbreak and spread of diseases has a long history. John Graunt [1] may have been the first to quantify causes of death systematically. However, the first mathematical model describing an infectious disease was proposed by Daniel Bernoulli [2] in 1760. Moreover, McKendrick and Kermack introduced a simple deterministic (compartmental) model describing the behavior of epidemic outbreaks in 1927 [3], known as the Susceptible–Infected–Recovered (SIR) model. Based on these works, a number of mathematical models have been developed and established under various conditions and situations [4,5]. In addition, some authors used the basic idea of the SIR model to formulate several mathematical models [6,7,8,9] to predict the behavior of COVID-19, which is an infectious disease caused by the SARS-CoV-2 virus. Various ISR models have been implemented to investigate the progress of COVID-19 in different countries [10,11,12,13,14], while other mathematical models have been developed to include different factors [15,16].

As usual, the proposed models in the literature are governed by systems of linear/nonlinear ordinary differential equations. There is no doubt that obtaining accurate solutions for such mathematical models leads to better interpretation and prediction of the pandemic future. However, checking the accuracy of the obtained solution has been ignored in some studies in the literature. For example, the authors of [7] considered the following nonlinear SIR/COVID-19 model:

$$\begin{array}{ccc}& & \frac{dR}{d\tau }=I\left(\tau \right),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & \frac{dI}{d\tau }=\sigma \left[1-R\left(\tau \right)-I\left(\tau \right)\right]I\left(\tau \right)-I\left(\tau \right),\hfill \end{array}$$

(2)

where $\tau =t/T$ (t represents the time in days and T is the time of transmission of the virus, which changes from 2–4 weeks). $I\left(t\right)$ denotes the infected individuals who are carrying the virus, $R\left(t\right)$ represents the recovered individuals, $S\left(t\right)$ stands for susceptible individuals with $S\left(t\right)=1-R\left(t\right)-I\left(t\right)$, and $\sigma$ is the physical contact number between susceptible and infected individuals.

At the initial time $t=0$ of the outbreak, the initial number of infected people and the initial recovered people are regularly zero. However, these initial numbers of infected and the recovered people may be assumed in the general forms $R\left(0\right)=A$ and $I\left(0\right)=B$, respectively, where A and B are the real constants. Hence, the current model is solved under the initial conditions (ICs):

$$R\left(0\right)=A,I\left(0\right)=B.$$

(3)

There are numerous analytical methods that can be effectively used to deal with nonlinear models (1)–(3), such as the differential transform method (DTM) [17], the homotopy analysis method (HAM) [18,19], the homotopy perturbation method (HPM) [20,21], and the Adomian decomposition method (ADM) [22].

However, each of these methods has its own difficulties when applied to solve nonlinear systems. So, the main incentive of this work is to introduce a simple analytical approach that directly obtains the desired accurate approximate solution. Although the authors of [7] applied HPM to solve nonlinear models (1–2) under ICs (3), there are some remarks on the accuracy of their approximate solution.

The main objective of this work is to obtain accurate approximate solution of the present nonlinear COVID-19 model. The proposed approach suggests combining the governing Equations (1) and (2) to formulate a new single nonlinear equation in $R\left(\tau \right)$. The new single equation in $R\left(\tau \right)$ will be solved using two degrees or orders of approximations. The first-order and the second-order approximate analytical solutions are denoted as FOAS and SOAS, respectively. It will be proven that the FOAS and SOAS exhibit a better accuracy when compared with the homotopy perturbation method (HPM) [7]. The present results will be validated via performing various comparisons with the numerical results using the Runge–Kutta method. In addition, the effectiveness and efficiency of the present approach over the HPM in the literature will be discussed. Furthermore, the impact of the contact number $\sigma$ (transmission rate of the virus) on the number of infected individuals and recovered individuals will be explained and interpreted. Moreover, the features and advantages of the current analysis will be addressed.

2. Reduced Model

This section indicates that models (1–3) can be reduced to a new form with a single unknown. The proposed approach is based on combining the two ordinary differential equations (ODEs) (1) and (2). This approach leads to a single 2nd-order ODE in only $R\left(\tau \right)$, which can be easily converted to a 1st-order nonlinear ODE in $R\left(\tau \right)$. Before continuing to the main target of this section, we rewrite systems (1–2) as follows

$$\begin{array}{ccc}& & R^{\prime }\left(\tau \right)=I\left(\tau \right),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & I^{\prime }\left(\tau \right)=\sigma \left[1-R\left(\tau \right)-I\left(\tau \right)\right]I\left(\tau \right)-I\left(\tau \right).\hfill \end{array}$$

(5)

The primes denote the derivative with respect to $\tau$. Differentiating Equation (4) once with respect to $\tau$ yields

$$R^{″}\left(\tau \right)=I^{\prime }\left(\tau \right).$$

(6)

Substituting Equations (4) and (6) into Equation (5) provides

$$R^{″}\left(\tau \right)=\sigma \left[1-R\left(\tau \right)-R^{\prime }\left(\tau \right)\right]R^{\prime }\left(\tau \right)-R^{\prime }\left(\tau \right).$$

(7)

which is a 2nd-order nonlinear ODE and it is subjected to the ICs:

$$R\left(0\right)=A,R^{\prime }\left(0\right)=B.$$

(8)

Assume that

$$\psi \left(\tau \right)=R\left(\tau \right)+R^{\prime }\left(\tau \right).$$

(9)

Then, Equation (7) becomes

$$\frac{{\psi }^{\prime }\left(\tau \right)}{1-\psi \left(\tau \right)}=\sigma R^{\prime }\left(\tau \right).$$

(10)

Integrating both sides with respect to $\tau$ implies

$$\psi \left(\tau \right)=1+c{e}^{-\sigma R\left(\tau \right)},$$

(11)

where c is a constant of integration to be determined later. From Equations (9) and (11), we obtain

$$R^{\prime }\left(\tau \right)=1-R\left(\tau \right)+c{e}^{-\sigma R\left(\tau \right)}.$$

(12)

Applying ICs (8) to (12) provides

$$c=e^{\sigma A}(B+A-1).$$

(13)

In the next section, the 1st-order nonlinear ODE (12) will be solved to derive the solution for $R\left(\tau \right)$ via two different analytical approximations.

3. Analytical Approximations

In this section, two different approximate solutions are determined for the reduced model. Once the solution of the reduced model is evaluated for $R\left(\tau \right)$, the approximate solution of the original systems (1)–(3) is established. In addition, the accuracy of such approximations will be validated in a subsequent section by performing several comparisons with the Runge–Kutta numerical method.

3.1. First-Order Approximate Solution (FOAS)

Expanding $e^{-\sigma R}$ and using its first-order approximation leads to

$$R^{\prime }\left(\tau \right)+(1+c\sigma )R\left(\tau \right)=c+1.$$

(14)

It should be noted that the constant c in Equation (14) does not have the same value in Equation (13). In addition, the constant c in Equation (14) can be directly determined through applying the ICs (8) and this provides

$$c=\frac{B+A-1}{1-A\sigma }.$$

(15)

The exact solution of the 1st-order linear ODE (14) can be obtained as

$$R\left(\tau \right)=A{e}^{-(1+c\sigma )\tau }+\frac{c+1}{c\sigma +1}\left[1-e^{-(1+c\sigma )\tau }\right].$$

(16)

Hence, we have

$$I\left(\tau \right)=\left(c+1-A(c\sigma +1)\right)e^{-(1+c\sigma )\tau },$$

(17)

by taking the derivative of Equation (16). By inserting the value of c in the expressions above and simplifying, we obtain

$$R\left(\tau \right)=A{e}^{-\left(\frac{\sigma (B-1)+1}{1-A\sigma }\right)\tau }+\frac{B+(1-\sigma )A}{\sigma (B-1)+1}\left[1-e^{-\left(\frac{\sigma (B-1)+1}{1-A\sigma }\right)\tau }\right],$$

(18)

and

$$I\left(\tau \right)=B{e}^{-\left(\frac{\sigma (B-1)+1}{1-A\sigma }\right)\tau }.$$

(19)

3.2. Second-Order Approximate Solution (SOAS)

Extracting the second-order approximation from the expansion of $e^{-\sigma R}$ and substituting it into Equation (12), we obtain

$$R^{\prime }\left(\tau \right)=c+1-(1+c\sigma )R\left(\tau \right)+\frac{c{\sigma }^2}{2}{R}^2.$$

(20)

Applying ICs (8) to Equation (20) provides c as

$$c=\frac{R^{\prime }\left(0\right)+R\left(0\right)-1}{1-\sigma R\left(0\right)+\frac{{\sigma }^2}{2}{R}^2\left(0\right)},$$

(21)

or

$$c=\frac{2(B+A-1)}{1+{(1-\sigma A)}^2}.$$

(22)

Equation (20) can be written as

$$R^{\prime }\left(\tau \right)=\frac{c{\sigma }^2}{2}\left[R^2-\frac{2(1+c\sigma )}{c{\sigma }^2}R+\frac{2(1+c)}{c{\sigma }^2}\right],$$

(23)

or

$$R^{\prime }\left(\tau \right)=\frac{c{\sigma }^2}{2}\left[{\left(R-\frac{1+c\sigma }{c{\sigma }^2}\right)}^2+\frac{2(1+c)}{c{\sigma }^2}-{\left(\frac{1+c\sigma }{c{\sigma }^2}\right)}^2\right],$$

(24)

i.e.,

$$R^{\prime }\left(\tau \right)=\frac{c{\sigma }^2}{2}\left[{\left(R\left(\tau \right)-\beta \right)}^2-{\rho }^2\right],$$

(25)

where

$$\beta =\frac{1+c\sigma }{c{\sigma }^2},{\rho }^2={\beta }^2-\frac{2(1+c)}{c{\sigma }^2}.$$

(26)

Using the separation of the variable method, Equation (25) can be rewritten as

$$\frac{dR\left(\tau \right)}{{\left(R\left(\tau \right)-\beta \right)}^2-{\rho }^2}=\frac{c{\sigma }^2}{2}d\tau .$$

(27)

Integrating both sides, we obtain

$$-\frac{1}{\rho }tan{\mathrm{h}}^{-1}\left(\frac{R\left(\tau \right)-\beta }{\rho }\right)+\frac{1}{\rho }tan{\mathrm{h}}^{-1}\left(\frac{R\left(0\right)-\beta }{\rho }\right)=\frac{c{\sigma }^2}{2}\tau .$$

(28)

i.e.,

$$-\frac{1}{\rho }tan{\mathrm{h}}^{-1}\left(\frac{R\left(\tau \right)-\beta }{\rho }\right)+\frac{1}{\rho }tan{\mathrm{h}}^{-1}\left(\frac{A-\beta }{\rho }\right)=\frac{c{\sigma }^2}{2}\tau ,$$

(29)

which can be written as

$$R\left(\tau \right)=\beta +\rho tan\mathrm{h}\left(tan{\mathrm{h}}^{-1}\left(\frac{A-\beta }{\rho }\right)-\frac{c\rho {\sigma }^2}{2}\tau \right).$$

(30)

Expanding the right side of the last equation provides

$$R\left(\tau \right)=\beta +\rho \left[\frac{A-\beta -\rho tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}{\rho -(A-\beta )tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}\right].$$

(31)

It can be seen from this solution that it satisfies IC $R\left(0\right)=A$. Further, Solution (31) can be verified by direct substitution into Equation (20). In addition, the solution of the infected cases can be evaluated by differentiating $R\left(\tau \right)$, thus

$$I\left(\tau \right)=\frac{c{\rho }^2{\sigma }^2\left({(A-\beta )}^2-{\rho }^2\right)\sec {\mathrm{h}}^2\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}{2{\left(\rho -(A-\beta )tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)\right)}^2}.$$

(32)

This expression can be simplified by employing the relation:

$$\frac{c{\sigma }^2}{2}\left[{\left(A-\beta \right)}^2-{\rho }^2\right]=B,$$

(33)

which can be obtained from (25) through inserting ICs $R\left(0\right)=A$ and $R^{\prime }\left(0\right)=B$. Hence, the soultion of the infected cases is provided in the final form:

$$I\left(\tau \right)=\frac{B{\rho }^2\sec {\mathrm{h}}^2\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}{{\left(\rho -(A-\beta )tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)\right)}^2},$$

(34)

where c is already provided in (22). At first glance of the solution, one can see that it satisfies the IC $I\left(0\right)=B$. Moreover, it should be noted that the above expressions for $R\left(\tau \right)$ and $I\left(\tau \right)$ are valid when $\rho$ is real.

4. Features and Behaviors

Normally, it is expected that the value of the infected individuals reduces over time, while the number of recovered individuals reaches a certain constant value at a prescribed interval of time, which may be large or short depending on the nature or data of each country. Equation (12) can be used directly to obtain an accurate prediction of the recovered individuals after a relatively or sufficiently large amount of time ${\tau }_{\infty }$ (say). After this time, i.e., $\tau >{\tau }_{\infty }$, no variation is expected in the number of the recovered individuals, i.e., $R^{\prime }\left(\tau \right)=0$ for all $\tau >{\tau }_{\infty }$. Mathematically, one can consider $R^{\prime }\left(\tau \right)\to 0$ as $\tau \to \infty$. Accordingly, the number of recovered individuals, $\tau \to \infty$, denoted by $R_{\infty }$, can be predicted through solving the transcendental equation:

$$1-R_{\infty }+c{e}^{-\sigma R_{\infty }}=0,$$

(35)

where c is already provided in Equation (13) in terms of $\sigma$, A, and B. The above transcendental equation can be solved numerically using software. Moreover, approximate values of $R_{\infty }$ can be estimated from the solutions obtained in the previous section. In this regard, FOAS (18) estimates $R_{\infty }$ as

$$R_{\infty }=\underset{\tau \to \infty }{lim}R\left(\tau \right)=\frac{B+(1-\sigma )A}{\sigma B+1-\sigma }.$$

(36)

It must be noted that the value $R_{\infty }$ in the last equation is restricted by

$$\frac{\sigma (B-1)+1}{1-\sigma A}>0,$$

(37)

so that solution $I\left(\tau \right)$ in (19) of the infected individuals satisfies the condition:

$$\underset{\tau \to \infty }{lim}I\left(\tau \right)=0.$$

(38)

The possible ranges of initial values for A and B, provided that Condition (37) is satisfied, are shown in Figure 1 and Figure 2 at $\sigma =0.9$ and $\sigma =1$, respectively. Figure 1 indicates that FOAS is valid in the entire ranges of A and B when $\sigma =0.9$. However, Figure 2 reveals that FOAS is valid in certain domains of A and B when $\sigma =1$.

On the other hand, SOAS determines $R_{\infty }$ as

$$\begin{array}{ccc}\hfill R_{\infty }& =& \beta +\rho \underset{\tau \to \infty }{lim}\left[\frac{A-\beta -\rho tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}{\rho -(A-\beta )tan\mathrm{h}\left(\frac{c\rho {\sigma }^2}{2}\tau \right)}\right]=\beta +\rho \hfill \\ & =& \frac{1+c\sigma }{c{\sigma }^2}-\sqrt{{\left(\frac{1+c\sigma }{c{\sigma }^2}\right)}^2-\frac{2(1+c)}{c{\sigma }^2}},\hfill \end{array}$$

(39)

where $c=\frac{2(B+A-1)}{1+{(1-\sigma A)}^2}<0$ for all A and B, such that $0<A+B<1$. Note that Expression (39) for $R_{\infty }$ is valid when $\rho$ is real. This requires discussing the domain of applicability for the involved parameters. From Equation (26), we note that $\rho$ is real if the following condition is satisfied:

$${\left(\frac{1+c\sigma }{c{\sigma }^2}\right)}^2-\frac{2(1+c)}{c{\sigma }^2}>0.$$

(40)

Inserting the constant $c=\frac{2(B+A-1)}{1+{(1-\sigma A)}^2}$ into this condition and performing simplifications, then

$${\left(\frac{(2\sigma B+2(1-\sigma )+{\sigma }^2{B}^2)}{2{\sigma }^2(1-A-B)}\right)}^2-\frac{1-2(A+B)-{(1-\sigma A)}^2}{{\sigma }^2(1-A-B)}>0.$$

(41)

Figure 3 and Figure 4 depict the applicable and inapplicable domains of A and B for the validity of SOAS at $\sigma =0.9$ and $\sigma =1$, respectively.

Regarding the accuracy of the current approximate solutions,

Table 1. Comparisons between the approximate values of $R_{\infty }$ (the expected final number of the recovered cases) using the present FOAS and SOAS with the numerical results at different values of contact number $\sigma$ (transmission rate) for initial infected individuals $A=0.01$ and initial recovered individuals $B=0.001$.

$\mathit{\sigma }$	${\mathit{R}}_{\mathit{\infty }}$ (FOAS)	${\mathit{R}}_{\mathit{\infty }}$ (SOAS)	${\mathit{R}}_{\mathit{\infty }}$ (Numerical)
0.1	0.0120977	0.0120969	0.0120969
0.2	0.0134663	0.0134619	0.0134619
0.3	0.0152205	0.0152059	0.0152059
0.4	0.0175497	0.0175096	0.0175097
0.5	0.0207921	0.0206876	0.0206879
0.6	0.0256158	0.0253348	0.0253362
0.7	0.0335505	0.0327069	0.0327130
0.8	0.0490385	0.0458424	0.0458764
0.9	0.0926606	0.0730530	0.0733310
1.0	1.0000000	0.1324565	0.1352252

shows several comparisons between FOAS, SOAS, and the numerical solution using MATHEMATICA for the values of $R_{\infty }$. It can be seen that the values of $R_{\infty }$ via FOAS are close to the numeric ones when $0.1\le \sigma \le 0.6$. In addition, SOAS provides better estimations for $R_{\infty }$, as follows. The results in the SOAS column are identical to the numerical results when $\sigma =0.1$, $\sigma =0.2$, and $\sigma =0.3$ (which can be viewed as exact answers at these cases of $\sigma$). Moreover, the values of SOAS agree with the numerical results to six decimal places when $\sigma =0.4$ and $\sigma =0.5$. For $\sigma =0.6$, SOAS’s value agrees with the numerical result to five decimal places. For $\sigma =0.7$ and $\sigma =0.8$, the present values agree with the numerical results to four decimal places. For $\sigma =0.9$ and $\sigma =1$, the current results agree with the numerical ones to three and two decimal places, respectively. It is important to note that if the exact solution of the nonlinear SIR model is not available, the numerical solution is considered as an optimal solution for the purpose of estimating the accuracy of the present results.

5. Numerical Results and Comparisons with Other Methods

In Ref. [7], the authors applied HPM to systems (1–3) to obtain the following approximate solutions:

$$\begin{array}{c}R\left(t\right)=A-B{e}^{-\tau }+B+\sigma B\left[-B\left(-\tau e^{-\tau }-e^{-\tau }\right)-\frac{1}{2}B{e}^{-2\tau }+B{e}^{-\tau }-\tau e^{-\tau }-e^{-\tau }\right]-\hfill \\ \sigma B\left(\frac{3}{2}B-1\right),\hfill \end{array}$$

(42)

$$\begin{array}{c}I\left(t\right)=B{e}^{-\tau }+e^{-\tau }\left[-\sigma B\left(B\tau -B{e}^{-\tau }-\tau \right)-\sigma B^2\right]-\frac{1}{2}\sigma B{e}^{-\tau }(4\sigma B^2-2\sigma B+2B)-\hfill \\ \frac{1}{2}\sigma B{e}^{-\tau }[4\sigma B^2\tau e^{-\tau }+6\sigma B^2{e}^{-\tau }-2\sigma B^2{e}^{-2\tau }-\sigma B^2{\tau }^2-2\sigma B^2{\tau }^{-}2\sigma B{e}^{-\tau }+2\sigma B\tau +\hfill \\ 2\sigma B{\tau }^2-4\sigma B\tau e^{-\tau }+2B\tau -2B{e}^{-\tau }-\sigma {\tau }^2].\hfill \end{array}$$

(43)

The comparisons between the present FOAS, HPM [7], and numerical solution (Runge–Kutta) are displayed in Figure 5 and Figure 6 for the instantaneous $R\left(\tau \right)$ and $I\left(\tau \right)$, respectively, for the initial recovered individuals $A=0.001$, initial infected individuals $B=0.01$, and transmission rate $\sigma =0.5$. The results in these figures show the agreement of the present FOAS with the numerical solution, while HPM [7] deviates; hence, it may need revision. Figure 7 and Figure 8 also indicate the superiority of the present FOAS over HPM [7]; however, the accuracy of FOAS is slightly affected by the small increase in $\sigma$ ($\sigma =0.7$). Figure 9 and Figure 10 confirm this conclusion that the increase in $\sigma$ ($\sigma =0.8$) leads to a slight deviation between the present FOAS and the numerical solution; however, the accuracy of FOAS is still better than HPM [7].

Although the FOAS provides accurate solutions in a short range of the transmission rate parameter $\sigma$, one can see in Figure 11, Figure 12, Figure 13 and Figure 14 that SOAS is in full agreement with the numerical solution, even for higher values of $\sigma$, $\sigma =0.9$ (Figure 11 and Figure 12) and $\sigma =1$ (Figure 13 and Figure 14). This leads to the conclusion that the present SOAS represents an accurate solution for the current COVID-19 model for all values of $\sigma$ in the range of $0<\sigma \le 1$ while HPM [7] does not. Finally, the variations in recovered individuals $R\left(\tau \right)$ and infected individuals $I\left(\tau \right)$ are plotted in Figure 15 and Figure 16 at various values of the transmission rate $\sigma$.

As the main point is to compare the solution between different methods, the absolute difference between the present approximations FOAS/SOAS and the Runge–Kutta method is plotted in Figure 17 for the recovered individuals, while the absolute difference between the HPM [7] and Runge-Kutta method is displayed in Figure 18. It can be seen from these two figures that the accuracy of the approximation FOAS is slightly better than HPM [7]. However, the approximation SOAS displayed a better accuracy than HPM [7]. This conclusion can also be confirmed or observed in Figure 19 and Figure 20.

6. Conclusions

In this paper, two different approximate solutions were obtained for a nonlinear COVID-19 model. The obtained approximate solutions were expressed explicitly in terms of exponential and hyperbolic functions. It was shown that the current two approximate solutions displayed many advantages over the HPM solution in the literature [7]. The numerical comparisons between our approximations and the Runge–Kutta method were declared. Moreover, it was proven that the present approximations were more accurate compared with HPM [7], where the Runge–Kutta numerical method was taken as a reference method. Furthermore, the values of the recovered individuals at a sufficiently large time were derived analytically through the two approximations that were obtained. The numerical results reveal that our approach is straightforward and efficient compared with the previous study [7].