Implementing a seventh-order linear multistep method in a predictor-corrector mode or block mode: which is more efficient for the general second order initial value problem

Abstract A Seventh-Order Linear Multistep Method (SOLMM) is developed and implemented in both predictor-corrector mode and block mode. The two approaches are compared by measuring their total number of function evaluations and CPU times. The stability property of the method is examined. This SOLMM is also compared with existing methods in the literature using standard numerical examples. AMS Subject Classification 65L05; 65L06


Introduction
Linear multistep methods (LMMs) of the form k j=0 α j y n+j = h 2 k j=0 β j f n+j , k ≥ 2, have been extensively applied to solve the special second order initial value problem (IVP) y = f (t, y), y(t 0 ) = y 0 , y (t 0 ) = y 0 , t [t 0 , t N ] on the discrete set of points t n = t 0 + nh, n = 0, ..., N, h = t N −t 0 N , (see Lambert and Watson (1976), Ramos and Vigo-Aguiar (2005), Ixaru and Berghe (2004). Despite the successful application of (1) to solving problems of the form (2), fewer methods of the form (1) have been proposed for solving the general second order IVP y = f (t, y, y ), y(t 0 ) = y 0 , y (t 0 ) = y 0 . ( Some of the methods available for directly solving (3) are due to Awoyemi (2001) and Ramos and Vigo-Aguiar (2006). These methods are generally implemented in a step-bystep fashion in a predictor-corrector mode.
In this paper, we construct the continuous form of (1) which has ability to generate several methods which are combined and implemented in block form to solve (3) directly (see Jator and Li (2009) and Jator (2012Jator ( , 2010Jator ( , 2007. http://www.springerplus.com/content/3/1/447 The paper is organized as follows. In Section 'SOLMM', we derive a continuous approximation which is used to obtain the discrete methods that are combined to form the block method. The analysis and computational aspects of the SOLMM is given in Section 'Implementation of the SOLMM'. Numerical examples are given in Section 'Numerical examples' to show the accuracy and efficiency of the method. Finally, the conclusion of the paper is discussed in Section 'Conclusion' .

Continuous form
On interval t n ≤ t ≤ t n+6 , the exact solution to (3) is approximated by the continuous form of the SOLMM whose first derivative is given by where α 0 (t), α 1 (t), and β j (t), j = 0, 1, 2 are continuous coefficients that are uniquely determined. We assume that y n+j is the numerical approximation to the analytical solution y(t n+j ), y n+j is an approximation to y (t n+j ), and f n+j = f (t n+j , y n+j , y n+j ), j = 0, 1, . . . , 6 is supplied by the differential equation. The coefficients of the method (4) are specified by the following theorem.
Theorem 1. In order to obtain the coefficients of the continuous method (4), a nine by nine system is solved with the aid of Mathematica by demanding that the following conditions are satisfied u t n+j = y n+j , j = 0, 1, u t n+j = f n+j , j = 0, 1, . . . , 6.
After some algebraic manipulations, the equivalent form (6) produces the coefficients of (4) whose first derivative is given by (7), where we define the matrix W as , http://www.springerplus.com/content/3/1/447 and W j is obtained by replacing the j th column of W by V ; P j (t) = t j , j = 0, . . . , 8 are basis functions, and V is a vector given by V = (y n , y n+1 , f n , f n+1 , . . . , f n+6 ) T . We note that T is the transpose.

Discrete by-products
The following methods which are used to construct the block form are obtained by evaluating (4) and (5) The derivatives are given by

Block form
The methods (8) and (9) are combined and expressed in the form and A 0 , A 1 , B 0 , and B 1 are matrices of dimension 12 whose entries denoted by α j = α i,j , β j = β i,j , i = 1, . . . , 12 are given by the coefficients of (8) and (9).

Definition 2. The block method (10) is said to be consistent if it has order at least one.
The block method (10) has order and error constant given by the vector p = 6 and C p+2 =

Linear stability of the SOLMM
The linear-stability of the SOLMM is discussed by applying the method to the test equation y = λy, where λ is expected to run through the (negative) eigenvalues of the Jacobian matrix ∂f ∂y (see Sommeijer (1993)). Letting q = λh 2 , it is easily shown that the application of (10) to the test equation yields where the matrix M(q) is the amplification matrix which determines the stability of the method. Sommeijer (1993)).

Implementation of the SOLMM
The SOLMM was implemented in both block mode and predictor-corrector mode using a written code in PERL programming language and executed on a laptop computer with AMD Quad-Core A10-4600M Processor, 8GB of RAM and Windows 8.1 OS. The total program running time was acceptable, as shown in Tables 1, 2 and 3. The computational time complexity and space complexity of the algorithms for both modes of SOLMM used for the examples in this paper are polynomial. Details of the block mode implementation is given in Jator (2012) and the predictor-corrector implementation is discussed next.

Predictor-corrector mode algorithm
The initial block was used to start predictor-corrector algorithm, after which the predictor (14) and corrector (15) were used in a step-by-step fashion to provide the numerical solution from the second block to the end of the interval.
Correctors. The last members of (8) and (9) are used as correctors.

Comparison of block mode and predictor-corrector mode
The SOLMM is implemented in both predictor-corrector and block modes. The two approaches are compared by measuring their total number of function evaluations (NFEs) and CPU times in seconds. The block mode implementation is shown to be superior to the predictor-corrector mode implementation in terms of accuracy and the number of function evaluations. However, the predictor-corrector mode implementation uses less time than the block implementation. Details of the numerical examples are displayed in Tables 1, 2 and 3.

Comparison of block method with other methods
The theoretical solution at t = 8 is y(8) = 2 8π sin (8) 0.279092789108058969. The errors in the solution were obtained at t = 8 using the SOLMM of order 7 and compared the the errors in (Vigo-Aguiar and Ramos 2006) which is based on the variable-step Falker method of order eight (VAR (8)) implemented in the predictor-corrector mode. The results given in Table 4 show that the SOLMM is more accurate than the method in (Vigo-Aguiar and Ramos 2006).
The maximum norm of the global error for the y-component is given in the form 10 −CD , where CD denotes the the number correct decimal digits at the endpoint (see (Sommeijer 1993)). This problem has also been solved in (Sommeijer 1993) using the eighth-order, eight-stage RKN (H8) method constructed by Hairer (1977). We have chosen to compare this method of order 8 with our method of order 7, because the orders of the methods are very close. The results obtained using the H8 are reproduced in Table 5 and compared  with the results given by our method. It is seen from Table 5 that our method performs better than those in (Sommeijer 1993) in terms of accuracy (smaller errors) and efficiency (smaller NFEs).

Conclusion
A SOLMM is proposed and implemented in both predictor-corrector and block modes. It is shown that the block mode algorithm is superior to the predictor-corrector mode algorithm in terms of accuracy and the number of function evaluations. However, the predictor-corrector mode implementation uses less time that the block implementation. the Details of the comparison of the numerical examples are displayed in Tables 1, 2, 3, 4 and 5. Our future research will be focus on developing a variable step version of the SOLMM in both modes.