J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 Journal of the Nigerian Society of Physical Sciences An Order Four Continuous Numerical Method for Solving General Second Order Ordinary Differential Equations F. O. Obarhua∗, O. J. Adegboro Department of Mathematical Sciences, The Federal University of Technology, Akure, Nigeria Abstract Continuous hybrid methods are now recognized as efficient numerical methods for problems whose solutions have finite domains or cannot be solved analytically. In this work, the continuous hybrid numerical method for the solution of general second order initial value problems of ordinary differential equations is considered. The method of collocation of the differential system arising from the approximate solution to the problem is adopted using the power series as a basis function. The method is zero stable, consistent, convergent. It is suitable for both non-stiff and mildly-stiff problems and results were found to compete favorably with the existing methods in terms of accuracy. DOI:10.46481/jnsps.2021.150 Keywords: Numerical Scheme, Continuous hybrid method, Zero stability, Linear and nonlinear Article History : Received: 05 January 2021 Received in revised form: 31 January 2020 Accepted for publication: 04 February 2021 Published: 27 February 2021 c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: B. J. Falaye 1. Introduction We consider the second order Ordinary Differential Equa- tion (ODEs) y′′ = f (x, y, y′), y(µ) = ω0, y ′(µ) = ω1 (1) Equation (1) occur virtually every areas of physical or biologi- cal process in connection with numerous problems that are en- countered in various aspects of everyday life. It is well con- ceived that this type of equation can either be solved directly or solved by reducing to system of first order differential equa- tions before applying different methods available to solve the resulting system of first order ODEs Chan et al. [1], Gho- lamtabar and Parandin [2]. ∗Corresponding author tel. no: +2348062488066 Email address: obarhuafo@futa.edu.ng (F. O. Obarhua ) Various linear multistep methods with different order of ac- curacy have been developed for the solution of 1 which varies from discrete linear multistep methods to continuous linear mul- tistep methods. Lambert and Watsan, [3] reported that linear multistep methods generally are more efficient in terms of ac- curacy with weak stability properties for a given number of evaluation per step and suffered the disadvantage of requiring additional starting values and special procedures for changing step length h. It is also good to note that continuous linear mul- tistep methods have advantages over the discrete methods in such a way that they give better error estimation, they provide a simplified form of coefficients for further evaluation at different points, and they provide solutions at all interior points within the interval of integration Saravi and Mirrajei, [4], Kayode and Awoyemi, [5], Golbabai and Arabshahi [6]. Among the first methods developed are first derivative methods that are imple- mented in predictor-corrector mode, and Taylor series expan- 42 Obarhua Adegboro / J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 43 sion are adopted to provide the starting values. The identified setbacks of the predictor-corrector methods are; they are very costly to implement and reduced order of accuracy of the pre- dictors. Recently, authors have proposed different methods of higher order differential equations to improve on the existing setbacks. Such improved methods are Kayode and Adeyeye, [7, 8] and Kayode and Obarhua, [9, 10]. They independently proposed linear multistep methods of higher order of accura- cies and same order of main predictors and the correctors and hence improved significantly the accuracies of the methods. This work proposed an accurate continuous numerical hy- brid method for direct solution of initial value problems of ODEs. The derived method is capable to handle stiff, mildly stiff, non- linear and engineering problems modeled as a second order ini- tial value problem of ODEs. 2. Derivation of the Method We define the general power series approximation solution in the form y(x) = (c+i)−1∑ j=0 a j x j (2) y′(x) = (c+i)−1∑ j=1 ja j x j−1 (3) y′′(x) = (c+i)−1∑ j=2 j( j − 1)a j x j−2 (4) Equating (4) with (1) gives f (x, y, y′) = (c+i)−1∑ j=2 j( j − 1)a j x j−2 (5) Equation (2) is interpolated at xn+i, i = 2, 5 2 and (5) is collo- cated at xn+c, c = 0(1)3. Therefore, interpolation and collocation equation at the selected grid and offstep points give rise to system of equations which can be express in matrix form 1 xn+2h x2n+2h x 3 n+2h x 4 n+2h x 5 n+2h 1 xn+rh x2n+rh x 3 n+rh x 4 n+rh x 5 n+rh 0 0 2 6xn 12x2n 20x 3 n 0 0 2 6xn+h 12x2n+h 20x 3 n+h 0 0 2 6xn+2h 12x2n+2h 20x 3 n+2h 0 0 2 6xn+3h 12x2n+3h 20x 3 n+3h   a0 a1 a2 a3 a4 a5  =  yn+2 yn+r fn fn+1 fn+2 fn+3  (6) Gaussian elimination method is then applied to solve equation (6) to obtain the unknown coefficients a′j s which is then substi- tuted into (2). Continuous system is obtained after some alge- braic simplifications. Applying transformation t = 1h (x − xn+k−1) , k = 3, t = (0, 1] in Obarhua [11], the continuous coefficients are obtained as fol- lows α2 = − ( −rh + th + 2h h(r − 2) ) αr = th h(r − 2) β0 = h5 360 ( −3t5 + 10t3 + 8t + 3tr4 − 24tr3 + 62tr2 − 56tr h3 ) β1 = − h5 120 ( −3t5 − 5t4 + 20t3 − 72t + 3tr4 − 19tr3 + 22tr2 + 44tr h3 ) β2 = h5 120 ( −3t5 − 10t4 + 10t3 + 60t2 + 48t + 3t4r4 h3 −14t4r3 + 2t4r2 + 4tr h3 ) (7) The first derivatives of equation (7) are α′2 = − 1 (r − 2) α′r = 1 (r − 2) β′0 = h5 360 ( −15t4 + 30t2 + 3r4 − 24r3 + 62r3 − 56r + 8 h3 ) β′1 = − h5 120 ( −15t4 − 20t3 + 60t2 + 3r4 − 19r3 + 22r2 + 44r − 72 h3 ) β′2 = h5 120 ( −15t4 − 40t3 + 30t2 + 120t + 12t3r4 − 56t2r3 + 8t3r2 h3 +4r + 48 h3 ) β′3 = − h5 360 ( −15t4 − 60t3 − 60t2 + 3r4 − 9r3 + 8t3r2 + 16t3r + 8 h3 ) (8) Evaluating equation (7) and (8) at t = 1 yield the discrete order continuous numerical scheme yn+3 = − 1 360(r − 2) ( −60h2r5 fn+2 − 360yn+2 − 630h 2 fn+2 − 30h2 fn − 3h 2r5 fn+3 + 9h 2r5 fn+2 − 9h 2r5 fn+1 + 3h 2r5 fn + 15h2r4 fn+3 − 20h 2r3 fn+3 − 180h 2r3 fn+1 + 90h 2r3 fn+2 − 180h2r2 fn + 110h 2r3 fn − 30h 2r4 fn + 75h 2r4 fn+1 − 60h2r4 fn+2 + 360ryn+2 + 291h 2r fn+2 + 360yn+r − 360h 2 fn+1 + 38h2r fn+3 + 444h 2r fn+1 + 127h 2r fn ) (9) 43 Obarhua Adegboro / J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 44 its first derivative is given as y′n+3 = − 1 360h(r − 2) ( 9h2r5 fn+2 + 9h 2r5 fn+1 + 3h 2r5 fn + 90h 2r3 fn+2 − 60h2r4 fn+2 + 110h 2r3 fn − 180h 2r3 fn+1 − 180h 2r2 fn − 30h2r4 fn + 75h 2r4 fn+1 − 3h 2r5 fn+3 − 20h 2r3 fn+3 + 15h2r4 fn+3 − 254h 2 fn+3 − 858h 2 fn+2 − 282h 2 fn+1 − 46h2 fn + 360yn+r + 405h 2r fn+2 + 405h 2r fn+1 + 135h 2r fn+3 + 135h2r fn − 360yn+2 ) (10) The values of r is taken in the interval r ∈ (2, 3) to obtain a particular discrete hybrid method. For the purpose of testing the properties of equation (9), the value of r is taken to 52 to have yn+3 = 2yn+ 52 −yn+2+ h2 384 (33 fn+3+83 fn+2−25 fn+1+5 fn)(11) with its first derivative given by hy′n+3 = 2yn+ 52 − 2yn+2 + h 5760 (2047 fn+3 +3069 fn+2 − 999 fn+1 + 203 fn) (12) 3. Implementation of the Method (11) In order to implement the implicit linear one-point discrete scheme (11) and its derivative (12), the symmetric explicit schemes and their derivatives are also developed by the same procedure for the evaluation of yn+3 and y′n+3 contained in fn+3 in the main scheme (11) and its derivative (12). yn+3 = 2yn+ 52 −yn+1 + h2 240 (66 fn+ 52 −10 fn+2 +5 fn+1− fn)(13) and its first derivative as hy′n+3 = −2hyn+ 52 + 2hyn+1 + h2 3600 ( −4094 fn+ 52 + 1920 fn+2 − 655 fn+1 + 129 fn ) (14) Other explicit schemes were also generated to evaluate other starting values using Taylor series expansions to evaluate the values for yn+i, y′n+i as yn+i = yn+( jh)y ′ n+ ( jh)2 2! fn+ ( jh)3 3! ( ∂ fn ∂xn + y′n ∂ fn ∂yn + fn ∂ fn ∂y′n ) +o(h4)(15) and y′n+i = y ′ n+( jh) fn+ ( jh)2 2! ( ∂ fn ∂xn + y′n ∂ fn ∂yn + fn ∂ fn ∂y′n ) +o(h4)(16) 4. Stability Analysis 4.1. Region of Absolute Stability In other to investigate the periodic stability properties of the numerical methods for solving the initial-value problem equa- tion 1 and the interval of periodicity, Lambert and Watson [3] introduced the following scalar test problem as y′′ = −λ2y (17) Based on the theory developed in Lambert and Watson [3], when multistep method l∑ j=0 α jyn+ j = h 2 l∑ j=0 θ j fn+ j, (18) is applied to the scalar test equation (17), a difference equation of the form l∑ i=0 (αi + H 2θi)yn+i = 0 (19) is obtained, where H = ph, h is the step length and yn is the computed approximation to y(x0 + nh), n = 0, 1, 2, . . . Then, we have following definitions. Definition. (See Konguetsof and Simos, [12]) Numerical method (19) has an interval of periodicity (0, H20 ), if ∀ H 2 ∈ (0, H20 ), Qi, i = 1(1)l satisfy |Q1| = |Q2| = 1, |Q j| ≤ 1, j = 3(1)l. (20) Definition. Following [3], a numerical method is P-stable if its interval of periodicity is (0, ∞). Therefore, we obtain the interval of periodicity of the new method, which is Equal to (0, -2.4) and the stability domain of the method is as shown in Figure 1. Figure 1: The stability domain of the new method 4.2. Order and Error Constant of the Method The method proposed by Lambert (1973) in Olanegan et al. [13] is adopted in this paper, with linear operator: k∑ j=0 α jyn+ j = h 2 k∑ j=0 β j fn+ j (21) 44 Obarhua Adegboro / J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 45 We associate the linear operator L with the scheme and defined as L{y(x), h} = k∑ j=0 [α jy(x + jh) − h 2β jy ′′(x + jh)] (22) Where α0 and β0 are both non-zero and y(x) is an arbitrary function which is continuous and differentiable on the interval [a, b]. Expanding the form y(x) and y′′(x) in Taylor series and comparing coefficients of h, we obtained ∆[y(x); h] = C0y(x) + C1hy ′(x) + · · · + C ph pyp(x) + C p+1h p+1yp+1(x) + C p+2h p+2yp+2(x) + · · · (23) The method (11) and its associate linear difference operator (13) are said to have order p if c0 = c1 = ··· = cp+1 = 0 and cp+2 , 0. The value cp+2 is called error constant. Therefore, in this paper, the method (11) is of order 4 and the error constant cp+2 = − 21 2560 or −8.2031 × 10−3 and its derivative (12) is of order 4 and error constant cp+2 = − 35 1536 or −2.2786 × 10 −2. 4.3. Consistency of the Method A numerical method is said to be consistent if the following conditions are satisfied 1. The order of the method must be greater or equal to 1, p ≥ 1. 2. ∑k j=0 α j = 0 3. ρ(r) = ρ′(r) = 0 4. ρ′′(r) = 2!σ(r) Where ρ(r) and σ(r) are the first and second characteristics polynomial of our method. according to Adesanya et al. [14], the first condition is a sufficient condition for a method to be consistent. Since our method is of order 4 then it is consistent. 4.4. Zero Stability Since |zi| = |0, 0, 1| ≤ 1 the method is zero stable. 4.5. Convergence A method is said to be convergent if and only if it is consis- tent and zero stable, hence our method is convergent. 5. Numerical Examples The method is applied to solve the following linear and non- linear second order initial value problems of ordinary differen- tial equations directly without reduction to system of first order equations. Problem 1: y′′ = y′, y(0) = 0; y′(0) = −1; h = 0.1 Theoretical solution: y(x) = 1 − ex This problem has been used in Kayode and Adeyeye [8] of or- der six to check the behavior of the methods. Table 1 shows the absolute errors of the methods for the purpose of comparison. The obtained results in the Table give the good performance of the proposed method. Problem 2: y′′ + 1001y′ + 1000y = 0, y(0) = 1; y′(0) = −1; h = 0.05 Theoretical solution: y(x) = e−x The Problem 2 was solved by Anake [15] of order 4. The new method was applied to solve it for the purpose of comparison. The results are shown in Table 2. Problem 3: y′′ = 100y′, y(0) = 1; y′(0) = −10; h = 0.01 Theoretical solution: y(x) = e−10x Table shows the absolute errors at different points of the in- tegration interval when h = 0.01was solved by Awari [16] of order five. The new method was applied to solve it for the pur- pose of comparison. The results show that the proposed method performed better than Awari [16]. Problem 4: y′′ − x(y′)2 = 0, y(0) = 1; y′(0) = 12 ; h = 0.003125 Theoretical solution:y(x) = 1 + 12 ln [ 2+x 2−x ] We have solved this problem with the new method and the re- sults have been compared with Kayode and Adeyeye [7] of or- der six shown in Table 5. Problem 5: Resonance Vibration of a Machine A stamping machine applies hammering forces on metal sheets by a die attached to the plunger moves vertically up and down by a fly wheel spinning at constant set speed. The constant ro- tational speed of the fly wheel makes the impact force on the sheet metal, and therefore the supporting base, intermittent and cyclic. The bearing base on which the metal sheet is situated has a mass, M = 2000kg. The force acting on the base fol- lows a function: F(t) = 2000 sin(10t), in which t − time in seconds. The base is supported by an elastic pad with an equiv- alent spring constant k = 2 × 105 N/m. Determine the differ- ential equation for the instantaneous position of the base y(t) if the base is initially depressed down by an amount 0.1m. Solution: The mass-spring system above is modeled as differ- ential equation as: The Bearing base mass = 2000kg Spring constant k = 2 × 105 N/m. Force (ma) on the metal sheet = m d 2 y dt2 = my ′′ i.e. ma = my′′ = 2000 sin(10t); where a = y′′Initial condi- tions on the system are y(t0) = y0; dy dt |t=o = y ′(t0) = y′0; t0 = 0, y′0 = 0.1 Therefore, the governing equation for the instantaneous posi- tion of the base y(t) is given by My′′ + ky = F(t); y(t0) = y0, y ′(t0) = y ′ 0 2000y′′ + 2 × 105y = 2000 sin 10t, y′(0) = 0, y(0) = 0.1 Theoretical solution: y(t) = 110 cos 10t + 1 200 sin 10t− t 20 cos 10t 6. Conclusion An order four continuous numerical method for solving gen- eral second order ordinary differential equations is proposed and applied to solve directly without reducing to system of first 45 Obarhua Adegboro / J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 46 Table 1 x 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Error in [8] 8.17(-7) 3.10(-6) 6.57(-6) 1.14(-5) 1.79(-5) 2.64(-5) 3.72(-5) 5.06(-5) 6.72(-5) Error in (11) 4.25(-8) 7.47(-8) 1.52(-7) 2.45(-7) 3.54(-7) 5.31(-7) 7.37(-7) 9.73(-7) 1.31(-7) Table 2 x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Error in [15] 0.10(-9) 0.20(-9) 0.28(-9) 0.34(-9) 0.39(-9) 0.43(-9) 0.45(-9) 0.4(-9) Error in (11) 2.00(-10) 3.15(-10) 2.74(-10) 5.44(-10) 7.53(-10) 2.76(-10) 1.18(-10) 1.76(-10) Table 3: Absolute errors at different points of the integration x 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Error in [16] 1.15(-7) 3.65(-7) 6.05(-7) 8.50(-7) 1.10(-6) 1.36(-6) 1.45(-6) 1.59(-6) Error in (11) 1.29(-8) 3.01(-8) 5.04(-8) 9.32(-10) 1.40(-7) 1.90(-7) 2.58(-7) 3.32(-7) Table 4: Absolute errors for Numerical example 4 x 0.0063 0.0094 0.0125 0.0188 0.0250 0.0313 0.0375 0.0437 0.0500 Error in (11) 0.00(0) 0.00(0) 2.81(-14) 2.36(-13) 8.73(-13) 1.91(-12) 2.97(-12) 5.21(-12) 7.55(-12) Table 5: Comparison of errors with [7] x Error in Kayode and Adeyeye [7] Error in the new Method (11) 0.0063 4.831(-11) 0.0000(00) 0.0094 3.382(-9) 0.0000(00) 0.0125 1.580(-8) 2.819(-14) 0.0156 4.333(-8) 1.709(-13) 0.0188 9.391(-8) 2.362(-13) Table 6: The new derived method was applied to solve this problem modeled as a second order (IVPs) and it was seen from the results in the Table that the method are useable in Engineering field. x Exact solution Computed solution Error 0.01 0.099404629653415691 0.099404630038381694 3.849660(-10) 0.02 0.097958005773976925 0.097958006644224049 8.702471(-10) 0.03 0.095207162458893865 0.095207165387981033 2.929087(-9) 0.04 0.091970827382988077 0.091970831862903016 4.479915(-9) 0.05 0.087961427477332363 0.087961431552930208 4.075598(-9) 0.06 0.082363909854646533 0.082363917430066838 7.575420(-9) 0.07 0.076833743309093400 0.076833753917924741 1.060883(-9) 0.08 0.069604876901833215 0.069604894375183718 1.747335(-9) 0.09 0.062811758617177721 0.062811776980930309 1.836375(-9) 0.10 0.055536073981512724 0.055536101603349465 2.762184(-8) order ordinary differential equations. The method is very flex- ible and easy to develop and may be applied to solve kinds of second order initial value problems as can be seen in the numer- ical examples. The method gives a high accuracy when com- pared the numerical results to the exact solution and a very good performance compared with existing methods in the literature. References [1] R. P. K. Chan, P. Leone, & A. Tsai, “Order Conditions and Symmetry for Two-Step Hybrid Methods”, Int. J. Comp. Math. 81 (2004) 1519. [2] S. Gholamtabar & N. Parandin, “Numerical Solutions of Second-Order Differential Equations by Adam Bashforth Method”, American J. of Engi. Research 3 (2014) 318. 46 Obarhua Adegboro / J. Nig. Soc. Phys. Sci. 3 (2021) 42–47 47 [3] J. D. Lambert & A. Watsan, “A Symmetric Multistep Method for Periodic Initial Value Problem”, Journal of the Institute of Mathematics and its Applications 18 (1976) 189. [4] M. F. Saravi & S.R. Mirrajei, “Numerical Solution of Linear Ordinary Differential Equations in Quantum Chemistry by Clenshaw Method”, World Acad. Sci. Eng. Technol. 49 (2009) 1051. [5] S. J. Kayode and D.O. Awoyemi, “A Multiple Derivative Collocation Method for Fifth Order Ordinary Differential Equations”, J. Math. stat. 6 (2010) 60. [6] A. Golbabai & M. M. Arabshahi, “Smart Numerov Generalized Alternat- ing Group Explicit (SNAGE-PR (2)) Scheme for 2-Point P.V. Problem”, World Appl. Sci. J. 8 (2010) 267. [7] S. J. Kayode & A. Adeyeye, “ A 3-Step Hybrid Method for Direct So- lution of Second Order Initial Value Problems”, Aust. J. of Basic and Applied Sciences 5 (2011) 2121. [8] S. J. Kayode & A. Adeyeye, “Two-Step Two-Point Hybrid Methods for General Second Order Differential Equations” Afri. J. of Math. & Compt. Sci. Res. 6 (2013) 191. [9] S. J. Kayode & F. O. Obarhua, “Continuous y− f unction Hybrid Methods for Direct Solution of Differential Equations”, Intern. J. of Diff. Eq. 12 (2013) 37. [10] S. J. Kayode & F. O. Obarhua, “3-Step y-function Hybrid Methods for Di- rect Numerical Integration of Second Order IVPs in ODEs”, Theo. Math. & Appl. 5 (2015) 39. [11] F. O. Obarhua, “An Efficient Two-Step Symmetric Hybrid Block Method for Solving Second Order Initial Value Problems of Ordinary Differential Equations”, Int. J. of Res. and Sci. Innovation 6 (2019) 200. [12] A. Konguetsof & T. E. Simos, “An Exponentially Fitted and Trigonomet- rically Fitted Method for the Numerical Solution of Periodic Initial-Value Problems”, Computers and Mathematics with Appl. 45 (2003) 547. [13] O. O. Olanegan, D. O. Awoyemi, B. G. Ogunware & F. O. Obarhua, “Con- tinuous double-hybrid point method for the solution of second order or- dinary differential equations”, Int. J. of Adv. Sci. & Tech. Res. 5 (2015) 549. [14] A. O. Adesanya, M. R. Odekunle & A. A. James, “Order Seven Hybrid Methods for the Solution of First Order Ordinary Differential Equations”, Canada Journal on Science and Engineering Mathematics 3 (2012) 154. [15] T. A. Anake, Continuous Implicit Hybrid One-Step Methods for Solutions of Initial Value Problems of General Second Order Ordinary Differen- tial Equations, Ph.D Thesis, Covenant University, Ota, Nigeria, 2011. [16] Y. S. Awari, “Derivation and Application of Six-Point Linear Multistep Numerical Method for Solution of Second Order Initial Value Problems”, IOSR Journal of Mathematics 7 (2013) 23. 47