A simple computational approach to the Susceptible-Infected-Recovered (SIR) epidemic model via the Laplace-Adomian Decomposition Method
aa r X i v : . [ q - b i o . P E ] J un A simple computational approach to the Susceptible-Infected-Recovered (SIR)epidemic model via the Laplace-Adomian Decomposition Method
Tiberiu Harko
1, 2, 3, ∗ and Man Kwong Mak Astronomical Observatory, 19 Ciresilor Street, 400487 Cluj-Napoca, Romania, Faculty of Physics, Babes-Bolyai University, 1 Kogalniceanu Street, 400084 Cluj-Napoca, Romania School of Physics, Sun Yat-Sen University, Xingang Road, 510275 Guangzhou, People’s Republic of China Departamento de Fisica, Facultad de Ciencias Naturales,Universidad de Atacama, Copayapu 485, Copiapo, Chile † The Susceptible-Infected-Recovered (SIR) epidemic model is extensively used for the study of thespread of infectious diseases. Even that the exact solution of the model can be obtained in an exactparametric form, in order to perform the comparison with the epidemiological data a simple buthighly accurate representation of the time evolution of the SIR compartments would be very useful.In the present paper we obtain a series representation of the solution of the SIR model by using theLaplace-Adomian Decomposition Method to solve the basic evolution equation of the model. Thesolutions are expressed in the form of infinite series. The series representations of the time evolutionof the SIR compartments are compared with the exact numerical solutions of the model. We findthat there is a good agreement between the Laplace-Adomian semianalytical solutions containingonly three terms, and the numerical results.
Keywords : Susceptible-Infected-Recovered (SIR) epidemic model; Laplace-Adomian Decompo-sition Method; series solution
Contents
I. Introduction II. Series solution of the SIR epidemic model
III. Conclusions References I. INTRODUCTION
The major Covid-19 coronavirus pandemic that began in 2019 in Wuhan, China [1], has again underlined thenecessity of understanding on both quantitative and qualitative levels the outbreak and spread of infectious diseases.This problem has been already studied a long time ago, John Graunt being the first scientist who attempted toevaluate in a methodical way the causes of fatalities during epidemics [2]. His analysis led to a theory that is stillwidely accepted by the modern epidemiologists. The first mathematician who did propose a mathematical modeldescribing the dynamics of an infectious disease was Daniel Bernoulli, who modeled the spread of Variola, which wasrampant at the time, in 1760 [3]. Moreover, in a later study Daniel Bernoulli did present the advantages of inoculation[4], the method first used to immunize an individual against smallpox with material taken from a patient.More recently an important simple (compartmental) deterministic model predicting the spread of epidemic out-breaks was proposed by McKendrick and Kermack in 1927 [5]. This mathematical epidemic model is called theSusceptible-Infected-Recovered (SIR) model, or the xyz model. In the SIR model it is assumed that in the presence ofan infectious disease a population with a fixed number can be compartmentalized into three groups (compartments), ∗ Electronic address: [email protected] † Electronic address: [email protected] denoted (S), (I) and (R), respectively. More precisely, the compartments used in the McKendrick-Kermack model aredefined as follows [5]:a) The (S) (susceptible) compartment consists of the individuals in the total population not yet infected with thedisease at time t , or those susceptible to the disease. The number of individuals in this compartment is denoted x ( t ).b) The (I) (infected) compartment consists of the individuals who have been infected with the epidemic disease, andwho can spread the disease to those in the susceptible category. The number of individuals in the (I) compartment isdenoted y ( t ).c) The (R) (recovered) compartment consists of the individuals who have been infected during the epidemics, anddid fully recover. The individuals in the (R) compartment, whose number is denoted by z ( t ), cannot be infectedagain, and they are not able to transmit the disease to others.From a mathematical point of view the SIR model is formulated in terms of a strongly nonlinear system of ordinarydifferential equations. The behavior of the SIR model essentially depends on two constants β and γ that give thetransition rates between compartments.The transition rate between the S (Susceptible) and I (infected) groups is denoted by β . It represents the contactrate, that is, the probability of being infected when the contact between a susceptible and an infectious individualtakes place [6–9].The transition rate between the I (Infected) and R (recovered) compartments is denoted by γ . γ gives the rate ofrecovery or death. If we denote by D the total duration of the infection, then γ = 1 /D , expressing the fact that anindividual experiences one recovery in D units of time [6–9]. Since in the SIR model both β and γ represent transitionrates (probabilities), their range of values is 0 ≤ β ≤ ≤ γ ≤
1, respectively.The SIR model has been investigated from both mathematical and epidemiological point of view in a large number ofstudies. Its exact solution has been obtained in a parametric form in [10]. The SIR model has been used extensively forthe simulation of the dynamics of the Covid-19 pandemic and its impact on society [11–21]. Despite the large numberof mathematical and numerical investigations of the SIR model, including the use of the Adomian decompositionmethod [22], of the variational iteration method [23], of the homotopy perturbation method [24], and of the differentialtransformation method [25], a simple but precise representation of the solution of the SIR model that could beefficiently used in epidemiological and populations studies is still missing.It is the purpose of this work to present a simple series solution of the SIR epidemic model. To achieve this goal wefirst derive the basic equation of the SIR model, which is represented by a first order nonlinear differential equationfor z ( t ). Even that this equation can be solved exactly, in a parametric form, we will investigate it by using theLaplace-Adomian Decomposition Method (LADM), which is an extension of the standard Adomian DecompositionMethod [26–28]. The ADM is a powerful mathematical method that permit to find semianalytical solutions of differenttypes of ordinary, partial and integral differential nonlinear equations. The ADM allows to obtain the solution of adifferential equation in the form of a series, with the terms of the series determined recursively by using the Adomianpolynomials [26, 27]. The Laplace-Adomian Method has been used to investigate nonlinear differential equations indifferent fields of science in [29–35]. By applying the Laplace-Adomian Method to the basic equation of the SIR modelits solution can be represented as a series that gives an excellent description of the numerical results.The present paper is organized as follows. The series solution of the SIR epidemic model is obtained in Section II,where its comparison with the exact numerical solution is also performed. We briefly conclude our results in Section III. II. SERIES SOLUTION OF THE SIR EPIDEMIC MODEL
In the present Section we will present a semianalytical solution of the SIR model represented in the form of a seriescontaining exponential terms. The basic equations of the SIR model are dxdt = − βx ( t ) y ( t ) , (1) dydt = βx ( t ) y ( t ) − γy ( t ) , (2)and dzdt = γy ( t ) , (3)respectively, where x ( t ) > y ( t ) > z ( t ) > ∀ t ≥
0. The system of equations (1)-(3) must be integrated withthe initial conditions x (0) = N ≥ y (0) = N ≥ z (0) = N ≥
0, respectively, where N i ∈ ℜ , i = 1 , , β , and the meanrecovery rate γ , which in the following are assumed to be positive constants.By adding Eqs. (1)–(3) we immediately obtain the differential equation, ddt [ x ( t ) + y ( t ) + z ( t )] = 0 , (4)which upon integration gives the first integral of the SIR system x ( t ) + y ( t ) + z ( t ) = N, ∀ t ≥ , (5)Hence it follows that in this epidemiological model the total size of the population N = P i =1 N i is, from a math-ematical point of view, an arbitrary positive integration constant. This result is consistent with our assumption ofconsidering only three compartments in a fixed population with N members. A. The basic evolution equation for the SIR model
We will now show that the SIR model can be equivalently formulated in terms of a first order nonlinear differentialequation. After differentiating Eq. (1) with respect to the time t we obtain the following second order differentialequation, dydt = − β " x ′′ x − (cid:18) x ′ x (cid:19) , (6)where in the following we denote by a prime the derivative with respect to time t . By inserting Eqs. (1) and (6) intoEq. (2), we obtain x ′′ x − (cid:18) x ′ x (cid:19) + γ x ′ x − βx ′ = 0 . (7)From Eqs. (1) and (3) we find the evolution equation for z , in terms of x as dzdt = − γβ (cid:18) x ′ x (cid:19) . (8)which can be integrated to give x = x e − ( β/γ ) z , (9)The value of the positive integration constant x can be obtained from the initial conditions of the SIR model as x = N e ( β/γ ) N . (10)With the use of Eq. (7), Eq. (6) becomes dydt = − x ′ + γβ x ′ x , (11)immediately giving y ( t ) = γβ ln x − x + C , (12)where C is an arbitrary integration constant, which can be determined from the initial conditions as C = N + N − γβ ln N . (13)After differentiating Eqs. (9) and (8) with respect to the time t , we obtain the basic second order differentialequation for z ( t ), describing the spread of the non-fatal disease in a given population, z ′′ = x βz ′ e − βγ z − γz ′ . (14)Eq. (14), a strongly nonlinear differential equation, is mathematically equivalent to the first order system of differ-ential equations Eqs. (1)–(3), respectively, giving the SIR epidemic model. However, a first order nonlinear differentialequation can also de obtained to describe the evolution of z . Eq. (14) can be rewritten as d zdt + γ dzdt = − γx ddt e − ( β/γ ) z , (15)and thus we obtain its first integral as given by dzdt + γz = − γx e − ( β/γ ) z + C, (16)where C is an arbitrary constant of integration. By estimating the above equation at t = 0 and with the help ofEq. (3) we obtain the value of the integration constant as C = γ ( N + N ) + γx e − βγ N = γN. (17) B. Series solution of the basic equation of the SIR model via the Laplace-Adomian method
The general solution of Eq. (16) can be obtained in a closed form, in an exact parametric representation. However,the parametric representation is not particularly useful from a computational point of view. That’s why in thefollowing we will obtain a series solution of Eq. (16). In order to do this we will use the Laplace-Adomian method.In the Laplace-Adomian method we first apply the Laplace transformation operator L x , defined as L x [ f ( x )]( s ) = R ∞ f ( x ) e − sx dx , to Eq. (16), thus obtaining L t (cid:20) dz ( t ) dt (cid:21) ( s ) + γ L t [ z ( t )]( s ) − L t [ C ]( s ) = − γx L t h e − ( β/γ ) z ( t ) i ( s ) . (18)By using the properties of the Laplace transform we immediately find( γ + s ) L t [ z ( t )]( s ) − γN + sz (0) s = − γx L t h e − ( β/γ ) z ( t ) i ( s ) . (19)With the use of the initial condition z (0) = N we obtain L t [ z ( t )]( s ) = γN + N ss ( γ + s ) − γx γ + s L t h e − ( β/γ ) z ( t ) i ( s ) . (20)We assume now that the function z ( t ) can be represented in the form of an infinite series, z ( t ) = ∞ X n =0 z n ( t ) , (21)where the terms z n ( t ) can be computed recursively. As for the nonlinear operator f ( z ( t )) = − γx e − ( β/γ ) z ( t ) , weassume that it can be decomposed as f ( z ( t )) = − γx e − ( β/γ ) z ( t ) = ∞ X n =0 A n ( t ) , (22)where the A n ’s are Adomian polynomials, defined generally for an arbitrary function f ( z ( t )) according to [28] A n = 1 n ! d n dλ n f ∞ X i =0 λ i z i !(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ =0 . (23)The first five Adomian polynomials can be obtained in the following form, A = f ( z ) , (24) A = z f ′ ( z ) , (25) A = z f ′ ( z ) + 12 z f ′′ ( z ) , (26) A = z f ′ ( z ) + z z f ′′ ( z ) + 16 z f ′′′ ( z ) , (27) A = z f ′ ( z ) + (cid:20) z + z z (cid:21) f ′′ ( z ) + 12! z z f ′′′ ( z ) + 14! z f (iv) ( z ) . (28)Substituting Eqs. (21) and (22) into Eq. (20) we obtain L t " ∞ X n =0 z n ( t ) = γN + N ss ( γ + s ) − γ + s L t [ ∞ X n =0 A n ] . (29)The matching of both sides of Eq. (29) gives the following iterative algorithm for the power series solution of thebasic evolution equation of the SIR model, Eq. (16), L t [ z ( t )] ( s ) = γN + N ss ( γ + s ) , (30) L t [ z ( t )] ( s ) = − γ + s L t [ A ( t )] ( s ) , (31) L t [ z ( t )] ( s ) = − s + γ L t [ A ( t )] ( s ) , (32) ... L t [ z k +1 ( t )] ( s ) = − s + γ L t [ A k ( t )] ( s ) . (33)By applying the inverse Laplace transformation to Eq. (30), we obtain the value of z as z ( t ) = L − t (cid:26)(cid:20) γN + N ss ( γ + s ) (cid:21) ( s ) (cid:27) ( t ) = N − ( N + N ) e − γt . (34)In the following, for computational simplicity, we will approximate z ( t ) by its first order series expansion, given by z ( t ) ≈ N + γ ( N + N ) t. (35)The first Adomian polynomial is given by A ( t ) = − γx e − ( β/γ ) z ( t ) , and thus, within the adopted approximation,we obtain A ( t ) = − γN e − β ( N + N ) t , (36)and z ( t ) = −L − t (cid:26) s + γ L t h γN e − β ( N + N ) t i ( s ) (cid:27) ( t ) = γN (cid:2) e − β ( N + N ) t − e − γt (cid:3) β ( N + N ) − γ , (37)respectively. The second Adomian polynomial is obtained as A ( t ) = − βγN e − [ γ +2 β ( N + N )] t (cid:2) e β ( N + N ) t − e γt (cid:3) β ( N + N ) − γ , (38)and thus we find z ( t ) = L − t (cid:26) s + γ L t [ A ( t )] ( s ) (cid:27) ( t ) = γN e − [ γ +2 β ( N + N )] t (cid:8) − β ( N + N ) e γt + [2 β ( N + N ) − γ ] e β ( N + N ) t + [ γ − β ( N + N )] e β ( N + N ) t (cid:9) ( N + N ) [ β ( N + N ) − γ ] [2 β ( N + N ) − γ ] . (39)Similarly, after computing the third Adomian polynomial, we obtain z ( t ) = γN γ − β ( N + N )] ( β [4 β ( N + N ) − γ ] e − β ( N + N ) t γ + 6 β ( N + N ) − βγ ( N + N ) + β e − [2 γ + β ( N + N )] t γ + β ( N + N ) +2 [ γ − β ( N + N )] e − [ γ + β ( N + N )] t ( N + N ) [2 β ( N + N ) − γ ] + [ γ − β ( N + N )] e − [ γ +2 β ( N + N )] t ( N + N ) − e − γt [ γ − β ( N + N )] [ γ + 2 β ( N + N )]( N + N ) [3 β ( N + N ) − γ ] [ γ + β ( N + N )] ) . (40)Finally, we will present the z ( t ) term in the Adomian series expansion, which is given by z ( t ) = 16 γN ( β [8 γ − β ( N + N )] e − β ( N + N ) t [ β ( N + N ) − γ ] [3 β ( N + N ) − γ ] [4 β ( N + N ) − γ ] + β e − [3 γ + β ( N + N )] t [ β ( N + N ) − γ ] [2 γ + β ( N + N )] +6 β h γ − β ( N + N ) i e − γ + β ( N + N )] t ( N + N ) [ β ( N + N ) − γ ] [ γ + β ( N + N )] [ γ + 2 β ( N + N )] +6 β e − [2 γ + β ( N + N )] t ( N + N ) [ β ( N + N ) − γ ] [2 β ( N + N ) − γ ] [ γ + β ( N + N )] + e − γt h − γ − β ( N + N ) − β γ ( N + N ) − βγ ( N + N ) i ( N + N ) [4 β ( N + N ) − γ ] [ γ + β ( N + N )] [ γ + 2 β ( N + N )] [2 γ + β ( N + N )] + h − γ + 12 β ( N + N ) − β γ ( N + N ) + 7 βγ ( N + N ) i e − [ γ +3 β ( N + N )] t ( N + N ) [ β ( N + N ) − γ ] [2 β ( N + N ) − γ ] +[3 γ + 6 β ( N + N )] e − [ γ + β ( N + N )] t ( N + N ) [3 β ( N + N ) − γ ] [ γ + β ( N + N )] − e − [ γ +2 β ( N + N )] t ( N + N ) [ β ( N + N ) − γ ] ) . (41)The higher order terms in the Laplace-Adomian power series representation of z ( t ) can be easily obtained with thehelp of some symbolic computation software. Hence, the solution of the SIR system can be obtained as x ( t ) = x e − βγ P ∞ i =0 z i ( t ) = N e ( β/γ ) ( N − P ∞ i =0 z i ( t ) ) , (42) y ( t ) = N + N + γβ ln x e − βγ P ni =0 z i ( t ) N − x e − βγ P ni =0 z i ( t ) = N + N + N − ∞ X i =0 z i ( t ) − N e ( β/γ ) ( N − P ∞ i =0 z i ( t ) ) , (43)and z ( t ) = ∞ X i =0 z i ( t ) , (44)respectively. C. Comparison with the exact numerical solution
In order to test the applicability of the power series solutions solving the SIR model we will compare the truncatedLaplace-Adomian series with the exact numerical solution. To perform the comparison we will consider only the first x ( t ) , y ( t ) , z ( t ) x ( t ) , y ( t ) , z ( t ) FIG. 1: Comparison of the exact numerical solution of the SIR model and of the Laplace-Adomian series expansion truncatedto three terms, for β = 0 . γ = 0 . β = 0 .
098 and γ = 0 .
164 (right figure). x ( t ) is represented bythe solid curve, y ( t ) by the dotted curve, and z ( t ) by the dashed curve. The initial conditions of the SIR system are x (0) = 25, y (0) = 15, and z (0) = 10.
40 t x ( t ) , y ( t ) , z ( t ) x ( t ) , y ( t ) , z ( t ) FIG. 2: Comparison of the exact numerical solution of the SIR model and of the Laplace-Adomian series expansion truncatedto three terms, for β = 0 .
119 and γ = 0 .
017 (left figure) and β = 0 . γ = 0 .
584 (right figure). x ( t ) is represented bythe solid curve, y ( t ) by the dotted curve, and z ( t ) by the dashed curve. The initial conditions of the SIR system are x (0) = 25, y (0) = 15, and z (0) = 10. three terms in the expansion of z , so that z ( t ) ≈ z ( t ) + z ( t ) + z ( t ), or, explicitly, z ( t ) ≈ N − ( N + N ) e − γt + γN (cid:2) e − β ( N + N ) t − e − γt (cid:3) β ( N + N ) − γ + γN e − [ γ +2 β ( N + N )] t × (cid:8) − β ( N + N ) e γt + [2 β ( N + N ) − γ ] e β ( N + N ) t + [ γ − β ( N + N )] e β ( N + N ) t (cid:9) ( N + N ) [ β ( N + N ) − γ ] [2 β ( N + N ) − γ ] . (45)We fix the initial conditions of the SIR system as x (0) = N = 25, y (0) = N = 15, and z (0) = N = 10, and we willvary the values of the parameters β and γ . The results of the comparison between the exact numerical solution andthe truncated series solutions are represented, for four distinct sets of values of ( β, γ ), in Figs. 1 and 2, respectively.As one can see from the Figures, the three terms truncated Laplace-Adomian series gives a very good description ofthe time dynamics and evolution of the SIR epidemiological model. If necessary, the precision of the semi-analyticalsolution can be easily increased by adding more terms in the series expansion. III. CONCLUSIONS
Even that the exact solution of the SIR model is known, a simple but highly accurate representation of the modelsolution would significantly simplify the comparison of the model predictions with the epidemiological data, alsoallowing an efficient determination of the model parameters β and γ from the data fitting. The exact solutionof the SIR model is given in a parametric form by x = x u , y = N + ( γ/β ) ln u − x u , z = − ( γ/β ) ln u , and t − t = R uu dξ/ [ ξ ( x β ln ξ − γξ − N )] [10], and it requires the numerical evaluation of the time integral for each step.In the present paper we have presented a series solution of the three compartment SIR epidemiological model, givingthe explicit time dependence of x ( t ), y ( t ) and z ( t ). Despite its mathematical simplicity, the semianalytical solutiondescribes very precisely the numerical behavior of the model for arbitrary values of β and γ . The time dependence ofthe compartments in the semianalytical solution is given by the sum of exponential terms, and such a representationmay facilitate the fitting with the statistical results. The series truncated to three terms only already give a verygood description of the numerical results for arbitrary values of β and γ , and of the initial conditions. In fact the twoterms approximation also gives a good approximation of the results of the numerical integration of the SIR model.Exact solutions of the epidemiological models are important because they can be used by epidemiologists to studythe evolution of infectious diseases in different environments, and to develop the best social strategies for their control.Hopefully the results of the present paper could also contribute to the better understanding of the dynamics of thepresent and of the future epidemics. [1] C. Huang, Y. Wang, X. Li, L. Ren, J. Zhao, Y. Hu et al., Clinical features of patients infected with 2019 novel coronavirusin Wuhan, China , Lancet , 497-506 (2020).[2] John Graunt,
Natural and political observations made upon the bills of mortality (1662).[3] D. Bernoulli,
Essai d’une nouvelle analyse de la mortalite causee par la petite verole et des avantages de l’inoculation pourla prevenir , Mem. Math. Phys. Acad. Roy. Sci., Paris, 1 (1760).[4] D. Bernoulli,
Reflexions sur les avantages de l’inoculation , Mercure de France, June issue, (1760).[5] W. O. Kermack and A. G. McKendrick,
Contribution to the mathematical theory of epidemics , Proc. Roy. Soc. Lond
A115 , 700-721 (1927).[6] F. Brauer and C. Castillo-Ch´avez,
Mathematical Models in Population Biology and Epidemiology , Springer, New York(2001).[7] J. D. Murray,
Mathematical Biology: I. An Introduction , Springer-Verlag, New York, Berlin, Heidelberg (2002).[8] D. J. Daley and J. Gani,
Epidemic Modeling: An Introduction , Cambridge University Press, Cambridge (2005).[9] F. Brauer, P. van den Driessche, and J. Wu (Editors),
Lecture Notes in Mathematical Epidemiology , Springer-Verlag,Berlin, Heidelberg (2008).[10] T. Harko, F. S. N. Lobo, and M. K. Mak,
Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemicmodel and of the SIR model with equal death and birth rates , Applied Mathematics and Computation , 184-194 (2014).[11] A. A. Toda,
Susceptible-Infected-Recovered (SIR) Dynamics of COVID-19 and Economic Impact , arXiv:2003.11221 [q-bio.PE] (2020).[12] J. Dehning, J. Zierenberg, F. P. Spitzner, M. Wibral, J. P. Neto, M. Wilczek, and V. Priesemann,
Inferring change pointsin the COVID-19 spreading reveals the effectiveness of interventions , arXiv:2004.01105 [q-bio.PE] (2020).[13] M. Cadoni,
How to reduce epidemic peaks keeping under control the time-span of the epidemic , Chaos, Solitons & Fractals , 109940 (2020).[14] J. N. Dhanwant and V. Ramanathan,
Forecasting COVID 19 growth in India using Susceptible-Infected-Recovered (SIR)model , arXiv:2004.00696 [q-bio.PE] (2020).[15] A. L. Bertozzi, E. Franco, G. Mohler, M. B. Short, and D. Sledge,
The challenges of modeling and forecasting the spreadof COVID-19 , arXiv:2004.04741 [q-bio.PE] (2020).[16] E. Franco,
A feedback SIR (fSIR) model highlights advantages and limitations of infection-based social distancing ,arXiv:2004.13216 [q-bio.PE] (2020).[17] S. Y. Lee, B. Lei, and B. K. Mallick,
Estimation of COVID-19 spread curves integrating global data and borrowing infor-mation , arXiv:2005.00662v3 [stat.AP] (2020).[18] S. A. Hojman and F. A. Asenjo,
Phenomenological dynamics of COVID-19 pandemic: meta-analysis for adjustmentparameters , arXiv:2005.08686v2 [physics.soc-ph] (2020).[19] B. K. Beare and A. A. Toda,
On the Emergence of a Power Law in the Distribution of COVID-19 Cases , arXiv:2004.12772v1[physics.soc-ph] (2020).[20] D. I. Ketcheson,
Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention ,arXiv:2004.08848v2 [math.OC] (2020).[21] D.-V. Anghel and I. T. A. Anghel,
Understanding the COVID19 infection curves – finding the right numbers , DOI:10.21203/rs.3.rs-32516/v1 (2020).[22] J. Biazar,
Solution of the epidemic model by Adomian decomposition method , Applied Mathematics and Computation ,1101-1106 (2006).[23] M. Rafei, H. Daniali and D. D. Ganji,
Variational iteration method for solving the epidemic model and the prey and predatorproblem , Applied Mathematics and Computation , 1701-1709 (2007).[24] M. Rafei, D. D. Ganji and H. Daniali,
Solution of the epidemic model by homotopy perturbation method , Applied Mathe-matics and Computation, [25] A.-M. Batiha and B. Batiha,
A new method for solving epidemic model , Australian J. of Basic and Applied Sciences ,3122-3126 (2011).[26] G. Adomian, A review of the decomposition method in applied mathematics , J. Math. Anal. Appl. , 501-544 (1988).[27] G. Adomian,
Solving Frontier Problems of Physics: the Decomposition Method , Kluwer, Dordrecht, (1994).[28] G. Adomian and R. Rach,
Modified Adomian Polynomials , Mathematical and Computer Modelling , 39-46 (1996).[29] A.-M. Wazwaz, A comparison between the variational iteration method and Adomian decomposition method , Journal ofComputational and Applied Mathematics , 129-136 (2007).[30] A.-M. Wazwaz, R. Rach, and J.-S. Duan,
Solving New Fourth-Order Emden-Fowler-Type Equations by the AdomianDecomposition Method , International Journal for Computational Methods in Engineering Science and Mechanics , 121-131 (2015).[31] J.-S. Duan, R. Rach, and A.-M. Wazwaz, A reliable algorithm for positive solutions of nonlinear boundary value problemsby the multistage Adomian decomposition method , Open Engineering , id.7 (2015).[32] M. K. Mak, C. S. Leung, and T. Harko, Computation of the general relativistic perihelion precession and of light deflectionvia the Laplace-Adomian Decomposition Method , Advances in High Energy Physics , 7093592 (2018).[33] H. O. Bakodah, A. Ebaid, and A.-M. Wazwaz,