1 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 Β© Global Society of Scientific Research and Researchers http://asrjetsjournal.org/ Numerical Study of Kermack-Mckendrik SIR Model to Predict the Outbreak of Ebola Virus Diseases Using Euler and Fourth Order Runge-Kutta Methods Md. Tareque Hossaina, Md. Musa Miahb*, Md. Babul Hossainc aDepartment of Textile Engineering, City University, Birulia, Savar, Dhaka, Bangladesh b,cDepartment of Mathematics, Mawlana Bhashani Science and Technology University, Santosh, Tangail-1902, Bangladesh aEmail: tareque.ms@gmail.com, bEmail: musa_ju69@yahoo.com, cEmail: babulhossainh@yahoo.com Abstract Mathematical Modeling has emerged as a vital tool for understanding the dynamics of the spread of many infectious diseases, one amongst is Ebola virus. The main focus of this paper is to model mathematically the transmission dynamics of Ebola virus. For this purpose we tend to use basic SIR model of Ebola Virus to predict the outbreak of the diseases. As we cannot fully solve the 3 basic equations of SIR model with a certain formula solution, we introduce Euler and fourth-order Runge-Kutta methods (RK4). These two proposed strategies are quite efficient and practically well suited for solving initial value problem (IVP) for ordinary differential equations (ODE).We discuss the numerical comparisons between Euler method and Runge-Kutta methods and also discuss regarding their performances with the actual data. The population that we used for this model had roughly a similar number of individuals as the number was living in Republic of Liberia during 2014. Keywords: Ebola; Outbreak; Evolution; Mathematical Modeling. 1. Introduction Mathematical models are a powerful tool for investigating human infectious diseases, such as Ebola virus, contributing to the understanding of the dynamics of disease and providing useful predictions about the potential transmission of a disease and the effectiveness of possible control measures, which can provide valuable information for public health policy makers. ------------------------------------------------------------------------ * Corresponding author. http://asrjetsjournal.org/ American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 2 The Ebola virus disease was first discovered in 1976 in the present Democratic Republic of Congo [1]. Since then, there have been many outbreaks; with the greatest was 2014 outbreak [2] which has spread through many countries. According to Bekoe (2015), from the first confirmed case recorded on 23 March 2014 which more than 18 months, at least 11,312 people have been reported died from the disease in six countries; Liberia, Sierra Leone, Guinea, Nigeria, Mali and US. Since the virus keeps spreading through contacts and the mortality rate of 0.7[3], it is needed to understand the patterns and epidemiology of the disease. By these conditions, mathematical model of the outbreaks of Ebola virus can be helpful as it is a platform for understanding the behavior of a dynamical system. The objectives are to understand better the mathematical dynamics of an infected population when an outbreak occurs. Another goal is to use mathematical modeling to examine and to analyze the viral dynamics of the Ebola virus. To model this outbreak, the systems of differential equations are used. To study the known data, several distinct models will be used and each model is different depends on the parameters acquired. I decided to model the Ebola Epidemics in Liberia in 2014 and compare their spread using an SIR model. 1.1 Objectives The purposes of this research are 1. To apply SIR model to predict the outbreaks of Ebola virus. 2. To determine the effect of the initial number of infectives of the population. 3. Compare with real data and fit with model. 1.2 Scope In the proposal study, we will only investigate the outbreaks of Ebola virus by applying the mathematical model. The data has been collected and it covers the area in the continent of Africa, on Liberia that was recorded by CDC [4]. The calculation will be done by using tools of MATLAB software. 2. Methodology 2.1 Formulation of SIR model The SIR model is used to illustrate the transfer of the epidemic through the interaction of the following three different variables: 𝑆𝑆 = Number of people that are susceptible to Ebola 𝐼𝐼 = Number of people infected with Ebola 𝑅𝑅 = Number of people recovered American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 3 from Ebola with total immunity It makes sense to assume that a fixed population of 𝑁𝑁 people, whereby there are no births and deaths by natural cause i.e. 𝑁𝑁 = 𝑆𝑆 + 𝐼𝐼 + 𝑅𝑅 [5] This is because the population is fixed and therefore, there are only three compartments in which the population may fit into. Thus, the total of the number of people susceptible infected and recovered in equivalent to the total population. The assumption that 𝑁𝑁 is fixed, with no births or deaths, makes sense given 60 days, although it is a simplification. These variables change over time, so we will define the variable 𝑑𝑑 = time in days. We will set 𝑑𝑑 = 0 at the start of August 2014. The model uses two parameters 𝛽𝛽 and 𝛾𝛾 with 𝛽𝛽, 𝛾𝛾 > 0. Given these parameters, the model uses 3 differential equations. The rate of change of the number of people susceptible to the disease over time 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = βˆ’ 𝛽𝛽𝐼𝐼𝑆𝑆 (1) The rate of change of the number of people recovered over time 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛾𝛾𝐼𝐼 (2) The rate of change of the number of people infected. 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛽𝛽𝐼𝐼𝑆𝑆 βˆ’ 𝛾𝛾𝐼𝐼 (3) Parameterization of the model In order to calculate 𝛽𝛽 (the rate of infection) and 𝛾𝛾 (the rate of recovery), it helps to define two more parameters. 𝐷𝐷 = Duration of disease for those recovered 𝑀𝑀 = Mortality rate for those who die per day (0.7 for Ebola) This leads to two further equations. American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 4 The rate at which the disease is spread 𝛾𝛾 = 1 𝐷𝐷 [6] (4) The infection rate of the disease 𝛽𝛽 = 𝑀𝑀 𝑑𝑑 [7] (5) 2.2 Transformation of Runge-Kutta Equations for SIR modeling RK4 is one of the classic methods for numerical integration of ODE models. Consider the following initial problem of ODE 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝑓𝑓(𝑑𝑑,𝑦𝑦) 𝑦𝑦(π‘‘π‘‘π‘œπ‘œ) = π‘¦π‘¦π‘œπ‘œ Where 𝑦𝑦(𝑑𝑑) is the unknown function (scalar or vector) which I would like to approximate. The Iterative formula of RK4 method for solving ODE is as follows π‘˜π‘˜1 = β„Žπ‘“π‘“(𝑑𝑑𝑛𝑛,𝑦𝑦𝑛𝑛) π‘˜π‘˜2 = β„Žπ‘“π‘“ �𝑑𝑑𝑛𝑛 + 1 2 β„Ž,𝑦𝑦𝑛𝑛 + 1 2 π‘˜π‘˜1οΏ½ π‘˜π‘˜3 = β„Žπ‘“π‘“ �𝑑𝑑𝑛𝑛 + 1 2 β„Ž,𝑦𝑦𝑛𝑛 + 1 2 π‘˜π‘˜2οΏ½ π‘˜π‘˜4 = β„Žπ‘“π‘“(𝑑𝑑𝑛𝑛 + β„Ž,𝑦𝑦𝑛𝑛 + π‘˜π‘˜3) 𝑦𝑦𝑛𝑛+1 = 𝑦𝑦𝑛𝑛 + 1 6 (π‘˜π‘˜1 + 2π‘˜π‘˜2 + 2π‘˜π‘˜3 + π‘˜π‘˜4) For simplicity, here we use the simplest SIR model to examine whether the RK4 method has been implemented correctly. The SIR model is defined as follows 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = βˆ’ 𝛽𝛽𝐼𝐼𝑆𝑆 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛽𝛽𝐼𝐼𝑆𝑆 βˆ’ 𝛾𝛾𝐼𝐼 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 5 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛾𝛾𝐼𝐼 where 𝑆𝑆(𝑑𝑑) is the number of susceptible people in the population at time 𝑑𝑑, 𝐼𝐼(𝑑𝑑) is the number of infectious people at time 𝑑𝑑, 𝑅𝑅(𝑑𝑑) is the number of recovered people at time 𝑑𝑑, 𝛽𝛽 is the transmission rate, 𝛾𝛾 represents the recovery rate, and 𝑁𝑁 = 𝑆𝑆(𝑑𝑑) + 𝐼𝐼(𝑑𝑑) + 𝑅𝑅(𝑑𝑑) is the fixed population. According to the general iterative formula, the iterative formulas for 𝑆𝑆(𝑑𝑑), 𝐼𝐼(𝑑𝑑) and 𝑅𝑅(𝑑𝑑) of SIR model can be written out 𝑆𝑆𝑛𝑛+1 = 𝑆𝑆𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) π‘˜π‘˜1𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛, 𝑆𝑆𝑛𝑛 , 𝐼𝐼𝑛𝑛) = βˆ’π›½π›½π‘†π‘†π‘›π‘›πΌπΌπ‘›π‘› π‘˜π‘˜2𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝑆𝑆𝑛𝑛 + π‘˜π‘˜1 π‘†π‘†βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜1 πΌπΌβˆ†π‘‘π‘‘ 2 οΏ½ = βˆ’π›½π›½ �𝑆𝑆𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) π‘˜π‘˜3𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ = βˆ’π›½π›½ �𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) π‘˜π‘˜4𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘, 𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘, 𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) = βˆ’π›½π›½(𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) 𝐼𝐼𝑛𝑛+1 = 𝐼𝐼𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) π‘˜π‘˜1𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛, 𝑆𝑆𝑛𝑛, 𝐼𝐼𝑛𝑛) = 𝛽𝛽𝑆𝑆𝑛𝑛𝐼𝐼𝑛𝑛 βˆ’ 𝛾𝛾𝐼𝐼𝑛𝑛 π‘˜π‘˜2𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝑆𝑆𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ = 𝛽𝛽 �𝑆𝑆𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ��𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜3𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 6 = 𝛽𝛽 �𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ��𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜4𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘, 𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘, 𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) = 𝛽𝛽(𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) βˆ’ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) 𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) π‘˜π‘˜1𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛, 𝐼𝐼𝑛𝑛) = 𝛾𝛾𝐼𝐼𝑛𝑛 π‘˜π‘˜2𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ = 𝛾𝛾 �𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜3𝑑𝑑 = 𝑓𝑓 �𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘ 2 , 𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ = 𝛾𝛾 �𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜4𝑑𝑑 = 𝑓𝑓(𝑑𝑑𝑛𝑛 + βˆ†π‘‘π‘‘, 𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) = 𝛾𝛾(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) 2.3 Transformation of Euler Equations for SIR modeling If we have a "slope formula," i.e., a way to calculate 𝑑𝑑𝑦𝑦/𝑑𝑑𝑑𝑑 at any point (𝑑𝑑,𝑦𝑦), then we can generate a sequence of 𝑦𝑦-values, 𝑦𝑦0 ,𝑦𝑦1,𝑦𝑦2,𝑦𝑦3 … … … … By starting from a given 𝑦𝑦0 and computing each rise as slope x run. That is, 𝑦𝑦𝑛𝑛+1 = 𝑦𝑦𝑛𝑛 + 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑛𝑛 βˆ†π‘‘π‘‘ where βˆ†π‘‘π‘‘ is a suitably small step size in the time domain. It really doesn't matter in this calculation if the slope formula happens to depend not just on t and y but on other variables, say x and z -- as long as we know how x and z are related to t and y. If x and z happen to be other dependent variables in a system of differential equations, I can generate values of x and z in the same way. Of course, for the SIR model, we want the dependent variable names to be 𝑆𝑆, 𝐼𝐼 π‘Žπ‘Žπ‘Žπ‘Žπ‘‘π‘‘ 𝑅𝑅. Thus we have three Euler formulas of the form 𝑆𝑆𝑛𝑛+1 = 𝑆𝑆𝑛𝑛 + 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑛𝑛 βˆ†π‘‘π‘‘ American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 7 𝐼𝐼𝑛𝑛+1 = 𝐼𝐼𝑛𝑛 + 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑛𝑛 βˆ†π‘‘π‘‘ 𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑛𝑛 βˆ†π‘‘π‘‘ More specifically, given the SIR equations, 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = βˆ’ 𝛽𝛽𝐼𝐼𝑆𝑆 ; 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛽𝛽𝐼𝐼𝑆𝑆 βˆ’ 𝛾𝛾𝐼𝐼; 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 = 𝛾𝛾𝐼𝐼 The Euler formulas become 𝑆𝑆𝑛𝑛+1 = 𝑆𝑆𝑛𝑛 βˆ’ π›½π›½πΌπΌπ‘›π‘›π‘†π‘†π‘›π‘›βˆ†π‘‘π‘‘ 𝐼𝐼𝑛𝑛+1 = 𝐼𝐼𝑛𝑛 + (𝛽𝛽𝐼𝐼𝑛𝑛𝑆𝑆𝑛𝑛 βˆ’ 𝛾𝛾𝐼𝐼𝑛𝑛)βˆ†π‘‘π‘‘ 𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + π›Ύπ›ΎπΌπΌπ‘›π‘›βˆ†π‘‘π‘‘ To calculate something from these formulas, we must have explicit values for 𝛽𝛽, 𝛾𝛾, 𝑆𝑆(0), 𝐼𝐼(0), 𝑅𝑅(0) π‘Žπ‘Žπ‘Žπ‘Žπ‘‘π‘‘ βˆ†π‘‘π‘‘ . 3. Application, Comparison and Result Discussion If we now take the example of the Ebola outbreak in Liberia 2014, we can assign the parameters with the following values. The total population of Liberia, N = 4294000 [8], and according to data from WHO, the number of people infected, I =391[9] and the number of people dead is 227 [9]. Seeing as R includes the number of people who have received permanent immunity, this includes those who have died as they have permanent immunity, in addition to those who have recovered with permanent immunity. Therefore, number of people recovered 𝑅𝑅 = 227 + (0.3 Γ— 391) = 344 We will now use this data to provide the parameters with the following values. 𝑁𝑁 = 4294000; 𝐼𝐼 = 391; 𝑅𝑅 = 344 Therefore, 𝑆𝑆 = 𝑁𝑁 βˆ’ 𝐼𝐼 + 𝑅𝑅 = 4294000 βˆ’ (344 + 391) = 4293265 The duration of the disease ranges from 2 to 21 days [10], therefore we could roughly estimate the duration of the disease at the midpoint, i.e.11 (approx.) days. 𝐷𝐷 = 11; 𝛾𝛾 = 1 11 = 0.09 According to WHO, the mortality rate of Ebola is 0.7 [3] and the number of people susceptible is 4293265. American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 8 Therefore, 𝛽𝛽 (the rate of infection) = 0.7 4293265 = 1.63 Γ— 10βˆ’7 In order to use the SIR model to predict the evolution of the disease, it would be helpful if we could solve the system of differential equations. Unfortunately, we cannot completely solve these equations with an explicit formula solution. Therefore, we will use numerical approaches. We will use Euler method and fourth order Runge-Kutta method to extract the solution. 3.1 Euler Method For each day, we will calculate the values of 𝑆𝑆, 𝐼𝐼 and 𝑅𝑅 using 1. 𝑆𝑆𝑛𝑛+1 = 𝑆𝑆𝑛𝑛 βˆ’ 𝛽𝛽𝐼𝐼𝑛𝑛𝑆𝑆𝑛𝑛 2. 𝐼𝐼𝑛𝑛+1 = 𝐼𝐼𝑛𝑛 + (𝛽𝛽𝐼𝐼𝑛𝑛𝑆𝑆𝑛𝑛 βˆ’ 𝛾𝛾𝐼𝐼𝑛𝑛) 3. 𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + 𝛾𝛾𝐼𝐼𝑛𝑛 We take the initial values as 𝑆𝑆0 = 4293265; 𝐼𝐼0 = 391; 𝑅𝑅0 = 344; 𝛾𝛾 = 0.09 ; 𝛽𝛽 = 1.63 Γ— 10βˆ’7; We will do this explicitly for the transition from t = 0 to t = 1. Using equations 1, 2 and 3, the following values for S, I and R can be calculated. 𝑆𝑆1 = 𝑆𝑆0 βˆ’ 𝛽𝛽𝐼𝐼0𝑆𝑆0 Γ— βˆ†π‘‘π‘‘ = 4293265 βˆ’ ( 1.63 Γ— 10βˆ’7 Γ— 391 Γ— 4293265) Γ— 1 = 4292991.377 β‰ˆ 4292991 𝐼𝐼1 = 𝐼𝐼0 + (𝛽𝛽𝐼𝐼0𝑆𝑆0 βˆ’ 𝛾𝛾𝐼𝐼0) Γ— βˆ†π‘‘π‘‘ = 391 + (1.63 Γ— 10βˆ’7 Γ— 391 Γ— 4293265 βˆ’ 0.09 Γ— 391) Γ— 1 = 629.4326582 β‰ˆ 630 𝑅𝑅1 = 𝑅𝑅0 + 𝛾𝛾𝐼𝐼0 = 344 + 0.09 Γ— 391 = 379.19 β‰ˆ 379 Here we use MATLAB to evaluate S, I, R over a two month period. Table-3.1 shows the result. American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 9 Table 3.1: Euler Method Day S I R S+I+R 01 4293265 391 344 4294000 02 4292991 630 379 4294000 03 4292550 1014 436 4294000 04 4291842 1630 528 4294000 05 4290700 2625 675 4294000 06 4288865 4225 910 4294000 07 4285912 6798 1290 4294000 08 4281162 10936 1902 4294000 09 4273530 17584 2886 4294000 10 4261283 28248 4469 4294000 11 4241662 45327 7011 4294000 12 4210323 72586 11091 4294000 13 4160509 115868 17623 4294000 14 4081931 184017 28052 4294000 15 3959494 289892 44614 4294000 16 3772399 450898 70703 4294000 17 3495141 687575 111284 4294000 18 3103424 1017410 173166 4294000 19 2588759 1440508 264733 4294000 20 1980912 1918710 394378 4294000 21 1361382 2365556 567062 4294000 22 836453 2677585 779962 4294000 23 471386 2801669 1020945 4294000 24 256117 2764788 1273095 4294000 25 140695 2631379 1521926 4294000 26 80349 2454901 1758750 4294000 27 48198 2266111 1979691 4294000 28 30394 2079964 2183642 4294000 29 20090 1903072 2370838 4294000 30 13858 1738028 2542114 4294000 31 9932 1585531 2698537 4294000 32 7365 1445400 2841235 4294000 33 5630 1317049 2971321 4294000 34 4422 1199723 3089855 4294000 35 3557 1092613 3197830 4294000 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 10 36 2923 994911 3296166 4294000 37 2449 905843 3385708 4294000 38 2088 824679 3467233 4294000 39 1807 750738 3541455 4294000 40 1586 683393 3609021 4294000 41 1409 622065 3670526 4294000 42 1266 566222 3726512 4294000 43 1150 515378 3777472 4294000 44 1053 469091 3823856 4294000 45 972 426953 3866075 4294000 46 905 388595 3904500 4294000 47 847 353679 3939474 4294000 48 799 321896 3971305 4294000 49 757 292968 4000275 4294000 50 720 266637 4026643 4294000 51 689 242671 4050640 4294000 52 662 220858 4072480 4294000 53 638 201004 4092358 4294000 54 617 182935 4110448 4294000 55 599 166489 4126912 4294000 56 583 151521 4141896 4294000 57 568 137899 4155533 4294000 58 555 125501 4167944 4294000 59 544 114217 4179239 4294000 60 533 103948 4189519 4294000 3.2 Runge-Kutta (RK4) Method For each day, we will calculate the values of 𝑆𝑆, 𝐼𝐼 and 𝑅𝑅 using π‘˜π‘˜1𝑑𝑑 = βˆ’π›½π›½π‘†π‘†π‘›π‘›πΌπΌπ‘›π‘› π‘˜π‘˜1𝑑𝑑 = 𝛽𝛽𝑆𝑆𝑛𝑛𝐼𝐼𝑛𝑛 βˆ’ 𝛾𝛾𝐼𝐼𝑛𝑛 π‘˜π‘˜1𝑑𝑑 = 𝛾𝛾𝐼𝐼𝑛𝑛 π‘˜π‘˜2𝑑𝑑 = βˆ’π›½π›½ �𝑆𝑆𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 11 π‘˜π‘˜2𝑑𝑑 = 𝛽𝛽 �𝑆𝑆𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜2𝑑𝑑 = 𝛾𝛾 �𝐼𝐼𝑛𝑛 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜3𝑑𝑑 = βˆ’π›½π›½ �𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) π‘˜π‘˜3𝑑𝑑 = 𝛽𝛽 �𝑆𝑆𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜3𝑑𝑑 = 𝛾𝛾 �𝐼𝐼𝑛𝑛 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ π‘˜π‘˜4𝑑𝑑 = βˆ’π›½π›½(𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) π‘˜π‘˜4𝑑𝑑 = 𝛽𝛽(𝑆𝑆𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) βˆ’ (𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) π‘˜π‘˜4𝑑𝑑 = 𝛾𝛾(𝐼𝐼𝑛𝑛 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) 𝑆𝑆𝑛𝑛+1 = 𝑆𝑆𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) 𝐼𝐼𝑛𝑛+1 = 𝐼𝐼𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) 𝑅𝑅𝑛𝑛+1 = 𝑅𝑅𝑛𝑛 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) We take the initial values as 𝑆𝑆0 = 4293265 ; 𝐼𝐼0 = 391; 𝑅𝑅0 = 344 ; 𝛾𝛾 = 0.09; 𝛽𝛽 = 1.63 Γ— 10βˆ’7 We will do this explicitly for the transition from t = 0 to t = 1. Using those equations the following values for S, I and R can be calculated. π‘˜π‘˜1𝑑𝑑 = βˆ’π›½π›½π‘†π‘†0𝐼𝐼0 = βˆ’1.63 Γ— 10βˆ’7 Γ— 4293265 Γ— 391 = βˆ’273.622658245000 π‘˜π‘˜1𝑑𝑑 = 𝛽𝛽𝑆𝑆0𝐼𝐼0 βˆ’ 𝛾𝛾𝐼𝐼0 = 1.63 Γ— 10βˆ’7 Γ— 4293265 Γ— 391 βˆ’ 0.09 Γ— 391 = 238.432658245000 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 12 π‘˜π‘˜1𝑑𝑑 = 𝛾𝛾𝐼𝐼0 = 0.09 Γ— 391 = 35.1900000000000 π‘˜π‘˜2𝑑𝑑 = βˆ’π›½π›½ �𝑆𝑆0 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼0 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) = βˆ’1.63 Γ— 10βˆ’7 Γ— οΏ½4293265 βˆ’ 273.622658245000 2 οΏ½ Γ— οΏ½391 + 238.432658245000 2 οΏ½ = βˆ’357.039129094785 π‘˜π‘˜2𝑑𝑑 = 𝛽𝛽 �𝑆𝑆0 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ��𝐼𝐼0 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼0 + π‘˜π‘˜1π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ = βˆ’1.63 Γ— 10βˆ’7 Γ— οΏ½4293265 + βˆ’273.622658245000 2 οΏ½ Γ— οΏ½391 + 238.432658245000 2 οΏ½ βˆ’ οΏ½391 + 238.432658245000 2 οΏ½ = 311.119659473760 π‘˜π‘˜2𝑑𝑑 = 𝛾𝛾 �𝐼𝐼0 + π‘˜π‘˜1 πΌπΌβˆ†π‘‘π‘‘ 2 οΏ½ = 0.09 Γ— οΏ½391 + 238.432658245000 2 οΏ½ =45.9194696210250 π‘˜π‘˜3𝑑𝑑 = βˆ’π›½π›½ �𝑆𝑆0 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 οΏ½ (𝐼𝐼0 + π‘˜π‘˜2π‘‘π‘‘βˆ†π‘‘π‘‘ 2 ) = βˆ’1.63 Γ— 10βˆ’7 Γ— οΏ½4293265 + βˆ’357.039129094785 2 οΏ½ οΏ½391 + 311.119659473760 2 οΏ½ = βˆ’382.467864374178 π‘˜π‘˜3𝑑𝑑 = 𝛽𝛽 �𝑆𝑆0 + π‘˜π‘˜2 π‘†π‘†βˆ†π‘‘π‘‘ 2 οΏ½ �𝐼𝐼0 + π‘˜π‘˜2 πΌπΌβˆ†π‘‘π‘‘ 2 οΏ½ βˆ’ �𝐼𝐼0 + π‘˜π‘˜2 πΌπΌβˆ†π‘‘π‘‘ 2 οΏ½ = 1.63 Γ— 10βˆ’7 Γ— οΏ½4293265 + βˆ’357.039129094785 2 οΏ½ Γ— οΏ½391 + 311.119659473760 2 οΏ½ βˆ’ οΏ½391 + 311.119659473760 2 οΏ½ = 333.277479697859 π‘˜π‘˜3𝑑𝑑 = 𝛾𝛾 �𝐼𝐼0 + π‘˜π‘˜2 πΌπΌβˆ†π‘‘π‘‘ 2 οΏ½ = 0.09 Γ— οΏ½391 + 311.119659473760 2 οΏ½ = 49.1903846763192 π‘˜π‘˜4𝑑𝑑 = βˆ’π›½π›½(𝑆𝑆0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 13 = βˆ’1.63 Γ— 10βˆ’7 Γ— (4293265 βˆ’ 382.467864374178) Γ— (391 + 333.277479697859) = βˆ’506.805816985307 π‘˜π‘˜4𝑑𝑑 = 𝛽𝛽(𝑆𝑆0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘)(𝐼𝐼0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) βˆ’ (𝐼𝐼0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) = 1.63 Γ— 10βˆ’7 Γ— (4293265 βˆ’ 382.467864374178) Γ— (391 + 333.277479697859) βˆ’ (391 + 333.277479697859) = 441.62084381250 π‘˜π‘˜4𝑑𝑑 = 𝛾𝛾(𝐼𝐼0 + π‘˜π‘˜3π‘‘π‘‘βˆ†π‘‘π‘‘) = 0.09 Γ— (391 + 333.277479697859) = 65.1849731728073 ∴ 𝑆𝑆1 = 𝑆𝑆0 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) = 4293265 + 1 6 (βˆ’273.622658245000 βˆ’ 2 Γ— 357.039129094785 βˆ’ 2 Γ— 382.467864374178 βˆ’ 506.805816985307) = 4292888.42625631 β‰ˆ 4292888 ∴ 𝐼𝐼1 = 𝐼𝐼0 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) = 391 + 1 6 (238.432658245000 + 311.119659473760 + 333.277479697859 + 441.62084381250) = 719.141296733456 β‰ˆ 719 ∴ 𝑅𝑅1 = 𝑅𝑅0 + βˆ†π‘‘π‘‘ 6 (π‘˜π‘˜1𝑑𝑑+2π‘˜π‘˜2𝑑𝑑+2π‘˜π‘˜3𝑑𝑑+π‘˜π‘˜4𝑑𝑑) = 344 + 1 6 (35.1900000000000 + 45.9194696210250 + 49.1903846763192 + 65.1849731728073) = 392.432446961249 β‰ˆ 393 Here we use MATLAB to evaluate S, I, R over a two month period.Table-3.2 shows the result. American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 14 Table 3.2 RK4 Method Day S I R S+I+R 01 4293265 391 344 4294000 02 4292888 719 393 4294000 03 4292196 1322 482 4294000 04 4290923 2432 645 4294000 05 4288583 4471 946 4294000 06 4284286 8214 1500 4294000 07 4276407 15077 2516 4294000 08 4261994 27626 4380 4294000 09 4235753 50457 7790 4294000 10 4188376 91623 14001 4294000 11 4104130 164649 25221 4294000 12 3958283 290515 45202 4294000 13 3717118 496967 79915 4294000 14 3346977 809182 137841 4294000 15 2838838 1226389 228773 4294000 16 2237427 1696343 360230 4294000 17 1636735 2124419 532846 4294000 18 1127150 2428029 738821 4294000 19 747957 2580728 965315 4294000 20 489487 2605010 1199503 4294000 21 321556 2540850 1431594 4294000 22 214459 2424226 1655315 4294000 23 146134 2280695 1867171 4294000 24 102037 2126431 2065532 4294000 25 73073 1971024 2249903 4294000 26 53657 1819893 2420450 4294000 27 40360 1675938 2577702 4294000 28 31056 1540566 2722378 4294000 29 24413 1414309 2855278 4294000 30 19575 1297197 2977228 4294000 31 15987 1188973 3089040 4294000 32 13279 1089223 3191496 4294000 33 11204 997456 3285340 4294000 34 9589 913148 3371263 4294000 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 15 35 8316 835770 3449914 4294000 36 7300 764806 3521894 4294000 37 6479 699764 3587757 4294000 38 5809 640176 3648015 4294000 39 5257 585604 3703139 4294000 40 4798 535640 3753562 4294000 41 4414 489905 3799681 4294000 42 4089 448050 3841861 4294000 43 3813 409750 3880437 4294000 44 3577 374709 3915714 4294000 45 3374 342652 3947974 4294000 46 3199 313328 3977473 4294000 47 3046 286506 4004448 4294000 48 2913 261974 4029113 4294000 49 2796 239538 4051666 4294000 50 2694 219019 4072287 4294000 51 2603 200255 4091142 4294000 52 2523 183095 4108382 4294000 53 2452 167404 4124144 4294000 54 2389 153057 4138554 4294000 55 2333 139937 4151730 4294000 56 2283 127941 4163776 4294000 57 2238 116972 4174790 4294000 58 2197 106943 4184860 4294000 59 2161 97773 4194066 4294000 60 2128 89389 4202483 4294000 Figure 3.1-3.3 shows the comparison between Euler and RK-4 method in the cases of infectives, recovered and susceptives respectively. Figure 3.1 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 16 Figure 3.2 Figure 3.3 For predict the evolution, we only discuss about the infectives. Each strategies show that initially, the number of individuals infected will increase steeply, however, over an extended period of time, the numbers eventually decrease. This happens simultaneously because the range of individuals recovered will increase because as those infected decreases, they're being transferred into the recovered class. The main reason is owing to the raised awareness of the disease leasing to additional medical support being given in order to assist combat the transmission of the disease. Furthermore, an increased awareness ends up in additional folks being conscious of strategies of protection. The steep increase within the beginning of the primary fifteen days is possibly to flow from to the good uncertainty that lied with Ebola allowing a greater rate of transmission. The peak of every graph illustrates the utmost number of individuals ever to be infected and once this point; there's a transition whereby the numbers decrease. 3.3 Comparing the model to actual data American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 17 However, in order for the model to be valid and allow informing government policy, it obviously needs to correspond fairly close to reality. The Table-3.3.1 below compares the data collected from the SIR model (using Euler and RK4) for the number of people infected and the real life data of the number of people infected. Table 3.3.1 Time(Days) I -Actual I- Euler I -RK4 Time(Days) I -Actual I- Euler I -RK4 1 391 391 391 35 1871 1092613 835770 2 486 630 719 37 2046 905843 699764 6 554 4225 8214 41 2081 622065 489905 10 599 28248 91623 45 2407 426953 342652 11 670 45327 164649 47 2710 353679 286506 13 786 115868 496967 51 3022 242671 200255 17 834 687575 2124419 53 3280 201004 167404 19 972 1440508 2580728 55 3458 166489 139937 20 1082 1918710 2605010 60 3696 103948 89389 26 1378 2454901 1819893 Using the value of Table-3.3.1 we can plot a graph using MATLAB which compares the model data to the actual data for the number of people infected (Figure-3.4). Figure 3.4 As the graph demonstrates, the real data does not correspond very well to the data received from the model. Although the actual data may seem to follow a straight like graph, this is untrue as it is only depicted in this American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 18 manner due to the limitations on the axis of the graph. The difference between the real life data and the data from the model is so vast that the straight line looks like a graph of 𝑦𝑦 = 0. Therefore, we decided to plot it separately (Figure-3.5). Figure 3.5 3.4 Curve Fitting with Actual Data The model has significantly overestimated the number of individuals who can become infected with Ebola fever. This is often because of the many limitations that the model presents. One of the main limitations includes the incorrect beta and gamma values that were calculated. Once fixing the beta and gamma values, we were able to find another gamma value that resulted in similar values to the real data. Here is that the graph to indicate this, with the suitable gamma value of 0.66. The Table-3.4.1 below compares the data collected from the fitted model for the number of people infected and the real life data of the number of people infected. Table-3.4.1 Time(Days) I -Actual I- Euler I -RK4 Time(Days) I -Actual I- Euler I -RK4 1 391 391 391 35 1871 1555 1419 2 486 408 407 37 2046 1679 1525 6 554 482 475 41 2081 1952 1758 10 599 569 555 45 2407 2261 2020 11 670 593 576 47 2710 2430 2163 13 786 644 623 51 3022 2796 2472 17 834 759 726 53 3280 2994 2639 19 972 824 783 55 3458 3201 2814 20 1082 858 813 60 3696 3757 3287 26 1378 1092 1019 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 19 Using the value of Table-3.4.1 we can plot a graph using MATLAB which compares curve fitted data with real data for the number of people infected (Figure-3.6). Figure 3.6 3.5 Result Discussion Although, this doesn't match the graph exactly, it shows a much better correlation of the number of individuals infected. Between Euler and RK4 technique, Rk4 provides better correlation as because Rk4 methods divide its step size more than Euler methods and the local truncation error is less than Euler methods. Therefore, in order to improve the model, many changes should be done, together with altering the gamma value. The value that we eventually used to alter the model, led to being in 2 decimal places. This goes to illustrate the necessary precision required as very little deviance will cause massive changes. This can be because; the gamma is calculated through extreme simplification, leaving great possibilities for more room for errors. Furthermore, it is tough to differentiate between the numbers of people who have died and also the numbers of people who have survived with permanent immunity as they both fall under the same class of being β€˜recovered’. It is also a fact that rate of recovery is faster that the time scale of birth and death. The infection rate and recovery rate plays a vital role in the model and in this method we find more accurate value of infection rate and recovery rate to fit the model with the actual data collected from CDC. 4. Conclusion The results obtained from modeling data will lead to completely different views and interpretations. This is due to the unequal distribution of knowledge across the globe whereby in countries like Liberia, there is little access to the statistics which makes it troublesome to form constructive predictions regarding the outbreak. Through our research, we have gained further insight into the uses of mathematical modeling so as to work out the outbreak of diseases similarly as evaluating its flaws. Having chosen Ebola as the diseases of concentration, as it is incredibly relevant to this situation in continent, it has enabled a practical understanding of its rate of American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 20 transmission. 5. Recommendations The Numerical methods give a better prediction if the data collected from an authentic source. This result can help to estimate future predictions of the disease and consequently, it will facilitate to determine practical components like the quantity of beds required in the hospitable, number of vaccination and reallocation costs etc. References [1] "Ebola Virus Disease". World Health Organization. N.p., 2017. Web. 29 Mar. 2017. [2] "2014-2016 Ebola Outbreak In West Africa| Ebola Hemorrhagic Fever | CDC". Cdc.gov. N.p., 2017. Web. 29 Mar. 2017. [3] "WHO Finds 70 Percent Ebola Mortality Rate". Aljazeera.com. N.p., 2017. Web. 29 Mar. 2017. [4] "Previous Case Counts| Ebola Hemorrhagic Fever | CDC". Cdc.gov. N.p., 2017. Web. 29 Mar. 2017. [5] Dolgoarshinnykh, R. G., & Lalley, S. P. (2003). Epidemic Modelling: SIRS Models (Doctoral dissertation, University of Chicago, Department of Statistics). [6] "Modelling Infectious Diseases." IB Maths Resources From British International School Phuket". ibmathsresources.com. N.p., 2017. Web. 29 Mar. 2017. [7] "The Spread Of Infectious Diseases." The British Medical Journal 2.1281 (1885): 108. Web. [8] Leone, S. Appendix: Additional Results and Technical Notes for the EbolaResponse Modeling Tool Additional Results. Population, 4, 3. [9] "Previous Case Counts| Ebola Hemorrhagic Fever | CDC". Cdc.gov. N.p., 2017. Web. 29 Mar. 2017. Retrieve date: 03 August 2014 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 37, No 1, pp 1-21 21 [10] Clinaero, Inc. "Ebola Incubation Period". eMedTV: Health Information Brought To Life. N.p., 2017. Web. 29 Mar. 2017. [11] "Kermack-Mckendrick Model – From Wolfram Mathworld". Mathworld.wolfram.com. N.p., 2017. Web. 29 Mar. 2017. [12] "The SIR Model For Spread Of Disease - Euler's Method For Systems". Mathematical Association of America. N.p., 2017. Web. 29 Mar. 2017. [13] Rahman, Prof. Dr. Md. Fazlur. Mathematical Modelling In Biology. 7th ed. ISBN-984-8759-19-0, 2015. Print. [14] Tsai, Tony. β€œTony Tsai.” RK4 Method for Solving SIR Model.N.p., n.d. Web. 16 Apr. 2017