L-Stable Block Hybrid Numerical Algorithm for First-Order Ordinary Differential Equations

In this work, a one-step L-stable Block Hybrid Multistep Method (BHMM) of order five was developed. The method is constructed for solving first order Ordinary Differential Equations with given initial conditions. Interpolation and collocation techniques, with power series as a basis function, are employed for the derivation of the continuous form of the hybrid methods. The discrete scheme and its second derivative are derived by evaluating at the specific grid and off-grid points to form the main and additional methods respectively. Both hybrid methods generated are composed in matrix form and implemented as a block method. The stability and convergence properties of BHMM are discussed and presented. The numerical results of BHMM have proven its efficiency when compared to some existing methods.


Introduction
Most Mathematical models are formulated using ordinary differential equations (ODEs) of different orders. However, mathematical models that are based on ODEs of order one occur in several engineering, applied sciences and economics problems. Some of the ODEs have been proven not to have closed-form solution hence the need to develop numerical methods to provide approximate solutions to these problems. Consider the Initial Value Problem of the form y (x) = f (x, y(x)), y(x 0 ) = y 0 Several numerical methods have been developed by different researchers to circumvent the Dahlquist's order Barrier Theo-rem [1 -3] that restricted the use of Linear Multistep Methods of high order to solve (1). Some of such researchers are Enright [4], Akinnukawe and Okunuga [5], Gear [6], Okunuga [7], Cash [8], Akinfenwa et al. [9,10], Adesanya et al. [11], Ijezie and Muka [12] to mention but a few. High-order A-stable and L-stable numerical methods are developed by incorporating off-step points, additional stages and/or employing higher differentiation of the solution (Lambert [13]). Generating approximate solutions for (1) in instances where they exist and are unique but are insoluble via analytical means is very important because they help in the analysis and validation of models in which they evolve from. One of the objectives of developing numerical schemes for solving (1) is to obtain methods with wider stability regions, better convergent rates and computational efficiency.
In this article, Continuous Hybrid Methods (CHM) of order 5 are derived through Interpolation and Collocation techniques (Onumanyi et al. [14]). Both methods incorporate off-step points and the second derivative of the scheme to produce the discrete hybrid methods which are composed and implemented as a block method (BHMM) to simultaneously produce approximation at nodal points. The BHMM has the advantage of being self-starting, L-stable in nature and possesses high accuracy because it implemented as a block method (see [15], [16], [17]).

Derivation of BHMM
This section describes the derivation of the k-step hybrid multistep method of the form: where h, n and v = 1 2 are the step size, grid index and offstep point respectively while α j , β j , j = 0, 1, ..., k and φ k are parameters to be determined uniquely. An approximate solution to (1) by the interpolating function where b j , j = 0, 1, ..., 4k + 1 are the unknown coefficients. Imposing condition for the construction of the proposed class of methods are: y (x n+r ) = f n+r , r = 0(v)k (5) y (x n+r ) = g n+r , r = k Equations (4) -(6) will lead to a system of 4k + 2 equations. These equations are solved simultaneously to obtain b j and the values of b j are substituted into (3) to form the Continuous Hybrid Method (CHM) expressed in the form where t = x−x n+k−1 h and α j (t), β j (t), j = 0(v)k and φ k (t) are the continuous coefficients. The main scheme is generated evaluating CHM (7) at x n+k while the additional scheme is generated from the second derivative of (7) at x n+v of the form The developed main and additional schemes are combined and implemented simultaneously as a block hybrid multistep method BHMM for the numerical integration of IVPs (1).
Following the steps discussed above, the Continuous Hybrid Method (CHM) for k=1 is equation (7) Interpolating (7) at x = x n+1 to generate the main method at k = 1 becomes Equation (8), the second derivative of (7) for k = 1 is Interpolating(8) at x = x n+1/2 , the additional method at k = 1 becomes The discrete hybrid methods (9a)and (9b) together forms the One-step Block Hybrid Multistep Methods (BHMM). The BHMM can be presented in a matrix block form as The 2 by 2 matrices A (0) , A (1) , B (0) , B (1) , C (1) , D (1) of the BHMM (9) are defined as follows

Order and Error Constant of the Method
Following Lambert [13] and Fatunla [18], a method was proposed for finding the order p and error constant W p+1 of the block method (9) by first expanding y−, f − and g− functions by Taylors series expansion about x and then comparing the coefficients of h. It is established from our calculation that one-step Block Hybrid Multistep Method have order and error constants as p = (5, 5) T and W p+1 = ( 13 44160 , 1 66240 ) T respectively where T is transpose.

Zero Stability
A numerical method is said to be zero-stable if the roots R j , j = 1, 2, . . . , N of the first characteristic polynomial ρ(R) satisfies |R j | ≤ 1, j = 1, . . . , N and those roots with |R j | = 1 is simple (see Lambert [3]). Applying the above conditions to the derived block method, the first characteristic polynomial ρ(R) = 0 is calculated as The BHMM is found to be zero-stable since ρ(R) = 0 satisfies |R j | ≤ 1, j = 1, 2.

Convergence
According to Henrici [19], a numerical method converges if it is consistent and zero-stable. Since BHMM (9) is of order 5 > 1, then it is consistent and we have established earlier that the method satisfies the conditions of zero-stability. Therefore, the block method (9) converges.

Stability of BHMM
Applying the BHMM to the test equation where Q(z) is the amplification matrix given by has eigenvalues (ζ 1 , ζ 2 ) = (0, ζ 2 ). The dominant eigenvalue ζ 2 is the stability function with real coefficient as 10.4348 + 4.17391z + 0.652174z 2 + 0.0434783z 3 10.4348 − 6.26087z + 1.69565z 2 − 0.26087z 3 + 0.0217391z 4 The stability function is used to plot the Region of Absolute Stability (RAS) of the BHMM (see Figure 1). The proposed method BHMM is L-stable since the RAS covers the entire left plane of the graph (A-stable) and the limit of the stability function ζ 2 is zero as z → ∞ 162

Numerical Results
The following problems are considered to examine the accuracy and computational efficiency of BHMM. The computations were carried out using MATHEMATICA 9.0 software. The following acronyms are used: where λ = 10 4 and the exact solution of the IVP is given as In Table 1, the numerical results obtained using BHMM with h = 10 −1 compares favourably with MHIRK method [10] with h = 10 −1 and is superior to that of Hojjati et al. [20] that used h = 10 −4 .  Table 2 for the numerical results.
Problem 3: Consider the non-linear system y 1 = −2y 1 + y 2 + 2sin x y 2 = 998y 1 − 999y 2 + 999(cos x − sin x) with initial conditions y 1 (0) = 2, y 2 (0) = 3, 0 ≤ x ≤ 10 with analytical solution given as The numerical results of Problem 3 are shown in Table 3 using step size h = 10 −3 . The derived method integrated the problem efficiently that the numerical solution is close to that of the analytical solution. with initial conditions 15] with exact solution given as Problem 4 was integrated using the step size of h = 10 −4 to aid in comparing with other methods in literature as shown in Tables 4 and 5. It is discovered that BHMM has better accuracy than the others compared with.