An-Najah National University Faculty of Graduate Studies Numerical Methods for Solving Lienard’s Equation By Hiyam Ahmad Aldalou Supervisor Prof. Dr. Naji Qatanani This Thesis is submitted in Partial Fulfillment of the Requirements for the Degree of Master of Mathematics, Faculty of Graduate Studies, An- Najah National University, Nablus, Palestine. 2020 II III Dedication To Mustafa,,, IV Acknowledgment I am heartily thankful to my supervisor, Prof. Dr. Naji Qatanani, whose encouragement, guidance and support from initial to the final level enabled me to develop and understanding the subject. My thanks and appreciation goes also to my thesis committee members Dr. Adnan Daraghmeh and Dr. Rania Wannan for their encouragement, support, interest and valuable hints. I acknowledge An-Najah National University for supporting this work, and I wish to pay my great appreciation to all respected teachers and staff in department of mathmatics. Lastly, I offer my regards and blessings to all of those who supported me during the completion of this thesis. V VI Table of Contents Dedication III Acknowledgement IV Declaration V Table of contents VI List of Tables VIII List of figures IX Abstract X Chapter 1: 1 1.1 Introduction 2 Chapter 2: Lienard’s Equation 8 2.1 Historical review of Lienard’s Equation 9 2.2 Fractional Lienard’s Equation 16 Chapter 3: Numerical method for solving Lienard’s equation 19 3.1 Variational Iteration Method (VIM) 20 3.1.1 The Variational Iteration Method (VIM) for Lienard equation 20 3.2 Variational Homotopy Perturbation Method (VHPM) 23 3.2.1 The Homotopy Perturbation Method 23 3.2.2 The Variational Homotopy Perturbation Method 24 3.2.3 The Variational Homotopy Perturbation Method (VHPM) for Lienard’s equation 26 3.3 Adomian Decomposition Method (ADM) 28 3.3.1 The Adomian Decomposition Method (ADM) for Solving Lienard Equation 31 Chapter 4: Numerical examples on Lienard equation 34 4.1 Example (1) 35 4.1.1 Solving (1) Using VIM 36 4.1.2 Solving (1) Using VHPM 41 VII 4.1.3 Solving (1) Using ADM 46 4.2 Example (2) 53 4.2.1 Solving (2) Using VIM 54 4.2.2 Solving (2) Using VHPM 60 4.2.3 Solving (2) Using ADM 66 Conclusion 73 References 74 Appendix: 81 (1) MAPLE code for solving example 1 81 (1.1) Using VIM 81 (1.2) Using VHPM 82 (1.3) Using ADM 83 (2) MAPLE code for solving example 1 84 (2.1) Using VIM 84 (2.2) Using VHPM 85 (2.3) Using ADM 86 ب الملخص VIII List of Tables Table 4.1: The exact and first approximate solutions using VIM for example (1). 38 Table 4.2: The exact and second approximate solutions using VIM for example (1). 39 Table 4.3: The exact and first approximate solutions using VHPM for example (1). 43 Table 4.4: The exact and second approximate solutions using VHPM for example (1). 44 Table 4.5: The exact and first approximate solutions using ADM for example (1). 50 Table 4.6: The exact and second approximate solutions using ADM for example (1). 51 Table 4.7: The exact and first approximate solutions using VIM for example (2). 57 Table 4.8: The exact and second approximate solutions using VIM for example (2). 58 Table 4.9: The exact and first approximate solutions using VHPM for example (2). 63 Table 4.10: The exact and second approximate solutions using VHPM for example (2). 64 Table 4.11: The exact and first approximate solutions using ADM for example (2). 70 Table 4.12: The exact and second approximate solutions using ADM for example (2). 71 IX List of Figures Figure 4.1: The exact and first approximate solutions using VIM for example (1). 39 Figure 4.2: The exact and second approximate solutions using VIM for example (1). 40 Figure 4.3: The errors of first two approximate solutions using VIM for example (1). 40 Figure 4.4: The exact and first approximate solutions using VHPM for example (1). 44 Figure 4.5: The exact and second approximate solutions using VHPM for example (1). 45 Figure 4.6: The errors of first two approximate solutions using VHPM for example (1). 45 Figure 4.7: The exact and first approximate solution using ADM for example (1). 51 Figure 4.8: The exact and second approximate solutions using ADM for example (1). 52 Figure 4.9: The errors of first two approximate solutions using ADM for example (1). 52 Figure 4.10: The exact and first approximate solutions using VIM for example (2). 58 Figure 4.11: The exact and second approximate solutions using VIM for example (2). 59 Figure 4.12: The errors of first two approximate solutions using VIM for example (2). 59 Figure 4.13: The exact and first approximate solutions using VHPM for example (2). 64 Figure 4.14: The exact and second approximate solutions using VHPM for example (2). 65 Figure 4.15: The errors of first two approximate solutions using VHPM for example (2). 65 Figure 4.16: The exact and first approximate solutions using ADM for example (2). 71 Figure 4.17: The exact and second approximate solutions using ADM for example (2). 72 Figure 4.18: The errors of first two approximate solutions using ADM for example (2). 72 X Numerical Methods for Solving Lienard’s Equation By Hiyam Ahmad Aldalou Supervisor Prof. Dr. Naji Qatanani Abstract In this thesis, we focus on the numerical handling on one of the most important equation in physics, namely, the Lienard equation. This equation has a wide range of applications in both science and technology. The numerical techniques used to approximate the solution of the Lienard equation are: the Variational Iteration Method (VIM), the Variational Homotopy Perturbation Method (VHPM) and the Adomian Decomposition Method (ADM). The efficiency and the accuracy of these numerical techniques have been proved by solving some numerical examples. Numerical results show clearly that the variational iteration method is one of the most powerful numerical technique for solving the Lienard equation in comparison with other numerical techniques. 1 Chapter One 2 Chapter 1 1.1 Introduction Nonlinear problems appearing in many physical phenomena, engineering and scientific applications are modeled with nonlinear differential equations. One of these equations is the Lienard's equation, which is a second order differential equation, named after the French physicist Alfred-Marie Lienard, in 1928. Lienard’s equation has the general form: 𝑦′′(𝑥) + 𝑓(𝑦)𝑦′(𝑥) + 𝑔(𝑦) = ℎ(𝑥) (1.1) The Lienard equation is one of the most important differential equations not only in mathematics but also in the natural science and technology there are a lot of examples which related to equation (1.1). Such systems are used to model electrical, mechanical, or biomedical systems, and in the literature, many systems are transformed into ones of Lienard type to support in investigations. A lot of chemical kinetics models can be formulated as Lienard type coupled nonlinear first-order rate equations in several variables. The first- order approximation for the Lienard system works well near a bifurcation point, with higher-order terms being required the further the system is from the bifurcation point. The dynamics of a scalar inflaton field with a 3 symmetric double-well potential can also be formulated mathematically as a Lienard system [15]. The Lienard’s equation applied in the dynamics of malarial transmission to guarantee a globally asymptotically stable positive equilibrium or a stable persistence of infection in terms of the rate of biting per female mosquito on uninfected hosts [30]. Different choices of 𝑓, 𝑔 and ℎ will produce different models of the Lienard’s equation. Such that, if: 𝑓(𝑦) = (𝑦2 − 1), (1.2) 𝑔(𝑦) = 𝑦, (1.3) ℎ(𝑥) = 0. (1.4) Where is a scalar parameter. Then equation (1.1) will be the well-known Van Der Pol equation, which gives considerable insight into the potential of a monolithic version of the circuit for optical communication functions including clock recovery and chaotic source applications. These appears when it used to treat the oscillator resonant tunneling diod-laser diode hybrid integrated circuit [35]. Some famous nonlinear oscillators such as the Van Der Pol equation, the Duffing oscillator and the Helmholtz oscillator belong to family of the Lienard equation. Moreover, the Lienard equation often appears as a traveling wave reduction of nonlinear partial differential equations. It is an 4 important problem to find subclasses of the Lienard equation which can be analytically solved, which used for the Fisher equation, the Burgers- Korteweg-de Vries equation and the Burgers-Huxley equation [21]. Many researchers have studied the Lienard equation, such as Harko who study a class of exact solutions of the Lienard equation in [15], LI study the periodic solution with singularities for the Lienard equation in [22], Gaiko study the applied geometry of the Lienard equation in [12], Aghajani with Moradifam have studied the Lienard plane in [2]. Here, we will study a special case of Lienard equation, if: 𝑓(𝑦) = 0, (1.5) 𝑔(𝑦) = 𝑎 𝑦(𝑥) + 𝑏 𝑦3(𝑥) + 𝑐 𝑦5(𝑥), (1.6) ℎ(𝑥) = 0. (1.7) Then we have: 𝑦′′(𝑥) + 𝑎 𝑦(𝑥) + 𝑏 𝑦3(𝑥) + 𝑐 𝑦5(𝑥) = 0, (1.8) where 𝑎, 𝑏 and 𝑐 are real coefficients. The Lienard equation of this form is well known. For example, some important equations (such as the Rangwala Rao equation, the Ablowitz equation, the Gerdjikov Ivanov equation, etc.) are proposed as integrable conditions when one studies the Backlund transformation and solitary wave solutions for the nonlinear Schrodinger equation. We observe that finding a 5 kind of solitary wave solution for these equations can be always reduced to finding a solution for the Lienard equation (1.8) [20]. Form of Lienard equation (1.8) has been considered by many authors. Feng [9] obtained an explicit exact solutions to equation (1.8), and applied the results to find some explicit exact solitary wave solutions to the nonlinear Schrodinger equation and the Pochhammer-Chree equation, Kong [20] discussed the explicit exact solution for equation (1.8) and its application, Jianhau [37] studied explicit exact solution for the Lienard equation and its applications and also obtained more explicit exact solitary wave solutions for the generalized Pochhammer-Chree equation given by seeking qualitatively the homoclinic and heteroclinic orbits for this class of Lienard equation, etc. Many researchers solved equation (1.8) using several numerical method. For example, the Differential Transform Method [26], the Block Pulse Functions Method [16], the Hybrid Heuristic Computation Method [23], the Residual Power Series Method [38], the Chebyshev Operational matrix [34], etc. In this study, we will propose and implement three numerical methods for solving the Lienard equation. These methods are: 1) The Variational Iteration Method (VIM): In 1998, the Lagrange multiplier method is modified into an iteration method that is called variational iteration method (VIM). This method is used to solve easily, effectively, and accurately a large class of non-linear 6 problems with approximations converging rapidly to accurate solutions. The variational iteration method changes the differential equation to a recurrence sequence of functions, where the limit of that sequence is considered as the solution of the partial differential equations. The main advantage of the method is that it can be applied directly to all types of nonlinear differential and integral equations, homogeneous or inhomogeneous, with constant or variable coefficients. See [4, 6, 18, 33]. 2) The Variational Homotopy Perturbation Method (VHPM): Researchers combine the VIM with Homotopy Perturbation Method. This combination gives a new method called the Variational Homotopy Perturbation Method (VHPM). The VHPM provides the solution in a rapid convergent series, and applied directly without any discretization, restrictive assumption, or transformation. See [10, 27, 31]. 3) The Adomian Decomposition Method (ADM): George Adomian established the Adomian’s Decomposition method (ADM) in the 1980s. The ADM is a powerful method for solving functional equations. In this method, special kinds of polynomials are used which called Adomian polynomials. The solution is given by an infinite series in which each term is easily obtained. The convergence of the decomposition series and the speed of the convergence was studied by several researchers. See [1, 5, 7, 17, 19, 25, 41] for more details. 7 This thesis organized as follows: Chapter one contains an overall review and discussion of the Lienard equation and its solvability. In chapter two, we propose three numerical methods for solving the Lienard equation, namely, the Variational Iteration Method (VIM), the Variational Homotopy Perturbation Method (VHPM) and the Adomian Decomposition Method (ADM). Some numerical examples on the Lienard equation are illustrated in chapter three. Moreover, a comparison between the analytical and numerical results is also presented. Finally, conclusions are drawn. 8 Chapter Two 9 Chapter 2 Lienard’s Equation 2.1 Historical Review of Lienard’s Equation Models involving Lienard’s equation and Lienard systems play very important role in science and technology. For example, during the development of radio and vacuum tubes, Lienard equation and Lienard system were intensely studied as they can used to model oscillating circuils. Also, the Lienard system has been shown to describe the operation of an optoelectronic circuit that uses a resonant tunneling diode to drive a laser diode to make an optoelectronic voltage controlled oscillator. The dynamic behaviors of Lienard equation and the Lienard system have been widely investigated due to their applications in many fields such as physics, mechanics and engineering technique fields. In such applications, it is important to know the existence of periodic solutions of the Lienard system. In applied science some parochial problems associated with the Lienard system can be found in [3] and [2]. The general form of the Lienard’s equation is: 𝑦′′(𝑥) + 𝑓(𝑦)𝑦′(𝑥) + 𝑔(𝑦) = ℎ(𝑥) (2.1) Different choices of the functions 𝑓(𝑦)and 𝑔(𝑦) gives different cases of the equation. For example, if 𝑓(𝑦)𝑦′ is the damping force, 𝑔(𝑦) is the restoring 10 force and ℎ(𝑥) is the external force, we get the damped pendulum equation. However, if we choose: 𝑓(𝑦) = (𝑦2 − 1), (2.2) 𝑔(𝑦) = 𝑦, (2.3) ℎ(𝑥) = 0. (2.4) We get a nonlinear model of electronic oscillation, which named Van Der Pol equation. The Van Der Pol equation is considered as basic model for oscillatory process for physics, biology, electronics, and neurology. Van Der Pol himself built a number of electronic circuits to model human heart using his equation [29]. Moreover, the Lienard equation can be used to model fluid mechanical phenomena. The linearly forced isotropic turbulence can be described in terms of a cubic Lienard equation with linear damping of the form [15]: 𝑦′′(𝑥) + [𝑎 𝑦(𝑥) + 𝑏]𝑦′(𝑥) + 𝑐 𝑦(𝑥) − 𝑦3(𝑥) + 𝑑 = 0 (2.5) Where 𝑎, 𝑏, 𝑐 and 𝑑 are constants, also naturally appear in the mathematical description of some important astrophysical phenomena. For example, the time-dependence 𝜑(𝑡) of the perturbations of the stationary solutions of spherically symmetric accretion processes can be described by a generalized Lienard type equation of the form: 𝜑′′ + 𝜖 𝑓(𝜑, 𝜑′)𝜑′ + 𝑉′(𝜑) = 0 (2.6) 11 Where 𝜖 is a small parameter, and 𝑉(𝜑) is potential of the system, with the prime indicating the derivative with respect to 𝜑 . A dynamical systems analysis of this Lienard equation reveals a saddle point in real time, with the implication that when the perturbation is extended into the nonlinear regime, instabilities will develop in accreting system. In physical view, the Lienard equation represents the generalization of the equation of damped oscillations 𝑦′′ + 𝛾 𝑦′ +𝜔2𝑦 = 0 (2.7) where 𝛾 and 𝜔2 are constant parameters, respectively. For 𝛾 = 0, we get equation of the linear harmonic oscillator, which represents one of the fundamental equations of both classical and quantum physics [15]. Models involving Lienard’s equation and Lienard systems play very important role in science and technology. For example, during the development of radio and vacuum tubes, Lienard equation and Lienard system were intensely as they can used to model oscillating circles. Another example of Lienard equation is given by the Duffing’s equation, which named after Georg Duffing [42], is a non-linear second- order differential equation used to model certain damped and driven oscillators. Duffing equation is given by: 𝑦′′ + 𝛿 𝑦′ + 𝛼 𝑦 + 𝛽 𝑦3 = 𝛾 cos(𝜔𝑡) (2.8) where 𝑦 = 𝑦(𝑡) is the displacement at time 𝑡, 𝑦′ = 𝑑𝑦 𝑑𝑡⁄ (i.e. velocity), and 𝑦′′ = 𝑑2𝑦 𝑑𝑦2⁄ (i.e. acceleration), The numbers 𝛿, 𝛼, 𝛽, 𝛾 and 𝜔 are 12 given constant. Equation (2.8) describes the motion (which corresponds to the case 𝛽 = 𝛿 = 0), in physical terms, it models for example a spring pendulum whose springs stiffness does not exactly obey Hooke’s law [42]. Also, the Lienard system has been shown to describe the operation of an optoelectronic circuit that uses a resonant tunneling diode in order to drive a laser diode that made an optoelectronic voltage controlled oscillator. The classical Lienard equation: 𝑦′′(𝑥) + 𝑓(𝑦)𝑦′(𝑥) + 𝑔(𝑦) = 0 (2.9) is equivalent to system (2.10) in the phase plane. And equivalent to system (2.11) in the Lienard plane [32]. { 𝑦′ = 𝑣 𝑣′ = −𝑓(𝑦)𝑣 − 𝑔(𝑦) (2.10) { 𝑦′ = 𝑣 − 𝐹(𝑦) 𝑣′ = −𝑔(𝑦) (2.11) Where 𝐹(𝑦) = ∫ 𝑓(𝑡)𝑑𝑡 𝑦 0 , and 𝐺(𝑦) = ∫ 𝑔(𝑡)𝑑𝑡 𝑦 0 . The boundedness of solutions of the Lienard equation was put in researcher consideration. As an example, Sugie in [36] states his result as: Theorem (2.1) [36]: Suppose that 𝑔(𝑦)𝐹(𝑦) ≥ 0 for all 𝑦 ∈ ℝ. If there are sequences {𝑦𝑛} tending to ∞ and {�̅�𝑛} tending to −∞ such that 𝐺(𝑦𝑛) → ∞ and 𝐺(�̅�𝑛) → ∞ as 𝑛 → ∞ Then all solutions of system (2.11) are bounded. 13 On the other hand, Hara and Yoneyama [14] state a proposition that gives sufficient condition include the existence of an unbounded solution of the Lienard equation. Proposition (2.1) [14]: Suppose that ∫ 𝑔(𝑦)𝑑𝑦 ∞ 0 < ∞ and sup 𝑦≥0 𝐹(𝑦) < ∞ Or ∫ 𝑔(𝑦)𝑑𝑦 −∞ 0 < ∞ and inf 𝑦≤0 𝐹(𝑦) > −∞ Then there exists an unbounded solution of system (2.11). Several researcher studied the existence of periodic solution of the Leinard’s equation. For example, Villari [40] who studied equation (2.9) and found: Theorem (2.2) [40]: Let 𝑓:ℝ → ℝ and 𝑔:ℝ → ℝ be such that: (1) 𝑓 is continuous and 𝑔 is locally Lipschitz continuous, (2) 𝑦 𝑔(𝑦) > 0 if 𝑦 ≠ 0, (3) ∫ 𝑔(𝑠)𝑑𝑠 ∞ 0 = +∞, (4) 𝑓(0) < 0, and there exist 𝛿1, 𝛿2 ∈ ℝ with 𝛿1 < 0 < 𝛿2 and 𝑓(𝑦) > 0 if 𝑦 ∉ [𝛿1, 𝛿2], (5) There exist 𝛼 ≤ 𝛿1 and 𝛽 ≥ 𝛿2 such that ∫ 𝑓(𝑠)𝑑𝑠 𝛽 𝛼 > 0. Then there exists a periodic solution of equation (2.9). 14 Figueiredo proved in [11] that equation (2.9) (or equivalently system (2.11)) has at least 𝑛 periodic solutions. But Zuo-Huan [43] studied equation (2.9), and put a theorem discussed the non-existence of periodic solutions for the equivalent system. Theorem (2.3) [43]: Suppose system (2.11) satisfies the following conditions: (1) 𝑓(𝑦) and 𝑔(𝑦) are locally Lipschitz, and 𝑦𝑔(𝑦) > 0, 𝑦 ≠ 0. (2) ∃ℎ > 0 and 𝑀 > 0 such that ℎ𝐺(𝑦) − 𝐹(𝑦) ≤ 𝑀, 0 < 𝑦 < +∞. (3) ∃𝛼 < 0, 𝛽 > 0 such that 𝐺(𝛽) > 1 2 (𝑀 + 1 ℎ ) 2 , 𝐺(𝛼) > 1 2 (𝑀 + 1 ℎ ) 2 , 𝑦𝐹(𝑦) > 0, 𝑦 ∈ (𝛼, 0) ∪ (0, 𝛽). Then there are no periodic solutions of system (2.11). The Lienard equation in its general form is very difficult to solve. Many researchers have studied special cases of the Lienard equation, in this study, we discuss this form: 𝑦′′(𝑥) + 𝑎𝑦(𝑥) + 𝑏𝑦3(𝑥) + 𝑐𝑦5(𝑥) = 0 (2.12) where 𝑎, 𝑏 and 𝑐 are real coefficients. 15 Several researchers have studied the exact solution of this special case. For example, Feng [9] investigated the exact solution of equation (2.12) with some conditions on the coefficients of the equation 𝑎 , 𝑏 and 𝑐. Theorem (2.4) [9]: (i) Suppose that 𝑎 < 0 , 𝑐 ≥ 0, then equation (2.12) admits exact solutions 𝑦(𝑥) = ±(− 4𝑎𝑑0 𝑒 −2√−𝑎𝑥 (𝑑0 𝑒 −2√−𝑎𝑥 + 𝑏 2) 2 − 4𝑎𝑐 3 ) 1/2 Where 𝑑0 is an arbitrary positive integration constant. (ii) In particular, (a) If 𝑎 < 0, 𝑏 > 0, 𝑐 ≥ 0 or 𝑎 < 0, 𝑏 ≤ 0, 𝑐 > 0, then equation (2.12) admits exact solutions 𝑦1(𝑥) = ± [ 4√ 3𝑎2 3𝑏2 − 16𝑐𝑎 sech2√−𝑎𝑥 2 + (−1 + √3𝑏 √3𝑏2 − 16𝑐𝑎 ) sech2√−𝑎𝑥 ] 1/2 (b) If 𝑎 < 0, 𝑏 > 0 and 3𝑏2 − 16 𝑎𝑐 = 0, then equation (2.12) admits exact solutions: 𝑦2(𝑥) = ± √ −2𝑎 𝑏 (1 + 𝑡𝑎𝑛ℎ√−𝑎𝑥) 𝑦3(𝑥) = ± √ −2𝑎 𝑏 (1 − 𝑡𝑎𝑛ℎ√−𝑎𝑥) 16 The generalization for Theorem (2.4) is in [20]. Since most of the nonlinear differential equations are difficult to solve analytically, these equations must be tackled numerically. We will study some of these numerical methods to approximate solutions for equation (2.12). 2.2 Fractional Lienard’s Equation Here we will talk about the Fractional Lienard equation, but first of all, what is the fractional? Fractional calculus concerns the generalization of differentiation and integration to non-integer (fractional) orders. The subject has a long mathematical history being discussed for the first time already in the correspondence of Leibniz with L’Hopital when this replied “What does 𝑑𝑛 𝑑𝑥𝑛 𝑓(𝑥) mean if 𝑛 = 1 2 ?” in 1695. Over the centuries many mathematicians have built up a large body of mathematical knowledge on fractional integrals and derivatives [24] [38]. The fractional calculus of variations has gained an importance due to its wide applications in different fields of sciences and engineering. This new field of applied mathematics is based on fractional calculus which has led to many break-through in classical and modern physics. It is describing the problems of life sciences, for instance polarization, electro-magnetic waves, viscoelasticity, electrode-electrolyte heat conduction, many others physical processes, finance, control theory, biomedical engineering, and biology. It 17 was observed that real world problems mainly dissipative systems are more adequately described by means of fractional operators than integer operators. The fractional calculus allowed the investigation of the nonlocal response of multiple phenomena, the fractional derivatives may be memory operators which usually represent dissipative effects or damage. The fractional derivative considers the history and nonlocal distributed effects of any physical system [8, 29, 39]. As a generalization of the Lienard’s equation (2.12) to the fractional case, we will study the following class of fractional Lienard’s equation of the form: 𝐷2𝛼𝑦(𝑥) + 𝑎𝑦(𝑥) + 𝑏𝑦3(𝑥) + 𝑐𝑦5(𝑥) = 0 , 1 2 < 𝛼 ≤ 1 (2.13) for 𝑥 > 0 , Subject to: 𝑦(0) = 𝑦0 , 𝐷 𝛼𝑦(0) = 𝑦1 (2.14) where 𝑎, 𝑏, 𝑐, 𝑦0 and 𝑦1 are constant. There are several mathematical definitions of the fractional derivatives, here we adopt the most usual used definitions: the Caputo and its reverse operator Riemann- Liouville. Definition (2.1) [38]: Let 𝑛 be the smallest integer greater than or equal to α. The Caputo fractional derivative of order 𝛼 > 0 is defined as: 18 𝐷𝛼𝑦(𝑥) = { 1 Г(𝑛 − 𝛼) ∫ (𝑥 − 𝑡)𝑛−𝛼−1𝑦(𝑛)(𝑡)𝑑𝑡, 𝑛 − 1 < 𝛼 < 𝑛, 𝑥 0 𝑦(𝑛)(𝑥), 𝛼 = 𝑛 ∈ ℕ. where Г(𝑛 − 𝛼) is the Gamma function. The power rule of the Caputo derivative is given as follows: Theorem (2.5) [38]: The Caputo fractional derivative of the power function is given by: 𝐷𝛼𝑦(𝑥) = { Г(𝑝 + 1) Г(𝑝 − 𝛼 + 1) 𝑥𝑝−𝛼 , 𝑛 − 1 < 𝛼 < 𝑛, 𝑝 > 𝑛 − 1 , 𝑝 ∈ ℝ 0 , 𝑛 − 1 < 𝛼 < 𝑛, 𝑝 ≤ 𝑛 − 1 , 𝑝 ∈ ℕ. Syam [38] investigated the exact solution of equation (2.13) subject to equation (2.14), using the Residul Power Series method. Then, 𝑦2 = −(𝑎𝑦0 + 𝑏𝑦0 3 + 𝑐𝑦0 5), (2.15) 𝑦3 = −(𝑎𝑦1 + 3𝑏𝑦0 2𝑦1 + 5𝑐𝑦0 4𝑦1). (2.16) 19 Chapter Three 20 Chapter 3 Numerical methods for solving Lienard’s equation 3.1 The Variational Iteration Method (VIM) To illustrate the basic concepts of the VIM, we consider the following nonlinear differential equation: 𝐿𝑢 + 𝑁𝑢 = 𝑔(𝑥), (3.1) Where 𝐿 is linear operator, 𝑁 is a nonlinear operator, and 𝑔(𝑥) is an inhomogeneous term. According to the VIM, we can construct a correction functional as follows: 𝑢𝑛+1(𝑥) = 𝑢𝑛(𝑥) + ∫ 𝜆(𝜏){𝐿𝑢𝑛 +𝑁�̃�𝑛 − 𝑔(𝜏)}𝑑𝜏 𝑥 0 , (3.2) where 𝜆(𝜏) is a general Lagrange multiplier which can be identify optimally via the variational theory, the subscript 𝑛 denotes the 𝑛 th-order approximation and �̃�𝑛 is considered as a restricted variation, i.e. �̃�𝑛 = 0. 3.1.1 The Variational Iteration Method (VIM) for Lienard equation We consider VIM for equation (2.12). For 𝑖 ≥ 0 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ 𝜆(𝜏){𝑦𝑖 ′′(𝜏) + 𝑎𝑦𝑖(𝜏) + 𝑏�̃�𝑖 3(𝜏) + 𝑐�̃�𝑖 5(𝜏)}𝑑𝜏 𝑥 0 (3.3) 21 where �̃�𝑖 is considered as restricted variational, i.e. 𝛿�̃�𝑖 = 0. To find optimal value of 𝜆(𝜏), we have 𝛿𝑦𝑖+1(𝑥) = 𝛿𝑦𝑖(𝑥) + 𝛿 ∫ 𝜆(𝜏){𝑦𝑖 ′′(𝜏) + 𝑎𝑦𝑖(𝜏) + 𝑏�̃�𝑖 3(𝜏) + 𝑐�̃�𝑖 5(𝜏)}𝑑𝜏 𝑥 0 (3.4) When 𝛿�̃�𝑖 = 0, we get: 𝛿𝑦𝑖+1(𝑥) = 𝛿𝑦𝑖(𝑥) + 𝛿∫ 𝜆(𝜏){𝑦𝑖 ′′(𝜏) + 𝑎𝑦𝑖(𝜏)}𝑑𝜏 = 0, 𝑥 0 (3.5) Which implies: 𝛿𝑦𝑖+1(𝑥) = 𝛿𝑦𝑖(𝑥) + 𝛿𝜆(𝜏)𝑦𝑖 ′(𝜏)|𝜏=𝑥 − 𝛿𝜆′(𝜏)𝑦𝑖 (𝜏)|𝜏=𝑥 + ∫ {𝜆′′(𝜏) + 𝑎𝜆 (𝜏)}𝑦𝑖(𝜏)𝑑𝜏 𝑥 0 = 0. (3.6) Therefore, the stationary conditions are obtained in the following form: 1 − 𝜆′(𝜏) = 0|𝜏=𝑥, 𝜆(𝜏) = 0|𝜏=𝑥 , 𝜆′′(𝜏) + 𝑎𝜆(𝜏) = 0|𝜏=𝑥 . (3.7) Using above stationary conditions, we can find several values for 𝜆(𝜏): 22 𝜆1(𝜏) = 𝜏 − 𝑥, 𝜆2(𝜏) = 1 √−𝑎 sinh (√−𝑎(𝜏 − 𝑥)) , 𝜆3(𝜏) = sin(𝜏) cos(𝑥) − cos(𝜏) sin(𝑥) = sin(𝜏 − 𝑥). (3.8) Substituting each of above values instead of 𝜆(𝜏) into equation (3.3), we can obtain a new iteration formula, i.e.: substituting 𝜆1(𝜏) obtain the iteration formula (3.9), 𝜆2(𝜏) obtain (3.10) and 𝜆3(𝜏) obtain (3.11). 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ (𝜏 − 𝑥){𝑦𝑖 ′′ + 𝑎𝑦𝑖 + 𝑏𝑦𝑖 3 + 𝑐𝑦𝑖 5}𝑑𝜏, 𝑥 0 (3.9) 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ 1 √−𝑎 sinh (√−𝑎(𝜏 − 𝑥)) {𝑦𝑖 ′′ + 𝑎𝑦𝑖 + 𝑏𝑦𝑖 3 𝑥 0 + 𝑐𝑦𝑖 5}𝑑𝜏, (3.10) 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ sin(𝜏 − 𝑥) {𝑦𝑖 ′′ + 𝑎𝑦𝑖 + 𝑏𝑦𝑖 3 + 𝑐𝑦𝑖 5}𝑑𝜏. 𝑥 0 (3.11) These iteration formulas can obtain a sequence which tends to the exact solution of the Lienard equation (2.12). The difference between them is in their convergence speed. There exists also various values for initial approximation for these iteration formulas. Such as: 𝑦0(𝑥) = 𝑦(0) + 𝑦 ′(0)𝑥, (3.12) 𝑦0(𝑥) = 𝑦(0) 𝑐𝑜𝑠(√−𝑎𝑥) + 𝑦 ′(0) 𝑠𝑖𝑛(√−𝑎𝑥). (3.13) 23 3.2 The Variational Homotopy Perturbation Method (VHPM) Since the variational homotopy perturbation method comes by combine the variational iteration method with homotopy perturbation method, we will include review for the homotopy perturbation method. 3.2.1 The Homotopy Perturbation Method Let 𝐿(𝑢) = 0 (3.14) a general equation, where 𝐿 is any integral or differential operator. We define a convex homotopy 𝐻(𝑢, 𝑝) by 𝐻(𝑢, 𝑝) = (1 − 𝑝)𝐹(𝑢) + 𝑝 𝐿(𝑢), (3.15) Where 𝐹(𝑢) is a functional operator with known solutions 𝑣0, which can be obtained easily. It is clear that, for 𝐻(𝑢, 𝑝) = 0 (3.16) We have 𝐻(𝑢, 0) = 𝐹(𝑢), 𝐻(𝑢, 1) = 𝐿(𝑢). (3.17) This shows that 𝐻(𝑢, 𝑝) continuously traces an implicitly defined curve from a starting point 𝐻(𝑣0, 0) to a solution function 𝐻(𝑓, 1). The embedding parameter monotonically increases from zero to unit as the trivial problem 𝐹(𝑢) = 0 is continuously deforms the original problem 𝐿(𝑢) = 0. The embedding parameter 𝑝 ∈ [0.1] can be considered as an expanding 24 parameter. The homotopy perturbation method uses the homotopy parameter 𝑝 as an expanding parameter to obtain 𝑢 =∑𝑝𝑖𝑢𝑖 = 𝑢0 + 𝑝 𝑢1 + 𝑝 2𝑢2 + 𝑝 3𝑢3 +⋯ . ∞ 𝑝=0 (3.18) If 𝑝 → 1, then equation (3.18) corresponds to equation (3.15) and becomes the approximate solution of the form 𝑓 = lim 𝑝→1 𝑢 = ∑𝑢𝑖 ∞ 𝑖=0 . (3.19) It is well known that the series (3.18) is convergent for most of the cases and also the rate of convergence is dependent on 𝐿(𝑢). 3.2.2 The Variational Homotopy Perturbation Method To convey the basic idea of the variational homotopy perturbation method, we consider the general differential equation: 𝐿𝑢 + 𝑁𝑢 = 𝑔(𝑥), (3.20) Where 𝐿 is a linear operator, 𝑁 a nonlinear operator, and 𝑔(𝑥) the forcing term. According to variational iteration method (VIM), with the correct functional (3.21). 𝑢𝑛+1(𝑥) = 𝑢𝑛(𝑥) + ∫ 𝜆(𝜏){𝐿𝑢𝑛 +𝑁�̃�𝑛 − 𝑔(𝜏)}𝑑𝜏 𝑥 0 , (3.21) 25 Where 𝜆(𝜏) is a general Lagrange multiplier which can be identify optimally via the variational theory, the subscript 𝑛 denotes the 𝑛 th-order approximation and �̃�𝑛 is considered as a restricted variation, i.e. �̃�𝑛 = 0. Now, we apply the homotopy perturbation method: ∑𝑝𝑖𝑢𝑖 = 𝑢0 + 𝑝∫𝜆(𝜏) {𝑁 (∑𝑝𝑖𝑢�̃� ∞ 𝑖=0 )}𝑑𝜏 𝑥 0 − 𝑝∫𝜆(𝜏)𝑔(𝜏)𝑑𝜏 𝑥 0 ∞ 𝑖=0 , (3.22) Which is the variational homotopy perturbation method and it is formulated by the coupling of variational iteration method and Adomian’s polynomials. The embedding parameter 𝑝 ∈ [0,1] can be considered as an expanding parameter. The homotopy perturbation method uses the homotopy parameter 𝑝 as an expanding parameter to obtain: 𝑓 =∑𝑝𝑖𝑢𝑖 = 𝑢0 + 𝑝𝑢1 + 𝑝 2𝑢2 +⋯ ∞ 𝑖=0 (3.23) If 𝑝 → 1, then equation (3.23) becomes the approximate solution of the form: 𝑢 = lim 𝑝→1 𝑓 = 𝑢0 + 𝑢1 + 𝑢2 +⋯ (3.24) A comparison of like power of 𝑝 gives solutions of various orders. 26 3.2.3 The Variational Homotopy Perturbation Method (VHPM) for Lienard’s equation In this section, we consider Lienard equation (2.12) 𝑦′′(𝑥) + 𝑎𝑦(𝑥) + 𝑏𝑦3(𝑥) + 𝑐𝑦5(𝑥) = 0 with initial conditions: 𝑦(0) = 𝐶1 , 𝑦′(0) = 𝐶2 (3.25) using these initial conditions, we choose: 𝑦0(𝑥) = 𝑦(0) + 𝑦 ′(0)𝑥, (3.26) the initial approximation solution of equation (2.12). For solving equation (2.12) via VHPM, we consider: 𝐿(𝑦) = 𝑦′′, (3.27) 𝑁(𝑦) = 𝑎𝑦 + 𝑏𝑦3 + 𝑐𝑦5, (3.28) Where 𝐿 is a linear and 𝑁 is a nonlinear operators. According to the variational iteration method, we can construct a correct functional as follows: Where �̃�𝑛is considered as a restricted variation. Making the above functional stationary, the Lagrange multiplier can be determined as 𝜆(𝜏) = 𝜏 − 𝑥 𝑦𝑛+1(𝑥) = 𝑦𝑛(𝑥) + ∫ 𝜆(𝜏){𝑦𝑛𝜏𝜏 + 𝑎�̃�𝑛 + 𝑏�̃�𝑛 3 + 𝑐�̃�𝑛 5}𝑑𝜏, 𝑥 0 (3.29) 27 which yields the following iteration formula: 𝑦𝑛+1(𝑥) = 𝑦𝑛(𝑥) + ∫ (𝜏 − 𝑥){𝑦𝑛𝜏𝜏 + 𝑎𝑦𝑛 + 𝑏𝑦𝑛 3 + 𝑐𝑦𝑛 5}𝑑𝜏, 𝑥 0 (3.30) Applying the variational homotopy perturbation method, we have: 𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 +⋯ = 𝐶1 + 𝐶2𝑥 + 𝑝∫ (𝜏 − 𝑥)𝑎(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯)𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥)𝑏(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 3𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥)𝑐(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 5𝑑𝜏 𝑥 0 . (3.31) Comparing the coefficient of like powers of 𝑝 we have: 𝑝0: 𝑦0(𝑥) = 𝐶1 + 𝐶2𝑥 , 𝑝1: 𝑦1(𝑥) = ∫(𝜏 − 𝑥)𝑎 𝑦0𝑑𝜏 𝑥 0 +∫(𝜏 − 𝑥) 𝑏 𝑦0 3 𝑑𝜏 + ∫(𝜏 − 𝑥) 𝑐 𝑦0 5𝑑𝜏 , 𝑥 0 𝑥 0 ⋮ (3.32) So we obtain the components which constitute 𝑦(𝑥), thus we will have: 𝑦(𝑥) = 𝑦0 + 𝑦1 + 𝑦2 +⋯. (3.33) 28 For later numerical computation, we let the expression: 𝜑𝑛 =∑𝑦𝑖(𝑥) 𝑛 𝑖=0 (3.34) to denote the n-term approximation to 𝑦(𝑥). 3.3 The Adomian Decomposition Method (ADM) For a clear overview of ADM, consider a differential equation 𝐹(𝑢(𝑡)) = 𝑔(𝑡), (3.35) Where 𝐹 represents a general nonlinear ordinary or partial differential operator including both linear and nonlinear terms. Linear terms are decomposed into 𝐿 + 𝑅, where 𝐿 is invertible and taken as the highest order derivative, while 𝑅 is the remainder of the linear operator. Thus equation (3.35) may be written as: 𝐿𝑢 + 𝑁𝑢 + 𝑅𝑢 = 𝑔, (3.36) Where 𝑁 represents the nonlinear terms. Solving for 𝐿𝑢, we get: 𝐿𝑢 = 𝑔 − 𝑁𝑢 − 𝑅𝑢. (3.37) Since 𝐿 is taken as the highest order derivative, if 𝐿 = 𝑑2 𝑑𝑥2⁄ , Adomian defines 𝐿−1as the two-fold integration operator from 0 to 𝑥, i.e. 29 𝐿−1(. ) = ∫ ∫ (. ) 𝑑𝑥 𝑥 0 𝑑𝑥 𝑥 0 . (3.38) Applying operator (3.38) on both sides of equation (3.37) we have: 𝐿−1𝐿 𝑢 = 𝐿−1𝑔 − 𝐿−1𝑁𝑢 − 𝐿−1𝑅𝑢. (3.39) The decomposition method represents the solution 𝑢(𝑥, 𝑡) as a series of the form: 𝑢(𝑥, 𝑡) = ∑𝑢𝑛(𝑥, 𝑡). ∞ 𝑛=0 (3.40) The nonlinear term 𝑁 is decomposed as: 𝑁(𝑢) = ∑𝐴𝑛 ∞ 𝑛=0 . (3.41) Rewrite equation (3.39) by substitute equation (3.40) and equation (3.41), we get: ∑𝑢𝑛(𝑥, 𝑡) ∞ 𝑛=0 = 𝛿0 + 𝐿 −1𝑔(𝑥) − 𝐿−1𝑅 ∑𝑢𝑛 ∞ 𝑛=0 − 𝐿−1∑𝐴𝑛 ∞ 𝑛=0 , (3.42) Where, 𝛿0 = { 𝑢(0), 𝑖𝑓 𝐿 = 𝑑 𝑑𝑥 , 𝑢(0) + 𝑥𝑢′(0), 𝑖𝑓 𝐿 = 𝑑2 𝑑𝑥2 , ⋮ 𝑢(0) + 𝑥𝑢′(0) + 𝑥2 2! 𝑢′′(0) + ⋯+ 𝑥𝑛 𝑛! 𝑢(𝑛)(0), 𝑖𝑓 𝐿 = 𝑑𝑛+1 𝑑𝑥𝑛+1 . (3.43) 30 Therefore, the zeroth component is usually taken to be all terms arise from the initial conditions, i.e., 𝑢0 = 𝛿0 + 𝐿 −1𝑔(𝑥), (3.44) The remaining components 𝑢𝑛, 𝑛 ≥ 1, can be completely determined such that each term is computed by using the previous term. Since 𝑢0 is known, 𝑢𝑛 = −𝐿 −1𝑅𝑢𝑛−1 − 𝐿 −1𝐴𝑛−1 , 𝑛 ≥ 1 (3.45) The nonlinear term 𝑁(𝑢) is decomposed into an infinite series of polynomials that has the form 𝑁(𝑢) = ∑𝐴𝑛(𝑢0, 𝑢1, ⋯ , 𝑢𝑛) ∞ 𝑛=0 . (3.46) The components 𝐴𝑛 are called the Adomian polynomials, these polynomials can be calculated for all forms of nonlinearity according to specific algorithms constructed by Adomian. For this specific nonlinearity, we use the general formula for 𝐴𝑛 polynomials as 𝐴𝑛 = 1 𝑛! 𝑑𝑛 𝑑𝜗𝑛 [𝑁 (∑𝜗𝑘 ∞ 𝑘=0 𝑢𝑘)] 𝜗=0 , 𝑛 = 0, 1, 2, …. (3.47) 31 The first few Adomian polynomials for the nonlinearity 𝑁(𝑢) are 𝐴0 = 𝑁(𝑢0), 𝐴1 = 𝑢1𝑁 ′(𝑢0), 𝐴2 = 𝑢2𝑁 ′(𝑢0) + 1 2! 𝑢1 2𝑁′′(𝑢2), 𝐴3 = 𝑢3𝑁 ′(𝑢0) + 𝑢1𝑢2𝑁 ′′(𝑢0) + 1 3! 𝑢1 3𝑁(3)(𝑢0), ⋮ (3.48) So, the solution for 𝑛-terms approximation is 𝜙𝑛 = ∑𝑢𝑘 𝑛−1 𝑘=0 , (3.49) Where 𝑢(𝑥, 𝑡) = lim 𝑛→∞ 𝜙𝑛(𝑥, 𝑡) = ∑𝑢𝑘(𝑥, 𝑡) ∞ 𝑘=0 . (3.50) 3.3.1 The Adomian Decomposition Method (ADM) for Solving Lienard Equation Now we will obtain analytical solution for the Lienard equation using the ADM. First we write the Lienard equation in the standard operator form 𝐿𝑦 = −(𝑎𝑦 + 𝑏𝑦3 + 𝑐𝑦5), (3.51) 32 where the notation 𝐿 = 𝑑2 𝑑𝑥2 symbolizes the linear operator. The inverse operator 𝐿−1 exists. Thus, applying the inverse operator 𝐿−1 to (3.51) yields 𝐿−1𝐿𝑦 = −𝐿−1(𝑎𝑦 + 𝑏𝑦3 + 𝑐𝑦5). (3.52) Therefore, it follows that 𝑦(𝑥) = 𝑦(0) + 𝑥𝑦′(0) − 𝐿−1(𝑎𝑦 + 𝑏𝑦3 + 𝑐𝑦5). (3.53) The solution 𝑦(𝑥) of equation (2.12) represents as a series using the decompose method, we get 𝑦(𝑥) = ∑𝑦𝑖(𝑥) ∞ 𝑖=0 . (3.54) This yields to 𝑦0 = 𝑦(0) + 𝑥 𝑦 ′(0), 𝑦𝑖 = −𝐿 −1[𝑎𝑦𝑖−1 + 𝑏 𝐴𝑖−1 + 𝑐𝐵𝑖−1], 𝑖 ≥ 1, (3.55) where 𝑁(𝑦) = ∑ 𝐴𝑖(𝑦0, 𝑦1, … , 𝑦𝑖) ∞ 𝑖=0 . And the components 𝐴𝑖 are the Adomian polynomials, and are given by the general formula 𝐴𝑖 = 1 𝑖! 𝑑𝑖 𝑑𝜗𝑖 [𝑁 (∑𝜗𝑘𝑦𝑘 ∞ 𝑘=0 )] 𝜗=0 , 𝑖 = 0, 1, 2, …. (3.56) The first few Adomian polynomial for the nonlinearity 𝑁(𝑦) are 33 𝐴0 = 𝑁(𝑦0), 𝐴1 = 𝑦1𝑁 ′(𝑦0), 𝐴2 = 𝑦2𝑁 ′(𝑦0) + 1 2! 𝑦1 2𝑁′′(𝑦2), 𝐴3 = 𝑦3𝑁 ′(𝑦0) + 𝑦1𝑦2𝑁 ′′(𝑦0) + 1 3! 𝑦1 3𝑁(3)(𝑦0), ⋮ (3.57) Finally an 𝑛-term approximate solution is given by 𝜙𝑛 = ∑𝑦𝑖 𝑛−1 𝑖=0 (3.58) And the exact solution is 𝑦(𝑥) = lim 𝑛→∞ 𝜙𝑛. (3.59) 34 Chapter Four 35 Chapter 4 Numerical Examples on Lienard Equation In this chapter we attempt to solve some numerical examples on Lienard equation using the aforementioned methods. The implementation is carried out by using MAPLE and MATLAB software. 3.1 Example (1): Consider the Lienard equation: 𝑦′′(𝑥) + 𝑎𝑦(𝑥) + 𝑏𝑦3(𝑥) + 𝑐𝑦5(𝑥) = 0 With the following initial conditions: 𝑦(0) = √ −2𝑎 𝑏 , 𝑦′(0) = − 𝑎 √−𝑎 𝑏√ −2𝑎 𝑏 . (4.1) The exact solution of 𝑦(𝑥) in a closed form is readily found to be [9]: 𝑦1(𝑥) = √ −2𝑎 𝑏 (1 + 𝑡𝑎𝑛ℎ√−𝑎𝑥) . (4.2) Now, we will present the numerical solutions by the pervious methods. 36 4.1.1 Solving Example (1) Using Variational Iteration Method (VIM): To solve example (1) by Variational Iteration Method (VIM), we will use the iteration formula (3.9) which is: 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ (𝜏 − 𝑥){𝑦𝑖 ′′ + 𝑎𝑦𝑖 + 𝑏𝑦𝑖 3 + 𝑐𝑦𝑖 5}𝑑𝜏, 𝑥 0 With initial approximation 𝑦0 𝑉𝐼𝑀(𝑥) = 𝑦(0) + 𝑦 ′(0)𝑥 , i.e. 𝑦0 𝑉𝐼𝑀(𝑥) = √ −2𝑎 𝑏 − 𝑎 √−𝑎 𝑏√ −2𝑎 𝑏 𝑥 , (4.3) For solving 𝑦1 𝑉𝐼𝑀(𝑥) we use: 𝑦1 𝑉𝐼𝑀(𝑥) = 𝑦0(𝑥) + ∫ (𝜏 − 𝑥){𝑎𝑦0(𝜏) + 𝑏𝑦0 3(𝜏) + 𝑐𝑦0 5(𝜏)}𝑑𝜏, 𝑥 0 (4.4) Then, using MAPLE software, we obtain: 𝑦1 𝑉𝐼𝑀(𝑥) = 1 1680 1 𝑏3√− 𝑎 𝑏 (𝑎 (5(−𝑎) 9 2𝑐𝑥7 + 70𝑎4𝑐𝑥6 + 420(−𝑎) 7 2𝑐𝑥5 + 21(−𝑎) 5 2𝑏2𝑥5 − 1400𝑎3𝑐𝑥4 + 210𝑎2𝑏2𝑥4 + 2800(−𝑎) 5 2𝑐𝑥3 + 700(−𝑎) 3 2𝑏2𝑥3 + 3360𝑎2𝑐𝑥2 − 840𝑎𝑏2𝑥2 − 840√−𝑎𝑏2𝑥 − 1680𝑏2)√2), (4.5) 37 𝑦2 𝑉𝐼𝑀(𝑥) = 1 51073776977564870207078400000 𝑎√2 𝑏13√− 𝑎 𝑏 × (− 675367942093378058004480000 (−𝑎) 9 2 𝑏12𝑥9 − 21280740407318695919616000000 (−𝑎) 3 2 𝑏12 𝑥3 + 175996022308001841254400000(−𝑎) 11 2 𝑏12𝑥11 + 12342829436244843633377820000(−𝑎) 5 2 𝑏12𝑥5 + 12723623726943460531200000 (−𝑎) 39 2 𝑐6 𝑥27 + 353775758063888160000000 (−𝑎) 41 2 𝑐6𝑥29 + 5328568738092600000000(−𝑎) 43 2 𝑐6𝑥31 + 37845764487531000000(−𝑎) 45 2 𝑐6𝑥33 + 95409490304700000 (−𝑎) 47 2 𝑐6𝑥35 −⋯). (4.6) And so on using the same iteration formula (3.9). Here, using MAPLE software, we will calculate the errors between the exact solution 𝑦1(𝑥) and 𝑦𝑖 𝑉𝐼𝑀(𝑥)(𝑖 = 1,2), where 𝑎 = −1, 𝑏 = 4 , 𝑐 = −3 and for 𝑥 from 0.1 to 1. Table (4.1) contains both the exact and the first approximate solutions using the VIM. Table (4.2) contains both the exact and the second approximate solutions using VIM. 38 Moreover, Figure (4.1) compares both the exact and the first approximate solutions using the VIM, Figure (4.2) compares both the exact and the second approximate solutions using the VIM, and Figure (4.3) compares the absolute errors for both the first and the second approximate solutions using VIM. Table 4.1: The exact and first approximate solutions using VIM for example (1). 𝑥 𝑦1(𝑥) 𝑦1 𝑉𝐼𝑀(𝑥) |𝑦1(𝑥) − 𝑦1 𝑉𝐼𝑀(𝑥)| 0.1 0.7415079210 0.7415070380 8.830100 𝑒 − 7 0.2 0.7737490935 0.7737361617 1.293180 𝑒 − 5 0.3 0.8035274145 0.8034712744 5.614010 𝑒 − 5 0.4 0.8306470255 0.8305098173 1.372082 𝑒 − 4 0.5 0.8550196365 0.8548093246 2.103119 𝑒 − 4 0.6 0.8766554530 0.8765317729 1.236801 𝑒 − 4 0.7 0.8956471895 0.8960927030 4.455135 𝑒 − 4 0.8 0.9121504180 0.9142153132 2.064895 𝑒 − 3 0.9 0.9263632845 0.9319897243 5.626439 𝑒 − 3 1.0 0.9385079000 0.9509376127 1.242971 𝑒 − 2 39 Figure 4.1: The exact and first approximate solutions using VIM for example (1). Table 4.2: The exact and second approximate solutions using VIM for example (1). 𝑥 𝑦1(𝑥) 𝑦2 𝑉𝐼𝑀(𝑥) |𝑦1(𝑥) − 𝑦2 𝑉𝐼𝑀(𝑥)| 0.1 0.7415079210 0.7415079213 3.000000 𝑒 − 10 0.2 0.7737490935 0.7737491112 1.770000 𝑒 − 8 0.3 0.8035274145 0.8035275618 1.473000 𝑒 − 7 0.4 0.8306470255 0.8306475597 5.342000 𝑒 − 7 0.5 0.8550196365 0.8550206925 1.056000 𝑒 − 6 0.6 0.8766554530 0.8766565406 1.087600 𝑒 − 6 0.7 0.8956471895 0.8956479590 7.695000 𝑒 − 7 0.8 0.9121504180 0.9121570218 6.603800 𝑒 − 6 0.9 0.9263632845 0.9264090448 4.576030 𝑒 − 5 1.0 0.9385079000 0.9387036507 1.957507 𝑒 − 4 40 Figure 4.2: The exact and second approximate solutions using VIM for example (1). Figure 4.3: The errors of first two approximate solutions using VIM for example (1). 41 4.1.2 Solving Example (1) Using Variational Homotopy Perturbation Method (VHPM): To solve example (1) by VHPM, we will use equation (3.31) which is: 𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 +⋯ = 𝑦(0) + 𝑦′(0)𝑥 + 𝑝∫ (𝜏 − 𝑥) 𝑎 (𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯)𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥) 𝑏 (𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 3𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥) 𝑐 (𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 5𝑑𝜏 𝑥 0 . (4.7) With initial approximation 𝑦0 𝑉𝐻𝑃𝑀(𝑥) = 𝑦(0) + 𝑦 ′(0) 𝑥 , i.e. 𝑦0 𝑉𝐻𝑃𝑀 = √ −2𝑎 𝑏 − 𝑎√−𝑎 𝑏√ −2𝑎 𝑏 𝑥 , (4.8) For solving 𝑦1 𝑉𝐻𝑃𝑀(𝑥) we use: 𝑦1(𝑥) = ∫(𝜏 − 𝑥) (𝑎𝑦0(𝜏) + 𝑏𝑦0 3(𝜏) + 𝑐𝑦0 5(𝜏))𝑑𝜏 𝑥 0 Then using MAPLE software, we have: (4.9) 42 𝑦1 𝑉𝐻𝑃𝑀 = − 1 1680 1 𝑏3√− 𝑎 𝑏 (𝑎2𝑥2 (5(−𝑎) 7 2𝑐𝑥5 − 70𝑐𝑎3𝑥4 + 420(−𝑎) 5 2𝑐𝑥3 + 21(−𝑎) 3 2𝑏2𝑥3 + 1400𝑎2𝑐𝑥2 − 210𝑎𝑏2𝑥2 + 2800(−𝑎) 3 2𝑐𝑥 + 700√−𝑎𝑏2𝑥 − 3360𝑐𝑎 + 840𝑏2)√2), (4.10) 𝑦2 𝑉𝐻𝑃𝑀 = 1 11531520 1 𝑏5√− 𝑎 𝑏 (𝑎3√2𝑥4(275(−𝑎) 13 2 𝑐2𝑥9 + 7150𝑎6𝑐2𝑥8 + 858(−𝑎) 11 2 𝑐2𝑥7 + 2106(−𝑎) 9 2 𝑏2𝑥7 − 629200𝑎5𝑐2𝑥6 + 46332𝑎4𝑏2𝑐𝑥6 + 3146000(−𝑎) 9 2𝑐2𝑥5 + 446160(−𝑎) 7 2𝑏2𝑐𝑥5 +⋯. (4.11) Via variational homotopy perturbation method we will use 𝜑𝑛 𝑉𝐻𝑃𝑀: 𝜑𝑛 𝑉𝐻𝑃𝑀 =∑𝑦𝑖 𝑉𝐻𝑃𝑀(𝑥) 𝑛 𝑖=0 (4.12) Again, using MAPLE software, we will compute the errors between the exact solution 𝑦1(𝑥) and 𝜑𝑛 𝑉𝐻𝑃𝑀(𝑥)(𝑛 = 1,2), where 𝑎 = −1, 𝑏 = 4, 𝑐 = −3, and 43 for 𝑥 from 0.1 to 1. Table (4.3) contains both the exact and the first approximate solutions using the VHPM. Table (4.4) contains both the exact and the second approximate solutions using VHPM. Moreover, Figure (4.4) compares both the exact and the first approximate solutions using the VHPM, Figure (4.5) compares both the exact and the second approximate solutions using the VHPM, and Figure (4.6) compares the absolute errors for both the first and the second approximate solutions using VHPM. Table 4.3: The exact and first approximate solutions using VHPM for example (1). 𝑥 𝑦1(𝑥) 𝜑1 𝑉𝐻𝑃𝑀(𝑥) |𝑦1(𝑥) − 𝜑1 𝑉𝐻𝑃𝑀(𝑥)| 0.1 0.7415079210 0.7415070380 8.830100 𝑒 − 7 0.2 0.7737490935 0.7737361617 1.293180 𝑒 − 5 0.3 0.8035274145 0.8034712744 5.614010 𝑒 − 5 0.4 0.8306470255 0.8305098173 1.372082 𝑒 − 4 0.5 0.8550196365 0.8548093246 2.103119 𝑒 − 4 0.6 0.8766554530 0.8765317729 1.236801 𝑒 − 4 0.7 0.8956471895 0.8960927030 4.455135 𝑒 − 4 0.8 0.9121504180 0.9142153131 2.064895 𝑒 − 3 0.9 0.9263632845 0.9319897243 5.626439 𝑒 − 3 1.0 0.9385079000 0.9509376127 1.242971 𝑒 − 2 44 Figure 4.4: The exact and first approximate solutions using VHPM for example (1). Table 4.4: The exact and second approximate solutions using VHPM for example (1). 𝑥 𝑦1(𝑥) 𝜑2 𝑉𝐻𝑃𝑀 (𝑥) |𝑦1(𝑥) − 𝜑2 𝑉𝐻𝑃𝑀(𝑥)| 0.1 0.7415079210 0.7415079206 4.00000000 𝑒 − 10 0.2 0.7737490935 0.7737490297 6.3800000 𝑒 − 8 0.3 0.8035274145 0.8035262803 1.1342000 𝑒 − 6 0.4 0.8306470255 0.8306380606 8.9649000 𝑒 − 6 0.5 0.8550196365 0.8549746275 4.5009000 𝑒 − 5 0.6 0.8766554530 0.8764872386 1.6821440 𝑒 − 4 0.7 0.8956471895 0.8951376804 5.0950910 𝑒 − 4 0.8 0.9121504180 0.9108348743 1.3155437 𝑒 − 3 0.9 0.9263632845 0.9233738123 2.9894722 𝑒 − 3 1.0 0.9385079000 0.9324041880 6.1037120 𝑒 − 3 45 Figure 4.5: The exact and second approximate solutions using VHPM for example (1). Figure 4.6: The errors of first two approximate solutions using VHPM for example (1). 46 4.1.3 Solving Example (1) Using Adomian Decomposition Method (ADM): To solve example (1) by ADM, we use the equation (3.56), which is: 𝐴𝑖 = 1 𝑖! 𝑑𝑖 𝑑𝜗𝑖 [𝑁 (∑𝜗𝑘𝑦𝑘 ∞ 𝑘=0 )] 𝜗=0 , 𝑖 = 0, 1, 2, …. (4.13) The Adomian polynomials 𝐴𝑖 are computed according to (3.55) with 𝑁(𝑦) = 𝑦𝑝 and this gives: 𝐴0 = 𝑦0 𝑝 , 𝐴1 = 𝑝 𝑦0 𝑝−1 𝑦1 , 𝐴2 = (𝑝 − 1)𝑝 𝑦0 𝑝−2 1 2! 𝑦1 2 + 𝑝 𝑦0 𝑝−1 𝑦2 , 𝐴3 = (𝑝 − 2)(𝑝 − 1)𝑝 𝑦0 𝑝−3 1 3! 𝑦1 3 + (𝑝 − 1)𝑝 𝑦0 𝑝−2 𝑦1 𝑦2 + 𝑝 𝑦0 𝑝−1 𝑦3 . (4.14) And so on the other polynomials can be constructed. The Adomian polynomials 𝐴𝑖 and 𝐵𝑖 are computed according to equation (3.56) with 𝑦3 and 𝑦5, respectively. Which implies: 47 𝐴0 = 𝑦0 3, 𝐴1 = 3 𝑦0 2 𝑦1, 𝐴2 = 3 𝑦0 𝑦1 2 + 3 𝑦0 2 𝑦2, 𝐴3 = 𝑦1 3 + 6𝑦0𝑦1𝑦2 + 3 𝑦0 2 𝑦3, 𝐴4 = 3 𝑦1 2 𝑦2 + 3 𝑦0 𝑦2 2 + 6 𝑦0 𝑦1 𝑦3 + 3 𝑦0 2 𝑦4, ⋮ 𝐵0 = 𝑦0 5, 𝐵1 = 5 𝑦0 4 𝑦1, 𝐵2 = 10 𝑦0 3 𝑦1 2 + 5 𝑦0 4 𝑦2, 𝐵3 = 10𝑦0 2𝑦1 3 + 20𝑦0 3𝑦1𝑦2 + 5𝑦0 4𝑦3, 𝐵4 = 5𝑦0𝑦1 4 + 30𝑦0 2𝑦1 2𝑦2 + 10𝑦0 3𝑦2 2 + 20 𝑦0 3𝑦1𝑦3 + 5𝑦0 4𝑦4, Then, using these formulas: 𝑦0 = 𝑦(0) + 𝑥 𝑦 ′(0), (4.15) 𝑦𝑖 𝐴𝐷𝑀 = −𝐿−1[𝑎𝑦𝑖−1 + 𝑏 𝐴𝑖−1 + 𝑐𝐵𝑖−1], 𝑖 ≥ 1, (4.16) we obtain the following: 𝑦0 𝐴𝐷𝑀 = √ −2𝑎 𝑏 − 𝑎√−𝑎 𝑏√ −2𝑎 𝑏 𝑥, (4.17) 𝑦1 𝐴𝐷𝑀 = −𝐿 −1[𝑎𝑦0 + 𝑏 𝐴0 + 𝑐𝐵0] = −𝐿 −1[𝑎𝑦0 + 𝑏𝑦0 3 + 𝑐𝑦0 5] Using MAPLE software we have: 48 𝑦1 𝐴𝐷𝑀 = 1 42 𝑐𝑎5(−𝑎)5/2𝑥7 𝑏5 (− 2𝑎 𝑏 ) 5/2 − 1 6 𝑐𝑎6𝑥6 𝑏4 (− 2𝑎 𝑏 ) 3/2 + 1 5 ( 1 4 𝑎3(−𝑎)3/2 𝑏2 (− 2𝑎 𝑏 ) 3/2 + 5 2 𝑐𝑎3(−𝑎)3/2 √− 2𝑎 𝑏 𝑏3 ) 𝑥5 + 1 4 ( 𝑎3 𝑏√− 2𝑎 𝑏 + 10 3 𝑐√− 2𝑎 𝑏 𝑎3 𝑏2 ) 𝑥4 + 1 3 ( 1 2 𝑎2√−𝑎 𝑏 √− 2𝑎 𝑏 + 3 2 √− 2𝑎 𝑏 𝑎 √−𝑎 + 5 2 𝑐 (− 2𝑎 𝑏 ) 3/2 𝑎 √−𝑎 𝑏 ) 𝑥3 + 1 2 (−𝑐 (− 2𝑎 𝑏 ) 5/2 − 𝑏 (− 2𝑎 𝑏 ) 3/2 − 𝑎√− 2𝑎 𝑏 )𝑥2, (4.18) 49 𝑦2 𝐴𝐷𝑀 = − 5 6552 𝑐2𝑎9(−𝑎) 9 2𝑥13 𝑏9 (− 2𝑎 𝑏 ) 9 2 + 5 504 𝑐2𝑎12𝑥12 𝑏8 (− 2𝑎 𝑏 ) 7 2 + 1 11 ( − 1 140 𝑎7(−𝑎) 7 2 𝑐 𝑏6 (− 2𝑎 𝑏 ) 7 2 − 1 14 𝑐2𝑎7(−𝑎) 7 2 𝑏7 (− 2𝑎 𝑏 ) 5 2 − 1 3 𝑐2 𝑎9(−𝑎) 3 2 𝑏7 (− 2𝑎 𝑏 ) 5 2 − 1 8 𝑐 𝑎4 ( 1 20 𝑎3(−𝑎) 3 2 𝑏2 (− 2𝑎 𝑏 ) 3 2 + 1 2 𝑐𝑎3(−𝑎) 3 2 √− 2𝑎 𝑏 𝑏3 ) 𝑏2 ) +⋯. (4.19) And so on. Via Adomian Decomposition Method ADM We will use 𝜙𝑛 𝐴𝐷𝑀: 𝜙𝑛 𝐴𝐷𝑀 = ∑𝑦𝑖 𝐴𝐷𝑀(𝑥) 𝑛−1 𝑖=0 (4.20) Using MAPLE software, we will compute the errors between the exact solution 𝑦1(𝑥) and 𝜙𝑛 𝐴𝐷𝑀(𝑥)(𝑛 = 1,2), where 𝑎 = −1, 𝑏 = 4, 𝑐 = −3, 50 and for 𝑥 from 0.1 to 1. Table (4.5) contains both the exact and the first approximate solutions using the ADM. Table (4.6) contains both the exact and the second approximate solutions using ADM. Moreover, Figure (4.7) compares both the exact and the first approximate solutions using the ADM, Figure (4.8) compares both the exact and the second approximate solutions using the ADM, and Figure (4.9) compares the absolute errors for both the first and the second approximate solutions using ADM. Table 4.5: The exact and first approximate solutions using ADM for example (1). 𝑥 𝑦1(𝑥) 𝜙2 𝐴𝐷𝑀(𝑥) |𝑦1(𝑥) − 𝜙2 𝐴𝐷𝑀(𝑥)| 0.1 0.7415079210 0.7415070380 8.830100 𝑒 − 7 0.2 0.7737490935 0.7737361617 1.293180 𝑒 − 5 0.3 0.8035274145 0.8034712744 5.614010 𝑒 − 5 0.4 0.8306470255 0.8305098173 1.372082 𝑒 − 4 0.5 0.8550196365 0.8548093246 2.103119 𝑒 − 4 0.6 0.8766554530 0.8765317729 1.236801 𝑒 − 4 0.7 0.8956471895 0.8960927030 4.455135 𝑒 − 4 0.8 0.9121504180 0.9142153131 2.064895 𝑒 − 3 0.9 0.9263632845 0.9319897243 5.626439 𝑒 − 3 1.0 0.9385079000 0.9509376127 1.242971 𝑒 − 2 51 Figure 4.7: The exact and first approximate solutions using ADM for example (1). Table 4.6: The exact and second approximate solutions using ADM for example (1). 𝑥 𝑦1(𝑥) 𝜙3 𝐴𝐷𝑀(𝑥) |𝑦1(𝑥) − 𝜙3 𝐴𝐷𝑀(𝑥)| 0.1 0.7415079210 0.7415079206 4.00000000 𝑒 − 10 0.2 0.7737490935 0.7737490297 6.3810000 𝑒 − 8 0.3 0.8035274145 0.8035262803 1.1342000 𝑒 − 6 0.4 0.8306470255 0.8306380606 8.9649000 𝑒 − 6 0.5 0.8550196365 0.8549746275 4.5009000 𝑒 − 5 0.6 0.8766554530 0.8764872386 1.6821440 𝑒 − 4 0.7 0.8956471895 0.8951376804 5.0950910 𝑒 − 4 0.8 0.9121504180 0.9108348743 1.3155437 𝑒 − 3 0.9 0.9263632845 0.9233738123 2.9894722 𝑒 − 3 1.0 0.9385079000 0.9324041880 6.1037120 𝑒 − 3 52 Figure 4.8: The exact and second approximate solutions using ADM for example (1). Figure 4.9: The errors of first two approximate solutions using ADM for example (1). 53 4.2 Example (2): Consider the Lienard equation 𝑦′′(𝑥) + 𝑎𝑦(𝑥) + 𝑏𝑦3(𝑥) + 𝑐𝑦5(𝑥) = 0 (4.21) With the following initial conditions: 𝑦(0) = √ 𝐾 2 + 𝐷 , 𝑦′(0) = 0 , (4.22) Where 𝐾 = 4 √ 3𝑎2 (3𝑏2 − 16𝑐𝑎) , 𝐷 = −1 + 𝑏√3 √(3𝑏2 − 16𝑐𝑎) . (4.23) The exact solution of 𝑦(𝑥) in a closed form is readily found to be [9]. 𝑦2(𝑥) = √ 𝐾 sech2(√−𝑎𝑥) 2 + 𝐷 sech2(√−𝑎𝑥) . (4.24) Now, we will present the numerical solutions by the same methods. 54 4.2.1 Solving Example (2) Using Variational Iteration Method (VIM): To solve example (2) by VIM, we will use the iteration formula (3.9) which is: 𝑦𝑖+1(𝑥) = 𝑦𝑖(𝑥) + ∫ (𝜏 − 𝑥){𝑦𝑖 ′′ + 𝑎𝑦𝑖 + 𝑏𝑦𝑖 3 + 𝑐𝑦𝑖 5}𝑑𝜏, 𝑥 0 with initial approximation 𝑦0 𝑉𝐼𝑀(𝑥) = 𝑦(0) + 𝑦 ′(0)𝑥 = √ 𝐾 2 + 𝐷 Then: 𝑦1(𝑥) = √ 𝐾 2 + 𝐷 +∫ (𝜏 − 𝑥){𝑎 √ 𝐾 2 + 𝐷 + 𝑏(√ 𝐾 2 + 𝐷 ) 3 𝑥 0 + 𝑐 (√ 𝐾 2 + 𝐷 ) 5 }𝑑𝜏 (4.25) Since 𝑦0 ′′ = 0. By MAPLE software, since 𝐾 = 4 √ 3𝑎2 (3𝑏2 − 16𝑐𝑎) , 𝐷 = −1 + 𝑏√3 √(3𝑏2 − 16𝑐𝑎) . (4.26) we can obtain: 55 𝑦1 𝑉𝐼𝑀(𝑥) = − 2 × 3 1 4 (√3𝑏2 − 16𝑐𝑎 + 𝑏√3) 2 × ( √ √− 𝑎2 −3𝑏2 + 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 + 𝑏√3 × (−6𝑏2 + 16𝑎𝑐 − 2√3 𝑏 √3𝑏2 − 16𝑐𝑎 + 3𝑎 𝑏2𝑥2 + 16 𝑎2𝑐𝑥2 + √3𝑎𝑏√3𝑏2 − 16𝑐𝑎 𝑥2 + 6𝑏2√3𝑏2 − 16𝑐𝑎√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑥2 − 32√3 𝑎𝑏𝑐 √− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑥2 + 6 √3 𝑏3√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑥2) ) , (4.27) 56 𝑦2 𝑉𝐼𝑀(𝑥) = 128 1155 3 1 4 (√3𝑏2 − 16𝑐𝑎 + 𝑏√3) 12 ( √ √− 𝑎2 −3𝑏2 + 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 √3 𝑏 + √3𝑏2 − 16𝑐𝑎 × (−26943840 𝑏12 𝑥2 √− 𝑎2 −3𝑏2 + 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 + 8981280 √3 𝑏11 √3𝑏2 − 16𝑐𝑎 + 26943840 𝑏12 − 4087480320 𝑎5 𝑏2𝑐5 + 302776320 𝑎6𝑐6 + 1437004800 𝑎6 𝑏4𝑐4𝑥4 − 1176215040 𝑎5 𝑏6𝑐3𝑥4 + 676589760 𝑎4𝑏8𝑐2𝑥4 − 170644320 𝑎3 𝑏10𝑐𝑥4 + 1916006400 𝑎5𝑏4𝑐4𝑥2 + 223534080 𝑎4 𝑏6𝑐3𝑥2 − 574801920 𝑎3𝑏8𝑐2𝑥2 + 161663040 𝑎2𝑏10𝑐 𝑥2 − 1561190400 𝑎7𝑏2𝑐5𝑥4 −⋯) ) (4.28) And so on. By MAPLE software we can find the errors between the exact solution 𝑦2(𝑥) and 𝑦𝑖 𝑉𝐼𝑀(𝑥)(𝑖 = 1,2), where 𝑎 = −1, 𝑏 = 4, 𝑐 = 3, and for 𝑥 from 0.1 to 1. Table (4.7) contains both the exact and the first approximate solutions using the VIM. Table (4.8) contains both the exact and the second approximate solutions using VIM. 57 Moreover, Figure (4.10) compares both the exact and the first approximate solutions using the VIM, Figure (4.11) compares both the exact and the second approximate solutions using the VIM, and Figure (4.12) compares the absolute errors for both the first and the second approximate solutions using VIM. Table 4.7: The exact and first approximate solutions using VIM for example (2). 𝑥 𝑦2(𝑥) 𝑦1 𝑉𝐼𝑀(𝑥) |𝑦2(𝑥) − 𝑦1 𝑉𝐼𝑀(𝑥)| 0.1 0.6398446064 0.6398241650 2.04220000 𝑒 − 5 0.2 0.6288354073 0.6285139014 3.21503000 𝑒 − 4 0.3 0.6112463738 0.6096634622 1.58292900 𝑒 − 3 0.4 0.5880909709 0.5832728472 4.81811100 𝑒 − 3 0.5 0.5605742637 0.5493420568 1.12322022 𝑒 − 2 0.6 0.5299502894 0.5078710904 2.20792054 𝑒 − 2 0.7 0.4974059662 0.4588599484 3.85460206 𝑒 − 2 0.8 0.4639841630 0.4023086306 6.16755290 𝑒 − 2 0.9 0.4305455844 0.3382171374 9.23284300 𝑒 − 2 1.0 0.3977613922 0.2665854682 1.31175920 𝑒 − 1 58 Figure 4.10: The exact and first approximate solutions using VIM for example (2). Table 4.8: The exact and second approximate solutions using VIM for example (2). 𝑥 𝑦2(𝑥) 𝑦2 𝑉𝐼𝑀(𝑥) |𝑦2(𝑥) − 𝑦2 𝑉𝐼𝑀(𝑥)| 0.1 0.6398446064 0.6398446508 4.465157742 𝑒 − 8 0.2 0.6288354073 0.6288381358 2.728060385 𝑒 − 6 0.3 0.6112463738 0.6112755461 2.917236392 𝑒 − 5 0.4 0.5880909709 0.5882412139 1.502428512 𝑒 − 4 0.5 0.5605742637 0.5610876183 5.133544401 𝑒 − 4 0.6 0.5299502894 0.5312931303 1.342840754 𝑒 − 3 0.7 0.4974059662 0.5003091367 2.903170498 𝑒 − 3 0.8 0.4639841630 0.4694137573 5.429594282 𝑒 − 3 0.9 0.4305455844 0.4395881052 9.042520650 𝑒 − 3 1.0 0.3977613922 0.4114279492 1.366655671 𝑒 − 2 59 Figure 4.11: The exact and second approximate solutions using VIM for example (2). Figure 4.12: The errors of first two approximate solutions using VIM for example (2). 60 4.2.2 Solving Example (2) Using Variational Homotopy Perturbation Method (VHPM): To solve example (2) by VHPM, we will use the equation (3.31) which is 𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 +⋯ = 𝑦(0) + 𝑦′(0)𝑥 + 𝑝∫ (𝜏 − 𝑥)𝑎(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯)𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥)𝑏(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 3𝑑𝜏 𝑥 0 + 𝑝∫ (𝜏 − 𝑥)𝑐(𝑦0 + 𝑝𝑦1 + 𝑝 2𝑦2 + 𝑝 3𝑦3 +⋯) 5𝑑𝜏 𝑥 0 . (4.29) we have: 𝑦0 𝑉𝐻𝑃𝑀 = 𝑦(0) + 𝑦 ′(0) 𝑥 = √ 𝐾 2 + 𝐷 , (4.30) and 𝑦1(𝑥) = ∫(𝜏 − 𝑥) (𝑎√ 𝐾 2 + 𝐷 + 𝑏(√ 𝐾 2 + 𝐷 ) 3𝑥 0 + 𝑐(√ 𝐾 2 + 𝐷 ) 5 )𝑑𝜏 (4.31) 61 since 𝐾 = 4 √ 3𝑎2 (3𝑏2 − 16𝑐𝑎) , 𝐷 = −1 + 𝑏√3 √(3𝑏2 − 16𝑐𝑎) . (4.32) using MAPLE software we obtain: 𝑦1 𝑉𝐻𝑃𝑀 = − 1 (√3𝑏2 − 16𝑐𝑎 + 𝑏√3) 2 ( 2 𝑥2 (3 𝑎 𝑏2 + 8𝑎2𝑐 + √3𝑏2 − 16𝑐𝑎 √3 𝑎 𝑏 + 6 √2√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑏3 − 32 √2 √− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑎𝑏𝑐 + 2√2 √3√− 𝑎2 −3𝑏2 + 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 𝑏2) × 2 1 4 √ √− 𝑎2 −3𝑏2 + 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 √3𝑏2 − 16𝑐𝑎 + 𝑏√3 ) , (4.33) 62 𝑦2 𝑉𝐻𝑃𝑀 = 2 3 1 (√3𝑏2 − 16𝑐𝑎 + 𝑏√3) 4( 𝑥 4(81 𝑎 𝑏2 − 96 𝑎2𝑏2𝑐 + 288 𝑎3𝑐2 + 27 √3𝑏2 − 16𝑎𝑐 √3 𝑎 𝑏3 + 40 √3𝑏2 − 16𝑐𝑎 √3 𝑎2𝑏 𝑐 + 72 √2√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑏5 − 96 √2√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑎 𝑏3𝑐 − 1536 √2√− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑎2𝑏 𝑐2 + 24 √2 √3 √3𝑏2 − 16𝑐𝑎 √− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑏4 + 32 √2 √3 √3𝑏2 − 16𝑐𝑎 √− 𝑎2 −3𝑏2 + 16𝑐𝑎 𝑎 𝑏2𝑐) (4.34) And so on. Again, in Variational Homotopy Perturbation method (VHPM) we will use 𝜑𝑛 𝑉𝐻𝑃𝑀 with the formula: 𝜑𝑛 𝑉𝐻𝑃𝑀 =∑𝑦𝑖 𝑉𝐻𝑃𝑀(𝑥) 𝑛 𝑖=0 (4.35) 63 And by MAPLE software we can find the errors between the exact solution 𝑦2(𝑥) and 𝜑𝑛 𝑉𝐻𝑃𝑀(𝑥)(𝑖 = 1,2), where 𝑎 = −1, 𝑏 = 4, 𝑐 = 3, and for 𝑥 from 0.1 to 1. Table (4.9) contains both the exact and the first approximate solutions using the VHPM. Table (4.10) contains both the exact and the second approximate solutions using VHPM. Moreover, Figure (4.13) compares both the exact and the first approximate solutions using the VHPM, Figure (4.14) compares both the exact and the second approximate solutions using the VHPM, and Figure (4.15) compares the absolute errors for both the first and the second approximate solutions using VHPM. Table 4.9: The exact and first approximate solutions using VHPM for example (2). 𝑥 𝑦2(𝑥) 𝜑1 𝑉𝐻𝑃𝑀(𝑥) |𝑦2(𝑥) − 𝜑1 𝑉𝐻𝑃𝑀(𝑥)| 0.1 0.6398446064 0.6398241651 2.0441200 𝑒 − 5 0.2 0.6288354073 0.6285139016 3.2150580 𝑒 − 4 0.3 0.6112463738 0.6096634623 1.5829115 𝑒 − 3 0.4 0.5880909709 0.5832728474 4.8181234 𝑒 − 3 0.5 0.5605742637 0.5493420567 1.1232207 𝑒 − 2 0.6 0.5299502894 0.5078710904 2.2079198 𝑒 − 2 0.7 0.4974059662 0.4588599484 3.8546017 𝑒 − 2 0.8 0.4639841630 0.4023086306 6.1675532 𝑒 − 2 0.9 0.4305455844 0.3382171372 9.2328447 𝑒 − 2 1.0 0.3977613922 0.26658546821 1.3117592 𝑒 − 1 64 Figure 4.13: The exact and first approximate solutions using VHPM for example (2). Table 4.10: The exact and second approximate solutions using VHPM for example (2). 𝑥 𝑦2(𝑥) 𝜑2 𝑉𝐻𝑃𝑀(𝑥) |𝑦2(𝑥) − 𝜑2 𝑉𝐻𝑃𝑀(𝑥)| 0.1 0.6398446064 0.6398447249 8.5000000 𝑒 − 8 0.2 0.6288354073 0.6288428621 7.3730000 𝑒 − 6 0.3 0.6112463738 0.6113288252 8.2378000 𝑒 − 5 0.4 0.5880909709 0.5885362164 4.4526200 𝑒 − 4 0.5 0.5605742637 0.5621920788 1.6177780 𝑒 − 3 0.6 0.5299502894 0.5345168963 4.5665420 𝑒 − 3 0.7 0.4974059662 0.5082245934 1.0818496 𝑒 − 2 0.8 0.4639841630 0.4865225358 2.2538370 𝑒 − 2 0.9 0.4305455844 0.4731115299 4.2565890 𝑒 − 2 1.0 0.3977613922 0.4721858226 7.4424390 𝑒 − 2 65 Figure 4.14: The exact and second approximate solutions using VHPM for example (2). Figure 4.15: The errors of first two approximate solutions using VHPM for example (2). 66 4.2.3 Solving Example (2) Using Adomian Decomposition Method (ADM): To solve example (2) by ADM, we apply equation (3.56) which is 𝐴𝑖 = 1 𝑖! 𝑑𝑖 𝑑𝜗𝑖 [𝑁 (∑𝜗𝑘𝑦𝑘 ∞ 𝑘=0 )] 𝜗=0 , 𝑖 = 0, 1, 2, …. (4.36) The Adomian polynomials 𝐴𝑖 are computed according to equation (3.55) with 𝑁(𝑦) = 𝑦𝑝 and this gives 𝐴0 = 𝑦0 𝑝 , 𝐴1 = 𝑝 𝑦0 𝑝−1 𝑦1 , 𝐴2 = (𝑝 − 1)𝑝 𝑦0 𝑝−2 1 2! 𝑦1 2 + 𝑝 𝑦0 𝑝−1 𝑦2 , 𝐴3 = (𝑝 − 2)(𝑝 − 1)𝑝 𝑦0 𝑝−3 1 3! 𝑦1 3 + (𝑝 − 1)𝑝 𝑦0 𝑝−2 𝑦1 𝑦2 + 𝑝 𝑦0 𝑝−1 𝑦3 . (4.37) The Adomian polynomials 𝐴𝑖 and 𝐵𝑖 are computed according to equation (3.56) with 𝑦3 and 𝑦5, respectively. Which is: 67 𝐴0 = 𝑦0 3, 𝐴1 = 3 𝑦0 2 𝑦1, 𝐴2 = 3 𝑦0 𝑦1 2 + 3 𝑦0 2 𝑦2, 𝐴3 = 𝑦1 3 + 6𝑦0𝑦1𝑦2 + 3 𝑦0 2 𝑦3, 𝐴4 = 3 𝑦1 2 𝑦2 + 3 𝑦0 𝑦2 2 + 6 𝑦0 𝑦1 𝑦3 + 3 𝑦0 2 𝑦4, ⋮ 𝐵0 = 𝑦0 5, 𝐵1 = 5 𝑦0 4 𝑦1, 𝐵2 = 10 𝑦0 3 𝑦1 2 + 5 𝑦0 4 𝑦2, 𝐵3 = 10𝑦0 2𝑦1 3 + 20𝑦0 3𝑦1𝑦2 + 5𝑦0 4𝑦3, 𝐵4 = 5𝑦0𝑦1 4 + 30𝑦0 2𝑦1 2𝑦2 + 10𝑦0 3𝑦2 2 + 20 𝑦0 3𝑦1𝑦3 + 5𝑦0 4𝑦4, Then, using these formulas: 𝑦0 = 𝑦(0) + 𝑦 ′(0) 𝑥, (4.38) 𝑦𝑖 𝐴𝐷𝑀 = −𝐿−1[𝑎𝑦𝑖−1 + 𝑏 𝐴𝑖−1 + 𝑐𝐵𝑖−1], 𝑖 ≥ 1, (4.39) we obtain the following: 𝑦0 𝐴𝐷𝑀 = √ 𝐾 2 + 𝐷 , (4.40) 𝑦1 𝐴𝐷𝑀 = −𝐿 −1[𝑎𝑦0 + 𝑏 𝐴0 + 𝑐𝐵0] = −𝐿−1 [𝑎√ 𝐾 2 + 𝐷 + 𝑏(√ 𝐾 2 + 𝐷 ) 3 + 𝑐 (√ 𝐾 2 + 𝐷 ) 5 ] (4.41) 68 Since 𝐾 = 4 √ 3𝑎2 (3𝑏2 − 16𝑐𝑎) , 𝐷 = −1 + 𝑏√3 √(3𝑏2 − 16𝑐𝑎) . (4.42) using MAPLE software we have: 𝑦1 𝐴𝐷𝑀 = − 1 (𝑏 √3 + √−16 𝑎 𝑐 + 3 𝑏2) 2 ( 2 𝑥2 ( 3 𝑎 𝑏2 + 16 𝑎2 𝑐 + 6 √−16 𝑎 𝑐 + 3 𝑏2 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 + √3 √−16 𝑎 𝑐 + 3 𝑏2 𝑎 𝑏 + 6 √3 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑏3 − 32 √3 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑎 𝑏 𝑐 ) × 3 1 4 √ √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 √−16 𝑎 𝑐 + 3 𝑏2 √3 𝑏 + √−16 𝑎 𝑐 + 3 𝑏2 ) , (4.43) 69 𝑦2 𝐴𝐷𝑀 = 1 3 1 (√3 𝑏 + √−16 𝑎 𝑐 + 3 𝑏2 ) 4 ( 𝑥4 (3 𝑎 𝑏2 + 112 𝑎2 + 18√−16 𝑎 𝑐 + 3 𝑏2 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑏2 + √3 √−16 𝑎 𝑐 + 3 𝑏2 𝑎 𝑏 + 18 √3 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑏3 − 96 √3√− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑎 𝑏 𝑐 ) × 3 1 4 √ √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 √−16 𝑎 𝑐 + 3 𝑏2 √3 𝑏 + √−16 𝑎 𝑐 + 3 𝑏2 (3 𝑎 𝑏2 + 16 𝑎2𝑐 + 6 √−16 𝑎 𝑐 + 3 𝑏2 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑏2 + √3 √−16 𝑎 𝑐 + 3 𝑏2 𝑎 𝑏 + 6 √3 √− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑏3 − 32 √3√− 𝑎2 16 𝑎 𝑐 − 3 𝑏2 𝑎 𝑏 𝑐) ) . (4.44) 70 And so on by the same manner. Another time, in Adomian Decomposition Method ADM we will use 𝜙𝑛 𝐴𝐷𝑀 with the formula: 𝜙𝑛 𝐴𝐷𝑀 = ∑𝑦𝑖 𝐴𝐷𝑀 𝑛−1 𝑖=0 (4.45) And by MAPLE software we can find the errors between the exact solution 𝑦2(𝑥) and 𝜙𝑛 𝐴𝐷𝑀(𝑥)(𝑖 = 1,2), where 𝑎 = −1, 𝑏 = 4, 𝑐 = 3, and for 𝑥 from 0.1 to 1. Table (4.11) contains both the exact and the first approximate solutions using the ADM. Table (4.12) contains both the exact and the second approximate solutions using ADM. Moreover, Figure (4.16) compares both the exact and the first approximate solutions using the ADM, Figure (4.17) compares both the exact and the second approximate solutions using the ADM, and Figure (4.18) compares the absolute errors for both the first and the second approximate solutions using ADM. Table 4.11: The exact and first approximate solutions using ADM for example (2). 𝑥 𝑦2(𝑥) 𝜙2 𝐴𝐷𝑀(𝑥) |𝑦2(𝑥) − 𝜙2 𝐴𝐷𝑀(𝑥)| 0.1 0.6398446064 0.6398241651 2.0441200 𝑒 − 5 0.2 0.6288354073 0.6285139016 3.2150580 𝑒 − 4 0.3 0.6112463738 0.6096634623 1.5829115 𝑒 − 3 0.4 0.5880909709 0.5832728474 4.8181234 𝑒 − 3 0.5 0.5605742637 0.5493420567 1.1232207 𝑒 − 2 0.6 0.5299502894 0.5078710904 2.2079198 𝑒 − 2 0.7 0.4974059662 0.4588599484 3.8546017 𝑒 − 2 0.8 0.4639841630 0.4023086306 6.1675532 𝑒 − 2 0.9 0.4305455844 0.3382171372 9.2328447 𝑒 − 2 1.0 0.3977613922 0.26658546821 1.3117592 𝑒 − 1 71 Figure 4.16: The exact and first approximate solutions using ADM for example (2). Table 4.12: The exact and second approximate solutions using ADM for example (2). 𝑥 𝑦2(𝑥) 𝜙3 𝐴𝐷𝑀(𝑥) |𝑦2(𝑥) − 𝜙3 𝐴𝐷𝑀(𝑥)| 0.1 0.6398446064 0.6398447249 8.5000000 𝑒 − 8 0.2 0.6288354073 0.6288428621 7.3730000 𝑒 − 6 0.3 0.6112463738 0.6113288252 8.2378000 𝑒 − 5 0.4 0.5880909709 0.5885362164 4.4526200 𝑒 − 4 0.5 0.5605742637 0.5621920788 1.6177780 𝑒 − 3 0.6 0.5299502894 0.5345168963 4.5665420 𝑒 − 3 0.7 0.4974059662 0.5082245934 1.0818496 𝑒 − 2 0.8 0.4639841630 0.4865225358 2.2538370 𝑒 − 2 0.9 0.4305455844 0.4731115299 4.2565890 𝑒 − 2 1.0 0.3977613922 0.4721858226 7.4424390 𝑒 − 2 72 Figure 4.17: The exact and second approximate solutions using ADM for example (2). Figure 4.18: The errors of the first two approximate solutions using ADM for example (2). 73 Conclusion We have solved the Lienard equation using several numerical methods. These include the Variational Iteration Method, the Variational Homotopy Perturbation Method, and the Adomian Decomposition Method. Numerical methods were implemented in the form of algorithms to solve some numerical examples. The results show clearly that the Variational Iteration method is more efficient than the counterparts. 74 References [1] A. Abdelrazec, "Adomian Decomposition Method: Convergence Analysis and Numerical Approximations", MSc Thesis (Mathematics), McMaster University, Hamilton, Ontario, Canada, 2008. [2] A. Aghajane and A. Moradifam, "Some Sufficient Conditions for the Intersection with the Vertical Isocline in the Lienard Plane", Applied Mathematics Letters, vol. 19, pp. 491-497, 2006. [3] A. Aghajani and A. Moradifam, "Oscillation of Solutions of Second- Order Nonlinear Differential Equations of Euler Type", Journal of Mathematical Analysis and Applications, vol. 326, pp. 1076-1089, 2007. [4] M. N. Ahmadabadi, F. M. Maalek and M. Arab, "Application of He's Variational Iteration Method for Lienard Equation", World Applied Sciences Journal, vol. 7, pp. 1077-1079, 2009. [5] M. AhmadAbadi and F. Ghaini, "An Adomian Decomposition Method for Solving Lienard Equation in General Form", ANZIAM Journal, vol. 51, pp. 302-308, 2009. [6] A. AL-Saif and T. Hattim, "Variational Iteration Method for Solving Some Models of Nonlinear Partial Differential Equations", International Journal of Pure and Applied Sciences and Technology, vol. 4, pp. 30-40, 2011. 75 [7] B. S. Attili, "The Adomian Decomposition Method for Solving the Boussinesq Equation Arising in Water Wave Propagation", Numerical Methods for Partial Differential Equations, vol. 22, pp. 1337-1347, 2006. [8] Z. Dahmani, A. Anber and Y. Bouraoui, "Analytic Approximate Aolutions to the Coupled Lotka-Volterra Equation with Fractional Derivatives", International Journal of Nonlinear Science, vol. 9, pp. 276-284, 2010. [9] Z. Feng, "On Explicit Exact Solutions for the Lienard Equation and its Application", Physics Letters A, vol. 293, pp. 50-56, 2002. [10] F. Fernandez, "On the Variational Homotopy Perturbation Method for Nonlinear Oscillators", Journal of Mathematical Physics, vol. 1, pp. 1- 3, 2012. [11] R. Figueiredo, "On the Existence of N-Periodic Solutions of Lienard’s Equation", Nonlinear Analysis: Theory Methods and Applications, vol. 7, pp. 483-499, 1983. [12] V. Gaiko, "The Applied Geometry of a General Lienard Polynomial System", Applied Mathematics Letters, vol. 25, pp. 2328-2331, 2012. [13] J. Gomez-Aguilar, L. Torres, J.Reyes, D. Baleanu, H. Yepez-Martinez and I. Sosa, "Fractional Liénard Type Model of a Pipeline within the Fractional Derivative without Singular Kernel", Advances in Difference Equations, vol. 2016, p. 173, 2016. 76 [14] T. Hara and T. Yoneyama, "On the Global Center of Generalized Lienard Equation and its Application to Stability Problems", Funkcialaj Ekvacioj, vol. 28, pp. 171-192, 1985. [15] T. Harko, F. S. Lobo and M. K. Mak, "A Class of Exact Solutions of the Lienard Type Ordinary Non-Linear Differential Equation", Journal of Engineering Mathematics, vol. 89, no. 1, pp. 193-205, 2014. [16] M. Heydari, M. Hooshmandasl and F. Ghaini, "A Good Approximate Solution for Lienard Equation in a Large Interval using Block Pulse Functions", Journal of Mathematical Extension, vol. 7, pp. 17-32, 2013. [17] M. Hosseini and H. Nasabzadeh, "On the Convergence of Adomian Decomposition Method", Applied Mathematics and Computation, vol. 182, pp. 536-543, 2006. [18] Y. Huang and H. Liu, "A New Modification of the Variational Iteration Method for Van Der Pol Equations", Applied Mathematical Modelling, vol. 37, pp. 8118-8130, 2013. [19] D. Kaya and S. El-Sayed, "A Numerical Implementation of the Decomposition Method for the Lienard Equation", Applied Mathematics and Computation, vol. 171, pp. 1095-1103, 2005. [20] D. Kong, "Explicit Exact Solutions for the Lienard Equation and its Applications", Physics Letters A, vol. 196, pp. 301-306, 1995. [21] N. Kudryashov and D. Sinelshchikov, "On the Criteria for Integrability of the Lienard Equation", Applied Mathematics Letters, vol. 57, pp. 114-120, 2016. 77 [22] S. LI, F. Liao and W. Xing, "Periodic Solutions for Lienard Differential Equation with Singularities", Electronic Journal of Differential Equations, vol. 2015, pp. 1-12, 2015. [23] S. Mailk, I. Qureshi, M. Amir and I. Haq, "Numerical Solution of Lienard Equation Using Hybrid Heuristic Computation", World Applied Sciences Journal, vol. 5, pp. 636-643, 2013. [24] A. B. Malinowska, T. Odzijewicz and D. Torres, "Advanced Methods in the Fractional Calculus of Variations", Springer Briefs in Applied Sciences and Technology, 2015. [25] Manjak N. H, A. K Mishra and Rakiya M. K .A, "Application of Adomian Decomposition Method for Solving a Class of Diffusion Equation", International Journal of Engineering Research & Technology, vol. 8, pp. 1213-1218, 2013. [26] M. Matinfar, S. Bahar and M. Ghasemi, "Solving the Lienard Equation by Differential Transform Method", World Journal of Modelling and Simulation, vol. 8, pp. 142-146, 2012. [27] M. Matinfar, M. Mahdavi and Z. Raeisy, "Exact and Numerical Solution of Lienard's Equation by the Variational Homotopy Perturbation Method", Journal of Information and Computing Science, vol. 6, pp. 73-80, 2011. [28] M. Matinfar, H. Hosseinzadeh and M. Ghanbari, "A Numerical Implementation of the Variational Iteration Method for the Lienard Equation", World Journal of Modelling and Simulation, vol. 4, pp. 205- 210, 2008. 78 [29] V. Mishra, S. Das, H. Jafari and S. Ongc, "Study of Fractional Order Van Der Pol Equation", Journal of King Saud University - Science, vol. 28, pp. 55-60, 2015. [30] H. Moreira, "Lienard-type Equations and the Epidemiology of Malaria", Ecological Modelling, vol. 60, pp. 139-150, 1992. [31] M. Noor and S. Mohyud-Din, "Variational Homotopy Perturbation Method for Solving Higher Dimensional Initial Boundary Value Problems", Mathematical Problems in Engineering, vol. 2008, pp. 11, 2008. [32] S. Perez-Gonzalez, J. Torregrosa and P. J. Torres, "Existence and Uniqueness of Limit Cycles for Generalized φ-Laplacian Lienard Equations", Journal of Mathematical Analysis and Applications, vol. 439, no. 2, pp. 745-765, 2016. [33] E. salehpour, H. Jafari and C. M. Khalique, "A Modified Variational Iteration Method for Solving Generalized Boussinesq Equation and Lienard Equation", International Journal of the Physical Sciences, vol. 6, pp. 5406-5411, 2011. [34] H. Singh, "Solution of Fractional Lienard Equation Using Chebyshev Operational Matrix Method", Nonlinear Science Letter A, vol. 8, pp. 397-404, 2017. [35] T. Slight, B. Romeira, L. Wang, J. Figueiredo, E. Wasige and C. Ironside, "A Lienard Oscillator Resonant Tunneling Diode-Laser Diode Hybrid Integrated Circuit", IEEE Journal of Quantum Electronics, vol. 44, pp. 1158-1163, 2008. 79 [36] J. Sugie, "On the Boundedness of Solutions of the Generalized Lienard Equation without the Signum Condition", Nonlinear Analysis: Theory Methods and Applications, vol. 11, pp. 1391-1397, 1987. [37] J. Sun, W. Wang and L. Wu, "A Note on “On Explicit Exact Solutions for the Lienard Equation and its Application” ", Physics Letters A, vol. 318, pp. 93-101, 2003. [38] M. Syam, "A Numerical Solution of Fractional Lienard’s Equation by Using the Residual Power Series Method", Mathematics, vol. 6, 2018. [39] V. Verma, A. Prakash, D. Kumar and J. Singh, "Numerical Study of Fractional Model of Multi-Dimensional Dispersive Partial Differential Equation", Journal of Ocean Engineering and Science, vol. 4, no. 4, pp. 338-351, 2019. [40] G. Villari, "On the existence of Periodic Solutions for Lienard’s Equation", Nonlinear Analysis: Theory Methods and Applications, vol. 7, pp. 71-78, 1983. [41] A. Wazwaz, "A New Algorithm for Calculating Adomian Polynomials for Nonlinear Operators", Applied Mathematics and Computation, vol. 111, pp. 33-51, 2000. [42] H. Zeghdoudi, L. Bouchahed and R. Dridi, "A Complete Classification of Lienard Equation", European Journal of Pure and Applied Mathematics, vol. 6, pp. 126-136, 2013. 80 [43] Z. Zuo-Huan, "On the Nonexistence of Periodic Solutions for Lienard Equations Nonlinear Analysis: Theory Methods and Applications, vol. 16, pp. 101-110, 1991. 81 Appendix (1) Maple code for solving example 1: (1.1) Using Variational Iteration Method > restart; > with(PDEtools): > a:=-1;b:=4;c:=-3; > y[0]:=sqrt(-2*a/b)-a*sqrt(-a)/(b*sqrt(- 2*a/b))*t; > f:= proc (i) options operator, arrow; subs(t=x,y[i])+int((t- x)*((diff(subs(x=t,y[i]),t,t)+a*subs(x=t,y[i])+b*s ubs(x=t,y[i])^3+c*subs(x=t,y[i])^5)), t = 0 .. x); end proc; > m:=1; > for l from 0 by 1 to m do y[l+1]:=simplify(f(l)) end do; > Y[ex]:=sqrt(-2*a/b*(1+tanh(sqrt(-a)*x))); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(x=r,y[m+1])- subs(x=r,Y[ex]))) end do; 82 (1.2) Using Variational Homotopy Perturbation Method: > restart; > with(PDEtools): > a:=-1;b:=4;c:=-3; > y[0]:=sqrt(-2*a/b)-a*sqrt(-a)/(b*sqrt(- 2*a/b))*t; > A:= proc (i) options operator, arrow; subs(x=t,coeff(expand(p*(a*(sum(p^n*y[n],n=0..i- 1))+b*(sum(p^n*y[n],n=0..i- 1))^3+c*(sum(p^n*y[n],n=0..i-1))^5)),p^i)); end proc; > Y:= proc (l) options operator, arrow; simplify(int((t-x)*A(l), t = 0 .. x)); end proc; > m:=2; > for l from 1 by 1 to m do y[l]:=simplify(Y(l)) end do; > Phi[m]:=add( y[l],l=0..m); > Y[ex]:=sqrt(-2*a/b*(1+tanh(sqrt(-a)*x))); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(t=x,x=r,Phi[m])- subs(x=r,Y[ex]))) end do; 83 (1.3) Using Adomian Decomption Method: > restart; > with(PDEtools): > a:=-1;b:=4;c:=-3; > y[0]:=sqrt(-2*a/b)-a*sqrt(-a)/(b*sqrt(- 2*a/b))*x; > A := i -> 1/i!*subs(theta=0,diff((sum(theta^k*y[k],k=0..i))^ 3,theta$i)); > B := i -> 1/i!*subs(theta=0,diff((sum(theta^k*y[k],k=0..i))^ 5,theta$i)); > A(0):=y[0]^3;B(0):=y[0]^5; > f:= i -> int(-subs(x=t,a*y[i-1]+b*A(i-1)+c*B(i- 1)),[t=0..xi, xi=0..x]); > m:=2; > for l from 1 by 1 to m do y[l]:=(f(l)) end do; > Phi[m]:=add(y[l],l=0..m-1); > Y[ex]:=sqrt(-2*a/b*(1+tanh(sqrt(-a)*x))); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(x=r,Phi[m])- subs(x=r,Y[ex]))) end do; 84 (2) Maple code for solving example 2: (2.1) Using Variational Iteration Method > restart; > with(PDEtools): > a:=-1;b:=4;c:=3; > k:=4*sqrt(3*a^2/(3*b^2-16*c*a)); > d:=-1+b*sqrt(3)/sqrt(3*b^2-16*c*a); > y[0]:=simplify(sqrt(k/(2+d))); > f:= proc (i) options operator, arrow; subs(t=x,y[i])+int((t- x)*((diff(subs(x=t,y[i]),t,t)+a*subs(x=t,y[i])+b*s ubs(x=t,y[i])^3+c*subs(x=t,y[i])^5)), t = 0 .. x); end proc; > m:=1; > for l from 0 by 1 to m do y[l+1]:=simplify(f(l)) end do; > Y[ex]:=sqrt((k*(sech(sqrt(- a)*x))^2)/(2+d*(sech(sqrt(-a)*x))^2)); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(x=r,y[m+1])- subs(x=r,Y[ex]))) end do; 85 (2.2) Using Variational Homotopy Perturbation Method: > restart; > with(PDEtools): > a:=-1;b:=4;c:=3; > k:=4*sqrt(3*a^2/(3*b^2-16*c*a)); > d:=-1+b*sqrt(3)/sqrt(3*b^2-16*c*a); > y[0]:=simplify(sqrt(k/(2+d))); > A:= proc (i) options operator, arrow; subs(x=t,coeff(expand(p*(a*(sum(p^n*y[n],n=0..i- 1))+b*(sum(p^n*y[n],n=0..i- 1))^3+c*(sum(p^n*y[n],n=0..i-1))^5)),p^i)); end proc; > Y:= proc (l) options operator, arrow; simplify(int((t-x)*A(l), t = 0 .. x)); end proc; > m:=2; > for l from 1 by 1 to m do y[l]:=simplify(Y(l)) end do; > Phi[m]:=add( y[l],l=0..m); > Y[ex]:=sqrt((k*(sech(sqrt(- a)*x))^2)/(2+d*(sech(sqrt(-a)*x))^2)); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(t=x,x=r,Phi[m])- subs(x=r,Y[ex]))) end do; 86 (2.3) Using Adomian Decomption Method: > restart; > with(PDEtools): > a:=-1;b:=4;c:=3; > k:=4*sqrt(3*a^2/(3*b^2-16*c*a)); > d:=-1+b*sqrt(3)/sqrt(3*b^2-16*c*a); > y[0]:=simplify(sqrt(k/(2+d))); > A := i -> 1/i!*subs(theta=0,diff((sum(theta^k*y[k],k=0..i))^ 3,theta$i)); > B := i -> 1/i!*subs(theta=0,diff((sum(theta^k*y[k],k=0..i))^ 5,theta$i)); > A(0):=y[0]^3;B(0):=y[0]^5; > f:= i -> int(-subs(x=t,a*y[i-1]+b*A(i-1)+c*B(i- 1)),[t=0..xi, xi=0..x]); > m:=2; > for l from 1 by 1 to m do y[l]:=(f(l)) end do; > Phi[m]:=add(y[l],l=0..m-1); > Y[ex]:=sqrt((k*(sech(sqrt(- a)*x))^2)/(2+d*(sech(sqrt(-a)*x))^2)); > for r from 0.1 by 0.1 to 1 do Err[r]:=abs(simplify(subs(x=r,Phi[m])- subs(x=r,Y[ex]))) end do; جاح الوطنّيةامعة النج كلية الدراسات العليا عددية لمعادلة ليناردحلول إعداد هيام أحمد الدلو إشراف د. ناجي قطناني أ. لّية بك الرياضياتفي قدمت هذه األطروحة استكمااًل لمتطلبات الحصول على درجة الماجستير س، فلسطين.الدراسات العليا في جامعة الّنجاح الوطنّية في نابل م2020 ب عددية لمعادلة ليناردحلول إعداد هيام أحمد الدلو إشراف د. ناجي قطناني أ. الملخص ألنها تستخدم في مجموعة واسعة من التطبيقات معادلة ليناردفي هذه األطروحة ركزنا على حل وعلم االحياء وغيرها.، والهندسة ،الفيزيائية رق طمنا بشكل أساسي على ركزنا اهتما ،اسيات التي نحتاجهابعد أن قدمنا بعض التعاريف واألس اضطراب المتغير طريقة و ،(VIMطريقة التكرار المتغير). هذه الطرق هي: معادلة ليناردلحل حسابية (.ADMوطريقة أدوميان التفككية )، (VHPMالهوموتوبي ) يث رب الخاصة بها. حمع خصائص التقا حسابيةياضي لهذه الطرق الحيث سيتم عرض اإلطار الر من خالل بعض األمثلة العددية. حسابيةتم توضيح كفاءة هذه الطرق السي لحل سابيةحي واحدة من أقوى التقنيات اله التكرار المتغيرتظهر النتائج العددية بوضوح أن طريقة األخرى بناًء على األمثلة المستخدمة. حسابيةال بالمقارنة مع التقنيات معادلة لينارد