PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 86(100) (2009), 75–96 DOI: 10.2298/PIM0900075N

NINE-STAGE MULTI-DERIVATIVE RUNGE–KUTTA METHOD OF ORDER 12

Truong Nguyen-Ba, Vladan Božić, Emmanuel Kengne, and Rémi Vaillancourt

Communicated by Stevan Pilipović

Abstract. A nine-stage multi-derivative Runge–Kutta method of order 12, called HBT(12)9, is constructed for solving nonstiff systems of first-order differential equations of the form y′ = f(x, y), y(x_{0}) = y_{0}. The method uses y′ and higher derivatives y^{(2)} to y^{(6)} as in Taylor methods and is combined with a 9-stage Runge–Kutta method. Forcing an expansion of the numerical solution to agree with a Taylor expansion of the true solution leads to order conditions which are reorganized into Vandermonde-type linear systems whose solutions are the coefficients of the method. The stepsize is controlled by means of the derivatives y^{(3)} to y^{(6)} . The new method has a larger interval of absolute stability than Dormand–Prince’s DP(8,7)13M and is superior to DP(8,7)13M and Taylor method of order 12 in solving several problems often used to test high-order ODE solvers on the basis of the number of steps, CPU time, maximum global error of position and energy. Numerical results show the benefits of adding high-order derivatives to Runge–Kutta methods.

1. Introduction

A Taylor method of order 6, denoted by T6, and a 9-stage Runge–Kutta method of order 7 are cast into a nine-stage multi-derivative Runge–Kutta method of order 12, named HBT(12)9 because it uses Hermite–Birkhoff interpolation polynomials and high-order derivatives, y^{(2)} to y^{(6)}, for solving nonstiff systems of first-order initial value problems of the form

(1)  y′ = f(x, y),  y(x_{0}) = y_{0},  where ′ = \frac{d}{dx}.

2000 Mathematics Subject Classification: Primary 65L06; Secondary 65D05, 65D30. Key words and phrases: general linear method for non-stiff ODE; Hermite–Birkhoff method;

Taylor method; maximum global error; number of function evaluations; CPU time.

Partially supported by the Natural Sciences and Engineering Research Council of Canada and the Centre de recherches mathématiques of the Université de Montréal.

The link between the two types of methods is that values at off-step points are obtained by means of predictors which use values of derivatives of different orders at the current step point. By construction, HBT(12)9 uses lower order derivatives than the traditional Taylor method of order 12, denoted by T12 [16].

Taylor methods have been an excellent choice in astronomical calculations [3] and sensitivity analysis of ODEs/DAEs [2], and in solving general problems [7] and validating solutions of ODEs/DAEs by means of interval analysis [14, 17]. Deprit and Zahar [9] proved that recurrent power series in Taylor methods achieve high accuracy, with less computing time and larger stepsize than other methods.

The high-order derivatives used by HBT(12)9 can be obtained by differentiating f(x, y(x)) in the right-hand side of equation (1). But this approach is useful only in theoretical studies because of the computational complexity of high-order derivatives. Following Steffensen’s pioneering work [27, 25], fast automatic differentiation (AD) techniques are used to compute sums, differences, products and powers of power series, to name but a few (see [3, 16], and references therein). Formulae for generating these high-order derivatives can be found in textbooks (see, for instance, [12, pp. 46–49]).

Forcing an expansion of the numerical solution to agree with a Taylor expansion of the true solution leads to order conditions which are reorganized into linear Vandermonde-type systems leading to a convenient matrix formulation to handle order conditions. The solutions of these systems are the coefficients of the formulae which make HBT(12)9. These coefficients, which are available from the authors, were obtained to 32 digits by Gaussian elimination with Matlab variable precision arithmetic (VPA) for increased accuracy at stringent tolerance and use in extended precision computation.

The C++ performances of HBT(12)9, Dormand–Prince DP(8,7)13M [24] and T12, were compared on several problems frequently used to test higher order ODE solvers. It is seen that, generally, HBT(12)9 requires fewer steps, uses less CPU time, and has higher accuracy than DP(8,7)13M and T12.

Section 2 introduces HBT(12)9. Order conditions are listed in Section 3. In Section 4, HBT(12)9 is represented in terms of Vandermonde-type systems. Section 5 considers the region of absolute stability of the method. Section 6 deals with the step control. Numerical results are presented in Section 7.

2. One-step HBT(12)9

HBT(12)9 requires eight predictors and an integration formula to perform the integration step from x_{n} to x_{n+1} .

Let F_{j} := f(x_{n} + c_{j} h_{n+1} , Y_{j}) and set Y_{1} = y_{n}. Then the predictors P_{l},

(2)  Y_{l} = y_{n} + h_{n+1} ∑_{j=1}^{l−1} a_{lj} F_{j} + ∑_{j=2}^{6} h_{n+1}^{j} γ_{lj} y_{n}^{(j)} ,  l = 2, 3, …, 9,

are obtained recursively by means of Hermite–Birkhoff polynomials of degree l + 4, to order 6 for l = 2, order 7 for l = 3, and order 8 for l = 4, …, 9, respectively.

A Hermite–Birkhoff polynomial of degree 12 is used as integration formula IF to obtain y_{n+1} to order 12,

(3)  y_{n+1} = y_{n} + h_{n+1} ∑_{j=1}^{9} b_{j} F_{j} + ∑_{j=2}^{6} h_{n+1}^{j} γ_{j} y_{n}^{(j)} .

One sees that the derivatives y_{n}^{(2)} to y_{n}^{(6)} are computed only once per step at x_{n}. The defining formulae of HBT(12)9 involve the usual Runge–Kutta parameters c_{i}, a_{ij} and b_{j} and the Taylor expansion parameters γ_{lj} .

3. Order conditions for HBT(12)9

We impose the following simplifying assumptions on HBT(12)9 (with γ_{i1} = 0):

∑_{i=j+1}^{9} b_{i}a_{ij} = b_{j} (1 − c_{j}), j = 2, …, 8, (4)

b_{2} = b_{3} = 0,  a_{i2} = 0,  i = 4, ⋯, 9,

∑_{j=1}^{i−1} a_{ij} c_{j}^{k} + k!γ_{i,k+1} = \frac{1}{k + 1} c_{i}^{k+1} , \{ \begin{matrix}i = 2, 3, …, 9, \\ k = 0, 1, …, 5, \end{matrix}

∑_{j=1}^{i−1} a_{ij} c_{j}^{6} + 6!γ_{i7} = \frac{1}{7} c_{i}^{7}, i = 3, …, 9,

∑_{j=1}^{i−1} a_{ij} c_{j}^{7} + 7!γ_{i8} = \frac{1}{8} c_{i}^{8}, i = 4, …, 9.

There remain seven sets of equations to be solved [5]:

∑_{i=1}^{9} b_{i}c_{i}^{k} + k!γ_{k+1} = \frac{1}{k + 1} , k = 0, 1, …, 11, (5)

b_{8}(1 − c_{8})a_{87}c_{7}^{6}(c_{7} − c_{4})(c_{7} − c_{5})(c_{7} − c_{6}) (6)

= 9! ( \frac{1}{11!} − \frac{11}{12!}) − 8! ( \frac{1}{10!} − \frac{10}{11!}) (c_{4} + c_{5} + c_{6})

+7! ( \frac{1}{9!} − \frac{9}{10!}) (c_{4}c_{5} + c_{4}c_{6} + c_{5}c_{6}) − 6! ( \frac{1}{8!} − \frac{8}{9!}) c_{4}c_{5}c_{6},

b_{7}(1 − c_{7})(c_{8} − c_{7})a_{76}c_{6}^{6}(c_{6} − c_{4})(c_{6} − c_{5}) (7)

= 8! ( \frac{c_{8}}{10!} − (1 + c_{8})\frac{10}{11!} + 10\frac{11}{12!})

−7! ( \frac{c_{8}}{9!} − (1 + c_{8})\frac{9}{10!} + 9\frac{10}{11!}) (c_{4} + c_{5}) + 6! ( \frac{c_{8}}{8!} − (1 + c_{8})\frac{8}{9!} + 8\frac{9}{10!}) c_{4}c_{5},

b_{8}(1 − c_{8})a_{87}c_{7}^{6}(c_{7} − c_{4})(c_{7} − c_{5}) + b_{8}(1 − c_{8})a_{86}c_{6}^{6}(c_{6} − c_{4})(c_{6} − c_{5}) (8)

+b_{7}(1 − c_{7})a_{76}c_{6}^{6}(c_{6} − c_{4})(c_{6} − c_{5})

= ( \frac{8!}{10!} − 10\frac{8!}{11!}) − (c_{4} + c_{5}) ( \frac{7!}{9!} − 9\frac{7!}{10!}) + c_{4}c_{5} ( \frac{6!}{8!} − 8\frac{6!}{9!}) ,

∑_{i=4}^{7} b_{i}(1 − c_{i})(c_{8} − c_{i})a_{i3} = 0, (9)

∑_{i=4}^{8} b_{i}(1 − c_{i})a_{i3} = 0, (10)

∑_{i=5}^{8} b_{i}(1 − c_{i}) ∑_{j=4}^{i−1} a_{ij} a_{j3} = 0. (11)

The left-hand side of equation (6) is the result of the following expression similar to the left-hand side of Eq. (335j) in Butcher [6, pp. 206]:

(12)  ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{6}(c_{k} − c_{4})(c_{k} − c_{5})(c_{k} − c_{6}).

It is known that many terms in an expression of this form vanish (see [5]).

Expression (12) can also be written in terms of both sides of equations given in Appendix 8:

∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{9} − [ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{8}](c_{4} + c_{5} + c_{6})

+ [ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{7}](c_{4}c_{5} + c_{4}c_{6} + c_{5}c_{6})

− [ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{6}]c_{4}c_{5}c_{6}

= 9! ( \frac{1}{11!} − \frac{11}{12!}) − 8! ( \frac{1}{10!} − \frac{10}{11!}) (c_{4} + c_{5} + c_{6})

+ 7! ( \frac{1}{9!} − \frac{9}{10!}) (c_{4}c_{5} + c_{4}c_{6} + c_{5}c_{6}) − 6! ( \frac{1}{8!} − \frac{8}{9!}) c_{4}c_{5}c_{6},

= 9!(Eq.(45) − 8!(Eq.(33) − Eq.(41))(c_{4} + c_{5} + c_{6})

+ 7!(Eq.(27) − Eq.(31))(c_{4}c_{5} + c_{4}c_{6} + c_{5}c_{6})

− 6!(Eq.(24) − Eq.(26))c_{4}c_{5}c_{6},

since

c_{k}^{6}(c_{k} −c_{4})(c_{k} −c_{5})(c_{k} −c_{6}) = c_{k}^{9} −c_{k}^{8}(c_{4} +c_{5} +c_{6})+c_{k}^{7}(c_{4}c_{5} +c_{4}c_{6} +c_{5}c_{6})+c_{k}^{6}c_{4}c_{5}c_{6}

and

∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{9} = 9! ( \frac{1}{11!} − \frac{11}{12!}) = 9![Eq.(45) − Eq.(61)],

NINE-STAGE MULTI-DERIVATIVE RUNGE–KUTTA METHOD OF ORDER 12 79

∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{8} = 8! ( \frac{1}{10!} − \frac{10}{11!}) = 8![Eq.(33) − Eq.(41)],
∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{7} = 7! ( \frac{1}{9!} − \frac{9}{10!}) = 7![Eq.(27) − Eq.(31)],
∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{6} = 6! ( \frac{1}{8!} − \frac{8}{9!}) = 6![Eq.(24) − Eq.(26)].

Similarly, equations (7) and (8) are the results of the following two equations written in terms of equations given in Appendix 8, respectively:

∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})(c_{8} − c_{i})a_{ik} c_{k}^{6}(c_{k} − c_{4})(c_{k} − c_{5})

= c_{8} ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}a_{ik} c_{k}^{8} − (1 + c_{8}) ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}a_{ik} c_{k}^{8} + ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}^{2}a_{ik} c_{k}^{8}

+ [c_{8} ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}a_{ik} c_{k}^{7} − (1 + c_{8})[ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}a_{ik} c_{k}^{7}] + ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}^{2}a_{ik} c_{k}^{7}](c_{4} + c_{5})

+ [c_{8} ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}a_{ik} c_{k}^{6} − (1 + c_{8})[ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}a_{ik} c_{k}^{6}] + ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}c_{i}^{2}a_{ik} c_{k}^{6}](c_{4}c_{5})

= 8! [ \frac{c_{8}}{10!} − (1 + c_{8})\frac{10}{11!} + 10\frac{11}{12!}] − 7! [ \frac{c_{8}}{9!} − (1 + c_{8})\frac{9}{10!} + 9\frac{10}{11!}] (c_{4} + c_{5})

+ 6! [ \frac{c_{8}}{8!} − (1 + c_{8})\frac{8}{9!} + 8\frac{9}{10!}] c_{4}c_{5}

= 8![c_{8}Eq.(33) − (1 + c_{8})Eq.(41) + Eq.(57)] − 7![c_{8}Eq.(27) − (1 + c_{8})Eq.(31) + Eq.(39)](c_{4} + c_{5}) + 6![c_{8}Eq.(24) − (1 + c_{8})Eq.(26) + Eq.(30)](c_{4}c_{5}),

∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{6}(c_{k} − c_{4})(c_{k} − c_{5})

= ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{8} − [ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{7}](c_{4} + c_{5})

+ [ ∑_{i=1}^{9} ∑_{k=1}^{i−1} b_{i}(1 − c_{i})a_{ik} c_{k}^{6}]c_{4}c_{5}

= ( \frac{8!}{10!} − 10\frac{8!}{11!}) − ( \frac{7!}{9!} − 9\frac{7!}{10!}) (c_{4} + c_{5}) + ( \frac{6!}{8!} − 8\frac{6!}{9!}) c_{4}c_{5},

= 8![Eq.(33) − Eq.(41)] − 7![Eq.(27) − Eq.(31)](c_{4} + c_{5}) + 6![Eq.(24) − Eq.(26)]c_{4}c_{5}.

The nine off-step points used in this paper are

\begin{aligned}c_{1} & = 0, \\ c_{2} & = 0.34503974134180500927399082300440, \\ c_{3} & = 0.39433113296206286774170379771931, \\ c_{4} & = 0.45066415195664327741909005453635, \\ c_{5} & = 0.57269051227003684445548969961237, \\ c_{6} & = 0.28748636281590601727973892774720, \\ c_{7} & = 0.71586314033754605556936212451546, \\ c_{8} & = 0.91039578463195369728566674893955, \\ c_{9} & = 1, \end{aligned} ; (13)

which are chosen as follows.

Integration formula (3) contains 10 free parameters (b_{j}, j = 1, 4, 7, 8, 9, and γ_{j}, j = 2, 3, …, 6) and three free abscissae (x + hc_{j}, j = 4, 7, 8) for a total of 13 free parameters, while x+hc_{1} = x and x+hc_{9} = x+h are fixed abscissae and the three parameters c_{2}, c_{3} and c_{6} are to be determined later. Thus the formula is of order 13 since the five off-step points, c_{j}, j = 1, 4, 7, 8, 9 are obtained by the algebraic approach to Gauss integration formulae found in [8, pp. 85–87] and [13].

A 3-point Gauss-type integration formula with a 6-fold preassigned abscissa ξ_{1} = 0 and simple preassigned abscissa, ξ_{3} = 1 is of highest order 8 if the second abscissa is ξ_{2} = (7/8)ξ_{3}. Applying this formula to our case, we take c_{3} = (7/8)c_{4} and then c_{2} = (7/8)c_{3}.

The procedure to find c_{6} of (13) is described below. The abscissa c_{5} is adjusted so that c_{6} is a suitable value between 0 and 1. Firstly, we write the following reduced equation

(14) b_{8}(1 − c_{8})a_{87}a_{76}c_{6}^{6}(c_{6} − c_{4})(c_{6} − c_{5})

= 8! ( \frac{1}{11!} − \frac{11}{12!}) − (c_{4} + c_{5})7! ( \frac{1}{10!} − \frac{10}{11!}) + c_{4}c_{5}6! ( \frac{1}{9!} − \frac{9}{10!})

which is the result of the following equation written in terms of equations given in Appendix 8:

∑_{i=1}^{9} ∑_{j=1}^{i−1} ∑_{k=1}^{j−1} b_{i}(1 − c_{i})a_{ij} a_{jk} c_{k}^{6}(c_{k} − c_{4})(c_{k} − c_{5})

= ∑_{i=1}^{9} ∑_{j=1}^{i−1} ∑_{k=1}^{j−1} b_{i}(1 − c_{i})a_{ij} a_{jk} c_{k}^{8} − [ ∑_{i=1}^{9} ∑_{j=1}^{i−1} ∑_{k=1}^{j−1} b_{i}(1 − c_{i})a_{ij} a_{jk} c_{k}^{7}](c_{4} + c_{5})

+ [ ∑_{i=1}^{9} ∑_{j=1}^{i−1} ∑_{k=1}^{j−1} b_{i}(1 − c_{i})a_{ij} a_{jk} c_{k}^{6}]c_{4}c_{5}

= ( \frac{8!}{11!} − 11\frac{8!}{12!}) − ( \frac{7!}{10!} − 10\frac{7!}{11!}) (c_{4} + c_{5}) + ( \frac{6!}{9!} − 9\frac{6!}{10!}) c_{4}c_{5}

= 8![Eq.(49) − Eq.(65)] − 7![Eq.(35) − Eq.(43)](c_{4} + c_{5})

+ 6![Eq.(28) − Eq.(32)]c_{4}c_{5}.

Next, we write

θ = c_{7}^{6}(c_{7} − c_{4})(c_{7} − c_{5})(c_{7} − c_{6})b_{7}(1 − c_{7})(c_{8} − c_{7}),

so that the product of the left-hand sides of (6) and (7) is the product of θ with the left-hand side of (14). We therefore have

(15) [9! ( \frac{1}{11!} − \frac{11}{12!}) − (c_{4} + c_{5} + c_{6})8! ( \frac{1}{10!} − \frac{10}{11!})

+(c_{4}c_{5} + c_{4}c_{6} + c_{5}c_{6})7! ( \frac{1}{9!} − \frac{9}{10!}) − c_{4}c_{5}c_{6}6! ( \frac{1}{8!} − \frac{8}{9!}) ]

× [8! ( \frac{c_{8}}{10!} − (1 + c_{8}) \frac{10}{11!} + 10 \frac{11}{12!}) − (c_{4} + c_{5})7! ( \frac{c_{8}}{9!} − (1 + c_{8}) \frac{9}{10!} + 9 \frac{10}{11!})

+c_{4}c_{5}6! ( \frac{c8}{8!} − (1 + c_{8})\frac{8}{9!} + 8\frac{9}{10!}) ]

= [8! ( \frac{1}{11!} − \frac{11}{12!}) − (c_{4} + c_{5})7! ( \frac{1}{10!} − \frac{10}{11!}) + c_{4}c_{5}6! ( \frac{1}{9!} − \frac{9}{10!}) ] θ.

Setting c_{i} equal to the values of (13) for all i except i = 6, we can calculate c_{6} such that (15) and the linear system (16) below for the integration formula are satisfied. System (16) needs to be satisfied since θ is a function of b_{7}.

Put simply, c_{6} is chosen such that condition (65) in Appendix 8 is met automatically when all the other order conditions are satisfied.

4. Matrix formulation of HBT(12)9

Of the many methods to construct the formulae which make HBT(12)9, we choose to express the coefficients as unknowns of linear systems built from the order conditions and solved, in particular, by Matlab. The Matlab colon (:) notation is used, say, 1 : 4 for 1, 2, 3, 4.

4.1. Integration formula IF. Let the 12-vector of the reordered coefficients of IF in (3), u^{1} = [b_{9} b_{8} b_{7} b_{6} b_{5} b_{4} b_{1} γ_{2} γ_{3} γ_{5} γ_{5} γ_{6}]^{T}, be the solution of the Vandermonde-type system of order conditions:

(16)  [ [ \frac{c_{10−j}^{i−1}}{(i − 1)!} ]_{i=1:12,j=1:6}  [ \begin{matrix}I_{6} \\ 0_{6×6} \end{matrix} ] ] u^{1} = [ \frac{1}{i!}]_{i=1:12} .

With the choice of c_{i}, i = 4, 5, …, 9, in (13), the leading error term of IF is of order 14:

[b_{9} \frac{c_{9}^{13}}{13!} + ⋯+ b_{5} \frac{c_{5}^{13}}{13!} + b_{4} \frac{c_{4}^{13}}{13!} − \frac{1}{14!}]h_{n+1}^{14} y_{n}^{(14)} .

4.2. Predictor P_{2}. Let u^{2} = [a_{21} γ_{22} γ_{23} γ_{24} γ_{25} γ_{26}]^{T} be the 6-vector of reordered coefficients of predictor P_{2} in (2) with l = 2. Then the ith component of u^{2}, u_{2}(i), satisfies the order condition

u_{2}(i) = \frac{c_{2}^{i}}{i!} , i = 1, 2, …, 6.

A truncated Taylor expansion of the right-hand side of (2) with l = 2 about x_{n} gives

∑_{j=0}^{6} \frac{c_{2}^{j}}{j!} h_{n+1}^{j} y_{n}^{(j)}

which implies that P_{2} is of order 6 with leading error term (c_{2}^{7}/7!) h_{n+1}^{7} y_{n}^{(7)}.

4.3. Predictor P_{3}. The 7-vector u^{3} = [a_{32} a_{31} γ_{32} γ_{33} γ_{34} γ_{35} γ_{36}]^{T} of the reordered coefficients of predictor P_{3} in (2) with l = 3 is the solution of the system of order conditions

[ [ \frac{c_{2}^{i−1}}{(i − 1)!} ]_{i=1:7} [ \begin{matrix}I_{6} \\ 0_{1×6} \end{matrix} ] ] u^{3} = [ \frac{c_{3}^{i}}{i!} ]_{i=1:7} .

A truncated Taylor expansion about x_{n} of the right-hand side of (2) with l = 3 gives

∑_{j=0}^{7} \frac{c_{3}^{j}}{j!} h_{n+1}^{j} y_{n}^{(j)}

which implies that P_{3} is of order 7 with leading error term (c_{3}^{8}/8!) h_{n+1}^{8} y_{n}^{(8)}.

4.4. Predictors P_{4} and P_{5}. The vector u^{l} = [a_{l3} a_{l2} a_{l1} γ_{l2} γ_{l3} γ_{l4} γ_{l5} γ_{l6}]^{T} of the eight reordered coefficients of predictors P_{4} and P_{5} in (2) with l = 4 and l = 5, respectively, are the solution of the system of order conditions

[ [ \frac{c_{l−j}^{i−1}}{(i − 1)!} ]_{i=1:8,j=1:2} [ \begin{matrix}I_{6} \\ 0_{2×6} \end{matrix} ] ] u^{l} = [ \frac{c_{l}^{i}}{i!} ]_{i=1:8} .

4.5. The coefficients a_{ij} of P_{i}, for i = 6, 7, 8 and j = 3, 4, 5. It is

numerically convenient first to solve for a_{87} and a_{76} from (6) and (7), and a_{86} from (8). Next, we solve for the nine coefficients a_{63}, a_{64}, a_{65}, a_{73}, a_{74}, a_{75}, a_{83}, a_{84}, a_{85} of predictors P_{6} to P_{8} simultaneously before solving for their other coefficients. These nine coefficients are solutions of the system of order conditions

(17) [ \begin{matrix}c_{5}^{6}/6! & c_{4}^{6}/6! & c_{3}^{6}/6! & 0 \\ c_{5}^{7}/7! & c_{4}^{7}/7! & c_{3}^{7}/7! & 0 \\ 0 & 0 & 0 & c_{5}^{6}/6! \\ 0 & 0 & 0 & c_{5}^{7}/7! \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & b_{6}(1 − c_{6}) & 0 \\ b_{6}(1 − c_{6})a_{53} & b_{6}(1 − c_{6})a_{43} & b_{8}(1 − c_{8})a_{86} + b_{7}(1 − c_{7})a_{76} & b_{7}(1 − c_{7})a_{53} \\ 0 & 0 & b_{6}(1 − c_{6})(c_{8} − c_{6}) & 0 \end{matrix}

\begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ c_{4}^{6}/6! & c_{3}^{6}/6! & 0 & 0 & 0 \\ c_{4}^{7}/7! & c_{3}^{7}/7! & 0 & 0 & 0 \\ 0 & 0 & c_{5}^{6}/6! & c_{4}^{6}/6! & c_{3}^{6}/6! \\ 0 & 0 & c_{5}^{7}/7! & c_{4}^{7}/7! & c_{3}^{7}/7! \\ 0 & b_{7}(1 − c_{7}) & 0 & 0 & b_{8}(1 − c_{8}) \\ b_{7}(1 − c_{7})a_{43} & b_{8}(1 − c_{8})a_{87} & b_{8}(1 − c_{8})a_{53} & b_{8}(1 − c_{8})a_{43} & 0 \\ 0 & b_{7}(1 − c_{7})(c_{8} − c_{7}) & 0 & 0 & 0 \end{matrix} ][ \begin{matrix}a_{65} \\ a_{64} \\ a_{63} \\ a_{75} \\ a_{74} \\ a_{73} \\ a_{85} \\ a_{84} \\ a_{83} \end{matrix} ]

= r,

where r = r(1 : 9) has components

r(1) = c_{6}^{7}/7!,  r(5) = c_{8}^{7}/7! − a_{87}c_{7}^{6}/6! − a_{86}c_{6}^{6}/6!,
r(2) = c_{6}^{8}/8!,  r(6) = c_{8}^{8}/8! − a_{87}c_{7}^{7}/7! − a_{86}c_{6}^{7}/7!,
r(3) = c_{7}^{7}/7! − a_{76}c_{6}^{6}/6!,  r(7) = −b_{4}(1 − c_{4})a_{43} − b_{5}(1 − c_{5})a_{53},
r(4) = c_{7}^{8}/8! − a_{76}c_{6}^{7}/7!,  r(8) = −b_{5}(1 − c_{5})a_{54}a_{43},

and

r(9) = −b_{4}(1 − c_{4})(c_{8} − c_{4})a_{43} − b_{5}(1 − c_{5})(c_{8} − c_{5})a_{53}.

The equations for r(7), r(8) and r(9) correspond to equations (10), (11) and (9), respectively.

4.6. Predictors P_{l}, l = 6, 7, 8. Since a_{l5}, a_{l4}, a_{l3} are already obtained from system (17), the remaining six unknown coefficients of predictor P_{l} in (2) with l = 6 are in the 6-vector of reordered coefficients, u^{l} = [a_{l1}, γ_{l2}, γ_{l3}, γ_{l4}, γ_{l5}, γ_{l6}]^{T}, whose ith component, u_{l}(i), satisfies the order condition

u_{l}(i) = \frac{c_{l}^{i}}{i!} − \frac{1}{(i − 1)!} ∑_{j=3}^{l−1} a_{lj} c_{j}^{i−1} , i = 1, 2, …, 6,

where a_{76} is obtained from (7) and a_{87} and a_{86} are obtained from (6) and (8) respectively.

4.7. Predictor P_{9}. The 12-vector of reordered coefficients of predictor P_{9} in (2) with l = 9, u^{9} = [a_{98} a_{97} a_{96} a_{95} a_{94} a_{93} a_{91} γ_{92} γ_{93} γ_{94} γ_{95} γ_{96}]^{T}, is the solution of the system of order conditions

[ [ \begin{matrix}[ \frac{c_{9−j}^{i−1}}{(i − 1)!} ] & i=1:6,j=1:6 \\ b_{9}I_{6} \end{matrix}  ] [ \begin{matrix}I_{6} \\ 0_{6×6} \end{matrix} ] ] u^{9} = r^{9}

where r^{9} = r_{9}(1 : 12) has components

r_{9}(i) = c_{9}^{i}/i!,  i = 1, 2, …, 6,
r_{9}(7) = b_{8}(1 − c_{8}),
r_{9}(8) = b_{7}(1 − c_{7}) − (b_{8}a_{87}),
r_{9}(9) = b_{6}(1 − c_{6}) − (b_{8}a_{86} + b_{7}a_{76}) ,
r_{9}(10) = b_{5}(1 − c_{5}) − (b_{8}a_{85} + b_{7}a_{75} + b_{6}a_{65}) ,
r_{9}(11) = b_{4}(1 − c_{4}) − (b_{8}a_{84} + b_{7}a_{74} + b_{6}a_{64} + b_{5}a_{54}) ,
r_{9}(12) = b_{3}(1 − c_{3}) − (b_{8}a_{83} + b_{7}a_{73} + b_{6}a_{63} + b_{5}a_{53} + b_{4}a_{43}) .

The equations for r_{9}(i), i = 7, 8, …, 12, correspond to (4).

5. Region of absolute stability

To obtain the region of absolute stability, R, of HBT(12)9, we apply the predictors P_{2}, P_{3}, …, P_{9} and the integration formula IF with constant step h to the linear test equation

y′ = λy, y_{0} = 1.

Thus we obtain

(18) Y_{l} = y_{n} + λh_{n+1} ∑_{j=1}^{l−1} a_{lj} Y_{j} + ∑_{j=2}^{6} (λh_{n+1} )^{j} γ_{lj} y_{n},  l = 2, 3, …, 9,

and

(19)  y_{n+1} = y_{n} + λh_{n+1} ∑_{j=1}^{9} b_{j} Y_{j} + ∑_{j=2}^{6} (λh_{n+1} )^{j} γ_{j} y_{n}.

If we replace Y_{l}, for l = 2, 3, …, 9, in (18) and (19) with the corresponding righthand sides of (18), then (19) reduces to the following first-order difference equation and corresponding linear characteristic equation:

−r_{s}y_{n} + y_{n+1} = 0, −r_{s} + r = 0,

respectively. The root, r_{s}, of the characteristic equation is

r_{s} = 1 + ∑_{j=1}^{14} s_{j} λ^{j} h^{j} ,

Figure 1. Top half of the region of absolute stability of HBT(12)9 (left). Value of k vs. order p for listed tolerance (right).

with coefficients

s_{1} = 1.0,  s_{2} = 5.00000000000000 e-01,
s_{3} = 1.66666666666666 e-01,  s_{4} = 4.16666666666666 e-02,
s_{5} = 8.33333333333332 e-03,  s_{6} = 1.38888888888888 e-03,
s_{7} = 1.98412698412697 e-04,  s_{8} = 2.48015873015878 e-05,
s_{9} = 2.75573192239835 e-06,  s_{10} = 2.75573192239858 e-07,
s_{11} = 2.50521083854427 e-08,  s_{12} = 2.08767569879261 e-09,
s_{13} = 5.02071324573546 e-12,  s_{14} = 2.99693530083494 e-11.

A complex number λh is in R if r_{s} satisfies the root condition: |r_{s}| ⩽ 1 (see [12, pp. 378–380]).

The root condition is used to find the region of absolute stability of HBT(12)9 whose top half is shown in grey in Fig. 5, with interval of absolute stability (α, 0) = (−5.40, 0). We note that HBT(12)9 has a larger interval of absolute stability than DP(8,7)13M, namely, 5.40 > 5.12.

6. Controlling stepsize

6.1. The principal error terms of HBT(12)9 and T12 and the stepsizes. It is known that the step control predictor of Runge–Kutta pairs of orders p and p − 1 or p − 2 is of order p − 1 [24], or p − 2 [10]. The error of the step control predictor is kept within tolerance TOL while integration is done by the Runge–Kutta formula of higher order. Similarly, in our case, T12 and T11 act as step control predictors. The errors of orders 12 and 11 control the stepsizes.

If the principal error term of HBT(12)9 is C_{PET} h_{HBT}^{13} , then to obtain the same error at each integration step we set

C_{PET} h_{HBT}^{13} = \frac{y^{(13)}}{13!} h_{T}^{13},

where h_{HBT} and h_{T} are the stepsizes of HBT(12)9 and T12, respectively. Then we have

(20)  h_{HBT} = [ \frac{y^{(13)}}{C_{PET} 13!}]^{1/13} h_{T} =: η h_{T}.

If C_{PET} < y^{(13)} /13! then η > 1. This result will be used to justify the value of the factor η in the stepsize formula (22) in the next subsection.

6.2. Step size control. The stepsize, h_{n+1} , of the Taylor method of order p can be chosen within tolerance TOL by the formula (see [16, 3])

(21) h_{n+1} = min\{ k(TOL, p − 1) ∥ \frac{y^{(p−1)}}{(p − 1)!} ∥_{∞}^{−1/(p−1)} , k(TOL, p) ∥ \frac{y^{(p)}}{p!} ∥_{∞}^{−1/p} \} ,

where k(TOL, p) is the solution of the equation k^{p+1} /(1 − k) = TOL (see Fig. 5 (right)).

Since HBT(12)9 does not use derivatives of order higher than six, to determine the stepsize we shall use the following formula (22)

h_{n+1} = η min \{ k(TOL, 11) [ \frac{∥y^{(3)}∥_{∞} /3!}{[∥y^{(5)}∥_{∞} /5!]^{2}}]^{1/7} , k(TOL, 12) [ \frac{∥y^{(4)}∥_{∞} /4!}{[∥y^{(6)}∥_{∞} /6!]^{2}}]^{1/8} \}

similar to error estimators found in [3]. The exponents, in the above formula, come from 1/7 = (11/7)(1/11) and 1/8 = (12/8)(1/12).

It was observed that HBT(12)9 solves the ODEs considered in this paper more efficiently with stepsize h_{n+1} obtained by (22) without rejected steps than by means of a step control predictor. In (22), η acts as control factor in the variable step algorithm. If η is set to 1.0 as the assumption C_{PET} = y^{(13)} /13!, the stepsize of HBT(12)9 is very conservative. In our tests, we have fixed η = 1.4.

7. Numerical results

The derivatives, y^{(2)} to y^{(6)}, are calculated at each integration step by known recurrence formulae (see, for example, [12, pp. 46–49], [16]).

Computations were performed in C++ on a Mac with a dual 2.5 GHz PowerPC G5 and 4 GB DDR SSRAM running under Mac OS X Version 10.4.8.

7.1. Numerical results related to the step control. Table 1 lists the number of steps (NS) and the maximum global error (GE) of HBT(12)9 and T12 related to the step control for the DETEST problems [11] of class A, B, and E over the time interval [0, 20] with set tolerance (TOL). Thus we can compare the step controls of HBT(12)9 and T12.

In Table 2, we compare HBT(12)9 with results for Taylor’s method of order 12 obtained by Lara’s program [16], denoted by T12L. The considered problems are Kepler’s, Hénon–Heiles’ and the equatorial main problems over the time interval [0, t_{f} ] Table 1. For some test problems of [11], time interval [0, 20] and LT = log_{10}(TOL), the table lists the number of steps (NS) and the maximum global error (GE) for HBT(12)9 (left column) and T12 (right column).

HBT(12)9 and T12

Problem LT  NS  GE
A1  -04  8  6 6.03e-03 1.86e-05
-07 10  9 3.56e-05 2.45e-08
-10 14 15 2.06e-08 3.11e-11
A3  -04 19 22 2.44e-03 7.54e-05
-07 31 36 6.22e-08 3.68e-07
-10 53 63 6.65e-09 3.89e-10
A4  -04  5  5 1.99e-07 2.39e-05
-07  8  7 1.24e-09 6.91e-08
-10 13 13 2.75e-12 5.56e-11
B1  -04 30 40 3.71e-03 9.61e-03
-07 53 68 7.10e-07 1.48e-05
-10 91 120 1.52e-09 7.48e-09
B3 -04 9 9 8.84e-03 2.86e-05
-07 11 13 7.23e-06 6.63e-08
-10 17 21 7.82e-09 9.47e-11
B4  -04 21 18 2.64e-06 1.23e-04
-07 35 31 3.96e-09 1.15e-07
-10 60 53 2.17e-12 2.90e-10
B5  -04 19 23 3.00e-06 6.81e-04
-07 32 39 2.09e-09 1.96e-07
-10 55 68 2.78e-12 3.06e-10
E1  -04 11 11 1.45e-05 3.28e-05
-07 18 17 1.86e-08 4.47e-08
-10 30 29 2.13e-11 5.48e-11
E2  -04 38 50 7.51e-05 1.30e-04
-07 68 84 5.87e-08 4.25e-07
-10 117 147 1.61e-11 8.51e-10
E3  -04 29 33 2.43e-07 6.35e-04
-07 49 57 2.87e-10 5.68e-06
-10 85 98 7.21e-13 3.45e-08
E4  -04  4  2 1.27e-05 3.09e-05
-07  5  4 5.38e-08 5.24e-08
-10  8  7 1.86e-11 4.01e-11
E5  -04  4  3 1.17e-05 2.70e-04
-07  6  6 1.89e-09 5.92e-07
-10 10 10 7.25e-13 8.30e-10

Table 2. For given LT = log_{10}(TOL), the table lists the number of steps (NS) and the maximum global energy error (MGEE) for HBT(12)9 (left column) and T12L (right column) on the problems on hand over the time interval [0, t_{f} ].

HBT(12)9 and T12L

Problem LT NS MGEE
D1  -04 35 44 2.21e-04 5.56e-04
t_{f} = 16π  -07 60 73 3.11e-07 1.02e-06
-10 103 122 2.39e-10 1.12e-09
D3  -04 60 83 2.93e-03 8.32e-04
t_{f} = 16π  -07 106 139 7.77e-08 1.24e-06
-10 186 235 4.05e-11 1.43e-09
D5  -04 111 167 8.95e-03 2.47e-03
t_{f} = 16π  -07 205 273 2.20e-07 3.31e-09
-10 362 461 2.97e-11 3.62e-10
Hénon–Heiles -04 51 66 1.82e-05 2.42e-04
t_{f} = 70  -07 86 108 1.92e-07 3.81e-07
-10 148 185 2.13e-10 1.70e-10
Equatorial -04 102 172 1.78e-02 7.24e-04
main prob. -07 179 289 1.12e-06 1.08e-06
t_{f} = 70  -10 319 489 1.24e-09 1.77e-09

The maximum global energy error (MGEE) was obtained from the maximum of the absolute value of the relative error H/H_{0} − 1 at every integration step where H and H_{0} are the values of the Hamiltonian at t_{n+1} and t_{0}, respectively.

The Hamiltonians of Kepler’s, Hénon–Heiles’ and the equatorial main problems are

H_{Kepler} = \frac{1}{2} (y_{3}^{2} + y_{4}^{2}) − (y_{1}^{2} + y_{2}^{2}) ^{−1/2} ,
H_{Hénon–Heiles} = \frac{1}{2} (X^{2} + Y^{2}) + \frac{1}{2} (x^{2} + y^{2}) + ϵy ( x^{2} − \frac{1}{3} y^{2}) ,
H_{eq. main prob.} = \frac{1}{2} ( P^{2} + \frac{Λ^{2}}{ρ^{2}} + Z^{2}) + \frac{μ}{r} + \frac{α^{2}J_{2}μP_{2}(u)}{r^{3}} ,

respectively, where, in H_{eq. main prob.} , u = z/r, r = \sqrt{ρ^{2} + z^{2}} and P_{2}(x) = (3x^{2} − 1)/2 is the Legendre polynomial of degree 2.

Tables 1 and 2 show that our stepsize control is reliable for the problems on hand and usually compares favorably with the step control of T12 and T12L.

7.2. Comparison based on CPU time. We compare the CPU time in seconds used by HBT(12)9, T12, and DP(8,7) in solving several problems. The maximum global error (MGE) is taken to be max_{n}{∥y_{n+1} − y(t_{n+1} )∥_{∞} } of the difference between the numerical and the analytic solutions at every integration step for Kepler’s problem. For the other problems, y(t_{n+1} ) is replaced by reference solutions obtained by DP(8,7)13M at stringent tolerance 5 × 10^{−14} . In Fig. 7.2,

HBT(12)9 ∘, T12 × and DP(8,7)13M ▷

Figure 2. CPU time in seconds (horizontal axis) versus log_{10}(MGE) (vertical axis) for the listed problems.

CPU time in seconds (horizontal axis) is plotted versus log_{10}(MGE) (vertical axis) for the above problems.

The CPU percentage efficiency gain (CPU PEG) is defined by formula (cf. Sharp [26]),

(CPU PEG) = 100 [ \begin{matrix}∑ & j & CPU & 2,j \\ ∑ & j & CPU & 1,j \end{matrix} − 1] ,

where CPU_{1,j} and CPU_{2,j} are the CPU time of methods 1 and 2, respectively, and j = − log_{10}(MGE). The CPU time was obtained from the curves which fit, in a least-squares sense, the data (log_{10}(MGE), log_{10}(CPU)) by means of Matlab’s polyfit. The CPU PEG of HBT(12)9 over DP(8,7)13M and T12 for the above problems are listed in the middle part of Table 3.

It is seen from Fig. 7.2 and Table 3 that, at stringent tolerance, HBT(12)9 compares favorably with both DP(8,7)13M and T12 on the basis of CPU time versus MGE and versus CPU PEG.

7.3. Comparison based on the number of steps. The number of step percentage efficiency gain (NS PEG)_{i} is defined by the formula

(NS PEG) = 100 [ \begin{matrix}∑ & j & NS & T,j \\ ∑ & j & NS & HBT ,j \end{matrix} − 1] ,

where NS_{T,j} and NS_{HBT ,j} are the number of steps used by methods T12 and HBT(12)9, respectively, to integrate from t_{0} to t_{f} , and j = − log_{10}(MGEE). The Table 3. CPU time and NS PEG of HBT(12)9 over DP(8,7)13M and T12 for the listed problems.

CPU PEG of HBT(12)9 over: NS PEG of HBT(12)9 over:
Problem  DP(8,7)13M  T12  DP(8,7)13M  T12
Kepler (e=0.1)  143%  78%  159%  52%
Kepler (e=0.3)  101%  89%  143%  72%
Kepler (e=0.5)  110%  133%  136%  88%
Kepler (e=0.7)  137%  112%  158%  89%
Kepler (e=0.9)  146%  84%  183%  71%
Kepler (e=0.99)  258%  114%  161%  86%
Hénon–Heiles  35%  57%  123%  27%
Eq. main prob.  5%  111%  33%  41%
B1  77%  57%
B5  67%  147%
E2  32%  120%
Arenstorf  130%  223%

HBT(12)9 ∘, T12 × and DP(8,7)13M ▷

Figure 3. Number of steps (horizontal axis) versus log_{10}(MGEE)

(vertical axis) for the problems on hand.

number of steps (NS) was obtained from the curves which fit, in the least squares sense, the data (log_{10}(MGEE), log_{10}(NS)).

In Fig. 7.3, the number of step (horizontal axis) is plotted versus log_{10}(MGEE) (vertical axis) for the methods and problems on hand. It is observed that HBT(12)9 performs better than T12 on the basis of the number of steps versus MGEE shown in Fig. 7.3 and the number of step percentage efficiency gain listed in the rightmost part of Table 3.

The numerical results show that a combination of high-order derivatives with a Runge–Kutta method achieves a high degree of accuracy. It is to be noted that HBT(12)9 uses six derivatives of y compared to twelve for T12.

8. Conclusion

A one-step 9-stage Hermite–Birkhoff–Taylor method of order 12, HBT(12)9, was constructed by solving Vandermonde-type systems satisfying Runge–Kuttatype order conditions. By construction, HBT(12)9 uses lower order derivatives than the traditional Taylor method of order 12. The stability region of HBT(12)9 has a remarkably good shape. The stepsize is controlled by a formula which uses y_{n}^{(4)} and y_{n}^{(6)}. On the basis of CPU time versus the maximum global error, and the number of steps versus the maximum global energy error, HBT(12)9 wins over DP(8,7)13M and T12 in solving several well-known test problems. HBT methods with six high derivatives y^{(1)} to y^{(6)} appear to be promising for ODEs in the light of the numerical results since methods of high order can be derived and implemented efficiently. Furthermore, since these methods use a small number of high order derivatives, they may be useful for high dimensional problems.

Acknowledgment

The anonymous referees are deeply thanked for their suggestions which improved the paper considerably. Martín Lara is also thanked for supplying his Fortran programs.

Appendix. The order conditions used in equations (6), (7), (8), and (14)

Some of the order conditions listed here are used in equations (6), (7), (8) and (14).

Order 1 to 7:

∑ b_{i} = 1,  ∑ b_{i}c_{i} + γ_{2} = \frac{1}{2} ∑ b_{i}c_{i}^{2} + 2!γ_{3} = \frac{1}{3}
∑ b_{i}c_{i}^{3} + 3!γ_{4} = \frac{1}{4} ∑ b_{i}c_{i}^{4} + 4!γ_{5} = \frac{1}{5}, ∑ b_{i}c_{i}^{5} + 5!γ_{6} = \frac{1}{6}
∑ b_{i}c_{i}^{6} + 6!γ_{7} = \frac{1}{7},

Order 8:

(23)  ∑ b_{i}c_{i}^{7} = \frac{1}{8}
(24)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{6}}{6!}] = \frac{1}{8!}

92 NGUYEN-BA, BOŽIĆ, KENGNE, AND VAILLANCOURT

Order 9:

(25)  ∑ b_{i}c_{i}^{8} = \frac{1}{9}
(26)  ∑ b_{i} \frac{c_{i}}{8} [∑ a_{ij} \frac{c_{j}^{6}}{6!}] = \frac{1}{9!}
(27)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{7}}{7!}] = \frac{1}{9!}
(28)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{9!}

Order 10:

(29)  ∑ b_{i}c_{i}^{9} = \frac{1}{10}
(30)  ∑ b_{i} \frac{c_{i}^{2}}{8 × 9} [∑ a_{ij} \frac{c_{j}^{6}}{6!}] = \frac{1}{10!}
(31)  ∑ b_{i} \frac{c_{i}}{9} [∑ a_{ij} \frac{c_{j}^{7}}{7!}] = \frac{1}{10!}
(32)  ∑ b_{i} \frac{c_{i}}{9} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{10!}
(33)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{8}}{8!}] = \frac{1}{10!}
(34)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}}{8} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{10!}
(35)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{10!}
(36)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{10!}

Order 11:

(37)  ∑ b_{i}c_{i}^{10} = \frac{1}{11}
(38)  ∑ b_{i} \frac{c_{i}^{3}}{8 × 9 × 10} [∑ a_{ij} \frac{c_{j}^{6}}{6!}] = \frac{1}{11!}
(39)  ∑ b_{i} \frac{c_{i}^{2}}{9 × 10} [∑ a_{ij} \frac{c_{j}^{7}}{7!}] = \frac{1}{11!}
(40)  ∑ b_{i} \frac{c_{i}^{2}}{9 × 10} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{11!}
(41)  ∑ b_{i} \frac{c_{i}}{10} [∑ a_{ij} \frac{c_{j}^{8}}{8!}] = \frac{1}{11!}
(42)  ∑ b_{i} \frac{c_{i}}{10} [∑ a_{ij} \frac{c_{j}}{8} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{11!}

(43)  ∑ b_{i} \frac{c_{i}}{10} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{11!}

(44)  ∑ b_{i} \frac{c_{i}}{10} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{11!}

(45)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{9}}{9!}] = \frac{1}{11!}

(46)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{2}}{8 × 9} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{11!}

(47)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}}{9} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{11!}

(48)  ∑ b_{i} [∑ a_{ij} \frac{c_{i}}{9} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{11!}

(49)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{8}}{8!} ) ] = \frac{1}{11!}

(50)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}}{8} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{11!}

(51)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{7}}{7!}) ) ] = \frac{1}{11!}

(52)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{l,m} \frac{c_{m}^{6}}{6!} ) ) ) ] = \frac{1}{11!}.

Order 12:

(53)  ∑_{i} b_{i}c_{i}^{11} = \frac{1}{12}

(54)  ∑ b_{i} \frac{c_{i}^{4}}{8 × 9 × 10 × 11} [∑ a_{ij} \frac{c_{j}^{6}}{6!}] = \frac{1}{12!}

(55)  ∑ b_{i} \frac{c_{i}^{3}}{9 × 10 × 11} [∑ a_{ij} \frac{c_{j}^{7}}{7!}] = \frac{1}{12!}

(56)  ∑ b_{i} \frac{c_{i}^{3}}{9 × 10 × 11} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{12!}

(57)  ∑ b_{i} \frac{c_{i}^{2}}{10 × 11} [∑ a_{ij} \frac{c_{j}^{8}}{8!}] = \frac{1}{12!}

(58)  ∑ b_{i} \frac{c_{i}^{2}}{10 × 11} [∑ a_{ij} \frac{c_{j}}{8} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{12!}

(59)  ∑ b_{i} \frac{c_{i}^{2}}{10 × 11} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{12!}

(60)  ∑ b_{i} \frac{c_{i}^{2}}{10 × 11} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}

(61)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} \frac{c_{j}^{9}}{9!}] = \frac{1}{12!}

(62)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} \frac{c_{j}}{8 × 9} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{12!}^{(63)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} \frac{c_{j}}{9} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{12!}}^{(64)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} \frac{c_{j}}{9} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}}^{(65)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{8}}{8!} ) ] = \frac{1}{12!}}^{(66)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}}{8} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}}^{(67)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{7}}{7!}) ) ] = \frac{1}{12!}}^{(68)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{lm} \frac{c_{m}^{6}}{6!} ) ) ) ] = \frac{1}{12!}}^{(69)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{10}}{10!}] = \frac{1}{12!}}^{(70)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{3}}{8 × 9 × 10} ( ∑ a_{jk} \frac{c_{k}^{6}}{6!} ) ] = \frac{1}{12!}}^{(71)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}^{2}}{9 × 10} ( ∑ a_{jk} \frac{c_{k}^{7}}{7!} ) ] = \frac{1}{12!}}^{(72)  ∑ b_{i} [∑ a_{ij} \frac{c_{i}^{2}}{9 × 10} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}}^{(73)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}}{10} ( ∑ a_{jk} \frac{c_{k}^{8}}{8!} ) ] = \frac{1}{12!}}^{(74)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}}{10} ( ∑ a_{jk} \frac{c_{k}}{8} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}}^{(75)  ∑ b_{i} [∑ a_{ij} \frac{c_{j}}{10} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{7}}{7!}) ) ] = \frac{1}{12!}}^{(76)  ∑ b_{i} [∑ a_{ij} \frac{c_{i}}{10} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{lm} \frac{c_{m}^{6}}{6!} ) ) ) ] = \frac{1}{12!}}^{(77)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{9}}{9!} ) ] = \frac{1}{12!}}^{(78)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}^{2}}{8 × 9} ( ∑ a_{kl} \frac{c_{l}^{6}}{6!}) ) ] = \frac{1}{12!}}^{(79)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} \frac{c_{k}}{9} ( ∑ a_{kl} \frac{c_{l}^{7}}{7!}) ) ] = \frac{1}{12!}}^{(80)  ∑ b_{i} \frac{c_{i}}{11} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{l,m} \frac{c_{m}^{6}}{6!} ) ) ) ] = \frac{1}{12!}}

(81)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}^{8}}{8!}) ) ] = \frac{1}{12!}

(82)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} \frac{c_{l}}{8} ( ∑ a_{l,m} \frac{c_{m}^{6}}{6!} ) ) ) ] = \frac{1}{12!}

(83)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{l,m} \frac{c_{m}^{7}}{7!} ) ) ) ] = \frac{1}{12!}

(84)  ∑ b_{i} [∑ a_{ij} ( ∑ a_{jk} ( ∑ a_{kl} ( ∑ a_{l,m} ( ∑ a_{m,n} \frac{c_{n}^{6}}{6!} ) ) ) ) ] = \frac{1}{12!}.

References

[1] R. F. Arenstorf, Periodic solutions of the restricted three-body problem representing analytic continuations of Keplerian elliptic motions, Amer. J. Math. 85 (1963), 27–35.

[2] R. Barrio, Sensitivity analysis of ODEs/DAEs using the Taylor series method, SIAM J. Sci. Comput. 27(6) (2006), 1929–1947.

[3] R. Barrio, F. Blesa and M. Lara, VSVO formulation of the Taylor method for the numerical solution of ODEs, Comput. Math. Appl. 50 (2005), 93–111.

[4] J. Berntsen and T. O. Espelid, Error estimation in automatic quadrature routines, ACM Trans. Math. Software, 17 (1991), 233–255.

[5] J. C. Butcher, On Runge–Kutta processes of high order. J. Aust. Math. Soc. 18 (1964), 50–64.

[6] J. C. Butcher, The Numerical Analysis of Ordinary Differential Equations, Wiley, London, (1987).

[7] G. F. Corliss and Y. F. Chang, Solving ordinary differential equations using Taylor series, ACM Trans. Math. Software, 8(2) (1982), 114–144.

[8] P. J. Davis and P. Rabinowitz, Numerical Integration, Blaisdell, Waltham MA, 1967.

[9] A. Deprit and R. M. W. Zahar, Numerical integration of an orbit and its concomitant variations, Z. Angew. Math. Phys. 17 (1966), 425–430.

[10] J. R. Dormand, M. E. A. El-Mikkawy and P. J. Prince, Families of Runge–Kutta–Nystrom formulae, IMA J. Numer. Anal. 7(1987), 423–430.

[11] W. H. Enright and J. D. Pryce, Two FORTRAN packages for assessing initial value methods, ACM Trans. Math. Software, 13 (1987), 1–27.

[12] E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems, Section III.8, Springer-Verlag, Berlin, 1993.

[13] M. A. Hajji and R. Vaillancourt, Matrix derivation of Gaussian quadratures, Scientific Proceedings of Riga Technical University, 29, Boundary Field Problems and Computer Simulation, 48th issue, (2006), 198–213.

[14] J. Hoefkens, M. Berz, and K. Makino, Computing validated solutions of implicit differential equations, Adv. Comput. Math. 19 (2003), 231–253.

[15] T. E. Hull, W. H. Enright, B. M. Fellen and A. E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972), 603–637.

[16] M. Lara, A. Elipe and M. Palacios, Automatic programming of recurrent power series, Math. Comput. Simulation 49 (1999), 351–362.

[17] N. S. Nedialkov, K. R. Jackson, and G. F. Corliss, Validated solutions of initial value problems for ordinary differential equations, Appl. Math. Comput. 105 (1999), 21–68.

[18] T. Nguyen-Ba, E. Kengne and R. Vaillancourt, One-step 4-stage Hermite–Birkhoff–Taylor ODE solver of order 12, Can. Appl. Math. Q. 16(1) (2008), 77–94.

[19] T. Nguyen-Ba, V. Božić, E. Kengne and R. Vaillancourt, One-step 9-stage Hermite–Birkhoff– Taylor ODE solver of order 10, J. Appl. Math. Comput., DOI 10.1007/s12190-008-0216-3.

[20] T. Nguyen-Ba, V. Božić and R. Vaillancourt, One-step 7-stage Hermite–Birkhoff–Taylor ODE solver of order 13, Int. J. Pure Appl. Math. 43(4) (2008), 569–592.

[21] T. Nguyen-Ba, V. Božić, A. Przybylo and R. Vaillancourt, One-step 9-stage Hermite– Birkhoff–Taylor ODE solver of order 11, University Scientific J., Telecommunications and Electronics Series, University of Technology and Life Sciences (UTP), Bydgoszcz, Poland, 11 (2008), 33–52.

[22] T. Nguyen-Ba, V. Božić, E. Kengne and R. Vaillancourt, One-step 4-stage Hermite–Birkhoff– Taylor ODE solver of order 14, Scientific Proceedings of Riga Technical University, 33 Boundary Field Problems and Computer Simulation, 49th issue, (2007), 26–41.

[23] R. Piessens, E. de Doncker-Kapenga, C. W. Überhuber and D. K. Kahaner, QUADPACK, A subroutine package for automatic integration, Springer Series in Comput. Math. 1, 1983.

[24] P. J. Prince and J. R. Dormand, High order embedded Runge–Kutta formulae, J. Comput. Appl. Math. 7(1) (1981), 67–75.

[25] E. Rabe, Determination and survey of periodic Trojan orbits in the restricted problem of three bodies, Astron. J. 66(9) (1961),500–513.

[26] P. W. Sharp, Numerical comparison of explicit Runge–Kutta pairs of orders four through eight, ACM Trans. Math. Software 17 (1991), 387–409.

[27] J. F. Steffensen, On the restricted problem of three bodies, Danske Vid. Selsk., Mat.-fys. Medd. 30(18) (1956), 17 pp.

Department of Mathematics and Statistics (Received 01 11 2007) University of Ottawa (Revised 24 11 2008 and 23 07 2009) Ottawa Ontario Canada (all the authors)

trnguyen@uottawa.ca
vladan.bozic@gmail.com
ekengne6@yahoo.fr
remi@uottawa.ca