A NEW NUMERICAL INTEGRATOR FOR THE SOLUTION OF GENERAL SECOND ORDER ORDINARY DIFFERENTIAL EQUATIONS

This paper considered the development of numerical integrator for the solution of second order initial value problems. The method was derived through the interpolation and collocation of the basis polynomial which is combination of power series and exponential function to derive a continuous linear multistep method. The method was implemented in block method which gave solution at a non overlapping interval. The method was found to be convergence and A stable. The efficiency of the method was tested on some numerical examples. AMS Subject Classification: 65L05, 65L06, 65D30


Introduction
This paper considers the numerical solution to general second order initial value problem of the form where x n is the initial point, y (x n ) is the solution at the initial point, f is continuous within the interval of integration and assumed to satisfied the existence and uniqueness theorem stated below Theorem 1. (see [1]) Let k = 0, 1, ..., (n − 1) (u and f are scalars).Let R be the region defined by the inequalities x 0 ≤ x ≤ x 0 + a, |s j − c j | ≤ b, j = 0, 1, ..., (n − 1) (a > 0.b > 0).Suppose the function f (x, s 0 , s 1 , ..., s n−1 ) is defined in ℜ and in addition: (a) f is non negative and non decreasing in each of x, s 0, s 1,..., s n−1 in ℜ; (b) f (x, c 0 , c 1 , ..., c n−1 ) > 0, for x 0 ≤ x ≤ x 0 + a, and (c) c k ≥ 0, k = 0, 1, ..., n − 1.
Then the initial value Problem 2 has a ubique solution in ℜ.
See [1] for the proof.
Equation (1) are convectionally solved by reducing to systems of first order ordinary differential equation before appropriate method is applied.This method has been reported not to be efficient because it increases the dimension of the resulting systems of first order by the order of the differential equation; hence it wastes both the computer and human effort.Other setbacks of the method of reduction are discussed by [2,3,4,5,6], consequently, direct method for the solution of (1) has arouse the interest of scholars.
Scholars have adopted method of interpolation and collocation of polynomial basis function to developed a continuous linear multistep method which can be implemented over an overlapping interval (predictor-corrector method) or non overlapping interval (block method).scholars reported that block method is better than the predictor corrector method because it is more efficient in terms of accuracy, cost of implementation, time effectiveness and the flexibility of the derived method; hence scholars has proposed different numerical scheme implemented in block method using different polynomial basis function ranging from power series, Langrange polynomial, Chebychev polynomial to mention but few ( [7,8,9,10]).
It was discovered that most of this method are not adequate for problems which are non linear, oscillatory and stiff problems due to poor stability properties exhibited by these methods.Definition 2. (see [1]) A problem is said to be stiff if the following condition are fulfilled: (a) No solution component is unstable or equivalently, no eigenvalue of the jacobian matrix has a real part which is at all large and positive and at least some component is very stable.That is atleast one eigenvalue has a negative part which is negative and large.
(b) The solution is slowly varying with respect to the negative real part of the eigenvalues.Definition 3. (see [11]) The initial value Problem 1 is considered to be stiff oscillatory if the eigenvalues {λ j = u j +iv j , j = 1(1) m} of the Jacobian J = ∂f ∂y possess the following properties It should be recalled that stiff initial value problems was first encountered in the study of motion of spring of varying stiffness, since then, most problem have been discovered to be stiff especially in kinematics, chemical reactions, process control and electrical circuit theory, see [12].Also discribed stiffness of ordinary differential equations in a manner that Problem 1 possess some stiffness if R e (λ i ) < 0,where λ is the eigenvalue of the problem.
In order to improve on the stability of k-step method, scholars introduced hybrid points; hence hybrid linear multistep method was proposed, [3].Reported that the introduction of hybrid points improved both the stability and accuracy of the method, [3,7].Proposed one step method for the solution of higher order ordinary differential equation and both concluded that the method gave better stability properties and more accurate than the k-step method.This method was able to circumvent the Dahlquist stability barrier.It should be reminded that the implicit Eulers method was first proposed to solve stiff problems and was concluded that schemes which possess the properties of Astability are suitable for stiff problems.Definition 4. (see [12]) A numerical method is said to be A-stable if the whole of the left-half plane {z : Re(z) ≤ 0} is contained in the region {z : |Re(z) ≤ 1|}, where R(z) is called the stability polynomial of the method.
In the quest for the method that gives better stability properties, [13] proposed an approximate solution which is a combination of power series and exponential function.It was discovered that the method gives an A-stable method irrespective of how the points are selected.In this paper, we combined power series and exponential function as our basis function to develop our method.

Methodology
We consider a polynomial basis function in the form where r and s are the number of interpolation and collocation points respectively.x j is the polynomial basis function.a ′ j s ∈ ℜ are constants to be determined on the partition ∆ Substituting the second derivatives of (3) into (1) gives Interpolation (3) at x n and collocating (4) at x n+s , s = 0 1 5 1 gives a system of non linear equation which can be written compactly in the form and Solving for the unknown constants a ′ j s in ( 5) and ( 6) using Guassian elimination method and substituting into (3) gives an hybrid continuous linear multistep method of the form where Evaluating (7) and it first derivative at t = 2 5 1 5 1, and 0 1 5 1 respectively, solving for the independent solution at the selected grid points gives a discrete block method in the form where A (0) = 5 × 5 identity matrix, when i = 0,

Implementation of the Method
In order to implement the method, we propose a prediction equation of the form Substituting ( 9) into (8) gives Hence, (10) represents our new method.

Order of the Block
Let the linear operator ℓ{y(x) : h} associated with the discrete block method (8) be defined as Expanding (8) in Taylor series and comparing the coefficient of h gives ℓ{y(x) : h} = C 0 y(x) + C 1 hy 1 (x) + ... + C p h p y p (x) + C p+1 h p+1 y p+1 (x) + ... .Definition 5.The linear operator ℓ and associated block formular are said to be of order p if C 0 = C 1 = ...C p = C p+1 = 0, C p+2 = 0. C p+2 is called the error constant and implies that the truncation error is given by t n+k = C p+2 h p+2 y p+2 (x) + 0(h p+3 ).
For our method, expanding (8) in Taylor series and comparing the coefficient of h gives the error constant as T .

Zero Stabilty of the Method
Definition 6.A block method is said to be zero stable if as h → 0, the roots, r j = 1(1)k of the first characteristics polynomials p(x) = 0 that is p(r) = det A (0) R k−1 = 0 satisfying |R| ≤ 1 must have multiplicity equal to unity.
For our method: Hence, the new block integration is zero-stable.

Consistency
The block integrator ( 8) is consistent since it has order p = 6.

Convergence
Definition 7. (see [3]) A method is said to be convergence if lim y n where y (x n ) is the approximate solution, y n is the analytical solution.
Theorem 8.The necessary and sufficient conditions that a continuous linear multistep method be convergentare that it must be consistent and zerostable.

Region of Absolute Stability
Definition 9. Region of absolute stability is a region in the complex Z−plane, where τ = λh.τ is defined as those values of Z such that the numerical solution of y ′′ = −λ 2 y and y ′ = −λy satisfies y j → 0 as j → ∞ for any initial value condition.We adopted the boundary locus method to determine the stability of our method.Substituting y ′′ = −λ 2 y, y ′ = −λy into (7) gives the stability region as shown in Figure 1.

Numerical Examples
Problem 10.We consider the mildly stiff problem given by This problem was solved by [14].The result is shown graphically and in Figure 2 and tablulated in Table 1.It should be noted that the maximum error was taken at y (10.002) = 4.5309220642 (−5) Problem 11.We consider the non linear second order Van der pol equation This problem is known not to have analytical solution and is highly stiff if µ = 1000.We state the condition for finding the solution to this problem.The solution to Problem 2 is shown in Figure 3 and tabulated in Tables 2  and 3.
Problem 13.We consider the non linear initial value problem This problem was solved by [5] where a five step method implemented in block method was proposed.The result is shown in Figure 4 and Table 4.
Exact solution y (x) = cos 2x + sin 2x.This problem was solved by [16] where a method fo order five implemented in block predictor-block corrector was proposed.The solution is shoem in Figure 5 and Table 4.It should be noted that the result is taken at x = 1.0, y (1) = 0.493150590278540.

Discussion of Result
We have considered four numerical examples in this paper.The result as shown graphically in order to confirm if they conform with the exact solution.It should be noted that the method give better results if the value of h ≤ 0.0032 for Problem 1. and Table 1 re-affirm the convergence of the method.∆ 4 is  3 also affirm the convergence of the method.Table 4 compared our results with the existing method and shown that our method has better stability and gives better approximation as is increasing.Table 5 shown that our method is better than the existing method we compared with at equally affirm the stability by varrying the value of h.

Conclusion
We have proposed a numerical method in this paper.The approximate solution is combination of power series and exponential function.The developed method has a better stability properties compared with the existing method as shown in the numerical example and compete favourably.It must be noted that the method proposed in this paper is not a self starting.

Figure 1 :
Figure 1: Stability region of the method

Table 1 :
Results of Problem 1

Table 2 :
Results of Problem 2

Table 3 :
Results of Problem 3

Table 4 :
Results of Problem 3Note that ∆ 1

Table 5 :
Results of Problem 4Problem 14.We consider a highly oscillatory test problem