# A Note on Solutions of the SIR Models of Epidemics Using HAM

Authors: M. Sajid; Z. Abbas; N. Ali; T. Javed
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

Recently, Awawdeh et al. (2009) discussed the solutions of SIR epidemics model using homotopy analysis method. This comment points out some crucial flaws in (Awawdeh et al. 2009). Particularly, results presented in Figure 1 of the (Awawdeh et al. 2009) do not represent the 20 term solution of the considered problem as stated. The present paper also provides a new approach for solving SIR epidemics model using homotopy analysis method. The new approach is based on dividing the entire domain into subintervals. In each subinterval the three-term HAM solution is sufficient for obtaining accurate and convergent results. The comparison of the obtained solution using new approach is made with the numerical results and found in excellent agreement.


---

## Body

## 1. Introduction

The authors in [1] have considered the classic SIR epidemic model for the homotopy analysis method (HAM) solution. The constant population in SIR model is divided into susceptible, infectious, and recovered classes [2, 3]. The expressions for the susceptible, infectious, and recovered population presented in [1] for five- and nine-term HAM solutions clearly indicate that, for t=0, the initial population size is 20, 15, and 10, respectively. However, Figure 1 illustrates a different population size. The same analytic expressions for the susceptible, infectious, and recovered populations for the same problem are presented in the literature with Adomian decomposition method [4], variational iteration method [5], homotopy perturbation method [6], and differential transform method [7]. The graphical results in all the studies [4–7] reveal that the results are valid for small values of time and beyond that these analytic expressions fail to give results that can be compared with the existing numerical solutions. However, the authors in [1] obtained same expressions as presented in [4–7], but the graphical results are different from all these studies which are not possible. Thus Figure 1 present in [1] is not found through 20-term HAM solution.Figure 1
Plot for susceptible populationS(t) when NS=499, NI=1, and NR=1.The objective of the present paper is to revisit the HAM solution for the nonlinear initial value problems considered by Awawdeh et al. [1] and to provide a HAM solution which agrees well with the numerical results already obtained for the same problem. For the accurate and convergent HAM solution valid for all ranges of time values, we propose a new approach. In this new approach the whole domain of the problem is divided into subintervals, and the HAM solution is evaluated separately in each subinterval. A three-term HAM approximation in each subinterval is sufficient for an accurate and convergent solution.


## 2. SIR Model and Solution

The classic SIR epidemic model is given by the following system of nonlinear ordinary differential equations [1]:
(1)dSdt=-βSI,dIdt=βSI-γI,dRdt=γI,
with initial conditions
(2)S(0)=NS,I(0)=NI,R(0)=NR.


### 2.1. Solution by Homotopy Analysis Method

For the details of the homotopy analysis method for the initial value problems considered in (1)–(2), we refer the readers to [1] and only present the mth order deformation equations with h=-1 which is of interest for us and is given by
(3)ddt[Sm(t)-χmSm-1(t)]=-dSm-1(t)dt-β∑k=0m-1Sk(t)Im-1-k(t),ddt[Im(t)-χmIm-1(t)]=-dIm-1(t)dt-γIm-1(t)+β∑k=0m-1Sk(t)Im-1-k(t),ddt[Rm(t)-χmRm-1(t)]=-dRm-1(t)dt+γIm-1(t),Sm(0)=0,Im(0)=0,Rm(0)=0,χm={0,m≤1,1,m>1.
The final solution is given by
(4)S=∑m=0MSm,I=∑m=0MIm,R=∑m=0MRm,
in which M represents the number of iterations.  Awawdeh et al. [1] solved the system of (3) directly to obtain the solution at different orders for m=1,2,3,…,M. Solutions obtained thus have problems with the accuracy and convergence for large values of t even if one evaluates solution for hundred iterations. In the proposed approach for improving the accuracy and rate of convergence instead of directly solving the system of (3), the entire domain of t is split into subintervals. The mechanism of the proposed approach is discussed in next subsection.


### 2.2. Domain Split and Solution by New Approach

In the new approach, we divide the domain0≤t≤L into subintervals each represented by Hi such that
(5)∑Hi=L,i=1,2,3,….
For fixed length subintervals Hi=iH,i=1,2,3,…. The problem in each subinterval takes the form
(6)ddt[Smi(t)-χmSm-1i(t)]=-dSm-1i(t)dt-β∑k=0m-1Ski(t)Im-1-ki(t),(7)ddt[Imi(t)-χmIm-1i(t)]=-dIm-1i(t)dt-γIm-1i(t)+β∑k=0m-1Ski(t)Im-1-ki(t),(8)ddt[Rmi(t)-χmRm-1i(t)]=-dRm-1i(t)dt+γIm-1i(t),(9)S01=NS,I01=NI,R01=NR.
The final solution in each subdomain is thus given by
(10)Si=∑m=0MSmi,Ii=∑m=0MImi,Ri=∑m=0MRmi.
Choosing i=1, the system of (6)–(9) is solved for three terms such that 0≤t≤H. The numerical values of this solution at the end point of the interval, that is, t=H, becomes the initial values for the second subinterval when i=2, and therefore
(11)S02=S1(H),I02=I1(H),R02=R1(H).
System of (6)–(8) is solved again for i=2 subject to initial conditions (11), and we obtain the initial conditions for the third subinterval. In this manner, the numerical values obtained at the final point of each subinterval are the values at the initial point of the next subinterval. Hence at the ith subinterval, the initial values are
(12)S0i=Si-1((i-1)H),I0i=Ii-1((i-1)H),(13)R0i=Ri-1((i-1)H).
Following this procedure, the system of (6)–(8) are subsequently solved from one subinterval to the next subinterval. The convergence and accuracy of the obtained solutions are ensured by keeping the length of subinterval small.


## 2.1. Solution by Homotopy Analysis Method

For the details of the homotopy analysis method for the initial value problems considered in (1)–(2), we refer the readers to [1] and only present the mth order deformation equations with h=-1 which is of interest for us and is given by
(3)ddt[Sm(t)-χmSm-1(t)]=-dSm-1(t)dt-β∑k=0m-1Sk(t)Im-1-k(t),ddt[Im(t)-χmIm-1(t)]=-dIm-1(t)dt-γIm-1(t)+β∑k=0m-1Sk(t)Im-1-k(t),ddt[Rm(t)-χmRm-1(t)]=-dRm-1(t)dt+γIm-1(t),Sm(0)=0,Im(0)=0,Rm(0)=0,χm={0,m≤1,1,m>1.
The final solution is given by
(4)S=∑m=0MSm,I=∑m=0MIm,R=∑m=0MRm,
in which M represents the number of iterations.  Awawdeh et al. [1] solved the system of (3) directly to obtain the solution at different orders for m=1,2,3,…,M. Solutions obtained thus have problems with the accuracy and convergence for large values of t even if one evaluates solution for hundred iterations. In the proposed approach for improving the accuracy and rate of convergence instead of directly solving the system of (3), the entire domain of t is split into subintervals. The mechanism of the proposed approach is discussed in next subsection.


## 2.2. Domain Split and Solution by New Approach

In the new approach, we divide the domain0≤t≤L into subintervals each represented by Hi such that
(5)∑Hi=L,i=1,2,3,….
For fixed length subintervals Hi=iH,i=1,2,3,…. The problem in each subinterval takes the form
(6)ddt[Smi(t)-χmSm-1i(t)]=-dSm-1i(t)dt-β∑k=0m-1Ski(t)Im-1-ki(t),(7)ddt[Imi(t)-χmIm-1i(t)]=-dIm-1i(t)dt-γIm-1i(t)+β∑k=0m-1Ski(t)Im-1-ki(t),(8)ddt[Rmi(t)-χmRm-1i(t)]=-dRm-1i(t)dt+γIm-1i(t),(9)S01=NS,I01=NI,R01=NR.
The final solution in each subdomain is thus given by
(10)Si=∑m=0MSmi,Ii=∑m=0MImi,Ri=∑m=0MRmi.
Choosing i=1, the system of (6)–(9) is solved for three terms such that 0≤t≤H. The numerical values of this solution at the end point of the interval, that is, t=H, becomes the initial values for the second subinterval when i=2, and therefore
(11)S02=S1(H),I02=I1(H),R02=R1(H).
System of (6)–(8) is solved again for i=2 subject to initial conditions (11), and we obtain the initial conditions for the third subinterval. In this manner, the numerical values obtained at the final point of each subinterval are the values at the initial point of the next subinterval. Hence at the ith subinterval, the initial values are
(12)S0i=Si-1((i-1)H),I0i=Ii-1((i-1)H),(13)R0i=Ri-1((i-1)H).
Following this procedure, the system of (6)–(8) are subsequently solved from one subinterval to the next subinterval. The convergence and accuracy of the obtained solutions are ensured by keeping the length of subinterval small.


## 3. Results and Discussion

The procedure discussed in Section2.2 is implemented using the symbolic software Mathematica. The system of (6)–(8) in each subinterval is solved by using the same set of initial values of susceptible, infectious, and recovered population as given in [1]. The Runge-Kutta method is used to solve the initial value problems given in (1)-(2) numerically. The obtained numerical results are shown by bullets in the graphs, the twentieth order HAM solution developed by the standard procedure in [1] is given by the dotted line, and the HAM solution with the new approach is displayed through the solid line in graphs. Figures 1, 2, and 3 present the comparison between the solution given in [1], the numerical solution, and solution with the new approach. Figure 1 elucidate that a twentieth order standard HAM solution (dotted line) is valid for 0≤t≤10. In this case the chosen values of parameters β and γ are small, and for large values the domain becomes further short for the accurate and convergent solution. Beyond this time, the population is going in the negative direction which is not physically possible. The similar observations are made for the infectious and recovered populations and are presented in Figures 2 and 3. Figure 4 is made to show the comparison of our new approach HAM solution with the numerical solution. It is evident from Figure 4 that the results obtained are in excellent agreement with the numerical solutions. It is important to mention here that  Figure 1 given in [1] could not be obtained by the HAM solution presented in [1] and some other solutions are given in Figure 1 of [1].Figure 2
Plot for infectious populationI(t) when NS=499, NI=1, and NR=1.Figure 3
Plot for recovered populationR(t) when NS=499, NI=1, and NR=1.Figure 4
Comparison of the new approach and numerical solutions whenNS=499, NI=1, and NR=1.


## 4. Conclusions

In this paper we present a new approach for solving the SIR epidemic model problem using homotopy analysis method. It is found that this new approach is an effective method for providing numerical and analytic closed form solutions of such problems. The comparison of the obtained results with numerical results is in excellent agreement.



---

*Source: 457072-2013-09-26.xml*