A STIFFLY STABLE SECOND DERIVATIVE BLOCK MULTISTEP FORMULA WITH CHEBYSHEV COLLOCATION POINTS FOR STIFF PROBLEMS

Abstract: Most block methods in the literature which are implemented in predictor-corrector mode, usually suffer some stability setbacks and this may hinder their implementation on some stiff problems. In this paper, we construct a stiffly stable block second derivative backward differentiation formula with Chebyshev collocation points that is self-starting and is capable of solving stiff problems. The method is applied in block form as a simultaneous numerical integrator over non-overlapping subintervals. The method is proven to possess stiffly stable, A0 stable and A(α) stable properties. Some numerical examples reveal that this class of methods is very promising and are suitable for solving stiff problems.

Engineering and Sciences are known to be stiff in nature.Most realistic models developed from these problems cannot be solved analytically.Given a system of ordinary differential equations of the form where y = (y 1 , y 2 , • • • , y s ) and η = (η 1 , η 2 , • • • , η s ).Let λ i be the eigenvalues of the s × s matrix A, (1) is said to be stiff if Re (λ i ) < 0, i = 1, 2, • • • , s, and Since the famous theorem of G. Dahlquist [9] ( also known as the Dahlquist barrier), numerical analyst sought for robust numerical methods that can cope with stiff problems and this class of methods are called A−stable methods.It is in this notion that Widlund [29], Bickart and Rubin [2] proposed that to derive numerical methods with high order capable of solving stiff problems, some stability conditions may be relaxed and new class of methods can be generated to circumvent the Dahlquist barrier.
Curtiss and Hirschfelder [8] introduced the backward differentiation formula for solving stiff equations, through which most stiff codes are based.A survey of methods on stiff problems can be found in the literature [14,15].Hybrid methods for initial value problems which involve a combination of one-step procedures and the Runge-Kutta procedures were introduced by Gear [13] to overcome such barrier and these methods proved to be convergent under suitable conditions of stability and consistency.
For fast and efficient simulation of applied problems, numerical methods known as block methods which are capable of obtaining numerical solutions at several points were presented by several researchers [1, 6, 7, 10, 12, 17, 19, 22 -28].
In this paper, the interpolating function where the polynomial basis function φ j (x) = x j , j = 0, 1, • • • is used to approximate the theoretical solution of the problem.
To do this, we construct a block Second Derivative Backward Differentiation Formula (SDBDF) with Chebyshev collocation points, where these nodes are also included in the collocation points as zeroes of the Shifted Chebyshev polynomials.The Chebyshev polynomials is utilized to reasonably spread the errors uniformly on the interval of integration because of its great importance to approximation theory.Also, for several methods on collocation methods for stiff problems see [1,22,23,27].A continuous multistep formula is generated from this procedure and evaluation at some points yields the block methods which helps to make the implementation procedure self starting.[3,17,25] This article is organized as follows: we begin at Section 2 with the theoretical procedure which involves the construction of the method, the block representation of the method and the self starting implementation of the method.In Section 3, a class of this method is derived using the approach described in Section 2. Some properties of the method are investigated and analyzed in Section 4. Finally, we show the algorithm's performance with a few numerical experiments with comparison to some related work.

Theoretical Procedure
Let us consider the stiff initial value problem, on the interval I = [x 0 , x N ], where y and f are assumed to be continuously differentiable and satisfy the conditions to guarantee the existence and uniqueness of solution of the initial value problem.
In order to obtain an approximation to the solution of (3), we propose a Second Derivative Backward Differentiation Formula (SDBDF) with Chebyshev collocation points for the solution of (3) defined by where y n+j ≈ y(x n + jh), y n+v j ≈ y(x n + v j h), f n+j ≡ f (x n + jh, y n+j ) and ) as the off-grid points in the proposed integration formula, the interpolating function for a k-step method is given as The interpolating function y(x) is imposed such that it coincides with the analytical solution at the points [11,20,21].This therefore leads to a (2k + 2) system of equations given by . . .
in matrix form.The system (6) is solved using Gaussian elimination to obtain the coefficients Hence, substituting a j in (5), a continuous second derivative backward differentiation formula with Chebyshev collocation points given as is obtained.The main multistep method is obtained from the continuous multistep method (7) on evaluation at x = x n+k , while its self-starting deficiency is addressed by differentiating (7) which gives to obtain additional methods evaluated at some given points.Equation ( 8) In order to write the block multistep formula as a finite difference scheme, some variables will be defined which shall be necessary for the block representation of the numerical methods.However, the representation in this paper is due to the method of Fatunla [12], Ehigie et al. [10].The 2k-dimensional vector Y m , Y m−1 , F m and G m have collocation points v i 's specified as, The block second derivative BDF with Chebyshev collocation points shall be conveniently represented by a matrix finite difference equation in the block form, where A,B,C and D are 2k × 2k square matrices, m = 1, 2, • • • shall represents the block number.Hence the coefficients A, B, C and D respectively take the forms, Remark 2.1.We note that these class of methods can be applied to systems of ordinary differential equations, where y : Hence a system of 2k × m algebraic equations is solved for to obtain numerical results at a particular block.

Self-Starting Implementation of the Method
To implement the block multistep method without predictors, the interval of integration, [a, b] is partitioned with N ∈ Z equal spacing h such that h = b−a N .Using (9), with x 0 = a, y 0 , n = 0 and m = 1, the first block numerical solution {y v 1 , y 1 , y v 2 , y 2 , • • • , y v k , y k } T are generated simultaneously over the subinterval [x 1 , x k ] of integration.To obtain the numerical solution for the second block for m = 2, n = k, with the previous information y k , This process is repeated for to obtain numerical solutions to (3) on the entire range of integration over the subintervals

A Second Derivative Backward Differentiation Formula with Chebyshev Collocation Points (SDBDFC2)
To derive the continuous SDBDF with Chebyshev collocation points for k = 2, we set k = 2 so that the interpolating function Hence, it is necessary to interpolate (10) and collocate y ′ (x) and y ′′ (x) at x = x n+2 .We obtain a system of equations represented in the matrix form Solving the system of equations, we obtain  (12) Evaluating (12) at x = x n+2 yields the main method, while differentiating (12) and evaluating at The scheme ( 13), ( 14), ( 15) and ( 16) are used together as a simultaneous numerical integrator of the problem (3) to yield Rearranging and representing the block method in matrix finite difference form, we obtain a block 2-step Second Derivative Backward Differentiation Formula with Chebyshev collocation points (SDBDFC2) given by where

Analysis of the Method
The analysis of the Block SDBDF with Chebyshev collocation points for k = 2 (SDBDFC2) ( 17) is presented in this section.Numerical Properties such as Order and Error constant, consistency, stability and convergence are investigated.
Order and Error Constant.Let the individual SDBDF with Chebyshev collocation points (SDBDFC2) be associated with the formula (18) where y(x) is an arbitrary smooth function on [a, b].Expanding (18) with Taylor series expansions of y(x + jh), y(x + v j h), y ′ (x + jh), y ′ (x + v j h) and y ′′ (x + kh), j = 0, v 1 , 1, v 2 , 2, ..., v k , k to obtain the expression where C i are vectors in the form, The SDBDF with Chebyshev collocation points (17) and the associated linear difference operator is said to be of order p if,

Definition 4.2
The term C p+1 is called the Error Constant (EC) and the local truncation error for the method is given by, Hence, from the SDBDFC2 ( 17), the following coefficients are extracted and defined as Using ( 16) -( 22) for q = 0, 1, 2, • • • , it is easily obtained that Hence these methods are of order p = [5, 5, 5, 5] T with error constants Consistency.Since the Block SDBDFC2 is of order p = 5 ≥ 1, therefore it is consistent.Henrici [16].

Convergence
Since the SDBDFC2 is consistent and zero-stable, we can safely assert the convergence of SDBDFC2.(Henrici [16])

Region of Absolute Stability of SDBDFC2
Solving (26) for ξ, we obtain the following roots of the characteristic equation: the Stability function as Solving (29), we obtain the stability region S as,

Definition 4.3: Widlund [29]
A block multistep method is called A(α)−Stable if the angular sector, is contained in the stability region.

Definition 4.4: Gear [14]
A block multistep method is stiffly stable if for some D > 0, it is verified that, and the method be accurate in region D < Rez < δ, |Imz| < θ for some δ > 0.

Test for L 0 -Stability
The SDBDFC2 is A 0 -stable and lim z→−∞ R(z) = 0, we say that SDBDFC2 is L 0 Stable.

Numerical Illustration
In this section, some stiff problems are solved using a Maple code with 20 digits on a digital computer to show the performance of the new method.The nonlinear problems is tackled using the Newton's iteration features of the Maple software.Computational parameters such as Accurate Digit (∆), rate of convergence as well as absolute errors are presented accordingly.Problem 5.1: Stability test of Chartier [6] The stability of the method is compared with L-Stable method of order 5 of Chartier [6] and the Backward Differentiation Formula of order 5 (BDF5) using the problem whose Jacobian matrix J has purely imaginary eigenvalues: with exact solution, It is noted that for any given value of parameter α, J is the matrix, 0 −α α 0 , with eigenvalues of J are i • α and −i • α.For α = 10, we present in Table 1 the accurate digits ∆ , which is defined as, for the following methods defined by acronyms: • M (5, r 5 )-Chartier [6] order at least p = 5 • BDF5-Gear [14] Backward Differentiation Formula of order p = 5 • SDBDFC2-Block Backward Differentiation Formula with Chebyshev Collocation points Remark 5.1: Numerical overflow is indicated by ∞.Table 1 shows that the SDBDFC2 gained some digits and behaves correctly over other methods.

Problem 5.2: A linear stiff problem
Also, the linear system of 3 first order ordinary differential equations solved by Akinfenwa [1], Brugnano and Trigiante [3] and Ramos and Garcia-Rubio [22] given by, on the interval 0 < x < 10 is solved with the newly derived Block Multistep method SDBDFC2.We compare the maximum absolute errors (|y(x) − y n |) on the interval 0 < x < 10 with the Adams Type Block method of Akinfenwa [1] of order p = 7 (ATBM7) and Generalized Backward Differentiation formula of Brugnano and Trigiante [3] (GBDF8) using step lengths h = 1 2 n •100 , n = 0, 1, 2, 3 and 4 for numerical solution of y(x).The order of the methods are also verified by calculating the rate of convergence with the formula where err h is the maximum absolute error at step length h.Also in the range 0 ≤ x ≤ 1, AbsErr(t f ) in [22] is obtained by the SDBDFC2 in comparison with the CBDF 5 of degree s = 5 in Ramos and Garcia-Rubio [22] and the following results are presented  2, it can be seen that the new method even though it is of order p = 5, performs better that the ATBM7 and the GBDF8, both of orders 7 and 8 respectively.Also, the rate of convergence of the SDBDFC2 conforms almost exactly with the order of our methods unlike the ATBM7 and GBDF8.Table 3 shows that the SDBDFC2 is comparable with the CBDF 5 in [22].Numerical results also show that the new method is consistent with order of the method as the step size decreases.

Problem 5.3: Parabolic Partial Differential Equation Cash [5]
A partial differential equation due to Cash [4] and Jator [17], is considered.The problem is a heat equation that describes the heat flow in an insulated thin rod with zero temperature at the edges of the form, with analytical solution, u(x, t) = e −π 2 κt sin πx + e −ω 2 π 2 κt sin ωπx.
First, a semi-discretization known as the central difference formula, is applied to the space variable for the function u(x, t), to obtain an initial value problem of the form, where A is a tridiagonal matrix given by, , and the initial values, Remark 5.3a: The eigenvalues of A given that κ = 1 are, which belong to the range −4 (∆x) 2 , 0 and hence for large N , the system becomes very stiff.Lambert [18].Remark 5.3b: Cash [5] also noted that as ω ≫ 1 the systems of differential equations exhibits some characteristics of some highly stiff models, hence only stiffly stable methods are expected to perform very well on this problem.We compare our numerical results for κ = 1 and ω = 1, 2, 3, 5 and 10 with the Cranck-Nicholson Scheme (C-N), method of Cash [5] 2.13a, b, c (CashABC) and Jator [17] for order p = 5 (HAMC5, HAML5, HAMS5) and numerical results are presented in Table 4. Numerical result reveals the consistency of our method as ω increases.
Problem 5.4: Cash [4].Lastly, the integration of the system using the problem whose Jacobian matrix J has imaginary eigenvalues given by For the case α = 1 and β = 15 with a fixed step size h = 0.25, Table 5 presents the results obtained by the SDBDFC2 in comparison with the results in Cash [4] for conventional second derivative linear multistep methods (C2MM) and the second derivative extended backward differentiation formulas (E2BD).
Table 5 shows clearly that the new method SDBDFC2 on implementation on Problem 5.4 compares favourably with the methods E2BD and C2MM as obtained in Cash [4].

Conclusion
In this paper, a procedure for deriving a continuous second derivative backward differentiation formula with Chebyshev collocation points has been introduced.Investigation of some numerical properties of a class of these methods shows that the methods are uniformly accurate of order p = 2k + 1 and further analysis reveals that the method possesses some stiff stability properties such as A 0 -stability and A(α)−stability with good region of absolute stability.An implementation procedure is also presented to demonstrate the self-starting properties of the method.Implementation on some stiff problems shows that a class of order p = 5 performs better than methods of order p = 8 in the literature, and could be competitive with some stiff codes for solving stiff problems.

Acknowledgement
The Authors are particularly grateful to Late Prof. A. B. Sofoluwe and Prof. S. N. Jator for their many constructive comments and suggestions.

Figure 1 :
Figure 1: Diagram showing the stiffly stable characteristics

Figure 2 :
Figure 2: RAS of SDBDF for k = 2, RAS near the origin and RAS showing the zeros and poles of R(z)

Table 1 :
Problem 5.1: Table of Accurate Digits ∆ with α = 10 for methods of Order p = 5

Table 3 :
Problem 5.2: Numerical Results in comparison with CBDF 5 in the range 0 ≤ x ≤ 1