On contractivity-preserving 2-and 3-step predictor-corrector series for ODEs

New optimal, contractivity-preserving (CP), d-derivative, 2and 3-step, predictor-corrector, Hermite– Birkhoff–Obrechkoff series methods, denoted by HBO(d, k, p), k = 2, 3, with nonnegative coefficients are constructed for solving nonstiff first-order initial value problems y′ = f (t, y), y(t0) = y0. The upper bounds pu of order p of HBO(d, k, p), k = 2, 3 methods are approximately 1.4 and 1.6 times the number of derivatives d, respectively. Their stability regions have generally a good shape and grow with decreasing p − d. Two selected CP HBO methods: 9-derivative 2-step HBO of order 13, denoted by HBO(9,2,13), which has maximum order 13 based on the CP conditions, and 8-derivative 3-step HBO of order 14, denoted by HBO(8,3,14), compare well with Adams–Cowell of order 13 in PECE mode, denoted by AC(13), in solving standard N-body problems over an interval of 1000 periods on the basis of the relative error of energy as a function of the CPU time. They also compare well with AC(13) in solving standard N-body problems on the basis of the growth of relative error of energy and 10000 periods of integration. The coefficients of selected HBO methods are listed in the appendix.

The Taylor series methods have been an excellent choice in astronomical calculations [3], numerical integration of ordinary differential equations (ODEs) and differential algebraic equations (DAEs) [1], sensitivity analysis of ODEs/DAEs [2], in solving general problems [5] and validating solutions of ODEs by means of interval analysis [18,11].
The main cost in solving ODEs by the Taylor method of order p (T(p)) lies in the repeated evaluation of the p Taylor coefficients of the functions involved.
Deprit and Zahar [6] showed that such recursive computation is very effective in achieving high accuracy, even with little computing time and large step sizes.
The current investigation and the results are offered as potentially useful additions to the contemporary repertory of numerical integrators.This paper explores an alternative way to improve the stability and the order of predictor-corrector methods.In our construction of HBO(d, k, p), we replace the forward Euler (FE) method, used by Gottlieb et al. and Huang [7,12] in establishing strong stability preserving (SSP) Runge-Kutta (RK) methods as convex combinations of FE methods, by rewriting HBO(d, k, p) as a convex combination of the special d-derivative extension of FE, which we denote by S(d): where the coefficients η m satisfy the inequality η m ≤ 1 m! .If equality holds, then S(d) reduces to the Taylor method of order d, T(d).The error in S(d) is of order ≥ 2 if there exists a smallest ∈ {2, 3, . . ., d} such that η < 1  ! .If S(d) is contractive in a given norm, then HBO(d, k, p) will be contractive as a convex combination of S(d) with modified step sizes.
The region of absolute stability of HBO(d, k, p) is derived under the assumption that two solutions, y and ỹ, to problem (1) are contractive: y(t + ∆t) − ỹ(t + ∆t) ≤ y(t) − ỹ(t) , ∀ ∆t ≥ 0. (4) We assume that there exists a maximum stepsize ∆t S(d) such that f satisfies a discrete analog of (4) when S(d) is employed with ∆t ≤ ∆t S(d) : Here y n and ỹn are two numerical solutions generated by S(d) with different neighbouring starting (or previous) values y 0 = y(t 0 ) and ỹ0 = ỹ(t 0 ).
By interpreting ỹ0 as a perturbation of y 0 due to numerical error, we see that contractivity implies that these errors do not grow as they are propagated.
We are interested in a higher-order HBO(d, k, p) that maintains the contractivity-preserving property for 0 ≤ ∆t ≤ ∆t max = c∆t S(d) whenever inequality (5) holds.Here c, called the CP coefficient, depends only on the numerical integration method but not on f .This definition of the CP coefficient of HBO(d, k, p) follows closely the definition of the SSP coefficient of RK (see [7]).
In [15], similar CP RK methods have been constructed and tested on DETEST problems [14].The aim of HBO(d, k, p) is to maintain the CP property (6) while achieving higher-order accuracy, perhaps with a modified time-step restriction, measured here with the CP coefficient c(HBO(d, k, p)): This coefficient describes the ratio of the maximal HBO(d, k, p) time step to the time step ∆t S(d) , for which condition (5) holds.Similar to Huang et al [13], we, first, compare the numerical performance of HBO (9,2,13), HBO (8,3,14) and Adams-Cowell of order 13 in PECE mode, denoted by AC (13), on Kepler circular orbit with eccentricity e = 0.3, e = 0.5 and e = 0.7 over an interval of 1000 periods on the basis of log 10 (EE) as a function of CPU time.It is seen that HBO(9,2,13) and HBO (8,3,14) win.Next, these two HBO methods compare well with AC (13) in solving eccentric Kepler orbit with eccentricity e = 0.3, e = 0.5 and e = 0.7 over an interval of 10000 periods on the basis of the growth of relative error of energy and long intervals of integration.

Predictor-corrector CP Hermite-Birkhoff-Obrechkoff methods
All HBO methods considered in this paper are CP, so the denomination "CP" will generally be omitted in what follows.
The k-step, predictor-corrector HBO(d, k, p) are constructed, as a subclass of multiderivative, multistep, multistage methods, by the following two formulae which perform integration from t n to t n+1 .Let ∆t denote the step size.The abscissa vector [c 1 , c 2 ] T defines the two off-step points t n + c j ∆t, j = 1, 2. In all cases c 1 = 0 and, by convention, c 0 1 = 1.Let F 1 = f n .A Hermite-Birkhoff (HB) polynomial is used as stage formula P 2 to obtain the stage value Y 2 to order p − 1, A HB polynomial is used as the integration formula to obtain y n+1 = Y 3 to order p, where F 2 := f (t n + c 2 ∆t, Y 2 ) denotes the stage derivative.Here the subscript BU refers to the Butcher form, while the subscript SO will be used later for the Shu-Osher form.
Formulae ( 8)-( 9) are the Butcher form of HBO(d, k, p).One sees that the derivatives y n , m = 2, 3, . . ., d, are computed only once per step at t = t n .The defining formulae of the 2-and 3-step HBO(d, k, p) involve the usual RK parameters c i , a i,j and b j and the Taylor expansion parameters γ , j,m , = 0, 1, . . ., k − 1.Thus we can represent HBO(d, k, p) by its coefficient scheme matrices.One can display the coefficient scheme (A, b, γ 0 , γ 1 , . . ., γ k−1 ) and the c i in the Butcher tableau and the 3

Order conditions for predictor-corrector HBO(d, k, p)
To derive the order conditions of HBO(d, k, p), we shall use the following expressions coming from the backsteps of the method.First, we need to satisfy the set of consistency conditions: Next, we have to solve the following equations for a 21 , c 2 , b 1 , b 2 : where the backstep parts, B 2 (j) and B( j), are defined by Similar to the generalization to HB methods [20] of the Shu-Osher form of RK methods [24], equations ( 8)-( 9) can be written in the form with consistency conditions: Form ( 19) is called the Shu-Osher form of HBO(d, k, p).By setting v i = α i1 and w i = β i1 , i = 2, 3, in (19), we have the modified Shu-Osher form of HBO(d, k, p): This form generalizes the modified Shu-Osher form for RK methods (see [8]).We can rearrange the stage values Y i , i = 2, 3, in (21) as follows: with consistency conditions: Thus we obtain the difference Y i − Y i , i = 2, 3, from (22) as follows: Provided all the coefficients of ( 22) are nonnegative, the following straightforward extension of a result presented in [9,12] holds.
Theorem 1.If f satisfies condition (4) of the S(d) method, then the k-step, predictor-corrector HBO(p) method (22) satisfies the CP property where • the feasible CP coefficient, c feasible , is the minimum of the following numbers: • the following conditions are imposed on δ ,i,m : under the assumption that all coefficients of ( 22) are nonnegative, with the convention that the ratios a/0 = +∞, and the ratios 0/0 are ignored.
Proof.The difference Y i − Y i of HBO(d, k, p) can be rewritten as a convex combination of the three terms on the right-hand side of (24).Thus, by the convexity of the norm • , we have The first two terms on the right-hand side of ( 27) have the following upper bounds: The third term on the right-hand side of (27) has the following upper bound: Because v i , k−1 =1 δ ,i,0 and α i j are nonnegative and We thus obtain It is to be noted here that each representation of ( 22) with coefficients v i , w i , α i j , β i j , and δ ,i,m will produce a feasible CP coefficient, c feasible , defined in Theorem 1 and a feasible HBO(d, k, p) in modified Shu-Osher form (22). What we really want is not merely a feasible method HBO(d, k, p) in Shu-Osher form but a best HBO(d, k, p).This question will be considered in Subsection 4.4.
To find the relation between the Shu-Osher coefficients and the Butcher coefficients, we can solve (29) for Y since I − α is invertible because the matrix α is strictly lower triangular, Comparing ( 35) with (33), we have the following relations between the generalized Shu-Osher coefficients and the Butcher coefficients, These relations allow a simple transformation of the vectors and matrices v, w, δ 0 , δ , = 1, 2, . . ., k − 1, β of a Shu-Osher form into v 0 , w 0 , γ 0 , γ , = 1, 2, . . ., k − 1, β 0 of a Butcher form and vice versa.

Canonical Shu-Osher form of HBO(d, k, p) in vector notation
To find the CP coefficient of HBO(d, k, p), it is useful to consider a particular modified Shu-Osher form (29) where the elements of the matrices α and β satisfy the ratios or, in vector notation, Substituting this relation into (37), we can solve for β r in terms of β 0 and r.First, we have Then, since I + rβ 0 is invertible, the coefficients of the Shu-Osher form (29) are given in terms of the coefficients of the Butcher form (33) by the expressions It is to be noted that using (37) and (47), β r can then be written as As in [8], we shall refer to the form given by the coefficients (42)-(47), as the canonical Shu-Osher form of HBO(d, k, p): which can be written solely in terms of the vectors and matrices v 0 , w 0 , γ 0 , γ , β 0 of the Butcher form (33), The canonical Shu-Osher form (49) and the Shu-Osher form (50), will allow us to formulate simply the optimization problem considered in Subsection 4.4.
Relations (42)-(47) will enable us to obtain simply the vectors and matrices of a canonical Shu-Osher form (49) from those of a Butcher form (33) and vice versa.
To simplify notation, in the following theorem, the ratio r = α 32 β 32 becomes a feasible CP coefficient of HBO(d, k, p).Hence, 1. r must satisfy the conditions: 2. conditions (26) are imposed on δ ,i,m , = 0, 1, . . ., k−1.To enhance the performance of the optimization software, conditions (26) can be rewritten as Therefore, the following slight modification of Theorem 1 holds.
In the optimization formulation with any feasible initial data, the ratio r becomes the variable r which satisfies the following nonlinear equation in the variables α 32 , r, β 32 , together with conditions (53) and (54).
We can now formulate the optimization problem, using the Shu-Osher form (50) which is written solely in terms of the vectors and matrices of the Butcher form.Hence, the problem of optimizing HBO(d, k, p) can be formulated as, subject to the componentwise inequalities together with conditions (53), (54) and order conditions ( 12) to ( 14) for order p.

Conditions to obtain CP predictor-corrector HBO(d, k, p) series
We consider the following formula obtained from formula (22) using the relations r = α i j β i j , r i0 = v i w i and r i = δ ,i,0 δ ,i,1 from (40) and (53), The first two terms and the third term on the right hand side of (62) form the Taylor part and Runge-Kutta part respectively of HBO methods.
From formula (62), we can form the maximum stepsize ∆t RK and ∆t T of Runge-Kutta part and Taylor part, respectively, When HBO(d, k, p) methods are applied to the test problem y = −y, ∆t RK and ∆t T of ( 63) and ( 64) follow approximately the conditions: where x min,T (d), the lower end of the interval of stability of T(d), Taylor method of order d, for all d, follows the formula given by [3].
To guarantee the property that the stability intervals of HBO(d, k, p) methods grow as order p increases, like Taylor methods, ∆t RK must exceed substantially ∆t T (and, hence, |x min,T (d)|) as follows: where the security factor Θ can be between 2 to 3, and, consequently, from (65) and ( 68 condition (54) and order conditions ( 12)-( 14) for order p.
Since HBO(d, k, p) contains many free parameters, the MATLAB Optimization Toolbox was used to search for the methods with largest c(HBO(d, k, p)) under the tolerance 10 −12 on the objective function c(HBO(d, k, p)), provided all the constraints are satisfied to tolerance 8 × 10 −14 .For a given d, Table 1 1.As an example, Fig. 1 depicts the upper bounds p u of CP HBO(d, 3, p) methods as an increasing function of d, number of derivatives and shows that the upper bound p u of order p seems to be linear with the number of derivatives d and a linear least square fit is given in Table 2 where the slope η = 1.6190 is expected, since an increase of order p by one requires two additional order conditions while an increase of d by one yields 6 additional variables.Here, for any positive d, lower bounds x min of the unscaled intervals of stability (x min , 0) of HBO(d, k, p u )) satisfy generally the relation |x min | ≥ 2 × c(HBO(d, k, p u )) from condition (7) and, hence, CP HBO(d, k, p)) are absolutely stable for d = 1 to infinity.This analysis suggests that, for large d, HBO(d, k, p) methods have order p large enough to take into account many problems which require a very high precision of the solution, similar to Taylor methods which have the slope η = 1.

Existence and stability properties
Since p u follows the equations shown in Table 2, when we slow down the growth of order p as d increases as follows: HBO(d, k, p) method has desirable stability properties considered in Subsection 5.3.Here, η max > 1 should not be much larger than 1.From Table 4, η max seems to be about 1.4 or less.

Regions of absolute stability of HBO(d, k, p)
To obtain the regions of absolute stability, R, of HBO(d, k, p), we apply the stage formula P 2 and the integration formula with constant step h = ∆t to the linear test equation Thus we obtain and It is seen that Y 2 in ( 73) is expressed only in terms of y n− , = 0, 1, . . ., k − 1.Then, y n+1 in ( 74) is expressed only in terms of y n− , = 0, 1, . . ., k − 1; thus we obtain the following second-order difference equation and associated linear characteristic equation: A complex number ĥ = λh is in R if the k roots of the characteristic equation satisfy the root condition |r s | ≤ 1, s = 1, 2, . . ., k, provided the multiple roots satisfy |r s | < 1.The scanning method used to find R is similar to the one used for Runge-Kutta methods (see [16, pp. 70 and 204]).As an example, the grey regions in Fig. 3 depict R in the half-plane Im( ĥ) > 0 for two selected HBO methods: HBO(9,2,13) and HBO (8,3,14).Here the lower ends of the intervals of stability, denoted by x min , of HBO(9,2,13) and HBO (8,3,14) are x min = −3.00 and x min = −2.08 respectively.

Stability properties of HBO methods for a given step number k
In this section, we analyze some stability properties and list the principal error term of CP HBO(d, k, p) methods with p = η d and η satisfying (72).Similar to the case of Taylor series methods, the use of high number d gives HBO series methods of high order.
We are now interested in the size of the stability regions of CP HBO(d, k, p) methods.The percentage stability gain (PSG) of method 2 over method 1 is defined by the formula (cf.Sharp [23]), where x min (d, 1) and x min (d, 2) are the lower end of the interval of stability of methods 1 and 2, respectively.As an example, Table 3   6.Regions of absolute stability and principal error terms of two selected methods: HBO(9,2,13) and HBO(8,3,14) Using the method described in Subsection 5.2, we obtain stability regions R in the half-plane Im( ĥ) > 0 for two selected HBO methods: HBO(9,2,13) and HBO (8,3,14).These stability regions R are depicted by the grey regions in Fig. 3. Figure 3: Grey unscaled regions of absolute stability, R, in the half-plane Im( ĥ) > 0 of two selected HBO methods with intervals of absolute stability (−3.00, 0) and (−2.08, 0).

Numerical results
Since HBO(d, k, p) is not a one-step method, we must provide not only an initial value, i.e., y 0 but also k − 1 additional starting value, i.e., y 1 , y 2 , . . ., y k−1 .The starting values for HBO(d, k, p) are calculated by the one-step, 4-stage, Hermite-Birkhoff-Taylor method of order d + 3 using y to y (d) with appropriate small step sizes [19].The d derivatives, y to y (d) , of the Taylor series are calculated at each integration step by known recurrence formulae (see, for example, [10, pp. 46-49], [17]).

CPU time of HBO
where E(t) is the energy at time t.
Our first result is a comparison of the relative energy error (EE(t)) as a function of CPU time of HBO (9,2,13), HBO (8,3,14) and AC(13) after a 1000 periods integration of a Hamiltonian system (Huang and Innanen [13]).For this comparison, we used Kepler's two-body problem with eccentricities of 0.3, 0.5 and 0.7 and an interval of integration of [0, 2000π].
Fig. 4 depicts the graph of log 10 (EE) (vertical axis) as a function of CPU time in seconds (horizontal axis) for HBO (9,2,13), HBO (8,3,14) and AC(13) after a 1000 periods integration of Kepler's two-body problem.It is seen, from Fig. 4, that HBO(9,2,13) and HBO (8,3,14) compare favorably with AC(13) at stringent tolerances.The CPU percentage efficiency gain (CPU PEG) of method 2 over method 1 is defined by the formula (cf.Sharp [23]), where CPU 1, j and CPU 2, j are the CPU time of methods 1 and 2, respectively, associated with a given problem P, and j = − log 10 (EE).The CPU time was obtained from the curves which fit, in a least-squares sense, the data (log 10 (EE), log 10 (CPU)) by means of MATLAB's polyfit.
For e = 0.3, e = 0.5 and e = 0.7, the relative energy errors of HBO(9,2,13) and HBO (8,3,14) are generally less than the relative energy error of AC( 13) across the interval of integration.These results are consistent with the CPU PEGs listed in Table 6 for two-body problem with e = 0.3, e = 0.5 and e = 0.7.It is seen, from Fig. 5, that HBO(9,2,13) and HBO (8,3,14) compare favorably with AC(13) on the basis of the growth of relative positional error and relative energy error over 10000 periods of integration.
We use linear least-squares to fit the power law C 1 t C 2 to the graphs of Fig. 5 to obtain C 1 and C 2 shown in Table 7 and Table 8.
The values of C 2 of the relative positional errors of HBO(9,2,13) and HBO (8,3,14) are generally less than one and compare favorably with the value of C 2 of AC (13), listed in Table 7.
The values of C 2 of the relative energy errors of HBO(9,2,13), HBO (8,3,14) and AC (13), listed in Table 8, are in good agreement with the expected asymptotic value of one for non-symplectic methods.

Conclusion
This paper explores an alternative way to improve the stability and the order of predictor-corrector methods.
We have used only d = 2, 3, . . ., 9, Taylor coefficients, with d < p to construct two series of optimal, contractivity-preserving (CP), d-derivative, 2-and 3-step, predictor-corrector, Hermite-Birkhoff-Obrechkoff methods of order p up to 13 and 16, respectively, with nonnegative coefficients.It can be shown that HBO(d, k, p)) are absolutely stable for d = 1 to infinity.The stability regions of HBO(d, k, p) have generally a good shape and grow with decreasing p−d.The lower ends x min of their intervals of stability are, in general, less than x min = −2 of the most stable member of the predictor-corrector Adams-Bashforth-Moulton method family.
We compare the numerical performance of HBO (9,2,13), HBO (8,3,14) and AC (13) in solving Kepler orbit over an interval of 1000 periods.It is seen that HBO(9,2,13), HBO (8,3,14) use less CPU time than AC(13) at stringent tolerances.On the basis of the growth of relative positional error and relative energy error over 10000 periods of integration, these HBO methods also compare positively with AC(13) in solving standard N-body problems.

4 . 6 .
), r must satisfy, Optimizing c of predictor-corrector HBO(d, k, p) series With the newly defined r in (69), the problem of optimizing HBO(d, k, p) can be formulated as,c(HBO(d, k, p)) = max r T ,(70)subject to component-wise inequalities (57)-(61) together with following conditions of HBO(d, k, p) series methods 5.1.Existence of HBO(d, k, p) series methods as a function of d derivatives HBO(d, k, p) coefficients were obtained using security factor Θ = 2 and formula |x min,T (d)| = 0.4 d + 1.4 which is larger than |x min,T (d)| in (67) for d > 0.

Figure 1 :
Figure 1: Upper bounds p u of order p of HBO(d, 3, p) methods as an increasing function of d, number of their derivatives.
lists bounds |x min | of unscaled intervals of stability (x min , 0) of CP HBO(d, k, p) method, k = 2, 3, with p = 1.2 d , d = 5, 6, . . ., 12 as a function of d and p.Here, on the basis of |x min | averages and percentage stability gains, |x min | of CP HBO(d, k, p) method, k = 2, 3, compare favorably with |x min | of T(d) with the order of Taylor method T(d) less than 87 % of the order of HBO method.It is also seen that HBO(d, 3, p) methods are generally more stable than HBO(d, 2, p) methods of the same order.The values |x min | of CP HBO(d, 3, p = 1.3 d ) method as a function of d are depicted in Fig. 2. Similar to Barrio et al. [3], we use linear least square fits of |x min | as a function of d, the number of derivatives of the method.These linear least square fits give different positive slopes ρ and σ of |x min | = ρ d + σ depending on the values of η.For a given formula of p as a function of d, Table 4 lists slopes ρ and σ of a linear least square fit of |x min | = ρ d + σ for the listed methods.It is seen that HBO(d, k, p = η d ) series methods have slopes ρ positive and non negligible up to η = 1.4 and compare favorably with T(d) methods on the basis of |x min | averages and percentage stability gains.These |x min | averages and percentage stability gains increase with decreasing η.

Figure 2 :Table 4 :
Figure 2: Unscaled bounds |x min | of the intervals of absolute stability of HBO(d, 3, p = 1.3 d ) methods as a function of d, number of their derivatives.

Table 1 :
lists the upper bound p u of p of HBO(d, k, p), k = 2, 3, which exist and x min (p u ) of order p u .Here, it is seen that x min (p u ) compare positively with x min = −2 of the most stable member of the predictor-corrector Adams-Bashforth-Moulton method family.For a given d, the table lists the upper bound p u of p of HBO(d, k, p), k = 2, 3, which exist and x min (p u ) of order p u .For a given k = 2, 3, Table2lists the equations of line p u = ηd + ζ of best fit for a set of ordered pairs (d, k, p u ), appearing in Table

Table 2 :
For a given k = 2, 3, the table lists the equations of line p u = ηd + ζ of best fit for a set of ordered pairs (d, k, p u ), appearing in Table1.

Table 3 :
Unscaled |x min | of HBO(d, k, p = 1.2 d ), k = 2, 3, series method and of Taylor series T(d) method as a function of d and p.The percentage stability gains of HBO(d, 2, p) and HBO(d, 3, p) methods over T(d) are more than 16 % and 33 % respectively.

Table 7 :
Values of C 1 and C 2 of power law C 1 t C 2 fitted to the graphs of log 10 (RE(t)) as a function of log 10 (t) for a 10000 periods integration of Kepler's two-body problem with e = 0.3, e = 0.5 and e = 0.7 respectively.

Table 8 :
Values of C 1 and C 2 of power law C 1 t C 2 fitted to the graphs of log 10 (EE(t)) as a function of log 10 (t)