An-Najah National University Faculty of Graduate Studies COMPARISON OF NUMERICAL METHODS FOR SOLVING SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS By Abdul Fattah Boshnaq Supervisor Prof. Naji Qatanani This Thesis is Submitted in Partial Fulfillment of the Requirements for the Degree of Master's in Computerized Mathematics, Faculty of Graduate Studies, An-Najah National . University, Nablus, Palestine 2024 II COMPARISON OF NUMERICAL METHODS FOR SOLVING SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS By Abdul Fattah Boshnaq This thesis was discussed on 13/31/2024, and it was approved. Prof. Naji Qatanani Supervisor Signature Dr. Saed Mallak External Examiner Signature Dr. Hadi Hamad Internal Examiner Signature III Dedication إلى من كان لي سندًا وعونًا طوال عمري، إلى الرجل األبرز في حياتي أبي العزيـز إلى القلب المعطاء والصدر الحاني أمي الحبيبة إلي رفيقة دربي وشريكة حياتي زوجتي الغالية إلى من شد الله بهم عضدي فكانوا خير معين إخــواني وأخواتي إلى كل من ساعدني ولو بحرف في حياتي الدراسية... إلى هؤالء جميعًا: أهديكم هذا العمل IV Acknowledgment أتوجه بالشكر إلى األساتذة الكرام الذين لم يبخلوا علينا بعلمهم، ولم يألوا جهدًا في سبيل المعرفة والعلم. ذللت الصعوبات من أجل إنجاز هذه الرسالة بجودة وكفاءة.والشكر موصول إلى إدارة الجامعة التي الممتحن و د. ناجي قطناني لجنة المناقشة الفاضلة المتمثلة في مشرفي أ. إلىالشكر والتقدير أوجهكما الداخلي د. هادي حمد والممتحن الخارجي أ. د. سائد مالك لكم مني جزيل الشكر والعرفان V Decaration I, the undersigned, declare that I have submitted the thesis entitled: COMPARISON OF NUMERICAL METHODS FOR SOLVING SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS I further declare that the work presented in this thesis, unless otherwise referenced, is the researcher's own work and has not been submitted elsewhere for any other degree or qualification. Abdul Fattah Boshnaq Student's Name: Abdul Boshnaq Signature: 31/10/2024 Date: VI Table of Contents Dedication ....................................................................................................................... III Acknowledgment ............................................................................................................ IV Decaration ........................................................................................................................ V Table of contents............................................................................................................. VI List of Tables ............................................................................................................... VIII List of Figures ................................................................................................................. IX List of Appendices ........................................................................................................... X Abstract .......................................................................................................................... XII Introduction....................................................................................................................... 1 Chapter One: Mathematical Preliminaries ........................................................................ 3 1.1 Systems of Linear Algebraic Equations .................................................................... 3 1.2 Systems of Linear Differential Equations ................................................................. 3 1.3 Initial Value Problem (IVP) ...................................................................................... 4 1.4 Existence and Uniqueness ......................................................................................... 5 1.5 Linear Dependence and Linear Independence........................................................... 8 1.6 Homogeneous Vector Differential Equations............................................................ 9 1.7 Eigenvalues and Eigenvectors ................................................................................. 10 1.8 Kind of Eigenvalues ................................................................................................ 10 1.8.1 Distinct real eigenvalues ....................................................................................... 10 1.8.2 Repeated eigenvalues ............................................................................................ 10 1.8.3 Complex Eigenvalues ............................................................................................ 11 1.9 Nonhomogeneous Vector Differential Equations ................................................. 11 1.10 The Method of Undetermined Coefficients .......................................................... 13 1.11 Variation of Parameters ......................................................................................... 13 1.12 Higher Order Initial Value Problems .................................................................... 14 Chapter Two: Numerical Methods for Solving System of Differential Equation .......... 16 2.1 The Adomian Decomposition Method .................................................................... 16 2.2 Runge-KuttaMethod ................................................................................................ 20 2.2.1 Euler’s method ...................................................................................................... 20 2.2.2 Runge-Kutta Methods )RKM) .............................................................................. 21 2.2.3 System of first order differential equations ........................................................... 24 2.3 Adams-Bashforth Method (ABM)........................................................................... 25 VII 2.4 Adams- Moulton Method (AMM) ........................................................................... 27 2.2 Predictor-Corrector Method (PCM) ........................................................................ 29 2.6 Error Analysis .......................................................................................................... 30 2.7 Stability and Convergence ....................................................................................... 31 Chapter Three: Numerical Examples and Results .......................................................... 32 3.1 The numerical approximation of system (3.1) using the fourth order Rung Kutta method ............................................................................................................................ 33 3.2 The numerical approximation of system (3.1) using the Adams-Bashforth Method ..................................................................................................................................................... 43 3.3 The numerical approximation of system (3.1) using the Adams-Moulton method ........................................................................................................................................ 45 3.4 The numerical approximation of system (3.1) using the Predictor-Corrector method ........................................................................................................................................ 47 3.5 The numerical approximation of system (3.3) using the fourth order Rung Kutta methods ........................................................................................................................... 50 3.6 The numerical approximation of system (3.3) using the Fourth Adams-Bashforth Method ............................................................................................................................ 53 3.7 The numerical approximation of system (3.3) using the Adams-Moulton method ........................................................................................................................................ 55 3.8 The numerical approximation of system (3.3) using the Predictor-Corrector method ........................................................................................................................................ 57 3.9 The numerical approximation of system (3.5) using the fourth order Rung Kutta methods ........................................................................................................................... 60 3.10 The numerical approximation of system (3.5) using the Fourth Adams-Bashforth Method ............................................................................................................................ 61 3.11 The numerical approximation of system (3.5) using the Adams-Moulton method ..................................................................................................................................................... 62 3.12 The numerical approximation of system (3.5) using the Predictor-Corrector method ..................................................................................................................................................... 63 Chapter Four: Comparison of Numerical Methods ........................................................ 65 4.1 Conclusions ............................................................................................................. 66 References....................................................................................................................... 67 Appendices ..................................................................................................................... 70 ب................................................................................................................................ الملخص VIII List of Tables Table 3.1(a): The 𝑘’s values for 𝑥(𝑡) ............................................................................. 40 Table 3.1 (b): The 𝑘’s values for 𝑦(𝑡) ............................................................................ 40 Table 3.2(a):The exact and numerical solution of 𝑥(𝑡) by applying RKM .................... 40 Table 3.2(b): The exact and numerical solution of 𝑦𝑡 by applying RKM ...................... 42 Table 3.3(a): The exact and numerical solution of applying Fourth Adams-Bashforth method at 𝑡 = 0.8 .......................................................................................... 45 Table 3.3(b): The exact and numerical solution of applying Fourth Adams-Moulton method at t = 0.8 ........................................................................................... 47 Table3.3(c): The exact and numerical solution of applying Predictor-Corrector method at 𝑡 = 0.8 ....................................................................................................... 49 Table 3.4(a): The 𝑘’s value for 𝑥(𝑡) ............................................................................... 50 Table 3.4(b): The 𝑘’s value for 𝑦(𝑡) .............................................................................. 50 Table 3.4(c): The 𝑘’s value for 𝑧(𝑡) ............................................................................... 50 IX List of Figures Figure 1.1: is convex set S and non-convex set D ........................................................... 6 Figure 3.1(a): The exact and numerical solutions of 𝑥𝑡 ................................................. 41 Figure 3.1(b): the absolute error resulting of applying RKM of 𝑥𝑡 ............................... 41 Figure 3.2(a): The exact and numerical solutions of 𝑦𝑡 ................................................. 42 Figure 3.2(b): absolute error resulting of applying RKM of 𝑦𝑡 ..................................... 43 Figure 3.3(a): The exact and numerical solutions of 𝑥𝑡 ................................................. 51 Figure 3.3(b): the absolute error resulting of applying RKM. ........................................ 51 Figure 3.4(a): The exact and numerical solutions 𝑦𝑡...................................................... 52 Figure 3.4(b): the absolute error resulting of applying RKM of 𝑦𝑡 ............................... 53 X List of Appendices Appendix A: Matlab code for example 3.1 ..................................................................... 70 Appendix B: Matlab code for example 3.2 ..................................................................... 82 Appendix C: Matlab code for example 3.3 ..................................................................... 97 Appendix D: Tables ...................................................................................................... 109 Table 3.5(a): The exact and numerical solution of 𝑥𝑡 by applying RKM .................... 109 Table3.5(b): The exact and numerical solution of 𝑦𝑡 by applying RKM. .................... 109 Table 3.5(c): The exact and numerical solution of 𝑧𝑡 by applying RKM .................... 109 Table 3.6(a): The exact and numerical solution of applying Fourth Adams-Bashforth method at 𝑡 = 0.8 ........................................................................................ 109 Table 3.6 (b): The exact and numerical solution of applying Fourth Adams-Moulton method at 𝑡 = 0.8 ........................................................................................ 109 Table 3.6(c): The exact and numerical solution of applying Predictor-Corrector method at 𝑡 = 0.8 ..................................................................................................... 109 Table 3.7(a): The 𝑘’s value for 𝑥(𝑡) ............................................................................. 110 Table 3.7(b): The 𝑘’s value for 𝑦(𝑡) ............................................................................ 110 Table 3.8(a): The exact and numerical solution of 𝑥𝑡 by applying RKM .................... 110 Table 3.8(b): The exact and numerical solution of 𝑦𝑡 by applying RKM. ................... 110 Table 3.9(a): The exact and numerical solution of applying Fourth Adams-Bashforth methodat 𝑡 = 0.8 ......................................................................................... 110 Table 3.9(b): The exact and numerical solution of applying Fourth Adams-Moulton method at 𝑡 = 0.8 ........................................................................................ 110 Table 3.9(c): The exact and numerical solution of applying Predictor-Corrector method at 𝑡 = 0.8 ..................................................................................................... 110 Table 4.1: The Exact and Numerical solutions at 𝑡 = 0.8 of system (3.1) ................... 111 Table 4.2: The Exact and numerical solutions of at 𝑡 = 0.8 of system (3.3) ............... 111 Table 4.3: The Exact and numerical solutions at 𝑡 = 0.8 of system (3.5).................... 111 Appendix E: Figures ..................................................................................................... 112 Figure 3.5(a): The exact and numerical solutions of 𝑧𝑡................................................ 112 Figure 3.5(b): the absolute error resulting of applying RKM of 𝑧𝑡 .............................. 112 Figure 3.6(a): the exact and numerical solution 𝑥𝑡...................................................... 113 Figure 3.6 (b): the absolute error resulting of applying RKM of 𝑥𝑡 ........................... 113 XI Figure3.7(a): the exact and numerical solutions 𝑦𝑡 ...................................................... 114 Figure 3.7(b): the absolute error resulting of applying RKM of 𝑦𝑡 ............................. 114 XII COMPARISON OF NUMERICAL METHODS FOR SOLVING SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS By Abdul Fattah Boshnaq Supervisor Prof. Naji Qatanani Abstract In this work, we shed the light on the numerical handling of systems of ordinary differential equations. These systems have wide range of applications in mathematical physics, chemistry, biology, stereology, heat conducing and engineering models. After introducing some important aspects of systems of differential equations including the solvability of homogeneous and non-homogeneous systems, we focus on the numerical techniques for solving systems of differential equations. Namely; one step and multistep methods. The one step methods include Euler and Runge-Kutta methods. The multistep methods involve Adams-Bashforth method, Adams-Moulton method and the Predictor-Corrector method. The mathematical framework of these numerical methods together with their convergence properties and their error bound associated with these methods will be presented. The proposed numerical methods will be illustrated by solving some numerical examples. Numerical results show clearly that the multistep methods are more efficient and give faster convergence than other methods. Keywords: Comparison, Numerical Methods, Solving Systems. 1 Introduction Many physical and technological problems in science and engineering are modeled mathematically by differential equations and systems of differential equations. For example, the mathematical systems involving several springs. Moreover, systems of ordinary differential equations are encountered in electro- chemistry, ecology, stereology, biology and various engineering applications [22, 29]. A standard class of problem for which considerable and software exist is the initial value problem for first order systems of ordinary differential equations [3, 4, 14]. There are many contributions from well-known mathematicians involved in the numerical approximations of the initial value problems in ordinary differential equations. One of these is the Euler’s method which is regarded as the most elementary approximation technique for solving initial value problems[17, 26]. The Runge-Kutta methods which are devised by the German mathematician Carl David Runge (1856-1827) and Martin Wilhelm Kutta (1867-1944). Also, the Adomian decomposition method developed between 1970 and 1990 [2]. This method is known as a semi-analytical method which is used by George Adomian to solve ordinary differential equations and nonlinear partial differential equations [12]. The multistep method known as the fourth order Adam-Bashforth method is devised by the English astronomer and mathematician John Couch Adams (1819 -1892) and the English mathematician Francis Bashforth (1819-1912) [5], then the forth order Adams-Moulton method is devised by the American astronomer Forest Moulton (1872- 1952) [20]. The combination of an explicit method to predict and an implicit to improve the prediction is called predictor- corrector method [19]. This thesis is organized as follows: In chapter one, we introduce some basic definitions and properties including the solvability of homogenous and non-homogenous systems of ordinary differential equations. In chapter two, we present the numerical methods for solving systems of linear ordinary differential equations. These methods are the one step and multistep methods. 2 In chapter three, we solve three numerical examples involving initial value problems for first order systems of ordinary differential equations of known exact solution. In chapter four, involves some comparisons between the proposed numerical methods. 3 Chapter One Mathematical Preliminaries In this chapter we introduce some important aspects for the systems of ordinary differential equations [1, 16, 21]. 1.1 Systems of Linear Algebraic Equations A set of 𝑛 simultaneous linear algebraic equations in 𝑛 variables can be written as 𝑎11𝑥1 + 𝑎12𝑥2 + … + 𝑎1𝑛𝑥𝑛 = 𝑏1 𝑎21𝑥1 + 𝑎22𝑥2 + … + 𝑎2𝑛𝑥𝑛 = 𝑏2 (1.1) ⋮ ⋮ ⋮ ⋮ 𝑎𝑛1𝑥1 + 𝑎𝑛2𝑥2 + … + 𝑎𝑛𝑛𝑥𝑛 = 𝑏𝑛 System (1.1) can be expressed in matrix form as 𝐴𝑋 = 𝑏 (1.2) where the 𝑛 × 𝑛 matrix 𝐴 and the 𝑛-dimensional vector 𝑏 are given, and the components of the 𝑛-dimensional vector 𝑋 are to be determined [31]. If the matrix 𝐴 is nonsingular that lead to know that, thedet 𝐴 is not zero, the 𝐴−1 exists and there is a unique solution founded by multiplying each side of equation (1.2) on the left by 𝐴−1. 𝑋 = 𝐴−1𝑏 (1.3) Otherwise, the solution of equation (1.2) infinitely many solutions and 𝐴−1 not exists [13]. 1.2 Systems of Linear Differential Equations This system of linear differential equations in 𝑛 unknowns has the general form: 𝑑𝑥1 𝑑𝑡 = 𝑎11(𝑡)𝑥1 + 𝑎12(𝑡)𝑥2 + … + 𝑎1𝑛(𝑡)𝑥𝑛 + 𝑓1(𝑡) 4 𝑑𝑥2 𝑑𝑡 = 𝑎21(𝑡)𝑥1 + 𝑎22(𝑡)𝑥2 + … + 𝑎2𝑛(𝑡)𝑥𝑛 + 𝑓2(𝑡) (1.4) ⋮ ⋮ ⋮ ⋮ ⋮ 𝑑𝑥𝑛 𝑑𝑡 = 𝑎𝑛1(𝑡)𝑥1 + 𝑎𝑛2(𝑡)𝑥2 + … + 𝑎𝑛𝑛(𝑡)𝑥𝑛 + 𝑓𝑛(𝑡) We assume that the coefficients 𝑎𝑖𝑗(𝑡)and 𝑓𝑖(𝑡) are continuous on a common interval 𝐼. When𝑓𝑖(𝑡) = 0 , 𝑖 = 1,2, … , 𝑛 , the linear system (1.4) is said to be homogeneous, otherwise, it is nonhomogeneous[8]. 1.3 Initial Value Problem (IVP) Let 𝑡0 denotes a number in an interval 𝐼 and 𝑋(𝑡0) = ( 𝑥1(𝑡0) 𝑥2(𝑡0) ⋮ 𝑥𝑛(𝑡0) ) and 𝑋0 = ( 𝛾1 𝛾2 ⋮ 𝛾𝑛 ) where the𝛾𝑖(𝑡), 𝑖 = 1, 2, … , 𝑛 are given constants,then the problem 𝑋′(𝑡) = 𝐴(𝑡)𝑋(𝑡) + 𝐹(𝑡) (1.5) subject to 𝑋(𝑡0) = 𝑋0 is called an initial value problem. System (1.4) can be written into the matrix form as: 𝑑 𝑑𝑡 ( 𝑥1 𝑥2 ⋮ 𝑥𝑛 ) = ( 𝑎11(𝑡) 𝑎12(𝑡) … 𝑎1𝑛(𝑡) 𝑎21(𝑡) 𝑎22(𝑡) … 𝑎2𝑛(𝑡) ⋮ ⋮ ⋱ ⋮ 𝑎𝑛1(𝑡) 𝑎𝑛2(𝑡) … 𝑎𝑛𝑛(𝑡) )( 𝑥1 𝑥2 ⋮ 𝑥𝑛 ) + ( 𝑓1(𝑡) 𝑓2(𝑡) ⋮ 𝑓𝑛(𝑡) ) where 𝑋 = ( 𝑥1 𝑥2 ⋮ 𝑥𝑛 ) , 𝐴(𝑡) = ( 𝑎11(𝑡) 𝑎12(𝑡) … 𝑎1𝑛(𝑡) 𝑎21(𝑡) 𝑎22(𝑡) … 𝑎2𝑛(𝑡) ⋮ ⋮ ⋱ ⋮ 𝑎𝑛1(𝑡) 𝑎𝑛2(𝑡) … 𝑎𝑛𝑛(𝑡) ) , 𝐹 = ( 𝑓1(𝑡) 𝑓2(𝑡) ⋮ 𝑓𝑛(𝑡) ) (1.6) 5 Moreover, in vector form, 𝑋′(𝑡) = 𝐴𝑋(𝑡) + 𝐹 (1.7) If the system is homogeneous, then system (1.7) reduces to 𝑋′(𝑡) = 𝐴𝑋(𝑡) (1.8) A solution vector of system (1.7) is equivalent to 𝑛 scalar equations 𝑥1 = 𝛷1(𝑡) , 𝑥2 = 𝛷2(𝑡), … , 𝑥𝑛 = 𝛷𝑛(𝑡)[23]. 1.4 Existence and Uniqueness To discuss the existence and the uniqueness of a solution for the initial value problem (1.5), we consider the following definitions and theorems. Definition 1.1[9]: A function 𝑓(𝑡 , 𝑥)is said to satisfy a Lipschitz condition in the variable 𝑥 on a set 𝑉 ⊂ ℝ2 if a constant 𝐿 > 0 exists with |𝑓(𝑡 , 𝑥1) − 𝑓(𝑡 , 𝑥2)| ≤ 𝐿 |𝑥1 − 𝑥2| whenever(𝑡 , 𝑥1) and (𝑡 , 𝑥2) are in 𝑉. The constant 𝐿 is called a Lipschitz constant for 𝑓. Definition 1.2 [9]: A set 𝐷 ⊂ ℝ2 is said to be convex set if, for any two point𝑥1, 𝑥2in𝐷, the line segment joining 𝑥1(𝑡) and 𝑥2(𝑡) lies entirely with in 𝐷. Formally, 𝐷 is convex if for all 𝑥1, 𝑥2 ∈ 𝐷 and for all 𝜃 in the interval [0,1] the point 𝜃𝑥1 + (1 − 𝜃)𝑥2 ∈ 𝐷 ∀ 𝜃 ∈ [0,1] The set 𝑆 in Figure (1.1) (a) is convex because the line segment joining any two point in 𝑆 lies entirely within 𝑆 The set 𝐷 in Figure (1.1) (b) is non-convex because there exist point 𝑎 and 𝑏 in 𝐷 such that line segment joining 𝑎 and 𝑏 dose not lie entirely within 𝐷 6 Figure (1.1) is convex set S and non-convex set D Theorem 1.1[9]: Suppose 𝑓(𝑡, 𝑥) is defined on a convex set 𝑉 ⊂ ℝ2. If a constant 𝐿 > 0 exists with | 𝜕𝑓 𝜕𝑥 (𝑡, 𝑥)| ≤ 𝐿 for all (𝑡, 𝑥) ∈ 𝑉 then 𝑓 satisfies a lipschitz condition on 𝑉 in the variable 𝑦 with a Lipschitz constant 𝐿. Proof: Holding 𝑡 constant and applying the Mean Value Theoremto the function 𝑓(𝑡, 𝑥), when 𝑥1 < 𝑥2, a number 𝛼 in (𝑥1, 𝑥2) exists with lim 𝑥2→𝑥1 𝑓(𝑡 , 𝑥2) − 𝑓(𝑡 , 𝑥1) 𝑥2 − 𝑥1 = 𝜕𝑓 𝜕𝑥 (𝑡, 𝑥 ) and |𝑓(𝑡 , 𝑥2) − 𝑓(𝑡 , 𝑥1)| = |𝑥2 − 𝑥1| | 𝜕𝑓 𝜕𝑥 (𝑡, 𝑥)| ≤ |𝑥2 − 𝑥1|𝐿. Thus, 𝑓 satisfies a Lipschitz condition on 𝑉 in the variable 𝑥 with a Lipschitz constant 𝐿[28]. Definition 1.3[9]: The initial value problem 𝑑𝑥 𝑑𝑡 = 𝑓(𝑡, 𝑥(𝑡)), 𝑎 ≤ 𝑡 ≤ 𝑏, 𝑥(𝑎) = 𝛼 7 is said to be a well-posed problem if: A unique solution, 𝑥(𝑡)to the problem exists and there exist constants 𝜖0 > 0 and 𝑘 > 0, such that for any 𝜖, in (0, 𝜖0) wherever δ(t) is continuous with |δ(t) | < 𝜖 for all t in [𝑎, 𝑏], and when |δ0| < 𝜖, the initial value problem 𝑑𝑦 𝑑𝑡 = 𝑓(𝑡, 𝑦) + δ(t), 𝑎 ≤ 𝑡 ≤ 𝑏, 𝑦(𝑎) = 𝛼 + δ0 has a unique solution 𝑦(𝑡) that satisfies | 𝑦(𝑡) − 𝑥(𝑡)| < 𝑘𝜖 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑡 𝑖𝑛 [𝑎, 𝑏] The second part of the definition says that small perturbations of the original problem and small perturbations of the initial condition have only small error effects on the approximated solution. Theorem 1.2[9]: Suppose that 𝑉 = {(𝑡, 𝑥): 𝑎 ≤ 𝑡 ≤ 𝑏 𝑎𝑛𝑑 − ∞ < 𝑥 < ∞} . If 𝑓(𝑡, 𝑥) is continuous and satisfies a lipschitz condition in the variable 𝑥 on the set 𝑉 then the initial value problem 𝑦′(𝑥) = 𝑓(𝑥, 𝑦), 𝑎 ≤ 𝑥 ≤ 𝑏, 𝑦(𝑎) = 𝑎 is well posed. We note, that this theorem ensures the uniqueness of the solution 𝑥(𝑡)for 𝑎 ≤ 𝑡 ≤ 𝑏. Theorem 1.3[10]: Suppose that 𝑓 and 𝜕𝑓 𝜕𝑥 its first partial derivative with respect to 𝑥, are continuous for 𝑡 in [𝑎, 𝑏] and for all 𝑥. Then the initial value problem 𝜕𝑥 𝜕𝑡 = 𝑓(𝑡, 𝑥(𝑡)), 𝑎 ≤ 𝑡 ≤ 𝑏, with initial condition 𝑥(𝑎) = 𝛼 has a unique solution 𝑥(𝑡) for 𝑎 ≤ 𝑡 ≤ 𝑏, and the problem is well posed. 8 Proof: the set 𝑉 = {(𝑡, 𝑥): 𝑎 ≤ 𝑡 ≤ 𝑏 𝑎𝑛𝑑 − ∞ < 𝑥 < ∞} is convex. Since 𝜕𝑓 𝜕𝑥 is continuous on [𝑎, 𝑏], then 𝜕𝑓 𝜕𝑥 is bounded. Therefore, there exists a real number 𝐿 > 0 such that | 𝜕𝑓 𝜕𝑥 | ≤ 𝐿. Thus, 𝑓 satisfies a Lipschitz condition on 𝑉 in the variable 𝑥. Also 𝑓 is continuous. It follows from theorem (1.2), that the IVP is well-posed [25]. 1.5 Linear Dependence and Linear Independence The set 𝑋 is linearly dependent on the interval 𝐼 if there exist real or complex constants𝑐1, 𝑐2, … , 𝑐𝑘at least one of which is nonzero, such that 𝑐1𝑋1 + 𝑐2𝑋2 + …+ 𝑐𝑘𝑋𝑘 = 0 (1.9) On the other hand, if 𝑐1 = 𝑐2 = ⋯ = 𝑐𝑘 = 0 , then 𝑋1, 𝑋2, … , 𝑋𝑘 are said to be linearly independent[8]. To consider the solution ofthe system of ordinary differential equations linearly dependence or independence we take the Wronskian. Theorem 1.4[31]: Criterion for linearly independent solution Let 𝑋1 = ( 𝑥11 𝑥21 ⋮ 𝑥𝑛1 ), 𝑋2 = ( 𝑥12 𝑥22 ⋮ 𝑥𝑛2 ), … ,𝑋𝑛 = ( 𝑥1𝑛 𝑥2𝑛 ⋮ 𝑥𝑛𝑛 ) be𝑛 solution vectors of the homogenous system (1.8) on an interval 𝐼. Then the set of solution vectors is linearly independent on interval𝐼 if and only if the Wronskian 𝑊(𝑋1, 𝑋2, … , 𝑋𝑛) = | 𝑥11 𝑥12 … 𝑥1𝑛 𝑥21 𝑥22 … 𝑥2𝑛 ⋮ ⋮ ⋱ ⋮ 𝑥𝑛1 𝑥𝑛2 … 𝑥𝑛𝑛 | ≠ 0 ( 1.10) for every 𝑡 in the interval𝐼. If 𝑋1, 𝑋2, … , 𝑋𝑛 are solution vectors of (1.8) and 𝑊(𝑋1, 𝑋2, … , 𝑋𝑛) ≠ 0 for some 𝑡0 in, then the solutions are linearly independent on the interval𝐼. 9 1.6 Homogeneous Vector Differential Equations Homogeneous vector differential equations are a class of differential equations where the function and its derivatives are vector-valued, and the equations are linear and homogeneous. They often appear in the context of systems of linear differential equations with vector solutions. The meaning of homogenous differential equations is the right side is zero, in the other words, it is in the form: 𝑋′(𝑡) = 𝐴𝑋(𝑡) To solve the system, find the eigenvalues and eigenvectors of the matrix 𝐴 , the eigenvalues 𝜆’s are found by solving the characteristic equation. 𝑑𝑒𝑡(𝐴 − 𝜆𝐼) = 0 (1.11) where𝐼 is the identity matrix. For each eigenvalue 𝜆, find the corresponding eigenvector 𝑉 by solving (𝐴 − 𝜆𝐼) = 0⃗ (1.12) the general solution to the system can be expressed in term of the eigenvalues and eigenvectors. If 𝜆𝑖 is an eigenvalue with corresponding eigenvector 𝑉𝑖, then the solution associated with eigenvalue𝜆𝑖 is 𝑋𝑖(𝑡) = 𝑒𝜆𝑖𝑡𝑉𝑖 (1.13) the general solution is a linear combination of these solutions 𝑋(𝑡) = 𝐶1𝑒 𝜆1𝑡𝑉1 + 𝐶2𝑒 𝜆2𝑡𝑉2 + ⋯+ 𝐶𝑛𝑒 𝜆𝑛𝑡𝑉𝑛 where 𝐶1, 𝐶2, … , 𝐶𝑛are constants determined by initial conditions. 10 1.7 Eigenvalues and Eigenvectors Suppose 𝐾𝑒𝜆𝑡 be a solution vector of the homogenous system (1.8), then after differentiation and rearranging we obtain (𝐴 − 𝜆𝐼)𝐾 = 0 The matrix equation (1.12) is equivalent to the simultaneous algebraic equations (𝑎11 − 𝜆)𝑘1 + 𝑎12𝑘2 + … + 𝑎1𝑛𝑘𝑛 = 0 𝑎21𝑘1 + (𝑎22 − 𝜆)𝑘2 + … + 𝑎2𝑛𝑘𝑛 = 0 ⋮ ⋮ ⋱ ⋮ 𝑎𝑛1𝑘1 + 𝑎𝑛2𝑘2 + … +(𝑎𝑛𝑛 − 𝜆)𝑘𝑛 = 0 the nontrivial solution need to found it 𝑑𝑒𝑡(𝐴 − 𝜆𝐼) = 0 1.8 Kind of Eigenvalues 1.8.1 Distinct real eigenvalues When an 𝑛 × 𝑛 matrix 𝐴 has 𝑛 distinct real eigenvalues𝜆1, 𝜆2, … , 𝜆𝑛, it implies, the matrix 𝐴 has 𝑛 linearly independent eigenvectors 𝐾1, 𝐾2, … , 𝐾𝑛 corresponding to the eigenvalues, the solutions of the system of differential equations given by 𝑋′ = 𝐴𝑋 can be expressed in the form 𝑋1(𝑡) = 𝐾1𝑒 𝜆1𝑡 , 𝑋2(𝑡) = 𝐾2𝑒 𝜆2𝑡 , … , 𝑋𝑛(𝑡) = 𝐾𝑛𝑒 𝜆𝑛𝑡 these solutions represent a fundamental set of solutions on the interval (−∞,∞) [31]. And since 𝐾1, 𝐾2, … , 𝐾𝑛 are linearly independent, the entire set {𝑋1, 𝑋2, … , 𝑋𝑛} forms a basis for the solution space of the system. 1.8.2 Repeated eigenvalues The eigenvalues 𝜆1, 𝜆2, … , 𝜆𝑛 of an 𝑛 × 𝑛 matrix 𝐴 my be not distinct, some of them my be repeated, for example 𝜆1 = 𝜆2 the second solution can be expressed as the form 𝑋2(𝑡) = 𝐾𝑡𝑒𝜆1𝑡 + 𝑃𝑒𝜆1𝑡 (1.14) 11 where 𝐾 = ( 𝑘1 𝑘2 ⋮ 𝑘𝑛 ) and 𝑃 = ( 𝑝1 𝑝2 ⋮ 𝑝𝑛 ). we substitute (1.14) into the system 𝑋′ = 𝐴𝑋 and simplify: (𝐴𝐾 − 𝜆1𝐾)𝑡𝑒𝜆1𝑡 + (𝐴𝑃 − 𝜆1𝑃 − 𝐾)𝑒𝜆1𝑡 = 0 since this last equation is to hold for all values of 𝑡, we must have (𝐴 − 𝜆1𝐼)𝐾 = 0 (1.15) (𝐴 − 𝜆1𝐼)𝑃 = 𝐾 (1.16) the first solution 𝑋1 is solving system (1.15) for the vector 𝐾, which the first solution 𝑋1 = 𝐾𝑒𝜆1𝑡, the second solution solving system(1.16) for the vector 𝑃 . 1.8.3 Complex Eigenvalues If 𝜆1 = 𝛼 + 𝑖𝛽and 𝜆2 = 𝛼 − 𝑖𝛽,𝛽 > 0,are complex eigenvalues of the coefficient matrix 𝐴 , the eigenvalue 𝜆2 is the complex conjugate of 𝜆1 . If 𝐾1 is the eigenvector corresponding to 𝜆1, then 𝐾2 corresponding to 𝜆2, will often be the complex conjugate of 𝐾1. Thus, both eigenvectors will have complex components. we defined 𝐁𝟏 = 1 2 (𝐾1 + 𝐾1) and 𝐁𝟐 = 𝑖 2 (−𝐾1 + 𝐾1) (1.17) The linearly independent solutions on (1.8) is 𝑋1(𝑡) = [𝐁𝟏 cos 𝛽𝑡 – 𝐁𝟐 sin 𝛽𝑡]𝑒𝛼𝑡 (1.18) 𝑋2(𝑡) = [𝐁𝟐 cos 𝛽𝑡 + 𝐁𝟏 sin 𝛽𝑡]𝑒𝛼𝑡 1.9 Nonhomogeneous Vector Differential Equations Nonhomogeneous vector differential equations are a generalization of homogeneous vector differential equations that include a non-zero forcing term or external input. 12 To solve system (1.7) we typically find the general solution by combining the solution to the corresponding homogeneous system and a particular solution to the nonhomogeneous system. To find a particular solution 𝑋𝑃(𝑡) to the nonhomogeneous system, there are various methods to find a particular solution, including: Method of Undetermined Coefficients and Variation of Parameters. Theorem 1.5 [31]: In the nonhomogeneous vector differential equations (1.5) if {𝑋1, 𝑋2, … , 𝑋𝑛} is a fundamental solution set of the corresponding homogenous equations (1.8), then every solution 𝑋(𝑡) to the nonhomogeneous equations can be expressed in the form 𝑋(𝑡) = 𝐶1𝑋1(𝑡) + 𝐶2𝑋2(𝑡) + ⋯+ 𝐶𝑛𝑋𝑛(𝑡) + 𝑋𝑃(𝑡) Where 𝑋𝑃(𝑡) is a particular solution. Proof: Consider the homogenous equation 𝑋′(𝑡) = 𝐴(𝑡)𝑋(𝑡) , the fundamental solution set {𝑋1, 𝑋2, … , 𝑋𝑛} consists of solutions that form a basis for the solution space of the homogenous equation, any solution 𝑋𝐶(𝑡) of the homogenous equation can be expressed as 𝑋𝐶(𝑡) = 𝐶1𝑋1(𝑡) + 𝐶2𝑋2(𝑡) + ⋯+ 𝐶𝑛𝑋𝑛(𝑡) where, 𝐶1, 𝐶2, … , 𝐶𝑛 are constants determined by initial conditions. Now, let 𝑋𝑃(𝑡) be a particular solution to the nonhomogeneous equation 𝑋𝑃 ′(𝑡) = 𝐴(𝑡)𝑋𝑃(𝑡) + 𝑏(𝑡) consider any solution 𝑋(𝑡)to the nonhomogeneous equation,we can express it as the sum of the complementary solution and a particular solution 𝑋(𝑡) = 𝑋𝐶(𝑡) + 𝑋𝑃(𝑡) 13 Then we can write 𝑋(𝑡) = 𝐶1𝑋1(𝑡) + 𝐶2𝑋2(𝑡) + ⋯+ 𝐶𝑛𝑋𝑛(𝑡) + 𝑋𝑃(𝑡) 1.10 The Method of Undetermined Coefficients This method is technique used to find a particular solution 𝑥𝑝 to non-homogenous equation (1.7) To choose a form for a particular solution 𝑋𝑝(𝑡) based on the form of (𝑡) , the guess should be similar to the nonhomogeneous term 𝐹(𝑡) , but with undetermined coefficients, the common forms;  Polynomial: if 𝐹(𝑡) is a polynomial of degree 𝑛, guess a polynomial of the same degree, if 𝐹(𝑡) = 4𝑡2 + 2𝑡 + 1 the suggested particular solution will be 𝑥𝑝(𝑡) = 𝐴𝑡2 + 𝐵𝑡 + 𝐶.  Exponential: if (𝑡) = 𝑒𝑟𝑡 , the suggested particular solution will be 𝑥𝑝(𝑡) = 𝐴𝑒𝑟𝑡.  Sine or Cosine: if 𝐹(𝑡) = sin (𝑟𝑡)or 𝐹(𝑡) = cos (𝑟𝑡), the suggested particular solution will be 𝑥𝑝(𝑡) = 𝐴 sin(𝑟𝑡) + 𝐵 cos (𝑟𝑡).  Combination: if 𝐹(𝑡) is a combination of function, combine the respective supposing, if 𝐹(𝑡) = 𝑒2𝑡sin (𝑡)the suggested particular solution will be 𝑥𝑝(𝑡) = 𝑒2𝑡(𝐴 sin(𝑟𝑡) + 𝐵 cos(𝑟𝑡)) 1.11 Variation of Parameters After finding the homogenous solution of a system, find a fundamental matrix solution 𝛷(𝑡), which consists of the solutions to the homogenous system if 𝑋1, 𝑋2, … , 𝑋𝑛 are the linearly independent solutions, then the fundamental matrix is 𝛷(𝑡) = [𝑋1, 𝑋2, … , 𝑋𝑛] assume a particular solution of the form 𝑋𝑃(𝑡) = 𝛷(𝑡)𝑈(𝑡) where𝑈(𝑡) is a vector of functions to be determined, differentiate 𝑋𝑃(𝑡) 𝑋𝑃 ′ (𝑡) = 𝛷(𝑡)𝑈′(𝑡) + 𝛷′(𝑡)𝑈(𝑡) 14 subsisting𝑋𝑃(𝑡) into the original equation 𝛷(𝑡)𝑈′(𝑡) + 𝛷′(𝑡)𝑈(𝑡) = 𝐴(𝑡)𝛷(𝑡)𝑈(𝑡) + 𝐹(𝑡) from the homogenous system 𝛷′(𝑡) = 𝐴(𝑡)𝛷(𝑡) thee equations simplifies to 𝛷(𝑡)𝑈′(𝑡) = 𝐹(𝑡) multiply both side by 𝛷−1(𝑡) 𝑈′(𝑡) = 𝛷−1(𝑡)𝐹(𝑡) integrate 𝑈′(𝑡) 𝑈(𝑡) = ∫𝛷−1(𝑡) 𝐹(𝑡)𝑑𝑡 substitute 𝑈(𝑡) back into the expression for 𝑋𝑃(𝑡) 𝑋𝑃(𝑡) = 𝛷(𝑡)∫𝛷−1(𝑡) 𝐹(𝑡)𝑑𝑡 the general solution to the non-homogenous system is 𝑋(𝑡) = 𝑋𝐶(𝑡) + 𝑋𝑃(𝑡) 1.12 Higher Order Initial Value Problems Higher order initial value problem involve differential equations of order greater than one. The higher order initial value problem typically takes the form 𝑦(𝑛)(𝑡) = 𝑓 (𝑡 , 𝑦(𝑡) , 𝑦′(𝑡), … , 𝑦(𝑛−1)(𝑡)) 𝑎 ≤ 𝑡 ≤ 𝑏 15 with initial conditions 𝑦(𝑡0) = 𝑦0 𝑦′(𝑡0) = 𝑦1 ⋮ 𝑦(𝑛−1)(𝑡0) = 𝑦𝑛−1 Methods of solution 1. Reduction of order: higher order equations can often be reduced to a system of first order equations. 2. Numerical methods: when analytical solutions are difficult to find, numerical methods such as Runge-Kutta methods, Multistep method are used. 3. Series solutions: in some cases, power series or Taylor series expansions can be used to approximate the solution. 4. Laplace transforms: this technique transforms the differential equation into an algebraic equation, which can then be solved for the Laplace transform of the function, followed by inverse transformation to find the solution. 16 Chapter Two Numerical Methods for Solving System of Differential Equation In this chapter, we will some of the numerical methods used for solving systems of ordinary differential equations, namely; one step and multistep methods. 2.1 The Adomian Decomposition Method The Adomian decomposition method(ADM), also known as the inverse operator method, is a mathematical effective approach for solving broad classes of linear and nonlinear mathematical physics equations with important applications in different fields of applied mathematics, engineering, physics and biology, it was proposed by George Adomian (1986, 1988, 1994)[1, 3,4]. The Adomian decomposition method involves the decomposing the unknown equation 𝐹(𝑡) = 𝑔(𝑡) into the sum of the components of different degree solution, and also to find the solution of each order, for these sum’s we want to approximate the true solution to a desired accuracy[30]. First, the whole equation is decomposed into several parts, which is linear and nonlinear, the linear part is divided into invertible and residual of the linear operator and other part is nonlinear [6,7]. The operator 𝐹 is decomposed as 𝐹 = 𝐿 + 𝑅 + 𝑁 (2.1) where𝐿 is linear part, 𝑅 is remainder term and 𝑁 is nonlinear term. therefore 𝐹𝑢(𝑡) = 𝑔(𝑡) = 𝐿𝑢 + 𝑅𝑢 + 𝑁𝑢 (2.2) where𝐿 is invertible we got 𝐿𝑢 = 𝑔(𝑡) − 𝑅𝑢 − 𝑁 we take suitable operator 𝐿 such that 𝐿−1 exist and take 𝐿−1 to both sidesgives 17 𝑢 = 𝐿−1𝑔 − 𝐿−1𝑅𝑢 − 𝐿−1𝑁𝑢 (2.3) we decompose the true solution 𝑛 into the sum of infinitely component 𝑢 = ∑ 𝑢𝑛 ∞ 𝑛=0 the nonlinear term evaluated to ∑ 𝐴𝑛 ∞ 𝑛=0 where 𝐴𝑛 are special polynomial to be discussed, so the equation ∑ 𝑢𝑛 ∞ 𝑛=0 = 𝐿−1𝑔 − 𝐿−1𝑅 ∑ 𝑢𝑛 ∞ 𝑛=0 − 𝐿−1 ∑ 𝐴𝑛 ∞ 𝑛=0 (2.4) we can write as 𝑢0 = 𝐿−1𝑔 𝑢1 = −𝐿−1𝑅𝑢0 − 𝐿−1𝐴0 𝑢2 = −𝐿−1𝑅𝑢1 − 𝐿−1𝐴1 𝑢𝑚+1 = −𝐿−1𝑅𝑢𝑚 − 𝐿−1𝐴𝑚 Now for solving system of linear and nonlinear ordinary differential equations consider the following system of ordinary differential equations [15]: 𝑦𝑖 ′(𝑥) = ∑𝑏𝑖𝑗(𝑥)𝑦𝑖 𝑛 𝑗=1 + 𝑁𝑖(𝑥, 𝑦1, 𝑦2, … , 𝑦𝑛) + 𝑔𝑖(𝑥) (2.5) withinitial conditions 𝑦𝑖(0) = 𝑐𝑖 𝑖 = 1,2, … , 𝑛 where 𝑏𝑖𝑗(𝑥), 𝑔𝑖(𝑥) ∈ [0, 𝑇] and 𝑁𝑖’s are nonlinear continuous functions of its argument. Integrating both side of equation (2.1) from 0 to 𝑥 and the using initial conditions, we got 𝑦𝑖(𝑥) = 18 𝑐𝑖 + ∫ 𝑔𝑖(𝑥) 𝑥 0 𝑑𝑥 + ∫ ∑𝑏𝑖𝑗(𝑥)𝑦𝑖 𝑛 𝑗=1 𝑥 0 𝑑𝑥 + ∫ 𝑁𝑖(𝑥, 𝑦1, 𝑦2, … , 𝑦𝑛) 𝑥 0 𝑑𝑥 (2.6) where 𝑖 = 1,2, … , 𝑛. The standard ADM yields the solution 𝑦𝑖(𝑥) by series 𝑦𝑖(𝑥) = ∑ 𝑦𝑖𝑚(𝑥) ∞ 𝑚=0 (2.7) and the nonlinear terms by an infinite series of Adomian polynomials 𝑁𝑖(𝑥, 𝑦1, 𝑦2, … , 𝑦𝑛) = ∑ 𝐴𝑖𝑚 ∞ 𝑚=0 (2.8) where 𝐴𝑖𝑚(𝑦10, 𝑦11, … , 𝑦1𝑚, 𝑦21, … , 𝑦2𝑚, … , 𝑦𝑛0, 𝑦𝑛1, … , 𝑦𝑛𝑚). These polynomials can be constructed using the general formula [2] 𝐴𝑖𝑚 = 1 𝑚! 𝑑𝑚 𝑑𝜆𝑚 [𝑁𝑖 (𝑥, ∑ 𝑦1𝑚 ∞ 𝑚=0 𝜆𝑚, … , ∑ 𝑦𝑛𝑚 ∞ 𝑚=0 𝜆𝑚)] 𝜆=0 = [ 1 𝑚! 𝑑𝑚 𝑑𝜆𝑚 𝑁𝑖 (𝑥, ∑ 𝑦1𝑚 ∞ 𝑚=0 𝜆𝑚, … , ∑ 𝑦𝑛𝑚 ∞ 𝑚=0 𝜆𝑚)] 𝜆=0 and the components 𝑦𝑖𝑚 , 𝑚 ≥ 0, can be determined in a recursive manner. In view of (2.6) - (2.8), the ADM defines the component 𝑦𝑖𝑚 , 𝑚 ≥ 0,by the following recursion relation : 𝑦𝑖0(𝑥) = 𝑐𝑖 + ∫ 𝑔𝑖(𝑥) 𝑥 0 𝑑𝑥 , 𝑖 = 1,2, … , 𝑛 . (2.9) 𝑦𝑖,𝑚+1(𝑥) = ∫ ∑𝑏𝑖𝑗(𝑥)𝑦𝑖𝑚 𝑛 𝑗=1 𝑥 0 𝑑𝑥 + ∫ 𝐴𝑖𝑚(𝑥) 𝑥 0 𝑑𝑥 , 𝑚 = 0 , 1 , … we approximate the solution 𝑦𝑖(𝑥) by the truncated series 19 𝑓𝑖𝑘(𝑥) = ∑ 𝑦𝑖𝑚(𝑥) 𝑘−1 𝑚=0 (2.10) and lim 𝑘→∞ 𝑓𝑖𝑘(𝑥) = 𝑦𝑖(𝑥), 𝑖 = 1,2, … , 𝑛 The advantages of Adomian decomposition method:  Simplicity: The method is straightforward to implement, requiring minimal computational effort for obtaining terms in the series.  No Linearization: It handles nonlinear problems directly without needing to linearize the equations, preserving the original system's characteristics.  Convergence: The series solution often converges quickly, providing accurate approximations even with a few terms.  Flexibility: Applicable to various types of differential equations, including initial and boundary value problems. The disadvantages of Adomian decomposition method:  Convergence Issues: While it often converges, some problems may exhibit slow or even non-convergence, depending on the nature of the nonlinearity.  Complexity for Higher Orders: As the order of the series increases, the calculations for higher-order Adomian polynomials can become complex and cumbersome.  Initial Condition Sensitivity: The method may be sensitive to initial conditions, which can affect solution accuracy.  Limited Applicability: It may not be suitable for all types of nonlinear equations, particularly those with discontinuities or singularities.  Error Estimation: Estimating the error of the approximation can be challenging, making it difficult to assess solution accuracy. 20 2.2 Runge-KuttaMethod The Runge-Kutta methods are a family of iterative techniques used for approximating the solutions of ordinary differential equations. 2.2.1 Euler’s method The Euler method is one of the simplest and most straightforward techniques for numerically solving ordinary differential equations (ODE’s). It's often used as an introductory method in numerical analysis due to its simplicity. The Euler method is designed to approximate the solution of an ODE with initial condition of the form: 𝑑𝑦 𝑑𝑡 = 𝑓(𝑥, 𝑦) 𝑦(𝑥0) = 𝑦0 (2.11) on some interval [𝑎, 𝑏]. We want to approximate 𝑦(𝑥) at the mesh point 𝑥𝑖 = 𝑎 + 𝑖ℎ with step size ℎ = (𝑏 − 𝑎)/𝑁, where 𝑖 = 0 , 1 , …𝑁 − 1. For each 𝑥𝑖 ,the Taylor expansion of 𝑦(𝑥) around 𝑥𝑖 is given by 𝑦(𝑥𝑖+1) = 𝑦(𝑥𝑖) + ℎ 𝑑𝑦 𝑑𝑥 | 𝑥𝑖 + ℎ2 2! 𝑑2𝑦 𝑑𝑥2 | 𝑥𝑖 + ⋯ assuming that 𝑦(𝑥) is twice continuously differentiable on [𝑎, 𝑏], we get 𝑦(𝑥𝑖+1) = 𝑦(𝑥𝑖 + ℎ) = 𝑦(𝑥𝑖) + ℎ𝑦′(𝑥𝑖) + ℎ2 2 𝑦′′(𝜖) (2.12) for some 𝜖 between 𝑥𝑖 and 𝑥𝑖 + ℎ we neglect the error term in equation (2.12) and using equation (2.11), we obtain the formula 𝑦(𝑥𝑖 + ℎ) ≈ 𝑦(𝑥𝑖) + ℎ𝑓(𝑥𝑖 , 𝑦𝑖) 21 if we denote 𝑦𝑖 ≈ 𝑦(𝑥𝑖) then 𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖) (2.13) equation (2.13) is known as Euler’s method. TheTaylor’s method of order 𝑘 given by the formula 𝑦𝑖+1 = 𝑦𝑖 + ℎ𝜑𝑘(𝑥𝑖 , 𝑦𝑖)(2.14) where 𝜑𝑘(𝑥𝑖 , 𝑦𝑖) = 𝑓(𝑥𝑖 , 𝑦𝑖) + ℎ 2! 𝑓′(𝑥𝑖 , 𝑦𝑖) + ⋯+ ℎ𝑘−1 𝑘! 𝑓(𝑘−1)(𝑥𝑖 , 𝑦𝑖) 2.2.2 Runge-Kutta Methods )RKM) Generalizing the Euler method involves enhancing its accuracy by allowing multiple evaluations of the derivative within a single step. This approach leads to the development of more sophisticated numerical methods for solving ordinary differential equations. The generally attributed to Runge(1895) [24], which him forming the basis of what we now call Runge-Kutta methds, Heun(1900)[18], he introduced what is now know as Heun’s method. A simple second-order Runge-Kutta method, Kutta (1901) [24] fully characterized fourth-order Runge-Kutta methods and proposed the first explicit methods of this order [11], Nyström (1925) Contributed methods for second-order differential equations and advanced the development of methods for first-order equations, Hutta (1956,1957)[27] Introduced sixth-order Runge-Kutta methods, extending the range of accuracy for numerical solutions. The Runge-Kutta methods provide a practical approach for numerically solving ordinary differential equations without requiring the differentiation of the function𝑓(𝑥 , 𝑦 ). The simplest of these methods, the Runge-Kutta method of order 2 (RK2), offers a balance between simplicity and accuracy. The formula for RK2 can be expressed as 𝑦𝑖+1 = 𝑦𝑖 + 𝛼𝑘1 + 𝛽𝑘2 (2.15) 22 where 𝑘1 = ℎ𝑓(𝑥𝑖 , 𝑦𝑖) (2.16) 𝑘2 = ℎ𝑓(𝑥𝑖 + 𝛼ℎ , 𝑦𝑖 + 𝛽𝑘1) (2.17) the common choice for the coefficients in the RK2 method is 𝛼 = 𝛽 = 1 2 thus, we can rewrite the 𝑘2 term as 𝑘2 = ℎ𝑓 (𝑥𝑖 + ℎ 2 , 𝑦𝑖 + 1 2 𝑘1) (2.18) this leads to the Runge-Kutta method of order 2, sometimes known as the modified Euler method 𝑦𝑖+1 = 𝑦𝑖 + ℎ 2 [𝑓(𝑥𝑖 , 𝑦𝑖) + 𝑓(𝑥𝑖 + ℎ , 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖))] (2.19) or 𝑦𝑖+1 = 𝑦𝑖 + ℎ 2 (𝑘1 + 𝑘2) (2.20) The fourth-order Runge-Kutta method (RK4) is indeed one of the most accurate and widely used numerical techniques for solving ordinary differential equations. To derive the RK4 formula, we need to use multiple evaluations of the function 𝑓 to capture more information about the functions behavior calculate 𝑘1: 𝑘1 = ℎ 𝑓(𝑥0 , 𝑦0) calculate 𝑘2: 𝑘2 = ℎ 𝑓 (𝑥0 + ℎ 2 , 𝑦0 + 1 2 𝑘1) 23 calculate 𝑘3: 𝑘3 = ℎ 𝑓 (𝑥0 + ℎ 2 , 𝑦0 + 1 2 𝑘2) calculate 𝑘4: 𝑘4 = ℎ 𝑓(𝑥0 + ℎ , 𝑦0 + 𝑘3) the RK4 formula is given by 𝑦1 = 𝑦0 + ℎ 6 (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4) The advantages of Runge-Kutta methods  Simplicity: It's easy to implement and understand, making it accessible for various applications.  Higher Accuracy: Higher-order Runge-Kutta methods provide improved accuracy compared to simpler methods like Euler's method, particularly for stiff equations.  Stability: It offers better stability properties for a wide range of problems, reducing numerical errors in solutions.  Flexibility: The method can be adapted for both ordinary and partial differential equations, as well as systems of equations.  No Need for Complex Formulations: Unlike some methods, it doesn’t require the derivation of complex formulas for each problem, streamlining the process. The disadvantages of Runge-Kutta methods  Computational Cost: Higher-order Runge-Kutta methods require more function evaluations per step, which can increase computational time and resource usage.  Stiff Equations: It can struggle with stiff differential equations, where alternative methods like implicit methods for example Backward Differentiation Formula may perform better.  Global Error: The method may accumulate significant global errors, especially when used with large step sizes, which can affect long-term accuracy.  Step Size Selection: Choosing an appropriate step size is crucial; too large can lead to instability and inaccuracies, while too small can lead to excessive computation. 24 2.2.3 System of first order differential equations Consider an initial value problem of the first order differential equation: 𝑑𝑥 𝑑𝑡 = 𝑓(𝑡, 𝑥, 𝑦) 𝑑𝑦 𝑑𝑡 = 𝑔(𝑡, 𝑥, 𝑦) with 𝑥(𝑡0) = 𝑥0 , 𝑦(𝑡0) = 𝑦0 and 𝑡0 ≤ 𝑡 ≤ 𝑡𝑛 , this system consists a single of ordinary differential equations, the solution domain is discredited such that 𝑡0 , 𝑡1 = 𝑡0 + ℎ , 𝑡2 = 𝑡0 + 2ℎ , 𝑡𝑛 = 𝑡0 + 𝑛ℎ , where ℎ is the step size of 𝑡. The solution that is obtained by the RK4 method is given by 𝑥𝑛+1(𝑡) = 𝑥𝑛(𝑡) + ℎ 6 (𝑓1 + 2𝑓2 + 2𝑓3 + 𝑓4) 𝑦𝑛+1(𝑡) = 𝑦𝑛(𝑡) + ℎ 6 (𝑔1 + 2𝑔2 + 2𝑔3 + 𝑔4) where 𝑓1 = ℎ 𝑓(𝑡𝑛, 𝑥𝑛, 𝑦𝑛) 𝑔1 = ℎ 𝑔(𝑡𝑛, 𝑥𝑛, 𝑦𝑛) 𝑓2 = ℎ 𝑓(𝑡𝑛 + ℎ 2 , 𝑥𝑛 + 𝑓1 2 , 𝑦𝑛 + 𝑔1 2 ) 𝑔2 = ℎ 𝑔(𝑡𝑛 + ℎ 2 , 𝑥𝑛 + 𝑓1 2 , 𝑦𝑛 + 𝑔1 2 ) 𝑓3 = ℎ 𝑓(𝑡𝑛 + ℎ 2 , 𝑥𝑛 + 𝑓2 2 , 𝑦𝑛 + 𝑔2 2 ) 𝑔3 = ℎ 𝑔(𝑡𝑛 + ℎ 2 , 𝑥𝑛 + 𝑓2 2 , 𝑦𝑛 + 𝑔2 2 ) 𝑓4 = ℎ 𝑓(𝑡𝑛 + ℎ, 𝑥𝑛 + 𝑓3 , 𝑦𝑛 + 𝑔3) 𝑔4 = ℎ 𝑔(𝑡𝑛 + ℎ, 𝑥𝑛 + 𝑓3 , 𝑦𝑛 + 𝑔3) 25 2.3 Adams-Bashforth Method (ABM) The Adams-Bashforth method is a family of explicit multistep methods used to solve ordinary differential equations. These methods use information from previous time steps to estimate the value at the next time step. If we integrating the solution of the initial value problem (2.13) between 𝑥𝑖 and 𝑥𝑖+1 , 𝑦(𝑥𝑖+1) − 𝑦(𝑥𝑖) = ∫ 𝑓(𝑥 , 𝑦(𝑥)) 𝑥𝑖+1 𝑥𝑖 𝑑𝑥 thus 𝑦(𝑥𝑖+1) = 𝑦(𝑥𝑖) + ∫ 𝑓(𝑥 , 𝑦(𝑥) ) 𝑥𝑖+1 𝑥𝑖 𝑑𝑥 (2.21) to carry out the integration in (2.21) we can use a finite-difference method to approximate 𝑓(𝑥 , 𝑦) at some of the data points 𝑥0, 𝑥1, . . . , 𝑥𝑖 .This will lead to the formula 𝑦𝑖+1 = 𝑦𝑖 + ∫ 𝑝(𝑥) 𝑥𝑖+1 𝑥𝑖 𝑑𝑥 (2.22) where 𝑦(𝑥𝑖) ≈ 𝑦(𝑥𝑖+1) and 𝑝(𝑥) ≈ 𝑓(𝑥 , 𝑦). The Newton backward-difference polynomial 𝑝𝑚−1(𝑥) based on 𝑚point (𝑥𝑖 , 𝑦(𝑥𝑖)) , (𝑥𝑖−1 , 𝑦(𝑥𝑖−1)), . . . , (𝑥𝑖−𝑚+1 , 𝑦(𝑥𝑖−𝑚+1)) can be expressed as 𝑝𝑚−1(𝑥) = 𝑦𝑖 + (𝑥 − 𝑥𝑖) ∇ 𝑦𝑖 + (𝑥 − 𝑥𝑖)(𝑥 − 𝑥𝑖−1) 2! ∇2𝑦𝑖 + ⋯ + (𝑥 − 𝑥𝑖)(𝑥 − 𝑥𝑖−1)… (𝑥 − 𝑥𝑖−𝑚+1) 𝑚! ∇𝑚𝑦𝑖 where∇𝑘𝑦𝑖 is the 𝑘𝑡ℎ backward difference of 𝑦. Substituting 𝑝𝑚−1(𝑥) into (2.11), we evaluate 𝑥𝑖+1 𝑝𝑚−1(𝑥𝑖+1) = 𝑦𝑖+1 26 this gives 𝑦𝑖+1 = 𝑦𝑖 + ℎ(𝑓(𝑥𝑖 , 𝑦𝑖) + ℎ 2! ∇𝑓(𝑥𝑖) + ℎ2 3! ∇2𝑓(𝑥𝑖) + ⋯+ ℎ𝑚−1 𝑚! ∇m−1𝑓(𝑥𝑖)) from this, we can derive the ABM for 𝑚 points 𝑦𝑖+1 = 𝑦𝑖 + ℎ 𝑚! ∑ 𝑏𝑗 𝑚−1 𝑗=0 𝑓(𝑥𝑖−𝑗 , 𝑦𝑖−𝑗) where𝑏𝑗 are the coefficients determined from the polynomial. The coefficients 𝑏𝑗 depend on the order 𝑚 and can be derived from the Newton polynomial, for 𝑚 = 1 𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖) for 𝑚 = 2 𝑦𝑖+1 = 𝑦𝑖 + ℎ 2 (3𝑓(𝑥𝑖 , 𝑦𝑖) − 𝑓(𝑥𝑖−1 , 𝑦𝑖−1)) for 𝑚 = 3 𝑦𝑖+1 = 𝑦𝑖 + ℎ 12 (23𝑓(𝑥𝑖 , 𝑦𝑖) − 16 𝑓(𝑥𝑖−1 , 𝑦𝑖−1) + 5𝑓(𝑥𝑖−2 , 𝑦𝑖−2)) The forth order Adams-Bashforth (𝑚 = 4): 𝑦𝑖+1 = 𝑦𝑖 + ℎ 24 (55 𝑓(𝑥𝑖 , 𝑦𝑖)– 59 𝑓(𝑥𝑖−1 , 𝑦𝑖−1) + 37𝑓(𝑥𝑖−2 , 𝑦𝑖−2) − 9𝑓(𝑥𝑖−3 , 𝑦𝑖−3)) The advantages of Adams-Bashforth methods  Efficiency: Adams-Bashforth methods are efficient for long-term integration since they use multiple previous points to estimate future values.  Explicit nature: Being explicit, these methods do not require solving systems of equations, simplifying their implementation. The disadvantages of Adams-Bashforth methods 27  Not self-starting: Adams-Bashforth methods need a set of initial values to start, often provided by a different method like Runge-Kutta.  Stability: Explicit methods can be less stable, especially for stiff differential equations, compared to implicit methods. 2.4 Adams- Moulton Method (AMM) The Adam-Moulton formulas are implicit multistep methods, to derive the second order Adams-Moulton formula, we first start with a first degree polynomial 𝑝1(𝑥) = 𝛼𝑥 + 𝛽, and we determine the coefficients 𝛼 and 𝛽 using the points (𝑥𝑛 , 𝑦𝑛) and (𝑥𝑛+1 , 𝑦𝑛+1). the polynomial 𝑝1(𝑥) must satisfy: 𝑝1(𝑥𝑛) = 𝑦𝑛and𝑝1(𝑥𝑛+1) = 𝑦𝑛+1 plugging in the values: at 𝑥𝑛, 𝑝1(𝑥𝑛) = 𝛼𝑥𝑛 + 𝛽 = 𝑦𝑛 at 𝑥𝑛+1, 𝑝1(𝑥𝑛+1) = 𝛼(𝑥𝑛 + ℎ ) + 𝛽 = 𝑦𝑛+1 where𝑥𝑛+1 = 𝑥𝑛 + ℎ this gives two equations: 𝛼𝑥𝑛 + 𝛽 = 𝑦𝑛 (2.23) 𝛼(𝑥𝑛 + ℎ ) + 𝛽 = 𝑦𝑛+1 (2.24) subtract equation (2.23) from equation (2.24) 𝛼(𝑥𝑛 + ℎ ) + 𝛽 – (𝛼𝑥𝑛 + 𝛽) = 𝑦𝑛+1 − 𝑦𝑛 simplifying this, we get 𝛼ℎ = 𝑦𝑛+1 − 𝑦𝑛 28 we have 𝛼 = 𝑦𝑛+1 − 𝑦𝑛 ℎ (2.25) substitute𝛼 from (2.25) into equation (2.23) to find 𝛽 𝑦𝑛+1 − 𝑦𝑛 ℎ 𝑥𝑛 + 𝛽 = 𝑦𝑛 rearranging gives 𝛽 = 𝑦𝑛 − 𝑦𝑛+1 − 𝑦𝑛 ℎ the Adams-Moulton method is gives by 𝑦𝑛+1 = 𝑦𝑛 + ℎ 2 (𝑓𝑛 + 𝑓𝑛+1) where𝑓𝑛 = 𝑓(𝑥𝑛 , 𝑦𝑛) and 𝑓𝑛+1 = 𝑓(𝑥𝑛+1 , 𝑦𝑛+1) the fourth-order Adams-Moulton (AM4) formula: 𝑦𝑖+1 = 𝑦𝑖 + ℎ 24 (9 𝑓(𝑥𝑖+1 , 𝑦𝑖+1) + 19 𝑓(𝑥𝑖 , 𝑦𝑖) − 5𝑓(𝑥𝑖−1 , 𝑦𝑖−1) + 𝑓(𝑥𝑖−2 , 𝑦𝑖−2))(2.26) The advantages of Adams-Moulton method  Increased stability: Adams-Moulton methods tend to be more stable than explicit methods, making them suitable for stiff differential equations.  Higher accuracy: These methods can be more accurate, particularly when higher- order variants are used. The disadvantages of Adams-Moulton method  Computational cost: The need to solve implicit equations at each step increases computational complexity and requires iterative solvers.  Not self-starting: Similar to Adams-Bashforth methods, Adams-Moulton methods need initial values from a single-step method for the starting values. 29 2.5 Predictor-Corrector Method (PCM) The Predictor-Corrector method is a class of numerical techniques used to solve ordinary differential equations. These methods combine two steps: a prediction step to estimate the next value and a correction step to refine this estimate. A prediction step is an explicit method that uses information from previous steps to predict the next value, it requires four evaluations from previous steps but only needs one new function evaluation per step, the explicit methods are derived using Newton’s backward-difference polynomial, which provides formulas like the Adams-Bashforth methods. A correction step is an implicit method that refines the predication using both current and previous values, it typically requires two function evaluations per step. To illustrate the derivation of an implicit formula consider the integral form of the solution to the ordinary differential equations 𝑦𝑖+1 = 𝑦𝑖 + ∫ 𝑓(𝑥, 𝑦) 𝑥𝑖+1 𝑥𝑖 𝑑𝑥 to approximate this integral, the trapezoidal rule can be used. The trapezoidal rule approximates the integral as follows ∫ 𝑓(𝑥, 𝑦) 𝑥𝑖+1 𝑥𝑖 𝑑𝑥 ≈ ℎ 2 [ 𝑓(𝑥𝑖 , 𝑦𝑖) + 𝑓(𝑥𝑖+1, 𝑦𝑖+1)] whereℎ = 𝑥𝑖+1 − 𝑥𝑖 is the step size substituting this approximation into the integral equation 𝑦𝑖+1 = 𝑦𝑖 + ℎ 2 [ 𝑓(𝑥𝑖, 𝑦𝑖) + 𝑓(𝑥𝑖+1, 𝑦𝑖+1)] The advantages of Predictor-Corrector method  Improved accuracy: predictor-corrector methods often achieve higher accuracy compared to using either method alone. 30  Flexibility: predictor-corrector methods can be adjusted based on the desired accuracy and stability.  Stability: the corrector step, being implicit, can improve stability, especially for differential equations. The disadvantages of Predictor-Corrector method  Computational cost: corrector step requires solving an implicit equation, which can increase computational complexity and effort.  Not self-starting: like many multistep methods, predictor-corrector methods need initial values provided by a single-step method. 2.6 Error Analysis Error analysis in the context of solving systems of linear differential equations involves understanding how numerical methods approximate solutions and how errors accumulate and affect these solutions. Euler’s method has local truncation error is 𝑂(ℎ2), this means the error made in a single step is proportional to the square of the step size (ℎ), also for the global truncation error is 𝑂(ℎ) because the global truncation error accumulates over 1 ℎ steps, leading to an overall error proportional to the step size (ℎ), for stability Euler’s method is conditionally stable, for stiff system, Euler’s method might exhibit numerical instability unless small step sizes are used. The Runge-Kutta methods has local truncation error is 𝑂(ℎ5), this indicates that the error is a single step is very small for a small (ℎ) for the global truncation error for 𝑂(ℎ4), because the error accumulates over 1 ℎ steps, leading to an overall error proportional to ℎ4. RK4is generally more stable compared to Euler’s method, especially for non-differential systems. It can be applied with larger step size without significant loss of accuracy. The Adams-Bashforth (Predictor) and Adams-Moulton (Corrector) methods are multistep methods that use previous values to predict and correct future values. For these methods, Local Truncation Error decreases with the order of the method. For example, Adams- Bashforth method of order 4 has Local Truncation Error 𝑂(ℎ5),the global truncation error for Adams-Bashforth methods is 𝑂(ℎ𝑝) where 𝑝 is the order of the method. Adams- 31 Moulton method also have a similar order of accuracy but typically require solving nonlinear equations for correction.The multistep methods can be more stable for certain problems but can be challenging to apply directly to differential systems. They generally require a good initial guess or values provided by a single-step method. 2.7 Stability and Convergence Stability: For differential systems, implicit methods (like Adams-Moulton) or implicit Runge-Kutta methods are preferred due to their better stability properties. Convergence: Ensuring that the numerical method converges to the true solution as the step size decreases is crucial. Methods of higher order generally converge faster, but for differential systems, implicit methods or specialized techniques might be required. 32 Chapter Three Numerical Examples and Results In this chapter, we will implement the numerical methods discussed in Chapter Two to solve several numerical examples. This implementation will involve using appropriate algorithms and MATLAB software. We will compare the exact solutions with the approximate solutions obtained through these methods. The results of this comparison will be presented both in tabular form and through graphical illustrations. Example 3.1 Consider the system of linear differential equations 𝑑𝑥 𝑑𝑡 = 2 𝑥 − 𝑦 (3.1) 𝑑𝑦 𝑑𝑡 = 𝑥 subject to the initial conditions 𝑥(0) = 6 (3.2) 𝑦(0) = 2 To find the exact solution of system (3.1), we have 𝑥(𝑡) = 𝑦′(𝑡) now, substitute 𝑥 = 𝑦′ into the equation 𝑥′ = 2 𝑥 − 𝑦 yields 𝑦′′ − 2𝑦′ + 𝑦 = 0 solving this equations gives 𝑦(𝑡) = (𝐶1 + 𝐶2𝑡)𝑒 𝑡 Now find 𝐶1 and 𝐶2 with initial condition gives 𝐶1 = 2 and 𝐶2 = 4 33 the exact solution to the system (3.1) is 𝑥(𝑡) = 6 𝑒𝑡 + 4 𝑡𝑒𝑡 𝑦(𝑡) = 2 𝑒𝑡 + 4 𝑡𝑒𝑡 3.1 The numerical approximation of system (3.1) using the fourth order Rung Kutta method The following algorithm implements the fourth order Rung Kutta method using the matlab software. Algorithm 3.1 In order to approximate the solution of the general initial value problem 𝑥′(𝑡) = 𝑓(𝑡, 𝑥) , 𝑎 ≤ 𝑡 ≤ 𝑏 and 𝑥(𝑎) = 𝑥0 We follow the algorithm 3.1 Set ℎ = 𝑎−𝑏 𝑛 ; 𝑡0 = 𝑎; 𝑘1 = 𝑓(𝑡𝑛, 𝑥𝑛); 𝑘2 = 𝑓(𝑡𝑛 + ℎ 2 , 𝑥𝑛 + ℎ 2 𝑘1); 𝑘3 = 𝑓 (𝑡𝑛 + ℎ 2 , 𝑥𝑛 + ℎ 2 𝑘2) ; 𝑘4 = 𝑓(𝑡𝑛 + ℎ , 𝑥𝑛 + ℎ𝑘3); 𝑥𝑛+1 = 𝑥𝑛 + ℎ 6 (𝑘1 + 2 𝑘2 + 2 𝑘3 + 𝑘4); 𝑡𝑘+1 = 𝑎 + 𝑘 + 1 ℎ ; 34 we repeated the code for𝑘 = 0,1, … , 𝑛 − 1 . After these steps the value 𝑥𝑘 is an approximated value of solution 𝑥(𝑡𝑘) of the differential at 𝑡𝑘. We illustrate the computational of 𝑥1 and 𝑦1with step size ℎ = 0.2 𝑓(𝑡, 𝑥, 𝑦) = 2 𝑥 − 𝑦 , 𝑔(𝑡, 𝑥, 𝑦) = 𝑥, 𝑡0 = 0 , 𝑥0 = 6 and 𝑦0 = 2 𝑘1,1 = ℎ 𝑓(𝑡0 , 𝑥0 , 𝑦0) = 0.2 𝑓(0 ,6 ,2) = 0.2 (10) = 2 𝑘2,1 = ℎ 𝑔(𝑡0 , 𝑥0 , 𝑦0) = 0.2 𝑔(0 ,6 ,2) = 0.2 (6) = 1.2 𝑘1,2 = ℎ 𝑓 (𝑡0 + ℎ 2 , 𝑥0 + 𝑘1,1 2 , 𝑦0 + 𝑘2,1 2 ) = 0.2 𝑓(0.1 , 6 + 1 , 2 + 0.6 ) = 0.2 𝑓(0.1 , 7 , 2.6 ) = 0.2 (11.4) = 2.28 𝑘2,2 = ℎ 𝑔 (𝑡0 + ℎ 2 , 𝑥0 + 𝑘1,1 2 , 𝑦0 + 𝑘2,1 2 ) = 0.2 𝑔(0.1 , 6 + 1 , 2 + 0.6 ) = 0.2 𝑔(0.1 , 7 , 2.6 ) = 0.2 (7) = 1.4 𝑘1,3 = ℎ 𝑓 (𝑡0 + ℎ 2 , 𝑥0 + 𝑘1,2 2 , 𝑦0 + 𝑘2,2 2 ) = 0.2 𝑓(0.1 , 6 + 1.14 , 2 + 0.7 ) = 0.2 𝑓(0.1 , 7.14 , 2.7 ) = 0.2 (11.58) = 2.316 𝑘2,3 = ℎ 𝑔 (𝑡0 + ℎ 2 , 𝑥0 + 𝑘1,2 2 , 𝑦0 + 𝑘2,2 2 ) = 0.2 𝑔(0.1 , 6 + 1.14 , 2 + 0.7 ) = 0.2 𝑔(0.1 , 7.14 , 2.7 ) = 0.2 (7.14) = 1.428 𝑘1,4 = ℎ 𝑓(𝑡0 + ℎ , 𝑥0 + 𝑘1,3 , 𝑦0 + 𝑘2,3) = 0.2 𝑓(0.2 , 6 + 2.316 , 2 + 1.428 ) = 0.2 𝑓(0.2 , 8.316 , 3.428 ) = 0.2 (13.204) = 2.6408 𝑘2,4 = ℎ 𝑔(𝑡0 + ℎ , 𝑥0 + 𝑘1,3 , 𝑦0 + 𝑘2,3) = 0.2 𝑔(0.2 , 6 + 2.316 , 2 + 1.428 ) = 0.2 𝑔(0.2 , 8.316 , 3.428 ) = 0.2 (8.316) = 1.6632 𝑥1 = 𝑥0 + 1 6 ( 𝑘1,1 + 2 𝑘1,2 + 2 𝑘1,3 + 𝑘1,4) 35 = 6 + 1 6 (2 + 2 (2.28) + 2 (2.316) + 2.6408) = 6 + 1 6 (13.8328) = 8.3054 𝑦1 = 𝑦0 + 1 6 ( 𝑘2,1 + 2 𝑘2,2 + 2 𝑘2,3 + 𝑘2,4) = 2 + 1 6 (1.2 + 2 (1.4) + 2 (1.428) + 1.6632) = 2 + 1 6 (8.5192) = 3.4198 These numbers give us the approximation 𝑥1 ≈ 𝑥(0.2) = 8.3054 and 𝑦1 ≈ 𝑦(0.2) = 3.4198 Now to find 𝑥(0.4) and 𝑦(0.4) we need to find 𝑡1 = 𝑡0 + ℎ = 0 + 0.2 = 0.2 For 𝑡1 = 0.2 we need to find 𝑘1,1 = ℎ 𝑓(𝑡1 , 𝑥1 , 𝑦1) = 0.2 𝑓(0.2 , 8.3054 , 3.4198) = 0.2 (13.191) = 2.6382 𝑘2,1 = ℎ 𝑔(𝑡1 , 𝑥1 , 𝑦1) = 0.2 𝑔(0.2 , 8.3054 , 3.4198) = 0.2 (8.3054) = 1.661 𝑘1,2 = ℎ 𝑓 (𝑡1 + ℎ 2 , 𝑥1 + 𝑘1,1 2 , 𝑦1 + 𝑘2,1 2 ) = 0.2 𝑓 (0.2 + 0.1 , 8.3054 + 2.6382 2 , 3.4198 + 1.661 2 ) = 0.2 𝑓(0.3 , 9.6245 , 4.2503 ) = 0.2 (14.9987) = 2.9997 𝑘2,2 = ℎ 𝑔 (𝑡1 + ℎ 2 , 𝑥1 + 𝑘1,1 2 , 𝑦1 + 𝑘2,1 2 ) = 0.2 𝑔 (0.2 + 0.1 , 8.3054 + 2.6382 2 , 3.4198 + 1.661 2 ) = 0.2 𝑔(0.3 , 9.6245 , 4.2503 ) = 0.2 (9.6245) = 1.9249 𝑘1,3 = ℎ 𝑓 (𝑡1 + ℎ 2 , 𝑥1 + 𝑘1,2 2 , 𝑦1 + 𝑘2,2 2 ) = 0.2 𝑓 (0.2 + 0.1 , 8.3054 + 2.9997 2 , 3.4198 + 1.9249 2 ) 36 = 0.2 𝑓(0.3 , 9.8052 , 4.3822 ) = 0.2 (15.2282) = 3.0456 𝑘2,3 = ℎ 𝑔 (𝑡1 + ℎ 2 , 𝑥1 + 𝑘1,2 2 , 𝑦1 + 𝑘2,2 2 ) = 0.2 𝑔 (0.2 + 0.1 , 8.3054 + 2.9997 2 , 3.4198 + 1.9249 2 ) = 0.2 𝑔(0.3 , 9.8052 , 4.3822 ) = 0.2 (9.8052) = 1.961 𝑘1,4 = ℎ 𝑓(𝑡1 + ℎ , 𝑥1 + 𝑘1,3 , 𝑦1 + 𝑘2,3) = 0.2 𝑓(0.2 + 0.2 , 8.3054 + 3.0456 , 3.4198 + 1.961 ) = 0.2 𝑓(0.4 , 11.351 , 5.3808 ) = 0.2 (17.3212) = 3.4642 𝑘2,4 = ℎ 𝑔(𝑡1 + ℎ , 𝑥1 + 𝑘1,3 , 𝑦1 + 𝑘2,3) = 0.2 𝑔(0.2 + 0.2 , 8.3054 + 3.0456 , 3.4198 + 1.961 ) = 0.2 𝑔(0.4 , 11.351 , 5.3808 ) = 0.2 (11.351) = 2.2702 𝑥2 = 𝑥1 + 1 6 ( 𝑘1,1 + 2 𝑘1,2 + 2 𝑘1,3 + 𝑘1,4) = 8.3054 + 1 6 (2.6382 + 2 (2.9997) + 2(3.0456) + 3.4642) = 8.3054 + 1 6 (18.193) = 11.3376 𝑦2 = 𝑦1 + 1 6 ( 𝑘2,1 + 2 𝑘2,2 + 2 𝑘2,3 + 𝑘2,4) = 3.4198 + 1 6 (1.661 + 2 (1.9249) + 2 (1.961) + 2.2702) = 3.4198 + 1 6 (11.703) = 5.3704 These numbers give us the approximation 𝑥2 ≈ 𝑥(0.4) = 11.3376 and 𝑦2 ≈ 𝑦(0.4) = 5.3704 Now to find 𝑥(0.6) and 𝑦(0.6) we need to find 𝑡2 = 𝑡1 + ℎ = 0.2 + 0.2 = 0.4 For 𝑡2 = 0.4 we need to find 37 𝑘1,1 = ℎ 𝑓(𝑡2 , 𝑥2 , 𝑦2) = 0.2 𝑓(0.4 , 11.3376 , 5.3704) = 0.2 (17.3047) = 3.4609 𝑘2,1 = ℎ 𝑔(𝑡2 , 𝑥2 , 𝑦2) = 0.2 𝑔(0.4 , 11.3376 , 5.3704) = 0.2 (11.3375) = 2.2675 𝑘1,2 = ℎ 𝑓 (𝑡2 + ℎ 2 , 𝑥2 + 𝑘1,1 2 , 𝑦2 + 𝑘2,1 2 ) = 0.2 𝑓 (0.4 + 0.1 , 11.3376 + 3.4609 2 , 5.3704 + 2.2675 2 ) = 0.2 𝑓(0.5 , 13.0679 , 6.504 ) = 0.2 (19.6318) = 3.9264 𝑘2,2 = ℎ 𝑔 (𝑡2 + ℎ 2 , 𝑥2 + 𝑘1,1 2 , 𝑦2 + 𝑘2,1 2 ) = 0.2 𝑔 (0.4 + 0.1 , 11.3376 + 3.4609 2 , 5.3704 + 2.2675 2 ) = 0.2 𝑔(0.5 , 13.0679 , 6.504 ) = 0.2 (13.0679) = 2.6136 𝑘1,3 = ℎ 𝑓 (𝑡2 + ℎ 2 , 𝑥2 + 𝑘1,2 2 , 𝑦2 + 𝑘2,2 2 ) = 0.2 𝑓 (0.4 + 0.1 , 11.3376 + 3.9264 2 , 5.3704 + 2.6136 2 ) = 0.2 𝑓(0.5 , 13.3006 , 6.677 ) = 0.2 (19.9242) = 3.9849 𝑘2,3 = ℎ 𝑔 (𝑡2 + ℎ 2 , 𝑥2 + 𝑘1,2 2 , 𝑦2 + 𝑘2,2 2 ) = 0.2 𝑔 (0.4 + 0.1 , 11.3376 + 3.9264 2 , 5.3704 + 2.6136 2 ) = 0.2 𝑔(0.5 , 13.3006 , 6.677 ) = 0.2 (13.3006) = 2.6601 𝑘1,4 = ℎ 𝑓(𝑡2 + ℎ , 𝑥2 + 𝑘1,3 , 𝑦2 + 𝑘2,3) = 0.2 𝑓(0.4 + 0.2 , 11.3376 + 3.9849 , 5.3704 + 2.6601 ) = 0.2 𝑓(0.6 , 15.3223 , 8.0304 ) = 0.2 (22.6142) = 4.5229 𝑘2,4 = ℎ 𝑔(𝑡2 + ℎ , 𝑥2 + 𝑘1,3 , 𝑦2 + 𝑘2,3) = 0.2 𝑔(0.4 + 0.2 , 11.3376 + 3.9849 , 5.3704 + 2.6601 ) 38 = 0.2 𝑔(0.6 , 15.3223 , 8.0304 ) = 0.2 (15.3223) = 3.0645 𝑥3 = 𝑥2 + 1 6 ( 𝑘1,1 + 2 𝑘1,2 + 2 𝑘1,3 + 𝑘1,4) = 11.3376 + 1 6 (3.4609 + 2 (3.9264) + 2 (3.9849) + 4.5229) = 11.3376 + 1 6 (23.8059) = 15.3054 𝑦3 = 𝑦2 + 1 6 ( 𝑘2,1 + 2 𝑘2,2 + 2 𝑘2,3 + 𝑘2,4) = 5.3704 + 1 6 (2.2675 + 2 (2.6136) + 2 (2.6601) + 3.0645) = 5.3704 + 1 6 (15.8791) = 8.017 These numbers give us the approximation 𝑥3 ≈ 𝑥(0.6) = 15.3054 and 𝑦3 ≈ 𝑦(0.6) = 8.017 Now to find 𝑥(0.8) and 𝑦(0.8) we need to find 𝑡3 = 𝑡2 + ℎ = 0.4 + 0.2 = 0.6 For 𝑡3 = 0.6 we need to find 𝑘1,1 = ℎ 𝑓(𝑡3 , 𝑥3 , 𝑦3) = 0.2 𝑓(0.6 , 15.3054 , 8.017) = 0.2 (22.5934) = 4.5187 𝑘2,1 = ℎ 𝑔(𝑡3 , 𝑥3 , 𝑦3) = 0.2 𝑔(0.6 , 15.3054 , 8.017) = 0.2 (15.3051) = 3.061 𝑘1,2 = ℎ 𝑓 (𝑡3 + ℎ 2 , 𝑥3 + 𝑘1,1 2 , 𝑦3 + 𝑘2,1 2 ) = 0.2 𝑓 (0.6 + 0.1 , 15.3051 + 4.5187 2 , 8.0168 + 3.061 2 ) = 0.2 𝑓(0.7 , 17.5644 , 9.5473 ) = 0.2 (25.5815) = 5.1164 𝑘2,2 = ℎ 𝑔 (𝑡3 + ℎ 2 , 𝑥3 + 𝑘1,1 2 , 𝑦3 + 𝑘2,1 2 ) = 0.2 𝑔 (0.6 + 0.1 , 15.3054 + 4.5187 2 , 8.017 + 3.061 2 ) = 0.2 𝑔(0.7 , 17.5644 , 9.5473 ) = 0.2 (17.5644) = 3.5129 39 𝑘1,3 = ℎ 𝑓 (𝑡3 + ℎ 2 , 𝑥3 + 𝑘1,2 2 , 𝑦3 + 𝑘2,2 2 ) = 0.2 𝑓 (0.6 + 0.1 , 15.3054 + 5.1164 2 , 8.017 + 3.5129 2 ) = 0.2 𝑓(0.7 , 17.8632 , 9.7732 ) = 0.2 (25.9532) = 5.1907 𝑘2,3 = ℎ 𝑔 (𝑡3 + ℎ 2 , 𝑥3 + 𝑘1,2 2 , 𝑦3 + 𝑘2,2 2 ) = 0.2 𝑔 (0.6 + 0.1 , 15.3054 + 5.1164 2 , 8.017 + 3.5129 2 ) = 0.2 𝑔(0.7 , 17.8632 , 9.7732 ) = 0.2 (17.8632) = 3.5727 𝑘1,4 = ℎ 𝑓(𝑡3 + ℎ , 𝑥3 + 𝑘1,3 , 𝑦3 + 𝑘2,3) = 0.2 𝑓(0.6 + 0.2 , 15.3054 + 5.1907 , 8.017 + 3.5727 ) = 0.2 𝑓(0.8 , 20.4957 , 11.5894 ) = 0.2 (29.402) = 5.8805 𝑘2,4 = ℎ 𝑔(𝑡3 + ℎ , 𝑥3 + 𝑘1,3 , 𝑦3 + 𝑘2,3) = 0.2 𝑔(0.6 + 0.2 , 15.3054 + 5.1907 , 8.017 + 3.5727 ) = 0.2 𝑔(0.8 , 20.4957 , 11.5894 ) = 0.2 (20.4957) = 4.0992 𝑥4 = 𝑥3 + 1 6 ( 𝑘1,1 + 2 𝑘1,2 + 2 𝑘1,3 + 𝑘1,4) = 15.3051 + 1 6 (4.5187 + 2(5.1164) + 2(5.1907) + 5.8805) = 15.3051 + 1 6 (31.0128) = 20.4744 𝑦4 = 𝑦3 + 1 6 ( 𝑘2,1 + 2 𝑘2,2 + 2 𝑘2,3 + 𝑘2,4) = 8.017 + 1 6 (3.061 + 2(3.5129) + 2(3.5727) + 4.0992) = 8.0168 + 1 6 (21.3309) = 11.5723 These numbers give us the approximation 𝑥4 ≈ 𝑥(0.8) = 20.4744 and 40 𝑦4 ≈ 𝑦(0.8) = 11.5723. Table 3.1(a) and 3.2(b), contain the 𝑘’s value when applying algorithm 3.1 for system 3.1 Table 3.1(a) The 𝑘’s values for 𝑥(𝑡) 𝑡 𝑘1 𝑘2 𝑘3 𝑘4 0 0.000000000 0.000000000 0.000000000 0.000000000 0.2 2.000000000 2.280000000 2.316000000 2.640800000 0.4 2.638231111 2.999764444 3.045623200 3.464241231 0.6 3.460913430 3.926414261 3.984932232 4.522923222 0.8 4.518724221 5.116422632 5.190741111 5.880512324 Table 3.1(b) The 𝑘’s values for 𝑦(𝑡) 𝑡 𝑘1 𝑘2 𝑘3 𝑘4 0 0.000000000 0.000000000 0.000000000 0.000000000 0.2 1.200000000 1.400000000 1.428000000 1.663200000 0.4 1.661011111 1.924936444 1.961641000 2.270222221 0.6 2.267512262 2.613614661 2.660116141 3.064526212 0.8 3.061613361 3.512941212 3.572711120 4.099261112 Table 3.2 (a), the exact and numerical solution for 𝑥(𝑡) when applying RKM on system (3.1) and showing the resulting error. Table 3.2(a) The exact and numerical solution of 𝑥(𝑡) by applying RKM 𝑡 Exact solution Approximate solution Absolute error 0 6.000000000 6.000000000 0.000000000 0.2 8.305511222 8.305444440 0.000072095 0.4 11.33784220 11.33761423 0.00018149 0.6 15.30571212 15.30542226 0.00034238 0.8 20.47492426 20.47446211 0.00057366 From table 3.2(a) we conclude the maximum error = 0.00057366 . Figure 3.1(a): Compares the exact solution 𝑥(𝑡) and approximate solution with time 𝑡. 41 Figure 3.1(a) The exact and numerical solutions of 𝑥(𝑡) Figure 3.1(b) the absolute error resulting of applying RKM of 𝑥(𝑡) 42 Figure 3.1(b): Shows the absolute error resulting of𝑥(𝑡). Table 3.2 (b), the exact and numerical solution for 𝑦(𝑡) when applying RKM on system (3.1) and showing the resulting error. Table 3.2(b) The exact and numerical solution of 𝑦(𝑡) by applying RKM t Exact solution Approximate solution Absolute error 0 2.00000000 2.000000000 0.000000000 0.2 3.419927723 3.419866666 0.000061057 0.4 5.370568912 5.370414373 0.000154539 0.6 8.017322722 8.017029717 0.000293005 0.8 11.57281283 11.57231957 0.00049326 From table 3.2(b) we conclude the maximum error = 0.00049326. Figure 3.2(a): Compares the exact solution 𝑦(𝑡) and approximate solution. Figure 3.2(a) The exact and numerical solutions of 𝑦(𝑡) 43 Figure 3.2(b) absolute error resulting of applying RKM of 𝑦(𝑡) Figure 3.2(b): Shows the absolute error resulting of 𝑦(𝑡).. 3.2 The numerical approximation of system (3.1) using the Adams-Bashforth Method The following algorithm implements the Adams-Bashforth Method using the matlab software. Algorithm 3.2 An algorithm for the Fourth-order Adams Bashforth method is given below Set ℎ = 𝑎−𝑏 𝑛 ; 𝑡0 = 𝑎; 𝑥𝑛 ′ = 𝑓(𝑡𝑛, 𝑥𝑛); 𝑥𝑛−1 ′ = 𝑓(𝑡𝑛−1, 𝑥𝑛−1); 𝑥𝑛−2 ′ = 𝑓(𝑡𝑛−2, 𝑥𝑛−2); 𝑥𝑛−3 ′ = 𝑓(𝑡𝑛−3, 𝑥𝑛−3); 44 𝑥𝑛+1 ∗ = 𝑥𝑛 + ℎ 24 ( 55𝑥𝑛 ′ − 59𝑥𝑛−1 ′ + 37𝑥𝑛−2 ′ − 9𝑥𝑛−3 ′ ); 𝑡𝑘+1 = 𝑎 + 𝑘 + 1 ℎ ; we repeated the code for𝑘 = 0, 1, … , 𝑛 − 1 . After these steps the value 𝑥𝑘 ∗ is an approximated value of solution 𝑥(𝑡𝑘) of the differential at 𝑡𝑘. The Adams-Bashforth Method to approximate 𝑥(0.8) and 𝑦(0.8) with ℎ = 0.2 in the system (3.1). 𝑥0 ′ = 𝑓(𝑡0 , 𝑥0 , 𝑦0) = 𝑓( 0 , 6, 2) = 10 𝑦0 ′ = 𝑔(𝑡0 , 𝑥0 , 𝑦0) = 𝑔( 0 , 6 , 2) = 6 𝑥1 ′ = 𝑓(𝑡1 , 𝑥1 , 𝑦1) = 𝑓( 0.2 , 8.30546666, 3.419866666) = 13.19106665 𝑦1 ′ = 𝑔(𝑡1 , 𝑥1 , 𝑦1) = 𝑔( 0.2 , 8.30546666, 3.419866666) = 8.30546666 𝑥2 ′ = 𝑓(𝑡2 , 𝑥2 , 𝑦2) = 𝑓( 0.4 , 11.33768621, 5.370414373) = 17.30495805 𝑦2 ′ = 𝑔(𝑡2 , 𝑥2 , 𝑦2) = 𝑔( 0.4 , 11.33768621, 5.370414373) = 11.33768621 𝑥3 ′ = 𝑓(𝑡3 , 𝑥3 , 𝑦3) = 𝑓( 0.6 , 15.30545554 , 8.017029717) = 22.59388136 𝑦3 ′ = 𝑔(𝑡3 , 𝑥3 , 𝑦3) = 𝑔( 0.6 , 15.30545554 , 8.017029717) = 15.30545554 𝑥4 ∗ = 𝑥3 + ℎ 24 ( 55𝑥3 ′ − 59𝑥2 ′ + 37𝑥1 ′ − 9𝑥0 ′ ) = 15.30545554 + 0.2 24 ( 55(22.59388136) − 59(17.30495805) + 37(13.19106665 ) − 9(10)) = 15.30545554 + 0.2 24 (619.7404164) = 15.30545554 + 5.16450347 = 20.46995901 45 𝑦4 ∗ = 𝑦3 + ℎ 24 ( 55𝑦3 ′ − 59𝑦2 ′ + 37𝑦1 ′ − 9𝑦0 ′) = 8.017029717 + 0.2 24 (55(15.30545554) − 59(11.33768621) + 37(8.30546666) − 9(6)) = 8.017029717 + 0.2 24 (426.1788347) = 8.017029717 + 3.551490289 = 11.56852001 Table 3.3(a), shows the exact and numerical solutions when using ABM on system (3.1) and showing the resulting error. Table 3.3(a) The exact and numerical solution of applying Fourth Adams-Bashforth method at 𝑡 = 0.8 Exact solution Approximate solution Absolute error 𝑥(𝑡) 20.47492426 20.46995901 0.00501753 𝑦(𝑡) 11.57281283 11.56852001 0.00429282 These results show the ABM is accuracy with exact solution. 3.3 The numerical approximation of system (3.1) using the Adams-Moulton method The following algorithm implements the Adams-Moulton methodusing the matlab software. Algorithm 3.3 An algorithm for the Fourth-order Adams-Moulton method is given below Set ℎ = 𝑎−𝑏 𝑛 ; 𝑥𝑛+1 ′ = 𝑓(𝑡𝑛+1, 𝑥𝑛+1 ∗ ); 𝑥𝑛+1 = 𝑥𝑛 + ℎ 24 ( 9𝑥𝑛+1 ′ + 19𝑥𝑛 ′ − 5𝑥𝑛−1 ′ + 𝑥𝑛−2 ′ ); 𝑡𝑘+1 = 𝑎 + 𝑘 + 1 ℎ ; we repeated the code for𝑘 = 0, 1, … , 𝑛 − 1 . 46 After these steps the value 𝑥𝑘 is an approximated value of solution 𝑥(𝑡𝑘) of the differential at 𝑡𝑘. Consider the RK4 method with 𝑥3 = 15.30545554, 𝑦3 = 8.017029717 and 𝑡4 = 0.8at system (3.1) we have 𝑥4 ′ = 𝑓(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗) = 𝑓(0.8 , 20.46995901 , 11.56852001) = 29.37139801 𝑦4 ′ = 𝑔(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗) = 𝑔(0.8 , 20.46995901 , 11.56852001) = 20.46995901 now 𝑦𝑛+1 = 𝑦𝑛 + ℎ 24 ( 9𝑓𝑛+1 + 19𝑓𝑛 − 5𝑓𝑛−1 + 𝑓𝑛−2) 𝑥4 = 𝑥3 + ℎ 24 ( 9𝑥4 ′ + 19𝑥3 ′ − 5𝑥2 ′ + 𝑥1 ′) = 15.30545554 + 0.2 24 ( 9(29.37139801) + 19(22.59388136) − 5(17.30495805) + (13.19106665)) = 15.30545554 + 0.2 24 (620.2926044) = 15.30545554 + 5.169105037 𝑥4 = 20.47456058 𝑦4 = 𝑦3 + ℎ 24 ( 9𝑦4 ′ + 19𝑦3 ′ − 5𝑦2 ′ + 𝑦1 ′) = 8.017029717 + 0.2 24 (9(20.46995901) + 19(15.30545554) − 5(11.33768621) + (8.30546666)) = 8.017029717 + 0.2 24 (426.650322) = 8.017029717 + 3.55541935 𝑦4 = 11.57244907 47 Table 3.3(b), the exact and numerical solution when applying AMM on system (3.1) and showing the resulting error. Table 3.3(b) The exact and numerical solution of applying Fourth Adams-Moulton method at 𝑡 = 0.8 Exact solution Approximate solution Absolute error 𝑥(𝑡) 20.47492426 20.47456058 0.00041596 𝑦(𝑡) 11.57281283 11.57244907 0.00036376 These results show the AMM is accuracy with exact solution. 3.4 The numerical approximation of system (3.1) using the Predictor-Corrector method The following algorithm implements the Predictor-Corrector method using the matlab software. Algorithm 3.4 Set ℎ = 𝑎−𝑏 𝑛 ; 𝑥𝑛+1 ^ = 𝑥𝑛 + ℎ 𝑓(𝑡𝑛 , 𝑥𝑛); 𝑥𝑛+1 = 𝑥𝑛 + ℎ 2 ( 𝑓(𝑡𝑛 , 𝑥𝑛) + 𝑓(𝑡𝑛+1 , 𝑥𝑛+1 ^ )) ; 𝑡𝑘+1 = 𝑎 + 𝑘 + 1 ℎ ; we repeated the code for𝑘 = 0,1, … , 𝑛 − 1. After these steps the value 𝑥𝑘 is an approximated value of solution 𝑥(𝑡𝑘) of the differential at 𝑡𝑘. Consider the RK4 method with 𝑥3 = 15.30545554, 𝑦3 = 8.017029717 and 𝑡4 = 0.8 at system (3.1) we have 𝑥𝑛+1 ^ = 𝑥𝑛 + ℎ 𝑓(𝑡𝑛 , 𝑥𝑛 , 𝑦𝑛) 𝑥𝑛+1 = 𝑥𝑛 + ℎ 2 ( 𝑓(𝑡𝑛 , 𝑥𝑛 , 𝑦𝑛) + 𝑓(𝑡𝑛+1 , 𝑥𝑛+1 ^ , 𝑦𝑛+1 ^ )) 48 now 𝑥4 ^ = 𝑥3 + ℎ 𝑓(𝑡3 , 𝑥3 , 𝑦3) = 15.30545554 + 0.2 𝑓(0.6 , 15.30545554 , 8.017029717) = 15.30545554 + 0.2 (22.59388136) = 15.30545554 + 4.518776272 = 19.82423181 𝑦4 ^ = 𝑦3 + ℎ 𝑔(𝑡3 , 𝑥3 , 𝑦3) = 8.017029717 + 0.2 𝑔(0.6 , 15.30545554 , 8.017029717) = 8.017029717 + 0.2 (15.30545554) = 8.017029717 + 3.061091108 = 11.07812083 𝑥4 = 𝑥3 + ℎ 2 ( 𝑓(𝑡3 , 𝑥3 , 𝑦3) + 𝑓(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) = 15.30545554 + 0.2 2 (𝑓(0.6 , 15.30545554 , 8.017029717) + 𝑓(0.8 , 19.82423181 , 11.07812083 )) = 15.30545554 + 0.2 2 ( 22.59388136 + 28.57034279) = 15.30545554 + 0.1( 51.16422415) = 20.42187796 𝑦4 = 𝑦3 + ℎ 2 ( 𝑔(𝑡3 , 𝑥3 , 𝑦3) + 𝑔(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) = 8.017029717 + 0.2 2 (𝑔(0.6 , 15.30545554 , 8.017029717) + 𝑔(0.8 , 19.82423181 , 11.07812083 )) = 8.017029717 + 0.1(15.30545554 + 19.82423181) = 8.017029717 + 0.1( 35.12968735) = 11.52999845 49 Table 3.3(c), the exact and numerical solution when applying PCM on system (3.1) and showing the resulting error. Table3.3(c) The exact and numerical solution of applying Predictor-Corrector method at 𝑡 = 0.8 Exact solution Approximate solution Absolute error 𝑥(𝑡) 20.47492426 20.42187796 0.05309858 𝑦(𝑡) 11.57281283 11.52999845 0.04281438 these results show the PCM is accuracy with exact solution. Example 3.2 Consider the system of linear differential equations 𝑑𝑥 𝑑𝑡 = 3𝑥 − 𝑦 − 𝑧 (3.3) 𝑑𝑦 𝑑𝑡 = 𝑥 + 𝑦 − 𝑧 𝑑𝑧 𝑑𝑡 = 𝑥 − 𝑦 + 𝑧 subject to the initial conditions 𝑥(0) = 6 (3.4) 𝑦(0) = 3 𝑧(0) = 4 to find the exact solution of System (3.3),we have det(𝐴 − 𝜆𝐼) = (1 − 𝜆)(𝜆 − 2)2 = 0. For 𝜆1 = 1 we obtain 𝐾1 = ( 1 1 1 ) for 𝜆2 = 2 we obtain 𝐾2 = ( 1 0 1 ) and 𝐾3 = ( 1 1 0 ) 50 then 𝑋(𝑡) = 𝑐1 ( 1 1 1 )𝑒𝑡 + 𝑐2 ( 1 0 1 ) 𝑒2𝑡 + 𝑐3 ( 1 1 0 ) 𝑒2𝑡. The exact solution of system (3.3) is 𝑥(𝑡) = 𝑒𝑡 + 5 𝑒2𝑡 𝑦(𝑡) = 𝑒𝑡 + 2 𝑒2𝑡 𝑧(𝑡) = 𝑒𝑡 + 3 𝑒2𝑡 3.5 The numerical approximation of system (3.3) using the fourth order Rung Kutta methods Table 3.4(a), 3.4(b) and 3.4(c), contain the 𝑘’s value when applying RKM for system 3.3 Table 3.4(a) The 𝑘’s value for 𝑥(𝑡) 𝑡 𝑘1 𝑘2 𝑘3 𝑘4 0 0.000000000 0.000000000 0.000000000 0.000000000 0.2 2.200000000 2.620000000 2.702000000 3.236400000 0.4 3.227746666 3.848868000 3.970649466 4.761776293 0.6 4.748900267 5.668843961 5.849849064 7.022603176 0.8 7.003435201 8.367680112 8.636884881 10.37728762 Table 3.4(b) The 𝑘’s value for 𝑦(𝑡) 𝑡 𝑘1 𝑘2 𝑘3 𝑘4 0 0.000000000 0.000000000 0.000000000 0.000000000 0.2 1.000000000 1.180000000 1.214000000 1.441200000 0.4 1.437666666 1.700772000 1.750950266 2.083816613 0.6 2.078578262 2.464457555 2.538649778 3.027801456 0.8 3.020026855 3.587590097 3.697458532 4.418108741 Table 3.4(c) The 𝑘’s value for 𝑧(𝑡) 𝑡 𝑘1 𝑘2 𝑘3 𝑘4 0 0.000000000 0.000000000 0.000000000 0.000000000 0.2 2.200000000 2.620000000 2.702000000 3.236400000 0.4 3.227746666 3.848868000 3.970649466 4.761776293 0.6 4.748900267 5.668843961 5.849849064 7.022603176 0.8 7.003435201 8.367680112 8.636884881 10.37728762 51 Table 3.5(a), -Appendix D- the exact and numerical solution for 𝑥(𝑡) when applying RKM on system (3.3) and showing the resulting error. From table 3.5(a) we conclude the maximum error = 0.006086350 . Figure 3.3(a): Compares the exact solution 𝑥(𝑡) and approximate solution with time 𝑡 . Figure 3.3(a) The exact and numerical solutions of 𝑥(𝑡) Figure 3.3(b) the absolute error resulting of applying RKM. 52 Figure 3.3(b): Shows the absolute error resulting of system 3.3. Table 3.5(b), -Appendix D- the exact and numerical solution for 𝑦(𝑡) when applying RKM on system (3.3) and showing the resulting error. From table 3.5(b) we conclude the maximum error = 0.002446607 Figure 3.4(a): Compares the exact solution 𝑦(𝑡) and approximate solution with time 𝑡 . Figure 3.4(a) The exact and numerical solutions 𝑦(𝑡) 53 Figure 3.4(b) the absolute error resulting of applying RKM of 𝑦(𝑡) Figure 3.4(b): Shows the absolute error resulting of 𝑦(𝑡). Table 3.5(c), -Appendix D- the exact and numerical solution for 𝑧(𝑡) when applying RKM on system (3.3) and showing the resulting error. From table 3.5(c) we conclude the maximum error = 0.003659851 Figure 3.5(a): -Appendix E- Compares the exact solution 𝑧(𝑡) and approximate solution with time 𝑡 . Figure 3.5(b): -Appendix E- Shows the absolute error resulting of 𝑧(𝑡). 3.6 The numerical approximation of system (3.3) using the Fourth Adams-Bashforth Method The Adams-Bashforth Method to approximate 𝑥(0.8), 𝑦(0.8) and 𝑧(0.8)with ℎ = 0.2 in the system (3.3). 𝑥0 ′ = 𝑓(𝑡0 , 𝑥0 , 𝑦0 , 𝑧0) = 𝑓( 0 , 6 , 3 , 4 ) = 11 𝑦0 ′ = 𝑔(𝑡0 , 𝑥0 , 𝑦0 , 𝑧0) = 𝑔( 0 , 6 , 3 , 4 ) = 5 54 𝑧0 ′ = 𝑞(𝑡0 , 𝑥0 , 𝑦0 , 𝑧0) = 𝑞( 0 , 6 , 3 , 4 ) = 7 𝑥1 ′ = 𝑓(𝑡1 , 𝑥1 , 𝑦1 , 𝑧1) = 𝑓( 0.2 , 8.680066666, 4.504866666, 5.6966) = 16.13873331 𝑦1 ′ = 𝑔(𝑡1 , 𝑥1 , 𝑦1 , 𝑧1) = 𝑔(0.2 , 8.680066666, 4.504866666, 5.6966) = 7.188333326 𝑧1 ′ = 𝑞(𝑡1 , 𝑥1 , 𝑦1 , 𝑧1) = 𝑞(0.2 , 8.680066666, 4.504866666, 5.6966) = 10.17179999 𝑥2 ′ = 𝑓(𝑡2 , 𝑥2 , 𝑦2 , 𝑧2) = 𝑓( 0.4 , 12.61815964, 5.942354635, 8.167622973) = 2374450131 𝑦2 ′ = 𝑔(𝑡2 , 𝑥2 , 𝑦2 , 𝑧2) = 𝑔( 0.4 , 12.61815964, 5.942354635, 8.167622973) = 10.3928913 𝑧2 ′ = 𝑞(𝑡2 , 𝑥2 , 𝑦2 , 𝑧2) = 𝑞( 0.4 , 12.61815964, 5.942354635, 8.167622973) = 14.84342798 𝑥3 ′ = 𝑓(𝑡3 , 𝑥3 , 𝑦3 , 𝑧3) = 𝑓( 0.6, 18.41964123, 8.461120366, 11.78062732) = 35.017176 𝑦3 ′ = 𝑔(𝑡3 , 𝑥3 , 𝑦3 , 𝑧3) = 𝑔( 0.6, 18.41964123, 8.461120366, 11.78062732) = 15.10013428 𝑧3 ′ = 𝑞(𝑡3 , 𝑥3 , 𝑦3 , 𝑧3) = 𝑞( 0.6, 18.41964123, 8.461120366, 11.78062732) = 21.73914818 𝑥4 ∗ = 𝑥3 + ℎ 24 ( 55𝑥3 ′ − 59𝑥2 ′ + 37𝑥1 ′ − 9𝑥0 ′ ) = 18.41964123 + 0.2 24 (55(35.017176 ) − 59(23.74450131) + 37(16.13873331 ) − 9(11)) 55 = 18.41964123 + 0.2 24 (1023.152235) = 18.41964123 + 8.526268625 = 26.94590986 𝑦4 ∗ = 𝑦3 + ℎ 24 ( 55𝑦3 ′ − 59𝑦2 ′ + 37𝑦1 ′ − 9𝑦0 ′) = 8.461120366 + 0.2 24 (55(15.10013428) − 59(10.3928913) + 37(7.188333326) − 9 (5)) = 8.461120366 + 0.2 24 (438.2951314) = 8.461120366 + 3.652459428 = 12.11357979 𝑧4 ∗ = 𝑧3 + ℎ 24 ( 55𝑧3 ′ − 59𝑧2 ′ + 37𝑧1 ′ − 9𝑧0 ′ ) = 11.78062732 + 0.2 24 (55(21.73914818) − 59(14.84342798) + 37(10.17179999) − 9(7)) = 11.78062732 + 0.2 24 (633.2474992) = 11.78062732 + 5.277062493 = 17.05768981 Table 3.6(a), -Appendix D- the exact and numerical solution when applying ABM on system (3.3) and showing the resulting error. These results show the ABM is accuracy with exact solution. 3.7 The numerical approximation of system (3.3) using the Adams-Moulton method Consider the RK4 method with 𝑥3 = 18.41964123 , 𝑦3 = 8.461120366 , 𝑧3 = 11.78062732and 𝑡4 = 0.8 at system (3.3) we have 𝑥4 ′ = 𝑓(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗ , 𝑧4 ∗) = 𝑓(0.8 , 26.94590986 , 12.11357979 , 17.05768981) = 51.66645998 56 𝑦4 ′ = 𝑔(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗ , 𝑧4 ∗) = 𝑔(0.8 , 26.94590986 , 12.11357979 , 17.05768981) = 22.00179984 𝑧4 ′ = 𝑞(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗ , 𝑧4 ∗) = 𝑞(0.8 , 26.94590986 , 12.11357979 , 17.05768981) = 31.89001988 now 𝑥𝑛+1 = 𝑥𝑛 + ℎ 24 ( 9𝑓𝑛+1 + 19𝑓𝑛 − 5𝑓𝑛−1 + 𝑓𝑛−2) 𝑥4 = 𝑥3 + ℎ 24 ( 9𝑥4 ′ + 19𝑥3 ′ − 5𝑥2 ′ + 𝑥1 ′) 𝑥4 = 18.41964123 + 0.2 24 (9(51.66645998) + 19(35.017176) − 5(23.74450131) + (16.13873331)) = 18.41964123 + 0.2 24 (1027.74071) = 18.41964123 + 8.564505917 𝑥4 = 26.98414715 𝑦4 = 𝑦3 + ℎ 24 ( 9𝑦4 ′ + 19𝑦3 ′ − 5𝑦2 ′ + 𝑦1 ′) 𝑦4 = 8.461120366 + 0.2 24 (9(22.00179984) + 19(15.10013428) − 5(10.3928913) + (7.188333326)) = 8.461120366 + 0.2 24 (440.142627) = 8.461120366 + 3.667855225 𝑦4 = 12.12897559 𝑧4 = 𝑧3 + ℎ 24 ( 9𝑧4 ′ + 19𝑧3 ′ − 5𝑧2 ′ + 𝑧1 ′) 57 𝑧4 = 11.78062732 + 0.2 24 (9(31.89001988) + 19(21.73914818) − 5(14.84342798) + (10.17179999)) = 11.78062732 + 0.2 24 (636.00865) = 11.78062732 + 5.300072083 𝑧4 = 17.0806994 Table 3.6 (b), -Appendix D- the exact and numerical solution when applying AMM on system (3.3) and showing the resulting error. These results show the AMM is accuracy with exact solution. 3.8 The numerical approximation of system (3.3) using the Predictor-Corrector method 𝑥4 = 𝑥3 + ℎ 2 ( 𝑓(𝑡3 , 𝑥3 , 𝑦3) + 𝑓(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) 𝑥4 ^ = 𝑥3 + ℎ 𝑓(𝑡3 , 𝑥3 , 𝑦3, 𝑧3) = 18.41964123 + 0.2 𝑓(0.6, 18.41964123, 8.461120366, 11.78062732) = 18.41964123 + 0.2 (35.017176) = 18.41964123 + 7.0034352 = 25.42307643 𝑦4 ^ = 𝑦3 + ℎ 𝑔(𝑡3 , 𝑥3 , 𝑦3, 𝑧3) = 8.461120366 + 0.2 𝑔(0.6, 18.41964123, 8.461120366, 11.78062732) = 8.461120366 + 0.2 (15.10013428) = 8.461120366 + 3.020026856 = 11.48114722 𝑧4 ^ = 𝑧3 + ℎ 𝑞(𝑡3 , 𝑥3 , 𝑦3, 𝑧3) = 11.78062732 + 0.2 𝑞(0.6, 18.41964123, 8.461120366, 11.78062732) 58 = 11.78062732 + 0.2 (21.73914818) = 11.78062732 + 4.347829636 = 16.12845696 𝑥4 = 𝑥3 + ℎ 2 ( 𝑓(𝑡3 , 𝑥3 , 𝑦3, 𝑧3) + 𝑓(𝑡4 , 𝑥4 ^ , 𝑦4 ^, 𝑧4 ^)) 𝑥4 = 18.41964123 + 0.2 2 (𝑓(0.6, 18.41964123, 8.461120366, 11.78062732) + 𝑓(0.8 , 25.42307643 , 11.48114722 , 16.12845696 )) = 18.41964123 + 0.2 2 (35.017176 + 48.65962511) = 18.41964123 + 0.1(83.67680111) = 18.41964123 + 8.367680111 = 26.78732134 𝑦4 = 𝑦3 + ℎ 2 ( 𝑔(𝑡3 , 𝑥3 , 𝑦3, 𝑧3) + 𝑔(𝑡4 , 𝑥4 ^ , 𝑦4 ^, 𝑧4 ^)) 𝑦4 = 8.461120366 + 0.2 2 (𝑔(0.6, 18.41964123, 8.461120366, 11.78062732) + 𝑔(0.8 , 25.42307643 , 11.48114722 , 16.12845696 )) = 8.461120366 + 0.1(15.10013428 + 20.77576669) = 8.461120366 + 0.1(35.87590097) = 8.461120366 + 3.587590097 = 12.04871046 𝑧4 = 11.78062732 + 0.2 2 (𝑞(0.6, 18.41964123, 8.461120366, 11.78062732) + 𝑞 (0.8 , 25.42307643 , 11.48114722 , 16.12845696 )) 59 = 11.78062732 + 0.1(21.73914818 + 30.07038617) = 11.78062732 + 0.1(51.80953435) = 11.78062732 + 5.180953435 = 16.96158076 Table 3.6(c), -Appendix D- the exact and numerical solution when applying PCM on system (3.3) and showing the resulting error. These results show the PCM is accuracy with exact solution. Example 3.3: Consider the system of linear differential equations 𝑑𝑥 𝑑𝑡 = 3𝑥 − 𝑦 + 4𝑒2𝑡 (3.5) 𝑑𝑦 𝑑𝑡 = −𝑥 + 3𝑦 + 4𝑒4𝑡 subject to the initial conditions 𝑥(0) = 1 (3.6) 𝑦(0) = 1 To find the exact solution of system (3.5), we have 𝜆 = 2 , 4 we obtain 𝛷(𝑡) = [−𝑒4𝑡 𝑒2𝑡 𝑒4𝑡 𝑒2𝑡] , 𝛷 −1(𝑡) = [ −1 2 𝑒−4𝑡 1 2 𝑒−4𝑡 1 2 𝑒−2𝑡 1 2 𝑒−2𝑡] 𝑋 = 𝛷𝛷−1(0)𝑋(0) + 𝛷 ∫ 𝛷−1𝐹 𝑡 0 𝑑𝑠 = 𝛷. ( 0 1 ) + 𝛷. (𝑒 −2𝑡 + 2𝑡 − 1 𝑒2𝑡 + 2𝑡 − 1 ) = ( 2 2 ) 𝑡𝑒2𝑡 + ( −1 1 ) 𝑒2𝑡 + ( −2 2 ) 𝑡𝑒4𝑡 + ( 2 0 ) 𝑒4𝑡 60 The exact solution to the system (3.5) is 𝑥(𝑡) = 2𝑡𝑒2𝑡 − 𝑒2𝑡 − 2𝑡𝑒4𝑡 + 2𝑒4𝑡 𝑦(𝑡) = 2𝑡𝑒2𝑡 + 𝑒2𝑡 + 2𝑡𝑒4𝑡 3.9 The numerical approximation of system (3.5) using the fourth order Rung Kutta methods The fourth order Runge-Kutta method to approximate 𝑥(0.8) and 𝑦(0.8) with ℎ = 0.2 in the system (3.5). Table 3.7 (a) and 3.7 (b), -Appendix D- contain the 𝑘’s value when applying algorithm 3.1 on system 3.5 Table 3.8(a), -Appendix D- the exact and numerical solution for 𝑥(𝑡) when applying RKM on system (3.5) and showing the resulting error. From table 3.8(a) we conclude the maximum error = 0.094200336 Figure 3.6(a): -Appendix E- Compares the exact solution 𝑥(𝑡) and approximate solution with time 𝑡 . Figure 3.6(b): -Appendix E- Shows the absolute error resulting of 𝑥(𝑡). Table 3.8(b), -Appendix D- the exact and numerical solution for 𝑦(𝑡) when applying RKM on system (3.5) and showing the resulting error. From table 3.8(b) we conclude the maximum error = 0.119853838 Figure 3.7(a): -Appendix E- Compares the exact solution 𝑦(𝑡) and approximate solution with time 𝑡. Figure 3.7(b): -Appendix E- Shows the absolute error resulting of 𝑦(𝑡). 61 3.10 The numerical approximation of system (3.5) using the Fourth Adams- Bashforth Method The Adams-Bashforth Method to approximate 𝑥(0.8) and 𝑦(0.8) with ℎ = 0.2 in the system (3.5). 𝑥0 ′ = 𝑓(𝑡0 , 𝑥0 , 𝑦0) = 𝑓( 0 , 1 , 1 ) = 6 𝑦0 ′ = 𝑔(𝑡0 , 𝑥0 , 𝑦0) = 𝑔( 0 , 1 , 1) = 6 𝑥1 ′ = 𝑓(𝑡1 , 𝑥1 , 𝑦1) = 𝑓(0.2 , 2.666086766, 2.977368648) = 10.98819044 𝑦1 ′ = 𝑔(𝑡1 , 𝑥1 , 𝑦1) = 𝑔(0.2 , 2.666086766, 2.977368648) = 15.16818289 𝑥2 ′ = 𝑓(𝑡2 , 𝑥2 , 𝑦2) = 𝑓(0.4, 5.502624064, 7.960421513) = 17.44961439 𝑦2 ′ = 𝑔(𝑡2 , 𝑥2 , 𝑦2) = 𝑔(0.4, 5.502624064, 7.960421513) = 38.19077018 𝑥3 ′ = 𝑓(𝑡3 , 𝑥3 , 𝑦3) = 𝑓(0.6 , 9.505252676, 20.49885426) = 21.29737145 𝑦3 ′ = 𝑔(𝑡3 , 𝑥3 , 𝑦3) = 𝑔(0.6 , 9.505252676, 20.49885426) = 96.0841564 𝑥4 ∗ = 𝑥3 + ℎ 24 ( 55𝑥3 ′ − 59𝑥2 ′ + 37𝑥1 ′ − 9𝑥0 ′ ) = 9.505252676 + 0.2 24 (55( 21.29737145 ) − 59(17.44961439) + 37(10.98819044) − 9(6)) = 9.505252676 + 0.2 24 (494.3912272) = 9.505252676 + 4.119926893 = 13.62517957 𝑦4 ∗ = 𝑦3 + ℎ 24 ( 55𝑦3 ′ − 59𝑦2 ′ + 37𝑦1 ′ − 9𝑦0 ′) = 20.49885426 + 0.2 24 (55(96.0841564) − 59(38.19077018) + 37(15.16818289) − 9(6)) 62 = 20.49885426 + 0.2 24 (3538.588187) = 20.49885426 − 29.48823489 = 49.98708915 Table 3.9(a), -Appendix D- the exact and numerical solution when applying ABM on system (3.5) and showing the resulting error. These results show the ABM is accuracy with exact solution. 3.11 The numerical approximation of system (3.5) using the Adams-Moulton method 𝑥4 ′ = 𝑓(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗) = 𝑓(0.8 , 13.62517957 , 49.98708915) = 10.70057926 𝑦4 ′ = 𝑔(𝑡4 , 𝑥4 ∗ , 𝑦4 ∗) = 𝑔(0.8 , 13.62517957 , 49.98708915) = 234.4662087 now 𝑦𝑛+1 = 𝑦𝑛 + ℎ 24 ( 9𝑓𝑛+1 + 19𝑓𝑛 − 5𝑓𝑛−1 + 𝑓𝑛−2) 𝑥4 = 𝑥3 + ℎ 24 ( 9𝑥4 ′ + 19𝑥3 ′ − 5𝑥2 ′ + 𝑥1 ′) = 9.505252676 + 0.2 24 ( 9(10.70057926) + 19(21.29737145) − 5(17.44961439) + (10.98819044)) = 9.505252676 + 0.2 24 (424.6953894) = 9.505252676 + 3.539128245 𝑥4 = 13.04438092 𝑦4 = 𝑦3 + ℎ 24 ( 9𝑦4 ′ + 19𝑦3 ′ − 5𝑦2 ′ + 𝑦1 ′) = 20.49885426 + 0.2 24 (9(234.4662087) + 19(96.0841564) − 5(38.19077018) + (15.16818289)) 63 = 20.49885426 + 0.2 24 (3760.006507) = 20.49885426 + 31.33338756 𝑦4 = 51.83224182 Table 3.9(b), -Appendix D- the exact and numerical solution when applying AMM on system (3.5) and showing the resulting error. These results show the AMM is accuracy with exact solution 3.12 The numerical approximation of system (3.5) using the Predictor-Corrector method 𝑥4 = 𝑥3 + ℎ 2 ( 𝑓(𝑡3 , 𝑥3 , 𝑦3) + 𝑓(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) and 𝑥4 ^ = 𝑥3 + ℎ 𝑓(𝑡3 , 𝑥3 , 𝑦3) = 9.505252676 + 0.2 𝑓(0.6 , 9.505252676 , 20.49885426) = 9.505252676 + 0.2(21.29737145) = 9.505252676 + 0.391535305 = 13.76472697 𝑦4 ^ = 𝑦3 + ℎ 𝑔(𝑡3 , 𝑥3 , 𝑦3) = 20.49885426 + 0.2 𝑔(0.6 , 9.505252676 , 20.49885426) = 20.49885426 + 0.2 (96.0841564) = 20.49885426 + 19.21683128 = 39.71568554 𝑥4 = 𝑥3 + ℎ 2 ( 𝑓(𝑡3 , 𝑥3 , 𝑦3) + 𝑓(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) = 9.505252676 + 0.2 2 (𝑓(0.6 , 9.505252676 , 20.49885426) + 𝑓(0.8 , 13.76472697 , 39.71568554 )) 64 = 9.505252676 + 0.2 2 ( 21.29737145 + 21.39062507 ) = 9.505252676 + 0.1(42.68799652) = 13.77405233 𝑦4 = 𝑦3 + ℎ 2 ( 𝑔(𝑡3 , 𝑥3 , 𝑦3) + 𝑔(𝑡4 , 𝑥4 ^ , 𝑦4 ^)) = 20.49885426 + 0.2 2 (𝑔(0.6 , 9.505252676 , 20.49885426) + 𝑔(0.8 , 13.76472697 , 39.71568554)) = 20.49885426 + 0.1( 96.0841564 + 203.5124504) = 20.49885426 + 0.1(299.5966068) = 50.45851494 Table 3.9(c), -Appendix D- the exact and numerical solution when applying PCM on system (3.5) and showing the resulting error. These results show the PCM is accuracy with exact solution. 65 Chapter Four Comparison of Numerical Methods In this chapter, we will compare the proposed numerical methods in the examples with exact solution and show the conclusions. In this chapter, we will carry some comparison between the proposed numerical methods, namely; Rung-Kutta Method, Adams-Bashforth Method, Adams-Moulton Method and Predictor-Corrector Method for solving the numerical examples presented in chapter three. For example 3.1 the exact and numerical solutions when applying Rung-Kutta Method, Adams-Bashforth Method, Adams-Moulton Method and Predictor-Corrector Method for solving system (3.1) and the error associated with their numerical methods are shown in Table 4.1 -Appendix D-. In fact, we see clearly that the Adams-Moulton Methods is more accurate with less absolute error than its counterparts. On the other hand, the Predictor-Corrector method requires less Cup-time than the other methods. For example 3.2 the exact and numerical solutions when applying Rung-Kutta Method, Adams-Bashforth Method, Adams-Moulton Method and Predictor-Corrector Method for solving system (3.3) and the error associated with their numerical methods are shown in Table 4.2 -Appendix D-. In fact, we see clearly that the Rung-Kutta Method is more accurate with less absolute error than its counterparts. On the other hand, the Predictor-Corrector method requires less Cup-time than the other methods. For example 3.3 the exact and numerical solutions when applying Rung-Kutta Method, Adams-Bashforth Method, Adams-Moulton Method and Predictor-Corrector Method for solving system (3.5) and the error associated with their numerical methods are shown in Table 4.3 -Appendix D-. 66 In fact, we see clearly that the Rung-Kutta Method is more accurate with less absolute error than its counterparts. On the other hand, the Adams-Moulton Method requires less Cup-time than the other methods. 4.1 Conclusions This work explores the numerical treatment of systems of ordinary differential equations (ODEs), which has wide rang of applications in mathematical physics, chemistry, biology, stereology, heat conduction, and engineering models. We implemented four numerical methods namely; Rung-Kutta Method, Adams- Bashforth Method, Adams-Moulton Method and Predictor-Corrector Method to solve systems of ordinary differential equations. Numerical results presented in Tables and Figures show clearly that the Rung-Kutta Method is one of the efficient method in comparison with its counterparts. 67 References [1] G. Adomian, Applications of Nonlinear Stochastic Systems Theory to Physics. Kulwer Academic, Dordrecht, 1989. [2] G. Adomian, A Review of The Decomposition Method in Applied Mathematics, J. Math. Anal. Appl. 135, 501-544, (1988). [3] G. Adomian, Nonlinear Stochastic Operator Equations. Academic Press. 1986. [4] G. Adomian, Solving Frontier Problems of Physics, The Decomposition Method. Kluwer Academic Publisher, Dordrecht, 1994. [5] F. Bashforth and J. Adams, An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid with an explanation of the method of integration employed in constructing the tables which give the theoretical forms of such drop, Cambridge University Press, Cambridge, 1883. [6] J. Biazar, E. Babolian and R. Islam, Solution of a System of Ordinary Differential Equations with Adomian Decomposition Method, Appl. Math. Comput., 147(3), 2004. [7] J. Biazar, Solution of a System Integral-Differential Equations by Adomian Decomposition Method. Appl. Math. Comput. 168. 2005. [8] W. Boyce, R. Diprima and D. Meade, Elementary Differential Equations and Boundary Value Problems, USA, 11th, (2017). [9] R. Burden, J. Faires and A. Burden, Numerical Analysis, 10th edition, Cengage Learning, 2014. [30] R. Burden and J. Faires, Numerical Methods, Brooks Cole, 2002. [31] J. Butcher, Numerical Method for Ordinary Differential Equations, 3rd ed., Wiley, UK, 2016. [32] J. Butcher, “Numerical Methods ForOrdinary Differential Equations in The 20thCentury”, Journal of Computational and Applied Mathematics, 125, 2000. 68 [33] S. Goode and S. Annin, Differential Equations and Linear Algebra, California State University, Fullerton, 4th ed., (2015), ppl.580-598. [34] H. Jafari, A. Borhanifar, and S. Karimi, New Solitary Wave Solution for the bad Boussinesq and good Boussinesq Equations, Numer. Methods Partial Differential Equations, 25(2009), pp.1231-1237. [15] H. Jafari and V. Daftardar-Gejji, Revised Adomian Decomposition Method for Solving Systems of Ordinary and Fractional Differential Equations, Applied Mathematics and Computation 181 (2006) pp.598-608. [16] H. Khalil, Nonlinear