IHJPAS. 36 (3) 2023 289 This work is licensed under a Creative Commons Attribution 4.0 International License *Corresponding Author: mahasssa@yahoo.com Abstract The aim of our study is to solve a nonlinear epidemic model, which is the COVID-19 epidemic model in Iraq, through the application of initial value problems in the current study. The model has been presented as a system of ordinary differential equations that has parameters that change with time. Two numerical simulation methods are proposed to solve this model as suitable methods for solving systems whose coefficients change over time. These methods are the Mean Monte Carlo Runge-Kutta method (MMC_RK) and the Mean Latin Hypercube Runge-Kutta method (MLH_RK). The results of numerical simulation methods are compared with the results of the numerical Runge-Kutta 4th order method (RK4) from 2021 to 2025 using the absolute error, which proves that the MLH_RK method is the best and closest to the expected values. The results have been discussed after being tabulated and represented graphically. Epidemic behavior for the next two years until 2025 has been projected using the proposed methods. Keywords: Epidemic models, Coronavirus (COVID-19) model, Runge-Kutta method, Simulation process, Numerical simulation techniques. 1. Introduction Coronavirus disease (COVID-19) is a new pandemic disease that began to appear in China, specifically in the city of Wuhan, and then spread rapidly around the world. This epidemic represented a global health problem, so some solutions and suggestions are discussed and studied to limit its spread as well as control it. Although, there are many researchers that have talked about this epidemiological model, studying its effects and the speed of its spread, as well as the number of people infected due to this epidemic and the extent of its impact in the future [1], The spread of the coronavirus began in 2019-2020, and this new disease is considered as a type of (SARS-COV) and has the fastest spread compared to other viruses [2]. Specifically, it quickly and suddenly spread around the world on 29th March 2020 [3, 4]. The World Health Organization declared that doi.org/10.30526/36.3.2945 Article history: Received 20 July 2022, Accepted 24 October 2022, Published in July 2023. Ibn Al-Haitham Journal for Pure and Applied Sciences Journal homepage: jih.uobaghdad.edu.iq Comparison of Some Numerical Simulation Techniques for COVID-19 Model in Iraq Mahdi A. Sabea Department of Mathematics, College of Education for Pure Sciences Ibn al-Haytham,University of Baghdad, Baghdad, Iraq. mahdiabd281@gmail.com Maha A. Mohammed* Department of Mathematics, College of Education for Pure Sciences Ibn al-Haytham,University of Baghdad, Baghdad, Iraq. mahasssa@yahoo.com https://creativecommons.org/licenses/by/4.0/ mailto:mahasssa@yahoo.com http://en.uobaghdad.edu.iq/?page_id=15060 http://en.uobaghdad.edu.iq/?page_id=15060 mailto:mahdiabd281@gmail.com mailto:mahdiabd281@gmail.com http://en.uobaghdad.edu.iq/?page_id=15060 http://en.uobaghdad.edu.iq/?page_id=15060 mailto:mahasssa@yahoo.com mailto:mahasssa@yahoo.com IHJPAS. 36(3)2023 290 this disease is a pandemic after spreading to 199 countries, as it first appeared in China in the city of Wuhan, and in early 2020, thousands of people died around the world [5]. The infection of this virus occurred unevenly in terms of damage, according to immunity and adherence to health prevention methods [5]. The spread of this virus has continuously led to great damage in most of the countries that have been affected. This has a significant impact on both global health and economic movement [6]. A large number of people have been infected with this epidemic as a result of mixing, interaction, and constant movement, just as some countries suffer from a lack of adherence to health prevention methods and a large number of them die. Some countries have adhered to the system of spacing, and A curfew was imposed for a specified period of time to limit the spread of the disease [7]. The epidemic model is the one that deals with a disease that spreads suddenly and rapidly and occupies large areas of society. These epidemiological models can be formulated as a system of nonlinear differential equations, and this epidemic model is considered a stochastic-deterministic model [8, 9]. Many researchers have studied analyzing the behavior of the epidemic and observing it over time to see whether it grows or settles in the population exposed to it. The SIR epidemic model was considered by the Banach contraction method, the Daftardar-Jafari method, and the Temimi-Ansari method [10]. LTAM was discussed for the first time to apply to the epidemic nonlinear model; this method is a merge of the Laplace transform with Tamimi and the Ansari iterative method [11]. The study's findings for the COVID-19 SIR epidemic model were designed with stochasticity in mind in order to eradicate the epidemic in Iraq [12]. To analyze the behavior of smoking habits, the Adomian decomposition method, variational iteration method, finite difference method, and Runge-Kutta method were used [13]. In addition, the Adomian decomposition method and variation iteration method were used for alcohol consumption [14]. Some researchers use the numerical Runge-Kutta method to analyze the behavior of obesity [15], smoking habits [13], influenza [16], and COVID-19 [17]. On the other hand, there are those who are interested in studying epidemiological models that have random coefficients and using suitable numerical simulation methods to solve them, such as [18-21]. The numerical method in this study is the Runge-Kutta method for the 4th order. It is considered one of the most efficient numerical methods for solving ordinary differential equations of different orders. High accuracy using this method can be achieved with order four to decrease the errors [22] and [9]. Four stages are considered in this method of order four [23]. Two numerical simulation methods are used, the first being the Mean Monte Carlo Runge-Kutta method (MMCRK), which is a reliable numerical simulation technique that merges two methods of different natures together in this study: the Monte Carlo simulation process (MC) and the Runge-Kutta method (RK) that are used to solve ordinary differential equations. Mahdi A. S. and Maha A. M. in 2019 discussed a modified numerical simulation technique for solving nonlinear epidemic models [19]. The second numerical simulation method is Mean Latin Hypercube Runge-Kutta (MLHRK). MLHRK is a numerical simulation method that is a mixture of the simulation method Latin Hypercube (LH) and the numerical method Runge-Kutta (RK). It is one of the reliable methods to solve such systems. Shatha J. M. and Maha A. M. (2021) studied the Mean Latin Hypercube Runge Kutta method to solve the influenza model. Emad T. G. and Maha A. M. discussed applying a suitable approximate- simulation technique to the COVID-19 epidemic model with random parameters in 2022; see [20] and [21]. IHJPAS. 36(3)2023 291 2. Mathematical Model of COVID-19 The mentioned model has been successfully used to study the people vaccinated against the COVID-19 epidemic in Iraq [24]. The population consists of five types of individuals S, V, A, I and R represent susceptible, vaccinated, asymptomatic, symptomatic, and the recovery respectively. These Individuals are functions according to time in the system. The governing equations for the epidemic under study are the nonlinear first-order ordinary differential equations [25]. 𝑆′(𝑑) = 𝑀 βˆ’ πœπ‘† βˆ’ 𝛼(1+𝛽𝐴)𝑆 𝑁 βˆ’ πœ‡π‘† + 𝛾𝑅, 𝑉′(𝑑) = πœπ‘† βˆ’ πœŒπ›Ό(1+𝛽𝐴)𝑉 𝑁 βˆ’ πœ‡π‘‰ , 𝐴′(𝑑) = 𝛼(1+𝛽𝐴)𝑆 𝑁 + πœŒπ›Ό(1+𝛽𝐴)𝑉 𝑁 βˆ’ 𝛿𝐴 βˆ’ πœ‡π΄, (1) 𝐼′(𝑑) = πœƒπ›Ώπ΄ βˆ’ 𝜎𝐼 βˆ’ πœ‡πΌ 𝑅′(𝑑) = (1 βˆ’ πœƒ)𝛿𝐴 + 𝜎𝐼 βˆ’ 𝛾𝑅 βˆ’ πœ‡π‘…, where Tables 1 and 2 represent subpopulation 𝑆, 𝑉, 𝐴, 𝐼, 𝑅 as the variables of model (1), and parameters 𝑀, 𝜏, 𝛼, 𝛽, πœ‡, 𝛾, 𝜌, 𝛿, πœƒ and 𝜎. System (1) have the initial conditions for the system obtained from the World Health Organization website for Iraq data [4] as the folwing: 𝑆(0) = 500, 𝑉(0) = 30 , 𝐴(0) = 36, 𝐼(0) = 24 and 𝑅(0) = 10, with the predicted parameters that are given in Table 2. Table 1 Variables of COVID-19 model [25]. Variable Definition 𝑺(𝒕) People who are not infected but are vulnerable not to have immunity 𝑽(𝒕) People vaccinated against coronavirus 𝑨(𝒕) People infected with the virus without showing any symptoms 𝑰(𝒕) Infected people have clear symptoms of infection in their bodies 𝑹(𝒕) People who have recovered from the coronavirus or who died because of coronavirus infection IHJPAS. 36(3)2023 292 Table 2 Parameters of COVID-19 model [25]. Parameter Definition Value Source 𝜢 The rate of transmission people who are infected with this virus 0.8883 [26] 𝜷 Correction factor for the rate of transmission of people without infection 0.45 [26] 𝝁 The rate of natural death 0.00003349 day [26] 𝜸 The rate of immunity 0.005 Assume 𝟏 βˆ’ 𝝆 The Vaccine efficacy and potency 0.8 Assume 𝟏/𝜹 The average period without symptoms of infection 7 days [27] 𝜽 The proportion of people does not show the effects of the symptoms of the virus, but it develops into a state of infection 0.2 [27] 𝟏 βˆ’ 𝜽 The proportion of individuals who recover but are asymptomatic 0.8 [27] 𝑴 The Birth rate in the community 1500/day Assume 𝝉 The Vaccination rate against the virus 0.01 day Assume 𝟏/𝝈 The average rate of people recovering from infection with the virus 10 days Assume 3. Numerical Method for Solving Covid-19 Model RK4 is one of the numerical iterative methods with high accuracy. The nonlinear system (1) for the COVID-19 model can be solved by RK4 with the mentioned initial conditions.The general form of RK4 is: 𝑦𝑖+1 = 𝑦𝑖 + β„Ž 6 (𝑙1 + 2𝑙2 + 2𝑙3 + 𝑙4), (2) where 𝑆𝑖+1 = 𝑓(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), = 𝑆𝑖 + 1 6 (𝑙𝑆1 + 2𝑙𝑆2 + 2𝑙𝑆3 + 𝑙𝑆4) βˆ— β„Ž, (3) 𝑉𝑖+1 = 𝑓(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), = 𝑉𝑖 + 1 6 (𝑙𝑉1 + 2𝑙𝑉2 + 2𝑙𝑉3 + 𝑙𝑉4) βˆ— β„Ž, (4) 𝐴𝑖+1 = 𝑓(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), = 𝐴𝑖 + 1 6 (𝑙𝐴1 + 2𝑙𝐴2 + 2𝑙𝐴3 + 𝑙𝐴4) βˆ— β„Ž, (5) 𝐼𝑖+1 = 𝑓(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), = 𝐼𝑖 + 1 6 (𝑙𝐼1 + 2𝑙𝐼2 + 2𝑙𝐼3 + 𝑙𝐼4) βˆ— β„Ž, (6) 𝑅𝑖+1 = 𝑓(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), = 𝑅𝑖 + 1 6 (𝑙𝑅1 + 2𝑙𝑅2 + 2𝑙𝑅3 + 𝑙𝑅4) βˆ— β„Ž, for all 𝑖 = 1,2, … , π‘š. (7) Now, the first iteration 𝑙𝑆1, 𝑙𝑉1, 𝑙𝐴1, 𝑙𝐼1 and 𝑙𝑅1 can be found as follows: 𝑙𝑆1 = 𝑓1(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), (8) IHJPAS. 36(3)2023 293 𝑙𝑉1 = 𝑓1(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), (9) 𝑙𝐴1 = 𝑓1(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), (10) 𝑙𝐼1 = 𝑓1(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ), (11) 𝑙𝑅1 = 𝑓1(𝑑𝑖 , 𝑆𝑖 , 𝑉𝑖 , 𝐴𝑖 , 𝐼𝑖 , 𝑅𝑖 ). for all 𝑖 = 0,1, … , π‘š. (12) Also, the second iteration 𝑙𝑆2, 𝑙𝑉2, 𝑙𝐼2, 𝑙𝑄2 and 𝑙𝑅2, in the same way, is given as follows: 𝑙𝑆2 = 𝑓2 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆1, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉1, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴1, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼1, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅1), (13) 𝑙𝑉2 = 𝑓2 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆1, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉1, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴1, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼1, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅1), (14) 𝑙𝐴2 = 𝑓2 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆1, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉1, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴1, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼1, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅1), (15) 𝑙𝐼2 = 𝑓2 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆1, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉1, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴1, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼1, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅1), (16) 𝑙𝑅2 = 𝑓2 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆1, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉1, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴1, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼1, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅1). for all 𝑖 = 0,1, … , π‘š. (17) To get the third iteration, 𝑙𝑆3, 𝑙𝑉3, 𝑙𝐴3, 𝑙𝐼3 and 𝑙𝑅3, follow the same way: 𝑙𝑆3 = 𝑓3 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆2, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉2, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴2, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼2, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅2), (18) 𝑙𝑉3 = 𝑓3 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆2, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉2, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴2, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼2, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅2), (19) 𝑙𝐴3 = 𝑓3 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆2, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉2, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴2, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼2, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅2), (20) 𝑙𝐼3 = 𝑓3 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆2, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉2, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴2, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼2, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅2), (21) 𝑙𝑅3 = 𝑓3 (𝑑𝑖 1 2 β„Ž, 𝑆𝑖 + 1 2 β„Ž βˆ— 𝑙𝑆2, 𝑉𝑖 + 1 2 β„Ž βˆ— 𝑙𝑉2, 𝐴𝑖 + 1 2 β„Ž βˆ— 𝑙𝐴2, 𝐼𝑖 + 1 2 β„Ž βˆ— 𝑙𝐼2, 𝑅𝑖 + 1 2 β„Ž βˆ— 𝑙𝑅2). for all 𝑖 = 0,1, … , π‘š. (22) To obtain the fourth iteration, 𝑙𝑆4, 𝑙𝑉4, 𝑙𝐴4, 𝑙𝐼4 and 𝑙𝑅4, follow this: 𝑙𝑆4 = 𝑓4(𝑑𝑖 + β„Ž, 𝑆𝑖 + β„Ž βˆ— 𝑙𝑆3, 𝑉𝑖 + β„Ž βˆ— 𝑙𝑉3, 𝐴𝑖 + β„Ž βˆ— 𝑙𝐴3, 𝐼𝑖 + β„Ž βˆ— 𝑙𝐼3, 𝑅𝑖 + β„Ž βˆ— 𝑙𝑅3), (23) 𝑙𝑉4 = 𝑓4(𝑑𝑖 + β„Ž, 𝑆𝑖 + β„Ž βˆ— 𝑙𝑆3, 𝑉𝑖 + β„Ž βˆ— 𝑙𝑉3, 𝐴𝑖 + β„Ž βˆ— 𝑙𝐴3, 𝐼𝑖 + β„Ž βˆ— 𝑙𝐼3, 𝑅𝑖 + β„Ž βˆ— 𝑙𝑅3), (24) 𝑙𝐴4 = 𝑓4(𝑑𝑖 + β„Ž, 𝑆𝑖 + β„Ž βˆ— 𝑙𝑆3, 𝑉𝑖 + β„Ž βˆ— 𝑙𝑉3, 𝐴𝑖 + β„Ž βˆ— 𝑙𝐴3, 𝐼𝑖 + β„Ž βˆ— 𝑙𝐼3, 𝑅𝑖 + β„Ž βˆ— 𝑙𝑅3), (25) 𝑙𝐼4 = 𝑓4(𝑑𝑖 + β„Ž, 𝑆𝑖 + β„Ž βˆ— 𝑙𝑆3, 𝑉𝑖 + β„Ž βˆ— 𝑙𝑉3, 𝐴𝑖 + β„Ž βˆ— 𝑙𝐴3, 𝐼𝑖 + β„Ž βˆ— 𝑙𝐼3, 𝑅𝑖 + β„Ž βˆ— 𝑙𝑅3), (26) 𝑙𝑅4 = 𝑓4(𝑑𝑖 + β„Ž, 𝑆𝑖 + β„Ž βˆ— 𝑙𝑆3, 𝑉𝑖 + β„Ž βˆ— 𝑙𝑉3, 𝐴𝑖 + β„Ž βˆ— 𝑙𝐴3, 𝐼𝑖 + β„Ž βˆ— 𝑙𝐼3, 𝑅𝑖 + β„Ž βˆ— 𝑙𝑅3). (27) for all 𝑖 = 0,1, … , π‘š. We substitute Eqs. (8), (13), (18) and (23) into Eq. (3) to find the numerical solutions of 𝑆𝑖. Also, we put Eqs. (9), (14), (19) and (24) into Eq. (4) to achieve the numerical solutions of 𝑉𝑖. In the same manner, we substitute Eqs. (10), (15), (20) and (25) into Eq. (5) to obtain the numerical solutions of 𝐴𝑖 and substite Eqs. (11), (16), (21) and (26) into Eq. (6) to find the numerical solutions of 𝐼𝑖. Finally, we substitute Eqs. (12), (17), (22) and (27) into Eq. (7) to get the numerical solutions of 𝑅𝑖, for all 𝑖 = 0,1, … , π‘š. IHJPAS. 36(3)2023 294 4. Numerical Simulation Methods for Solving Covid-19 Model 4.1 Mean Monte Carlo Range-Kutta (MMC_RK) Mean Monte Carlo Runge-Kutta (MMC_RK) is an efficient numerical simulation method for solving such mathematical models under study. This method consists of mixing two different methods, one numerical, the Runge-Kutta method (RK), with the Monte Carlo simulation process (MC). MC estimates the model coefficients that are random variables, while RK is used to solve the model numerically. The average of the last RK iteration resulting from each MC repetition is considered the estimated approximate solution for the model under study. Because the natural epidemic models have randomness in their coefficients, the MMC_RK numerical simulation process is considered as more suitable method than the classical methods like RK4 that solve models depending on fixed parameters. The MMC_RK method is explained in Fig 1 and implemented by using MATLAB software. More details are shown in [19]. MMC_RK Procedure Figure 1: Steps of MMC_RK process, [19]. Step 5 Compute the mean of the n solutions in Step 4 (MMC_RK solutions) Step 4 Repeat Step 2 & 3 for n times Step 3 solve the system using RK4 method with m iterations and the last iteration becomes the final RK4 solution Step 2 select randomly one value of each parameter and substitute it into the system Step 1 Simulate the parameters by MC sample for n times at once IHJPAS. 36(3)2023 295 4.2 Mean Latin Hypercube Runge-Kutta (MLH_RK) Mean Latin Hypercube Runge-Kutta (MLH_RK) is a numerical simulation method that is a mixture of a numerical method called Latin Hypercube (LH) and a numerical method called Runge- Kutta (RK). This method is implemented via the Matlab program, and more details are shown in [20]. LH estimates the coefficients of the model, which are considered random variables; MC estimates the model coefficients; and RK is used to solve the model numerically. The average of the last RK iteration results with each LH repetition is considered the estimated approximate solution for the model under study. Everything we talked about before the method MMC_RK applies to the method MLH_RK. In addition to that, the MLH_RK method is more accurate and faster than the MMC_RK method because it simulates all model parameters at once. It also saves time and effort. Fig 2 explains the MLHRK in greater detail. . MLH_RK Procedure Figure 2: Steps of MLH_RK process, [20]. Step 5 Step 5: Calculate the mean of final solutions from step 4, as a solution system, called MLH_RK Step 4 Step 4: Repeat steps 1 and 2 for n- timesl Step 3 Step 3: Solve the system m-times iterations numerically by RK . The last iterative result is the final solution Step 2 Step 2: For each random parameter, one value is specified and replaced in the system Step 1 Step1: All model parameters have been simulated by LHS times for n times at ones IHJPAS. 36(3)2023 296 5. Results and Discussion Numerical and numerical simulation solutions of the nonlinear of COVID-19 epidemic model in Iraq have been discussed and analyzed in this section. In this study, the real step size has been used such that β„Ž=0.02 in a week, (52 weeks in a year, the data of the COVID-19 epidemic is taken from each week; therefore, in order to change the weeks to a months, the real step size is calculated as β„Ž = 1 52 β‰ˆ 0.02) and β„Ž=0.08 in a month, (12 month in a year, the data of the COVID-19 epidemic is taken from each month; therefore, in order to change the months to a year, the real step size is calculated as β„Ž = 1 12 β‰ˆ 0.08). Table 3. Numerical and numerical simulation results of COVID-19 model from 2021 to 2022 Model Variables Step Size, h (weekly &monthly) Number of Iteration, π’Ž (weekly & Monthly) RK4 (1 years) MMC_RK 100 rep. (1 years) MLH_RK 100 rep. (1 years) 𝑺(𝒕) 0.08 (monthly) 150 111.2080 112.0276 112.02304 0.02 (weekly) 600 111.6410 111.9771 111.9769 𝑽(𝒕) 0.08 (monthly) 150 46.2424 45.3405 47.1184 0.02 (weekly) 600 46.2938 45.3136 47.0938 𝑨(𝒕) 0.08 (monthly) 150 198.0505 198.4761 198.4017 0.02 (weekly) 600 197.6970 198.0655 197.9236 𝑰(𝒕) 0.08 (monthly) 150 38.2352 37.9597 38.6755 0.02 (weekly) 600 36.0473 35.8729 36.1012 𝑹(𝒕) 0.08 (monthly) 150 205.2406 205.7047 205.4156 0.02 (weekly) 600 205.6102 206.3006 205.3439 Table 4. Expected numerical simulation results of COVID-19 model from 2021 to 2025 Model Variables Step Size, h (weekly &monthly) Number of Iteration, π’Ž (weekly &monthly) RK4 (4 years) MMC_RK 100 rep. (4 years) MLH_RK 100 rep. (4years) 𝑺(𝒕) 0.08 (monthly) 600 45.3077 45.7295 45.1721 0.02 (weekly) 2400 45.4914 45.7972 45.3906 𝑽(𝒕) 0.08 (monthly) 600 34.9403 35.5509 35.3254 0.02 (weekly) 2400 34.8710 35.6280 35.1763 𝑨(𝒕) 0.08 (monthly) 600 13.5280 13.7776 13.2342 0.02 (weekly) 2400 13.4236 13.7370 13.4793 𝑰(𝒕) 0.08 (monthly) 600 9.3148 9.3368 9.3201 0.02 (weekly) 2400 9.2985 9.3130 9.2956 𝑹(𝒕) 0.08 (monthly) 600 500.1596 500.4997 499.8646 0.02 (weekly) 2400 500.4212 500.7103 500.4103 IHJPAS. 36(3)2023 297 Table 5. Prediction intervals (5th percentile, 95th percentile) for MMC_RK and MLH_RK solutions MMC_RK from 2021 to 2025 (𝒕 ≀ πŸ’πŸ–) Subpopulation (100 repetitions) (1000 repetitions) 𝑺(𝒕) (34.4889 , 56.9510) (33.7740 , 58.7899) 𝑽(𝒕) (22.5928 , 57.1467) (22.0591 , 55.9641) 𝑨(𝒕) (7.3331 , 21.2286) (7.8735 , 20.6545) 𝑰(𝒕) (5.0933 , 13.9210) (5.5697 , 14.6953) 𝑹(𝒕) (447.8911 , 628.1435) (449.8386 , 628.6205) MLH_RK from 2021 to 2025 (𝒕 ≀ πŸ’πŸ–) Subpopulation (100 repetitions) (1000 repetitions) 𝑺(𝒕) (35.4277 , 59.8350) (33.5215 , 57.8323) 𝑽(𝒕) (26.3581 , 56.7235) (24.0661 , 56.6104) 𝑨(𝒕) (6.4966 , 17.3059) (6.5912 , 17.6509) 𝑰(𝒕) (4.2058 , 12.3463) (4.7953 , 12.5404) 𝑹(𝒕) (414.8612 , 578.4889) (410.9803 , 570.2245) Table 6. Absolute error for MMC_RK and MLH_RK with RK4 from 2021 to 2022. Model Variables Step Size, h (weekly & monthly) Iterations (months) MMC_RK MLH_RK 𝑺(𝒕) 0.08 (monthly) 12 0.8196 0.81504 0.02 (weakly) 0.3361 0.3359 𝑽(𝒕) 0.08 (monthly) 12 0.9019 0.8760 0.02 (weekly) 0.9802 0.8000 𝑨(𝒕) 0.08 (monthly) 12 0.4256 0.3512 0.02 (weekly) 0.3685 0.2266 𝑰(𝒕) 0.08 (monthly) 12 0.2755 0.04403 0.02 (weekly) 0.1744 0.0539 𝑹(𝒕) 0.08 (monthly) 12 0.4641 0.1750 0.02 (weekly) 0.6904 0.2663 IHJPAS. 36(3)2023 298 Table 7. Absolute error πΈβ„Ž = |π‘‡β„Ž βˆ’ π‘‡β„Ž 2 | for MMC_RK and MLH_RK from 2021 to 2025 Model Variables Step Size, h (weekly &monthly) Number of Iterations 𝑬𝒉_𝑴𝑴π‘ͺ_𝑹𝑲 100 rep. (4 years) 𝑬𝒉_𝑴𝑳𝑯_𝑹𝑲 100 rep. (4 years) 𝑺(𝒕) 0.08 monthly 600 0.045088 0.107386 0.08/2 1200 0.022554 0.111144 0.08/4 2400 0.019704 0.103381 𝑽(𝒕) 0.08 monthly 600 0.088598 0.109317 0.08/2 1200 0.104216 0.039814 0.08/4 2400 0.073518 0.029403 𝑨(𝒕) 0.08 monthly 600 0.027061 0.157000 0.08/2 1200 0.013543 0.048095 0.08/4 2400 0.010381 0.010124 𝑰(𝒕) 0.08 monthly 600 0.011903 0.010215 0.08/2 1200 0.011887 0.014358 0.08/4 2400 0.011402 0.009302 𝑹(𝒕) 0.08 monthly 600 0.332002 0.645725 0.08/2 1200 0.121379 0.099979 0.08/4 2400 0.109360 0.083376 Table 3 shows the numerical simulation solutions of MMC_RK and MLH_RK for one year (2021–2022) with 100 repetitions in the interval [0, 1] and real step sizes of 0.02, 0.08, and 0.12 weekly and monthly, respectively. Table 4 also includes the results of the numerical simulations MMC_RK and MLH_RK with 1000 for the societal groups in the monthly interval [0, 48] for the future time until 2025. The expected approximate MMC_RK and MLH_RK results have been calculated in the prediction intervals that contain the minimum bound (5th percentile) and maximum bound (95th percentile) for MMC_RK and MLH_RK results in the future until 2025. MMC_RK and MLH_RK results must be inside the predicted intervals; see Table 5. . To compare the two numerical simulation methods, MMC_RK and MLH_RK, with the numerical RK4 method, the absolute error for,,, and is calculated as a criterion for comparison from 2021 to 2022. In Table 6, it is clear that the MLH_RK method is more accurate and reliable than the MMC_RK method since the absolute error of the MLH_RK results is less than the MMC_RK results for all subpopulations of the studied model. Table 7 shows the numerical stability of the numerical methods used in the current study. Table 7 also explains how the numerical simulation method (MLH_RK) becomes more convergent as the step size is reduced, resulting in a smaller error between the proposed method and the actual method in step size. To prove the convergence of the methods used, see [28]. IHJPAS. 36(3)2023 299 Figure 3. Comparison of numerical and numerical simulation solutions by using RK4, MMC_RK and MLH_RK of (a) 𝑆(𝑑), (b) 𝑉(𝑑), (c) 𝐴(𝑑), (d) 𝐼(𝑑) and (e) 𝑅(𝑑) from 2021 to 2025 when β„Ž = 0.08 monthly in Iraq. IHJPAS. 36(3)2023 300 Figure 4. Comparison numerical simulation MMC_RK and MLH_RK results with numerical RK4 results for (a) 𝑆(𝑑), (b) 𝑉(𝑑), (c) 𝐴(𝑑), (d) 𝐼(𝑑) and (e) 𝑅(𝑑) from 2021 to 2025 when β„Ž = 0.04 monthly in Iraq. IHJPAS. 36(3)2023 301 Figure 5. The numerical simulation results MMC_RK and MLH_RK of (a) 𝑆(𝑑), (b) 𝑉(𝑑), (c) 𝐴(𝑑), (d) 𝐼(𝑑) and (e) 𝑅(𝑑) from 20201to 2025 when β„Ž = 0.02 weekly in iraq. IHJPAS. 36(3)2023 302 Figure 6. Curve fitting for Infected people 𝐼(𝑑) with observing data from 2021 to 2022 when β„Ž = 0.08 monthly in Iraq. Fig 3 shows the curves of the methods used in our study for 48 months, the interval 2021– 2025, in which all groups of society are shown according to the impact of the virus on them. Fig 3(a), which represents a group of people who are not infected with an epidemic but are susceptible to infection for this class for all RK4, MMC_RK, and MLH_RK methods, is used with step size weekly and monthly through 4 years. There is a sudden drop in the curve of this group for all the methods used in the study RK4, MMC_RK, and MLH_RK after the 33rd month, then it rises again after the 35th month to stabilize during the last months of our study. Because of continuos mixing between the people and the lack of commitment, people in the society can not continous profect their health properly. Therefore, we can observe the sudden descent in the curve as a result of a large number of infections, particularly between the 25th and 35th month of the study period, to settle down after the end of 2025. Fig 3(b) shows a curved group of vaccinated people against COVID-19 in Iraq, where we show that there is a simple rise in the curve for all RK4, MMC_RK, and MLH_RK methods used in the study from the beginning of the vaccine until the 25th month. The large number of people wanting to get the vaccine as a result of rising health and cultural awareness has resulted in a dramatic increase, particularly between the 27th and 35th month, after which the curve stabilized until the end of the study period. Fig 3(c) depicts the subpopulation as asymptomatic with the epidemic but not showing symptoms of infection.There is a simple rise in the first study months, bringing the curve of this group to its highest level between the 15th and 30th month. Because of following health prevention methods, the instructions of the World Health Organization, and the demand for vaccination against IHJPAS. 36(3)2023 303 the pandemic, the curve begins to decline gradually until the 40th month and settles at the beginning of the 40th month until the end of the year 2025. Fig 3(d) represents the number of people infected with the epidemic.The curve rose dramatically from the 15th month to the highest level in the 25th month, and then fell dramatically from the 33rd month to the 40th month to settle as a result of increased health awareness and following health guidelines, and more people are convinced to get the vaccine. Fig 3(e) represents the group of people who have been cured of the disease, as they have been removed from the list of positive cases. There is a gradual rise of the curve of this group until the 20th month, then it comes down in the 25th month, then it rises dramatically to reach its highest level in the 40th month due to following World Health Organization guidelines, so more people became fully vaccinated, then the curve is still at the same high level until the end of the study period. Fig 4 shows the convergence of the methods that are used, where the result that has the smallest step size has the greatest convergence. The convergence between the numerical simulation methods and the numerical method Runge-Kutta (RK4), which represented the exact solution, in Figure 4, is greater than the convergence in Fig 3 because the step size in Fig 4 is smaller. Also, we notice that the MLH_RK method is more convergent than MMC_RK to RK4. Fig 5 explains the similarities in the behavior of the solution for the proposed numerical simulations using the MMC_RK and MLH_RK methods despite the change in the step size (weekly with 100 repetitions) for each subpopulation , ,, , and for COVID-19 model for four years from 2021 to 2025. Finally, Fig 6 shows the curve fitting of real observed data from the World Health Organization, utilizing the concept of linear regression for the class of people who suffer from infection during the 12 months from 2021 to 2022 using the Magic Plot program. 6. Conclusion The epidemiological COVID-19 model in our current study has been solved for 48 months, from 2021 to 2022. Many methods are employed, including the numerical RK4 method and the other two numerical simulation methods MMC_RK and MLH_RK, to solve this model. In this study, we found convergence in the results for all methods that we used, but the MLH_RK method is the closest to the numerical method (RK4); the curve of MLH_RK results is more convergent with the curve of RK4 results than the curve of MMC_RK results. Therefore, it is considered the most efficient approximation method for solving this model. Also, we found that the epidemic model gives an impression of the impact of this virus on society. The results showed that the category of people not infected with the epidemic began to decrease during the study interval, while the category associated with vaccinated people increased. There is a gradual increase in the category of infected people who do not show symptoms because education does not adhere to health prevention methods, as well as social distancing.Then it decreases in the last months of the study period. When the category of infected people is clear and their symptoms are clear, there is a gradual increase in all methods (RK4, MMC_RK, and MLH_RK), then it decreases in the final months of the study interval due to a lack of people taking the virus vaccine. IHJPAS. 36(3)2023 304 Finally, in the category of people who have been cured of the disease, there is a clear increase for this class for all methods used for the study. Epidemic behavior for the next four years until 2025 has been projected using the proposed methods. The results of this study indicate a decrease in the epidemic rate in the future.. References 1. Coronavirus world meter website. [cited 13 July 2022; Available from: https://www.worldometers.info/coronavirus. 2. Rong Xinmiao; Yang Liu; Chu Huidi ; Fan Meng, Effect of delay in diagnosis on transmission of COVID-19. Math Biosci Eng, 2020. 17(3), 2725-2740. DOI: 10.3934/mbe.2020149. 3. National Health Commission of the People’s Republic of China, Situation report of the pneumonia cases caused by the novel coronavirus 2020. [cited 20 July 2022; Available from: http://www.nhc.gov.cn/xcs/yqfkdt/202002/3db09278e3034f289841300ed09bd0e1.shtml. 4. World Health Organization, Clinical management of severe acute respiratory infection when infection is suspected, 2020,. 5. World Health organization, coronavirus disease (COVID-19). [cited 13 July 2022; Available from: https://www.who.int/emergencies/overview. 6. Med., J.T. Pneumonia of unknown etiology in Wuhan, China: potential for international spread via commercial air travel. 2020; Available from: https://www.bing.com/search?q=Pneumonia+of+unknown+etiology+in+Wuhan%2C+China %3A+potential+for+international+spread+via+commercial+air+travel&cvid=ea862b7fe7bb4 9d490d0cc092b284831&aqs=edge..69i57.1194j0j1&pglt=299&FORM=ANNTA1&PC=U53 1. 7. Lin, Q., S. Zhao.; D. Gao.; Y. Lou.; S. Yang; SS. Musa, A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action. International journal of infectious diseases, 2020. 93: p. 211-216. DOI:10.1016/j.ijid.2020.02.058 8. Diekmann, O. ; J.A.P. Heesterbeek, Mathematical epidemiology of infectious diseases: model building, analysis and interpretation. Vol. 5. 2000: John Wiley & Sons. DOI:10.1016/j.ijid.2020.02.058 9. Brauer, F., C. Castillo-Chavez, ;C. Castillo-Chavez, Mathematical models in population biology and epidemiology. Vol. 2. 2012: Springer. 10. Abed, S.M.; M.A. AL-Jawary, Efficient iterative methods for solving the sir epidemic model. Iraqi Journal of Science, 2021: p. 613-622. DOI: 10.24996/ijs.2021.62.2.27. 11. Huisen, R.W.; S.H. Abd Almjeed; A.S. Mohammed, A Reliable Iterative Transform Method for Solving an Epidemic Model. Iraqi Journal of Science, 2021: p. 4839-4846. 12. Kareem, A.M.; S.N. Al-Azzawi, A stochastic differential equations model for the spread of coronavirus COVID-19): the case of Iraq. Iraqi Journal of Science, 2021: p. 1025-1035. DOI: 10.24996/ijs.2021.62.3.31. 13. Sabaa, M.A.; Mohammed, M.A., Approximate Solutions of Nonlinear Smoking Habit Model. Iraqi Journal of Science, 2020: p. 435-443. DOI: 10.24996/ijs.2020.61.2.23. 14. Sabaa, M.A.; Mohammed, M.A.; Abd Almjeed, S.H., Approximate solutions for alcohol consumption model in spain. Ibn AL-Haitham Journal For Pure and Applied Science, 2019. 32(3): p. 153-164. DOI: 10.30526/32.3.2292. 15. Mohammed, M. A.; Noor, N. F. M.; Siri, Z.; Ibrahim, A. I. N., Numerical solution for weight reduction model due to health campaigns in Spain. in AIP Conference Proceedings. 2015. AIP Publishing LLC. DOI:10.1063/1.4932414. 16. Mohammed, S.J.; M.A. Mohammed. Runge-kutta Numerical Method for Solving Nonlinear Influenza Model. in Journal of Physics: Conference Series. 2021. IOP Publishing. DOI:10.1088/1742-6596/1879/3/032040. http://www.worldometers.info/coronavirus http://www.nhc.gov.cn/xcs/yqfkdt/202002/3db09278e3034f289841300ed09bd0e1.shtml http://www.who.int/emergencies/overview http://www.bing.com/search?q=Pneumonia+of+unknown+etiology+in+Wuhan%2C+China%3A+potential+for+international+spread+via+commercial+air+travel&cvid=ea862b7fe7bb49d490d0cc092b284831&aqs=edge..69i57.1194j0j1&pglt=299&FORM=ANNTA1&PC=U531 http://www.bing.com/search?q=Pneumonia+of+unknown+etiology+in+Wuhan%2C+China%3A+potential+for+international+spread+via+commercial+air+travel&cvid=ea862b7fe7bb49d490d0cc092b284831&aqs=edge..69i57.1194j0j1&pglt=299&FORM=ANNTA1&PC=U531 http://www.bing.com/search?q=Pneumonia+of+unknown+etiology+in+Wuhan%2C+China%3A+potential+for+international+spread+via+commercial+air+travel&cvid=ea862b7fe7bb49d490d0cc092b284831&aqs=edge..69i57.1194j0j1&pglt=299&FORM=ANNTA1&PC=U531 http://www.bing.com/search?q=Pneumonia+of+unknown+etiology+in+Wuhan%2C+China%3A+potential+for+international+spread+via+commercial+air+travel&cvid=ea862b7fe7bb49d490d0cc092b284831&aqs=edge..69i57.1194j0j1&pglt=299&FORM=ANNTA1&PC=U531 https://doi.org/10.1016/j.ijid.2020.02.058 https://doi.org/10.24996/ijs.2021.62.2.27 https://doi.org/10.24996/ijs.2021.62.3.31 IHJPAS. 36(3)2023 305 17. Ghadeer, E.T.; M.A. Mohammed, Solving Nonlinear COVID-19 Mathematical Model Using a Reliable Numerical Method. Ibn AL-Haitham Journal For Pure and Applied Sciences, 2022. 35(2): p. 97-107. DOI: 10.30526/35.2.2818. 18. Mohammed, M. A., Ibrahim, A. I. N., Siri, Z., Noor, N. F. M., Mean Monte Carlo finite difference method for random sampling of a nonlinear epidemic system. Sociological Methods & Research, 2019. 48(1): p. 34-61. DOI: 10.1177/0049124116672683. 19. Mahdi, A.S., Modified Numerical Simulation Technique for Solving Nonlinear Epidemic Models. 2019, University of Baghdad: Baghdad, Iraq. 20. Mohammed, S.J.; M.A. Mohammed, Mean Latin Hypercube Runge-Kutta Method to Solve the Influenza Model. Iraqi Journal of Science, 2022. 63(3): p. 1158-1177. 21. Ghadeer, E.T.; M.A. Mohammed, Applying a suitable approximate-simulation technique of an epidemic model with random parameters. International Journal of Nonlinear Analysis and Applications, 2022. 22. Tam, C.K., Computational aeroacoustics-Issues and methods. AIAA journal, 1995. 33(10): p. 1788-1796. 23. Hu, F., M.Y. Hussaini; J. Manthey, Low-dissipation and low-dispersion Runge–Kutta schemes for computational acoustics. Journal of computational physics, 1996. 124(1): p. 177-191. DOI:10.1006/jcph.1996.0052. 24. Coronavirus world meter website, Iraq Population. [cited 13 July 2022; Available from: https://www.worldometers.info/world-population/iraq-population/. 25. Yang, B., Z. Yu; Y. Cai, The impact of vaccination on the spread of COVID-19: Studying by a mathematical model. Physica A: Statistical Mechanics and its Applications, 2022. 590: p. 126717. DOI:10.1016/j.physa.2021.126717. 26. Nadim, S.S.; I. Ghosh; J. Chattopadhyay, Short-term predictions and prevention strategies for COVID-19: a model-based study. Applied mathematics and computation, 2021. 404: p. 126251. DOI:10.1016/j.amc.2021.126251. 27. Serhani, M.; H. Labbardi, Mathematical modeling of COVID-19 spreading with asymptomatic infected and interacting peoples. Journal of Applied Mathematics and Computing, 2021. 66(1): p. 1-20. DOI:10.1007/s12190-020-01421-9. 28. Friedman, A.; McLeod, B., Blow-up of positive solutions of semi-linear heat equations, Indiana University Mathematics Journal, 1985. 34: p. 425-447. http://www.worldometers.info/world-population/iraq-population/