Implementing a Type of Block Predictor-corrector Mode for Solving General Second Order Ordinary Differential Equations

The paper is geared towards implementing a type of block predictor-corrector mode capable of integratinggeneral second order ordinary differential equations using variable step size. This technique will be carried out on nonstiff problems. The mode which emanated from Milne’s estimate has many computation advantages such as changing and designing a suitable step size, correcting to convergence, error control/minimization with better accuracy compare to other methods with fixed step size. Moreover, the approach will adopt the estimates of the principal local truncation error on a pair of explicit (predictor) and implicit (corrector) Adams family which are implemented in P(CE) m mode. Numerical examples are given to examine the efficiency of the method and compared with subsisting methods.


INTRODUCTION
Rising from the advent of computing machines and programming languages, the numerical solution of Initial Value Problems (IVPs) for Ordinary Differential Equations (ODEs) has been the topic to explore by numerical analysts, mostly procedures for the numerical solution of the general second-order ODEs of the form as seen in Ken et al. (2011): (1) The general solution to (1) can be coded as: where, the step size is h, α j = 1, α j , i = 1, … j, β j , are unknown constants which are uniquely specified such that the formula is of order j as discussed in Akinfenwa et al. (2013).
We assume that f ∈ R is sufficiently differentiable on x ∈ [a, b] and satisfies a global Lipchitz condition, i.e., there is a constant L ≥ 0 such that:  1) assured the existence and uniqueness defined on x ∈ [a, b] as discussed in Lambert (1973) and Xie and Tian (2014).where, a and b are finite and y (I) [y (i)  1 , y (i) 2 ,…, y (i)  n ] T for i = 0(1)3 and f = [f 1 , f 1 ,…, f n ] T , Again, Weierstrass approximation theorem stands as a justification for (1).See (Jain et al., 2007) for details.
Eq. ( 1) arises from many physical phenomena in a broad compass of applications.Largely in the field science and engineering such as in the electric circuits, damped and undamped mass-spring systems, forced oscillations and other areas of practical applications as introduced by Majid et al. (2012).Authors such as Anake et al. (2012), Ismail et al. (2009) and Majid and Suleiman (2009) indicated that (1) can be reduced to an equivalent first-order system of two times the dimension and evaluated utilizing the existing one-step method like Runge-Kutta method or block multistep method.This method has been described to increase the dimension of the problem, computational effect and very difficult.Block multistep methods are one of the numerical methods which have been suggested by several researchers, (Adesanya et al., 2012;James et al., 2013;Ken et al., 2011;Majid et al., 2012;Majid and Suleiman, 2009;Zarina et al., 2007).The common block methods used to solve the problems can be arranged into categories as one-step block method and multistep block method.Nevertheless, scholars have proposed an alternative method to solve at once (1) as discussed in Adesanya et al. (2012Adesanya et al. ( , 2013)), Anake et al. (2012), Ehigie et al. (2011), Ismail et al. (2009), Ken et al. (2011) and Majid et al. (2012).
However, (Adesanya et al., 2012;Ehigie et al., 2011;Ismail et al., 2009;Ken et al., 2011) proposed block multistep methods which were employed in predictor-corrector mode.Block multistep methods have the vantage of evaluating simultaneously at all points with the integration interval, thereby reducing the computational burden when evaluation is required at more than one point within the grid.Again, Taylor series expansion is used to provide the initial values in order to compute the corrector.
Researchers such as Adesanya et al. (2012), Ehigie et al. (2011), Ismail et al. (2009) and Ken et al. (2011), innovated block predictor-corrector method in which at each practical application of the method, the method was only intended to predict and correct the results generated.In this study, the motivation is stemmed by the fact that there are very few work been done in solving nonstiff ODEs using block predictor-corrector mode, effort will be geared towards developing a type of block predictor-corrector mode using variable step size technique otherwise called Milne's estimate.This method have several benefits as earlier stated in the abstract.
Thus, from the above definition a block method has the vantage that in each application, the solution is approximated at more than one point simultaneously.The number of points depends on the manner of construction of the block method.Therefore applying these methods can give quicker and faster solutions to the problem which can be managed to produce a desired accuracy.See (Majid and Suleiman, 2007;Mehrkanoon et al., 2010).
The block algorithm proposed in this study is based on interpolation and collocation.The continuous representation of the algorithm generates a main discrete collocation method to render the approximate solution Y n+i to the solution of (1) at points x n+i , i = 1, …, k as in Akinfenwa et al. (2013).The main aim of this study is to introduce a type of block predictorcorrector mode using variable steps size technique for mathematically integrating general second order ODEs directly.

DERIVATION OF THE METHOD
Adopting (Akinfenwa et al., 2013) in this section, the target is to derive the principal block predictorcorrector mode of the form (2). We move ahead by seeking an estimate of the exact solution y(x) by assuming a continuous solution Y(x) of the form: , m i are unknown coefficients and ϑ i (x) are polynomial basis functions of degree q+k -1, where q is the number of interpolation points and the collocation points k are respectively chosen to satisfy q = j and k≥1.The integer j≥1 denotes the step number of the method.We thus construct a j-step block method block method with ϑ i(x) = Ә ә by imposing the following conditions: where, y n+i is he approximation for the exact solution y(x n+i ), f n+i = f(x n+i , y n+i ) n is the grid index and x n+i = x n +ih .It should be noted that Eq. ( 4) and ( 5) leads to a system of q+1 equations of the AX = B. where, Solving Eq. ( 6) using Mathematica, we get the coefficients of m i and substituting the values of m i into (4) and after some algebraic computation, the block predictor-corrector mode is obtain as: where, α i and β i are continuous coefficients.

ANALYSIS OF SOME THEORETICAL PROPERTIES
Order of accuracy of the method: Conforming to Akinfenwa et al. (2013) and Lambert (1973), we specify the associated linear multistep method (7) and the difference operator as: [ ] Presuming that y(x) is sufficiently and continuously differentiable on an interval [a, b] and that y(x) has as many higher derivatives as demanded then, we write the conditions in (8) as a Taylor series expression of y(x n+i ) and f(x n+i ) ≡ y" (x n+i ) as: Substituting ( 8) and ( 9) into ( 7) we obtain the following expression: Hence, we noticed that the Predictor-Corrector Mode (P(CE) m ) of ( 7) has order p, if C p+1 , p = 0, 1, 2, …, I = 1, 2,…, j, are given as follows: , ...
Therefore, the method (7) has order p≥1 and error constants given by the vector, C p+2 ≠ 0.
Concurring with Lambert (1973), we say that the method (2) has order p if: Therefore, C p+2 is the error constant and C p+2 h p+2 y (p+2) (X n ) is the principal local truncation error at the point x n .

Stability analysis of the method:
To examine the method for stability, ( 7) is normalised and composed as a block method given by the matrix finite difference equations as presented in Akinfenwa et al. (2013), Ken et al. (2011), Mohammed et al. (2013) and Awari (2013). where, The matrices A (0) , A (1) , B (0) , B (1) are r by r matrices with real entries while Y m , Y m-1 , F m , F m-1 are r-vectors specified above.
Following (Ken et al., 2011;Lambert, 1973), we stick to the boundary locus method to decide on the region of absolute stability of the block predictorcorrector mode and to obtain the roots of absolute stability.Substituting the test equation y' = -λy and ℎ = h 2 λ 2 into the block method ( 12) to obtain: Replacing h = 0 in (13), we obtain all the roots of the derived equation to be r≤1.Therefore, according to Lambert (1973), the predictor-corrector mode is absolutely stable.
So, as considered in Adesanya et al. (2013), Lambert (1973) and Awari (2013), the boundary of the region of absolute stability can be obtained by filling (7) into:

IMPLEMENTATION OF THE METHOD
Embracing (Faires and Burden, 2012;Lambert, 1973), afterward this is implemented in the P(EC) m mode then it becomes very pertinent if the explicit (predictor) and the implicit (corrector) methods are individually of the same order and this prerequisite makes it essential for the step number of the explicit (predictor) method to be greater than that of the implicit (corrector) method.Consequently, the mode P(EC) m can be formally decided as follows for m = 1, 2, ….: P(EC) m : Note that as m→∞, the result of calculating with the above mode will incline to those given by the mode of correcting to convergence.
Moreover, predictor-corrector pair based on (1) can be applied.The mode P(EC) m specified by ( 15), where h 2 is the step size.Since the predictor and corrector both have the same order p, Milne's device is applicable and relevant.
According to Dormand (1996) and Lambert (1973), Milne's device suggests that it is possible to estimate the principal local truncation error of the explicitimplicit (predictor-corrector mode) method without estimating higher derivatives of y(x).Assume that p = p * , where p * and p represents the order of the explicit (predictor) and implicit (corrector) method with the same order.Now for a method of order p, the principal local truncation errors can be well defined as: where, W n+j and C n+j are called the predicted and corrected approximations given by method of order p while C * p+2 and C p+2 are independent of h.Neglecting terms of degree p+3 and above, it is easy to make estimates of the principal local truncation error of the mode as: Noting the fact that C p+2 ≠ C * p+2 and W n+j ≠ C n+j .Furthermore, the estimate of the principal local truncation error (18) is used to determine whether to accept the results of the current step or to reconstruct the step with a smaller step size.The step is accepted based on a test as prescribed by (18) as in Ascher and Petzold (1998).Equation ( 18) is the convergence criteria otherwise called Milne's estimate for correcting to convergence Furthermore, Eq. ( 18) ensures the convergence criterion of the mode during the test evaluation.

NUMERICAL EXPERIMENTS
The performance of the block predictor-corrector mode (P(CE) m ) was carried out on nonstiff and mildly stiff problems.For problem tested 5.1 and 5.2, the following tolerances (convergence criteria)10 -2 , 10 -4 , 10 -6 , 10 -8 , 10 -10 and 10 -12 were used to compare the performance of the newly proposed method with other existing methods as in James et al. (2013) and Ken et al. (2011).The exact solution is given by: The exact solution is given by: The first tested problem 5.1 to be discussed was extracted from James et al. (2013).Moreover, four steps continuous method for the solution of y" = f(x, y, y') was developed and implemented using fixed step size.Thus, the newly proposed method is formulated to solve nonstiff using variable step size technique.
Secondly, tested Problem 5.2 was extracted from Ken et al. (2011).However, the block methods for special second order ODEs was designed and executed using fixed step size technique.Moreover, the implementation of the explicit 2-point 1-block method and implementation of the explicit 3-point 1-block method was carried out using linear difference operator as well as their comparison.The newly proposed method belongs to the family of Adams otherwise called Milne's estimate and was created to solve nonstiff.ODEs.
simplifying (14) within [0 0 , 180 0 ].Accordingly, the boundary of the region of absolute stability rests on the real axis.Fig.1is free hand drawing.

Fig. 1 :
Fig. 1: Showing the region of absolute stability of the block predictor-corrector mode, since the root of the stability polynomial is r≤1

Table 1 :
(James et al., 2013)cal results ofJames et al. (2013)and newly proposed method (BPC) for solving problem 5.1(James et al., 2013)newly proposed method (BPC) Table Comparing the numerical results of Ken et al. (2011) and newly proposed method (BPC) for solving problem 5.2 TOL MTH