Numerical Solution of Higher Order Ordinary Differential Equations by Direct Block Code

: Problem statement: This study is concerned with the development of a code based on 2-point block method for solving higher order Initial Value Problems (IVPs) of Ordinary Differential Equations (ODEs) directly. Approach: The block method was developed based on numerical integration and using interpolation approach which is similarly as Adams Moulton type. Furthermore, the proposed method is derived in order to solve higher order ODEs in a single code using variable step size and implemented in a predictor corrector mode. This block method will act as simultaneous numerical integrator by computing the numerical solution at two steps simultaneously. Results: The numerical results for the direct block method were superior compared to the existing block method. Conclusion: It is clearly proved that the code is able to produce good results for solving higher order ODEs.


INTRODUCTION
The mathematical formulation of physical phenomena in science and engineering often leads to IVPs of ODEs. This type of problem can be formulated either in terms of first order ODEs or higher order ODEs. For instance, this application will often used in beam theory, electric circuits, control theory, mechanical system and celestial mechanics. Thus, this study will concern on solving directly higher order nonstiff IVPs of ODEs of the form: (1) The conventional methods of solving higher order ODEs will reduce such problems to a system of first order equations. This approach is cumbersome and will increase computational time as well as consume a lot of human effort. Thus, several researchers have concerned themselves with study to solve Eq. 1 directly such as Awoyemi (2003); Majid (2004); Awoyemi (2005); Majid and Suleiman (2006); Jator (2010); Jain et al. (1977) Kayode and Awoyemi (2010). Awoyemi (2005) has proposed a multiderivative collocation method for direct solution of fourth order IVPs of ODEs while Majid and Suleiman (2006) have introduced a direct integration implicit variable step method for solving higher order systems of ODEs. Majid (2004) has developed the 2-point block method for solving first and second order ODEs using variable step size. It was noted that Jator (2010) had used the application of a self starting linear multistep method for solving second order IVPs directly.
Block methods for numerical solution of higher order ODEs have been proposed by several researchers such as Chu and Hamilton (1987); Fatunla (1991) and Jator (2010). Chu and Hamilton (1987) have proposed multi-block methods for parallel solution of ODEs and Fatunla (1991) has presented a zero stable block method for second order ODEs. The uniqueness of block method is that in each application, the solution value will be computed simultaneously at several distinct points. There are several existence numerical methods for handling higher order ODEs directly but those methods will compute the numerical solutions at one point sequentially. Henceforth, we need a method that can give faster solution of the problem.
In this study, we are going to extend the study done in Majid (2004) by implemented the 2-point block method to solve ODEs up to order five in a single code.

MATERIALS AND METHODS
Formulation of the method: In 2-point block method, the closed finite interval [a,b] is divided into a series of blocks that contained the interpolation points involved in the derivation of the block method.
According to Fig. 1, the method will generate the numerical solution at two points simultaneously. The values of y n+1 and y n+2 at the points x n+1 and x n+2 respectively are simultaneously computed in a block with step size h using the same back values which is the values at the point x n , x n-1 and x n-2 with step size rh. The formulae of 2-point block method are derived by integrating Eq. 1 d times as follows Eq. 2: Let x n+v = x n + vh, where v = 1 or 2.
n v n n n n v n n n (2) which leads to the general formula below: where, g is the number of times which Eq. 3 is integrated over the corresponding interval. Lagrange interpolation polynomial is used to approximate the function of ( ) and dx = hds will be substitute into Eq. 3. By taking d = 5 in Eq. 1, the approximate solution of y n+1 and y n+2 will be obtained by integrating Eq. 1 once, twice, thrice, four times and five times over the interval [ ] n n 1 x , x + and [ ] n n 2 x , x + respectively. Finally, this integral will be evaluated using MAPLE and the corrector formulae in terms of r will be obtained. The same approaches were employed in the derivation of the predictor formulae for ( ) iv y, y , y , y , y ′ ′′ ′′′ at the points x n+1 and x n+2 respectively and the interpolation points involved are ( ) ( ) Hence, it will produce the predictor formulae in terms of q and r. For the sake of simplification, the general corrector formula of 2-point block method is developed in the manner shown in Eq. 4: g v, j β in Eq. 4 stands for the coefficients of the formulae and tabulated in Table 1-3 for r = 1, r = 2 and r = 0.5.  The order of this developed method is calculated in a block form as proposed by Fatunla (1991). The 2point block method for ODEs can be written in a matrix differentiation Eq. 5 below: where, α, β, λ and ξ are the coefficients of the 2-point block method. By applying the formulae for the constants C q in Fatunla (1991), the order and error constant of the method will be obtained and the formulae are defined by Eq. 6: The method is of order p if C q = 0, q = 0(1) p + d-1, C p+d ≠ 0. Thus, by implementing this approach to the 2-point block method, we found that the predictor is of order four and corrector is of order five. The error constant for the corrector formulae when r = 1 will be in matrix form as shown in Table 4. Implementation of the method: A single code of the PECE scheme has been implemented with variable step size to study the computational time and human effort saving in using a direct integration method. The developed code starts by calculating the initial step size then finding the initial points in the starting block of the method. In order to evaluate the initial three starting points, the Euler method was adopted in the code as a generator of the method. Hence, the Euler method will be used only once at the beginning of the code. Once the points for starting block are calculated, then the block method can be applied until the end of interval.
A test for checking the end of the interval is made in order to reach the end of the interval precisely and it will be functional at each step of integration. The strategy is by limiting the choices of the next step size to half, double or remains constant as the previous step size. At each step of the integration, if the approximated solutions fulfilled the desired accuracy, therefore the step is called as successive step.
Hence, the choices for the next step size will be doubled or constant which specified by step size controller. Otherwise the step is called failure step and the next step size becomes half.
The possible ratios for the next constant step size are (r = 1, q = 1), (r = 1, q = 2) and (r = 1, q = 0.5). At each doubled step size the ratios are (r = 0.5, q = 0.5) and (r = 0.5, q = 0.25). In the case of failure step size, the ratio is (r = 2, q = 2). A step failure happens due to the Local Truncation Error (LTE) exceeding the given tolerance. This corrector formulae show that, the code consist the formula of y, y , y , y ′ ′′ ′′′ and y (iv) . Hence, the developed algorithm could be use for solving first order up to fifth order problem of ODEs directly in a single code. The algorithm was written in C language. The approximation value of y n+1 and y n+2 are using predictor-corrector mode, (P k E) (C k+1 E) u where P k and C k+1 indicate the predictor of order k and corrector of order k+1 respectively and E indicate the evaluation of the function. In the code, the corrector will be iterated until it is converge and the convergence test employed was Eq. 7: where, u is the number of iterations. The LTE will be obtained by comparing the absolute difference of the corrector formula derived of order k with the similar corrector formula of order k+1 at the point x n+2 . If the LTE<TOL, hence the successful step achieved and the next step size will be obtained using the step size increment formula as follows Eq. 8: Where: δ = A safety factor k = The order of the corrector formulae while h new and h old are the step size for current and previous block respectively Finally, the errors in the code are calculated by the difference between an estimation value y n+2 resulting from using 2-point block method and the exact solution. The formula was defined as Eq. 9: where, (y i ) t is the t-th component of the approximation y. A = 1, B = 0 referred to the absolute error test, while A = 1, B = 1 will be used for mixed error test and A = 0, B = 1 correspond to the relative error test.

RESULTS
The results of numerical tests will be presented in order to illustrate the performance of 2-point block method. The following notations are used in the table: Solution: x 2 y(x) e x . = +

DISCUSSION
We are going to compare the numerical results obtained by solving the tested problems using a direct integration approach with the results using a reduction approach. It is clearly shown that the 2PDVS has superiority in terms of computational time and total number of steps taken especially as the tolerance getting smaller. This indicates the major advantage of direct integration method compared to reduce to dth order equation to d sets of first order equations. Table 5-8 show the maximum error of 2PRVS is better and comparable compared to 2PDVS.    Nevertheless, 2PDVS is still able to obtain the desired accuracy within the given tolerance. Furthermore, 2PDVS is economical to be implemented than 2PRVS since the number of functions to be evaluated by 2PDVS is much lesser than 2PRVS especially at finer tolerance. Therefore, by solving higher order ODEs directly will reduce the computational cost and will inevitable effects on the execution time.

CONCLUSION
In this study, we have shown the efficiency of the developed two point block method presented as in the form of Adams Moulton method using variable step size which is suitable for solving higher order ODEs directly.