Abstract

A four-step Continuous Block Hybrid Method (CBHM) with four non-step points of order (9, 9, 10, 9, 9, 9, 9, 9)^T is proposed for the direct solution of the 1st order Initial Value Problems (IVPs). The main method and additional methods are obtained from the same continuous scheme derived via interpolation and collocation procedures. The stability properties of the methods are discussed and the stability region shown. The methods are then applied in block form as simultaneous numerical integrators over non-overlapping intervals. Numerical results obtained using the proposed block form shows that it is attractive for solutions of stiff problems.

---

INTRODUCTION

Consider the IVP:

(1)

We seek a solution in the range a≤x≤b where, a and b are finite and we assume that f satisfies the conditions which guarantee that the problem has a unique continuously differentiable solution which we shall indicate by Y (x). Consider the sequence of points {x_n} defined by x_n = a+nh, n = 0, 1, 2,…, b-a/h where, the parameter h≥0 is a constant step-size. An essential property of the majority of computational methods for the solution of Eq. 1 is that of discretization that is we seek an approximate solution, not on the continuous interval a≤x≤b but on the discrete point set {x_n}. The k-step Linear Multistep Method (LMM) for the solution of Eq. 1 is generally written as:

(2)

Which has 2k+1 unknown α’s and β’s and therefore can be of order 2k where, k is the step number. But according to Dahlquist (1963), the order of Eq. 2 cannot exceed k+1 (k is odd) or k+2 (k is even) for the method to be stable. Several researchers such as Lambert (1973), Gear (1965), Gragg and Stetter (1964), Butcher (1965) and Brugnano and Trigiante (1998) proposed modified forms of Eq. 2 which were shown to overcome the Dahlquist barrier theorem. These methods known as Hybrid methods were obtained by incorporating at off-step points in the derivation process. We define a k-step continuous hybrid formula to be of the type:

(3)

Where, h is the stepsize, v is the number of off-points, α_k = 1, α_j (x), β_j (x), βη_j (x) are continuous coefficients which are uniquely determined from the derivation process. Hybrid methods were initially proposed to overcome the Dahlquist barrier theorem for solving 1st order IVPs. Hybrid methods have been considered by Gragg and Stetter (1964), Lambert (1973) and Kohfeld and Thompson (1967) and are of the form:

(4)

Which were shown to be of order up to 2k+2. Gupta (1978) and Jator (2010) noted that the design of algorithms for Eq. 4 is more tedious due to the occurrence of f_n+1-r in Eq. 4 which increases the number of predictors needed to implement the method. Donelson and Hansen (1971) proposed cyclic composite multistep methods which were also shown to circumvent the Dahlquist barrier theorem. The methods were applied in a predictor-corrector mode to scalar 1st order IVPs by combining different correctors and writing them as a single matrix difference equation. It was reported by Donelson and Hansen (1971) that stability was achieved with higher order methods without additional function evaluations however, an extra amount of programming was required. The Hybrid method proposed in this study is self-starting and implemented without the use of predictors. The derivation of the method is based on interpolation and collocation (Lie and Norsett, 1989; Atkinson, 1989; Onumanyi et al., 1994; Gladwell et al., 1980).

Block methods were first introduced by Milne (1953) for use only as a means of obtaining starting values for predictor-corrector algorithms and has since then been developed by several researchers (Rosser, 1967; Sarafyan, 1965; Shampine and Watts, 1969) for general use. We emphasize that the hybrid methods are developed for general use, not only as a means of obtaining starting values for predictor-corrector algorithms.

To this end, the continuous representation generates a main discrete hybrid method and four additional methods which are combined and implemented as a block method which simultaneously generates approximations y_n+4-jr to the exact solution, y(x_n+4-jr), j = 0, 1, 2, ..., 8 and r = 1/2. Without loss of generality, y_n+9/2, y_n+5, y_n+11/2, y_n+6, y_n+13/2, y_n+7, y_n+15/2, y_n+8 is obtained in the next block using y_n+4 as the starting value. Thus, the order of the algorithm is maintained. In this study, the aim is to generate a four step Continuous hybrid block method with four off-step points and to demonstrate the efficiency in its implementation on stiff and non stiff ODEs.

Derivation of the method: In this study, the objective is to derive the main Hybrid method of the form:

(5)

α_j, β_j, βη_j are coefficients and η_j is chosen from the interval (0, 4). In order to obtain Eq. 5, we proceed by seeking an approximation of the exact solution y(x) by assuming a continuous solution Y(x) of the form:

(6)

Such that xε[x_o, X], b_j are unknown coefficients and φ_j (x) are polynomial basis function of degree p+q-1 where, the number of interpolation points p and the number of distinct collocation points q are respectively chosen to satisfy 1≤p<k and q>0. The integer k≥1 denotes the step number of the method. We thus construct a k-step continuous multistep method with φ_j (x) = x^j, j = 1, 2, ..., 9, η_i = {1/2, 3/2, 5/2, 7/2}, i = 1, 2, 3, 4, p = 2, q = 8, k = 4 by imposing the following conditions:

(7)

(8)

For r∈(0, 1) then y_n+4-jr is the numerical solution for the exact solution y(x_n+4-jr), f_n+4-jr = f(x_n+4-jr, y_n+4-jγ) and n is the grid index. It should be noted that Eq. 7 and 8 leads to a system of p+q equations which must be solved to obtain the coefficient b_j. The four step Continuous hybrid method is obtained by substituting these values of b_j into Eq. 6. After some algebraic computation, the method yields the expression in the form (Eq. 3) as:

(9)

Which is used to generate the main discrete Hybrid method at the interpolation point x = x_n+4. And we obtain the additional methods by interpolating Eq. 9 at the points x = {x_n, x_n+1/2, x_n+1, x_n+3/2, x_n+2, x_n+5/2, x_n+7/2}. The coefficients of the CBHM are shown in Table 1-3. The method is implemented simultaneously by applying the CBHM to generate y_n+1/2, y_n+1, y_n+3/2, y_n+2, y_n+5/2, y_n+3, y_n+7/2, y_n+4 for the solution of Eq. 1 at discrete block points x_n+1/2, x_n+1, x_n+3/2, x_n+2, x_n+5/2, x_n+3, x_n+7/2, x_n+4 n = 0, 4, ..., N-4 on the partition [a, b].

Summary: On the partition I_N: a<x₀<x₁<...<x_N-1<x_N = b, h = b-a/N, n = 0, 1, 2, ..., N-1. In step 1, choose N for k = 4, the number of blocks II = N/4, using Eq. 13, n = 0, = 0, the values (y_n+1/2, y_n+1, y_n+3/2, y_n+2, y_n+5/2, y_n+3, y_n+7/2, y_n+4)^T are generated simultaneously over the subinterval [x₀, x₄] as y₀ is known from the IVP Eq. 1.

Table 1:	Coefficients βj/2 j = 0, ..., 2k for the method in Eq. 9 evaluated at x = x n+j/2, j = 0, ..., 2k, j≠6

 

Table 2:	Absolute error for Block hybrid method in Eq. 14 with four off-steps at the end point T = 10

 

 

Table 3:	Absolute error for Block hybrid method Eq. 14 with four off-steps at the end point T = 10 and ε = 10-4

 

In step 2, for n = 4, j = 1, (y_n+9/2, y_n+5, y_n+11/2, y_n+6, y_n+13/2, y_n+7, y_n+15/2, y_n+8)^T are obtained over the sub-interval [x₄, x₈] since, y₄ is known from the 1st block. In step 3, the process is continued for n = 8, ..., N-4, = 2,... II to obtain approximate solutions to Eq. 1 on sub-intervals [x₀, x₄], ..., [x_N-4, x_N]. We note that for linear problems, we solve Eq. 1 directly from the start with Gaussian elimination using partial pivoting and for non-linear problems, we use a modified Newton-Raphson method.

Order of accuracy and local truncation error: Following Fatunla (1991) and Lambert (1973), we define the local truncation error associated with Eq. 5 to be the linear difference operator:

(10)

Assuming that y(x) is sufficiently differentiable, we can write the terms in Eq. 10 as a Taylor series expression of y (x_n+j), f (x_n+j) = y’ (x_n+j) and f (x_n+ηj) = y’ (x_n+ηj) as:

(11)

Substituting Eq. 11 into equations in Eq. 10, we obtain the expression:

(12)

Where, the constant C_m, m = 0, 1, 2, ... are given as follows:

(13)

According to Henrici (1962), we say that the method in Eq. 5 has a maximal order of accuracy m if L [y (x); h] = O (h^m+1) hence, C₀ = C₁, C₂, ..., C_m = 0, C_m+1 ≠0 therefore, C_m+1 is the error constant and C_m+1 h^m+1 y^m+1 (x_n)+O (h^m+2), the principal local truncation error at the point x_n.

It is established from the calculations that the continuous block hybrid methods shown in Table 2 have high order and error constants.

Analysis of CBHM: The methods shown in Table 2-4 can be represented by a matrix finite difference equation in block form:

(14)

Where:

 

 

Where, j = 0, 1, 2, ... and n is the grid index.

Zero stability: The zero stability of the methods in Eq. 14 are determined as the limit as h→0, the difference system Eq. 14 tends:

(15)

Which is normalized to obtain the 1st characteristics polynomial ρ (R) given by:

(16)

Where A₀ is the identity matrix in I and A_l = A. Following Eq. 6, the block method (Eq. 14) is zero stable since, ρ (R) = 0 and R_j≤1, j = 0, 1, ..., 8. Thus for this method one of the root R₈ = 1 which is the principal root hence, the method is zero stable.

Consistency and convergence: Since each of the method in Table 1 has order m>1, the CBHM is consistent and by Henrici (1962). Convergence = zero stability+ consistency. Hence, the CBHM is convergent.

Stability analysis: The stability properties of the block formula (Eq. 14) is discussed according to Butcher (1965) and are determined through the application to the test equation:

(17)

which yield the stability polynomial:

(18)

From Eq. 18, we obtain the usual property of A-stability which requires that for all z = hλεC¯ and Re (z) < 0, P (t, z) must have a dominant eigenvalue t₈ such that |t₈|<1. From the analysis, we have that the eigenvalues {t_1-8} = {0, 0, 0, 0, 0,..., 0, 0, t₈) and the dominant eigenvalue, t₈ is a function of z given by:

(19)

Clearly from Eq. 19, Re (z)<0, t₈<1. Hence, the block method in Eq. 14 is A-stable since, its region of absolute stability contains the left half-plane {z C| Re (z)<0}.

Numerical experiment: This study deals with some numerical experiments, executed in Matlab language with double precision arithmetic which illustrate the result derived in the previous sections.

Example 1: Consider the stiffly linear problem:

The eigenvalues of the systems are λ₁ = -10000 and λ₂ = -1 with exact solution:

Example 2: Consider the linear SPP (Xiao et al., 2000):

 

Table 4:	Absolute error for Block hybrid method in Eq. 14 with four off-steps at the end point T = 10

 

 

Table 5:	Absolute error for Block hybrid method in Eq.14 with four off-steps at the end point T = 1

 

 

Table 6:	Absolute error for Block hybrid method in Eq. 14 with four off-steps at the end point T = 2

 

With the exact solution:

It is easy to see that the exact solution y₁ (t, ε) and y₁ (t, ε) approximate to the constants 2/2 ε+1 (ε y₁ (0)+y₂ (0)) and 1/2ε+1 (ε y₁ (0)+y₂ (0)), respectively after a short time.

Example 3: Consider the Stiffly nonlinear Kap’s problem:

Take ε = 10^-10, the smaller ε is the more serious, the stiffness of the system. The exact solution is:

Example 4: The last experiment is the system:

This problem has been extensively studied by Shampine (1977) and reported that the system is stiff when β = 21, γ = 10 and the Jacobian has eigenvalues -10±βi and -γ. Its exact solution is:

We compute the absolute error at the end point for T = 1 and 2 as shown in Table 5 and 6, respectively.

CONCLUSION

A continuous hybrid formula with four off-step points has been proposed and implemented as a self starting method for stiff ordinary differential equation. The strength of the method lies in the additional methods generated from the main method which are simultaneously implemented in block form without the need for predictors. The good convergent and stability properties of the method makes it attractive for numerical solution of stiff problems. The accuracy of the block method have been demonstrated on both linear and non-linear problems as shown in Table 4-6. Although for example 1, the methods showed poor accuracy for h = 0.1 but the accuracy improved greatly as the step size reduces.

ACKNOWLEDGEMENTS

We are very grateful to the State Key Laboratory of High-end Server and Storage Technology. Number: 009HSSA08 and Fundamental Research Funds for the Central Universities. Numbers: HEUCFT1007, HEUCF100 607 China for the their support in this research. We are also grateful to the referees whose useful suggestions greatly improve the quality of this research.

How to cite this article:

N.M. Yao, S.N. Jator and O.A. Akinfenwa. A Linear Multistep Hybrid Methods with Continuous Coefficient for Solving Stiff Ordinary Differential Equation.
DOI: https://doi.org/10.36478/jmmstat.2011.47.53
URL: https://www.makhillpublications.co/view-article/1994-5388/jmmstat.2011.47.53