Numer Algor (2010) 53:153–170 DOI 10.1007/s11075-009-9285-0 ORIGINAL PAPER

Trees and numerical methods for ordinary differential equations J. C. Butcher

Received: 21 November 2008 / Accepted: 4 March 2009 / Published online: 14 March 2009 © Springer Science + Business Media, LLC 2009

Abstract This paper presents a review of the role played by trees in the theory of Runge–Kutta methods. The use of trees is in contrast to early publications on numerical methods, in which a deceptively simpler approach was used. This earlier approach is not only non-rigorous, but also incorrect. It is now known, for example, that methods can have different orders when applied to a single equation and when applied to a system of equations; the earlier approach cannot show this. Trees have a central role in the theory of Runge–Kutta methods and they also have applications to more general methods, involving multiple values and multiple stages. Keywords Runge-Kutta methods · Trees · Order conditions · Taylor expansions

1 Introduction Trees now play a central role in the theory of numerical methods for ordinary differential equations, and this paper surveys some aspects of this role. In Section 2 trees are introduced along with the formulations of autonomous and non-autonomous differential equations, as single first order equations and as first order systems. This is followed in turn by Section 3, where the theory of order of numerical methods is presented using each of two approaches. The first approach is traditional and dates from the work of Runge; the essential idea is to develop order conditions for the scalar problem y = f (x, y). There has always been an underlying, unstated, assumption that a method of some

J. C. Butcher (B) The University of Auckland, Auckland, New Zealand e-mail: [email protected]

154

Numer Algor (2010) 53:153–170

specific order for a scalar problem will have this same order for a system of differential equations. This assumption is, however, incorrect as we see by a comparison with the second approach. Here only autonomous problems are considered, for formal simplicity, but the problem is allowed to have arbitrary dimension. This leads to a formulation of the order conditions in terms of trees and elementary weights, which are functions of trees. In Section 4, fourth order methods are discussed and this leads to the construction of a method which is fifth order under the first approach, but only fourth order under the second approach. Finally, in Section 5, extensions of the theory are introduced. This leads to the more general effective order concept for Runge–Kutta methods. It also leads to a tree-based approach to investigations of the order of multivalue methods.

2 Trees, differential equations and Runge–Kutta methods 2.1 Trees Throughout this paper, “Trees” will mean “Rooted Trees”. Both of these are simple ideas: a (non-rooted) tree is simply a connected graph with no cycles, and a (rooted) tree differs only in that one of its vertices is designated to be the root. The convention adopted in this paper is to indicate the root by placing it at the bottom of a tree diagram. To illustrate the idea, the rooted trees with up to four vertices are shown in Fig. 1. Familiar applications, in genealogies and the representation of nested arithmetic operations, are illustrated in Fig. 2. For non-commutative arithmetic operations, such as subtraction and division, the order in which the children, representing operands, of a particular parent, representing an operation, are drawn has a significance. However, for the applications of trees specific to this paper, the order in which subtrees are written has no significance. It is convenient to define τ as the tree and to write t = [t1 t2 · · · tν ] as the tree formed by adjoining the roots of the trees t1 , t2 , . . . , tν to a new vertex which becomes the root of t. If t1 , t2 , . . . , tm are distinct, and occur respectively k1 , k2 , . . . , km times amongst t1 , t2 , . . . , tν , then we will also write   km , t = t1k1 t2k2 · · · tm m where we note that i=1 ki = ν. This notation for representing trees makes it possible to define various functions on the set of all trees, which we will write

Fig. 1 Trees with up to 4 vertices

Numer Algor (2010) 53:153–170

155

Fig. 2 Some applications of trees

as T, recursively. In particular, the order r(t), which is simply the number of vertices in t, is generated by the rules r(τ ) = 1, r([t1 t2 · · · tν ]) = 1 +

ν 

r(ti ),

i=1

and the symmetry σ (t), which is the order of the automorphism group of t, by the recursion σ (τ ) = 1. m    km σ t1k1 t2k2 · · · tm ki !σ (ti )ki . = i=1

The so-called “density” γ (t) is defined by γ (τ ) = 1, γ ([t1 t2 · · · tν ]) = r([t1 t2 · · · tν ])

ν 

γ (ti ).

i=1

These functions are easy to compute for trees of order up to 4 and they are given, along with their notation, in Table 1. 2.2 Differential equations Throughout the paper, “differential equations” are initial value problems such as y (x) = f (y(x)),

y(x0 ) = y0 ,

f : RN → RN

(1)

Table 1 Trees of up to order 4 with their corresponding notation and some recursively defined functions t notation r(t) σ (t) γ (t)

τ

[τ ]

[τ 2 ]

[[τ ]]

[τ 3 ]

[τ [τ ]]

[[τ 2 ]]

[[[τ ]]]

1 1 1

2 1 2

3 2 3

3 1 6

4 6 4

4 1 8

4 2 12

4 1 24

156

Numer Algor (2010) 53:153–170

or y (x) = f (x, y(x)), y(x0 ) = y0 , f : R × R N → R N .

(2)

It looks deceptively easy to assume that N = 1, but this does not capture all the flavour of the problem, even for (2). If we make no such assumption, then (1) is simpler than (2) and really just as general, because we can always add an additional differential equation whose solution is x. 2.3 Runge–Kutta methods In terms of the autonomous initial value problem, (1), the aim of a Runge– Kutta method is to calculate approximations y1 ≈ y(x1 ), y2 ≈ y(x2 ), . . . . The basic example is the simple Euler method yn = yn−1 + hf (yn−1 ),

h = xn − xn−1 ,

for n = 1, 2, . . . .This can be made more accurate (as h → 0) by using a Runge– Kutta method based on the mid-point or the trapezoidal rule quadrature formula:

1 yn = yn−1 + hf yn−1 + hf (yn−1 ) , 2 1 1 yn = yn−1 + hf (yn−1 ) + hf (yn−1 + hf (yn−1 )) . 2 2 These methods, from Runge’s 1895 paper [13], are “second order”, because the error introduced in a single step behaves like O(h3 ). At a specific output point the error is O(h2 ). It is recognised that for accurate calculations, high order is better than low order and Runge’s work was followed by contributions by Heun [9], Kutta [11], Nyström [12], Hutˇa [10] and others, pushing the available order up to 4, 5, 6 and higher. Our aim now will be to formulate a typical step in an “s stage method” and then find criteria for this method to have a specific order. In carrying out a step we evaluate s stage values Y1 ,

Y2 ,

...,

Ys

F1 ,

F2 ,

...,

Fs ,

and s stage derivatives

where Fi = f (Yi ). Each Yi is found as a linear combination of the F j, added on to y0 : Y i = y0 + h

s  j=1

aij F j ≈ y(x0 + ci h)

Numer Algor (2010) 53:153–170

157

and the approximation at x1 = x0 + h is determined by s  y1 = y0 + h b i Fi ≈ y(x0 + h). i=1

We represent the method by a tableau: c1 a11 c2 a21 .. .. . . cs as1 b1

a12 · · · a1s a22 · · · a2s .. .. . . as2 · · · ass b2 · · · bs

or, if the method is “explicit”, by a lower-triangular tableau 0 c2 a21 c3 a31 a32 .. .. .. . . . . . . cs as1 as2 · · · as,s−1 b 1 b 2 · · · b s−1 b s The use of tableaux to represent specific methods can be illustrated by the three examples that have been introduced: the Euler method, the mid-point rule method and the trapezoidal rule method respectively. 0 1

0 1 1 2 2 0 1

0 1 1 1 1 2 2

3 Order of Runge–Kutta methods 3.1 The first approach to the order of Runge–Kutta methods We introduce two approaches to finding the order of a Runge–Kutta method. The first of these has a distinguished history and was used by Runge, Heun, Kutta, Nyström, Huˇta and others to obtain methods up to order 6. The starting point is the first order differential Eq. 2 with N = 1 and the first aim is to find the first few terms in the Taylor series about x0 . Hence we calculate the second and third derivatives: y = fx + f f y , y = fxx + 2 f fxy + f 2 f yy + f f y2 . Once we have the Taylor series for y(x0 + h) we can also find the series for a numerical approximation and a comparison of these gives a number of “order conditions”. For order p write the number of conditions derived in this way as Cp .

158

Numer Algor (2010) 53:153–170

The number of parameters in an explicit Runge–Kutta method with s stages is s(s + 1)/2. It might be expected that s should be chosen so that there are at least as many parameters as conditions. We show this comparison in Table 2 with the value of s chosen as the fewest number of stages satisfying this assumption, for each order p. Surprisingly, order 6 can be achieved with s = 7, where we note that this corresponds to solving 31 equations in 28 unknowns. Even more surprisingly, the first published method with order 6 [10], derived to satisfy these 31 conditions, actually satisfies the 37 order 6 conditions to be introduced in the second approach analysis. This was shown in [2]. 3.2 The second approach This description of the second approach uses the ideas introduced in [1]. We make two changes from the first approach; these are 1. Remove x from f (x, y) so that the problem reverts to the autonomous formulation (1). 2. Regard y as a vector-valued function so that, in general, N > 1. The first change seems to make everything simpler, because fx , fxy , and other terms involving partial derivatives with respect to x, just disappear. The second change seems to make everything more complicated, because terms such as f 2 f y f yy are replaced by more than one term when N > 1, even though they are identical in one dimension. We are now faced with complicated combinations of linear and multilinear operators and we need to find a convenient way of writing these. It is possible to work in terms of tensors, but it is preferable to use Fréchet derivatives. Written this way, the first, second and third derivatives are as follows, where the one-dimensional counterparts are also shown for comparison. y = f,

y = f

y = fx + f f y ,

y = f  f

y = fxx + 2 f fxy + f 2 f yy + f f y2 ,

y = f  (f, f) + f  f  f

In these formulae, f = f (y(x)), f  = f  (y(x)), f  = f  (y(x)). Furthermore, in notations such as f  (f, f), each of the two f vectors is one of the operands of Table 2 Number of order constraints with N = 1 compared with numbers of free parameters

p

Cp

s

s(s+1) 2

1 2 3 4 5 6

1 2 4 8 16 31

1 2 3 4 6 8

1 3 6 10 21 36

Numer Algor (2010) 53:153–170

159

the bi-linear operator f  , with a similar remark applying to other multi-linear operators f  , f (4) , . . . . If we count the number of order conditions using this second formulation, denoted by Dp , and include Cp for comparison, we obtain Table 3. Note that C7 , C8 , . . . are missing from this table, but their values and those of higher members of the C sequence do not seem to have any significant applications. In the second approach, we again list the formulae for the derivatives; but this time we go as far as order 4. Note the tree-like structures: y = f, 

f f f f f 



y = f f, f

y = f  (f, f)

f f f

 

+ f f f, y(4) = f  (f, f, f)

f

f

(3)

f f



f

+ 3f  (f, f  f)

f f

f 

 

+ f f (f, f) + f f f f

f

f f  f

f f f f

The plan now is to formulate Runge–Kutta methods in a systematic way and to explore the order conditions. This will all be done using the second approach. Later, we will to come back to the first approach and show why it is not completely satisfactory. Before we go on, note that the numbers Dp , is just the number of rooted trees with no more than p vertices. 3.3 Elementary differentials To find the Taylor expansion of the exact solution we need the higher derivatives of y evaluated at x0 . As we have seen, this becomes increasingly complicated the further we go and accordingly we will make use of the systematic pattern illustrated in (3). As in the low order examples, the various

Table 3 Comparison of order condition enumerations

p

Cp

Dp

1 2 3 4 5 6 7 8

1 2 4 8 16 31

1 2 4 8 17 37 85 200

160

Numer Algor (2010) 53:153–170

terms have a structure related to rooted-trees. These expressions are known as “elementary differentials”, and can be defined recursively as follows: F(τ ) = f, F([t1 t2 . . . tν ]) = f (n) (F(t1 ), F(t2 ), . . . , F(tν )). Note that this definition of elementary differentials is slightly different from standard terminology (see for example [5]) in that if F(t) is the usual elementary differential, then F(t) = F(t)(y0 ). 3.4 Taylor series for exact solution The formal Taylor expansion of the solution at x0 + h is y(x0 + h) = y0 +

 t∈T

hr(t) F(t). σ (t)γ (t)

(4)

To justify this result, consider the equation η = y0 1 + hL( f ◦ η),

(5)

where η is a continuous function on [0, 1], 1 is defined by 1(ξ ) = 1 for ξ ∈ [0, 1], L is a linear operator and it is assumed that f satisfies a Lipschitz condition. In the case where L is defined by  ξ φ ξ dξ , (6) L(φ)(ξ ) = 0

Equation (5) becomes the Picard integral equation with solution η(ξ ) = y(x0 + hξ ). For given L, define (t) by the recursion (τ ) = L(1), ([t1 , t3 , . . . , tν ]) = L( (t1 ) · (t2 ) · · · · · (tν )), where · denotes pointwise product. Because f satisfies a Lipschitz condition, (5) can be solved by functional iteration, if |h| is sufficiently small. Define η[0] = y0 1 and, for k = 1, 2, . . . , define η[k] as the polynomial approximation to η, found by truncating to hk terms the value of y0 1 + hL( f ◦ η[k−1] ). The first few approximations are η[1] = y0 1 + h (τ )F(τ ), η[2] = y0 1 + h (τ )F(τ ) + h2 ([τ ])F([τ ]), η[3] = y0 1 + h (τ )F(τ ) + h2 ([τ ])F([τ ])

1  2   2  +h3 τ F τ + ([[τ ]])F([[τ ]]) 2

Numer Algor (2010) 53:153–170

161

and in general η[k] = y0 1 +



hr(t) (t)F(t)/σ (t).

r(t)≤k

By taking the limit as k → ∞, we arrive at the formal Taylor expansion for the solution to (5). To derive (4), evaluate the formal series for η(ξ ) at ξ = 1 in the case (6). The last step is to show that (t)(ξ ) = ξ r(t) /γ (t), and this follows from (τ )(ξ ) = ξ and the recursion, when t = [t1 , t2 , . . . , tν ], 

ξ

(t) =

ξ

r(t)−1

dξ

0

 ν

 γ (tk ) = ξ r(t) γ (t).

k=1

This expression is 1/γ (t), when ξ = 1. Our aim will now be to find a corresponding formula for the result computed by one step of a Runge-Kutta method. By comparing these formulae, term by term, we will obtain conditions for a specific order of accuracy. 3.5 Taylor expansion for numerical approximation We need to evaluate various expressions which depend on the tableau for a particular method. These are known as “elementary weights”. For a given tree,

i (t) will denote the elementary weight associated with stage number i and

(t) will denote the elementary weight associated with the output. These are defined together by the recursions

i (τ ) = ci ,

(τ ) =

s 

i = 1, 2, . . . , s, b i,

i=1

i ([t1 t2 . . . tν ]) =

s 

aij

j=1

([t1 t2 . . . tν ]) =

s 

ν 

j(tk ),

i = 1, 2, . . . , s,

k=1

bi

i=1

ν 

i (tk ).

k=1

Expressions for (t) up to order 4 are shown in Table 4. The formal Taylor expansion of the computed solution at x0 + h is y1 = y0 +

 hr(t)

(t)F(t)(y0 ). σ (t) t∈T

162

Numer Algor (2010) 53:153–170

Table 4 Expressions for (t) t

(t)



bi



b i ci





b i ci2



b i aij c j

b i ci3



b i ci aij c j



b i aij c2j



b i aij a jk ck

To justify this result, we again make use of (5) but in this case φ is the set of functions on {1, 2, . . . , s, s + 1}, and L is defined by  (Lφ)(i) = aijφ( j), i = 1, 2, . . . , s, j=1

(Lφ)(s + 1) =



b i φ(i).

i=1

In this case, the output from a single step of a Runge–Kutta method will be y1 = η(s + 1) and the Taylor expansion of y1 will be  y1 = y0 + hr(t) (t)(s + 1)F(t)/σ (t). t∈T

Starting from

(τ )(i) = L(1)(i) =

⎧ s  ⎪ ⎪ ⎪ aij = ci , ⎪ ⎨

i ≤ s,

⎪ ⎪ ⎪ ⎪ ⎩

i = s + 1,

j=1 s 

b j,

j=1

it can be verified by induction that (t)(i) = i (t),

i ≤ s,

(t)(s + 1) = (t). 3.6 Order conditions To match the two Taylor series y(x0 + h) = y0 +

 t∈T

y1 = y0 +

hr(t) F(t)(y0 ), σ (t)γ (t)

 hr(t)

(t)F(t)(y0 ), σ (t) t∈T

up to h p terms we need to ensure that

(t) =

1 , γ (t)

Numer Algor (2010) 53:153–170

163

for all trees such that r(t) ≤ p. These are the “order conditions” and are at the heart of the study of Runge– Kutta methods. In the next section we will consider some families of fourth and fifth order methods.

4 Fourth and fifth order methods 4.1 Fourth order methods The definitive work on fourth order methods is Kutta’s 1901 paper [11]. Kutta found several families of methods; in the present paper only one of these families, based on Simpson’s rule, will be considered in detail. The starting point is the tableau: 0 1 2 1 2 1

1 2 1 − a32 2 1 − a42 − a43 1 6

a32 a42

a43

2 − b3 b3 3

1 6

Because the underlying quadrature formula has order 4, some order conditions are automatically satisfied; namely those for the trees

leaving, still to be considered, those associated with the trees

It is necessary, in fourth order methods with 4 stages, that



b j 1 − c j − i<4 b i aij a4 j = b −1 4

(7)

and it also implies that the conditions for three further trees become automatically satisfied. These are

164

Numer Algor (2010) 53:153–170

With (7) substituted into the tableau, the method becomes. 0 1 1 2 2 1 1 − a32 a32 2 2 1 6b 3 a32 −1 2−3b 3 −6b 3 a32 3b 3 2 1 − b3 b3 6 3 The remaining order condition is 

b i aijcj =

. 1 6

1 , 8

which reduces to b 3 (1 − a32 ) =

1 . 6

Substitute a32 = 1 − 1/6b 3 , and the method becomes 0 1 2 1 2 1

1 2 1 1 1 1− − 6b 3 2 6b 3 6b 3 − 2 3 − 9b 3 1 2 − b3 6 3

. 3b 3 b3

1 6

The value of b 3 = 0 is at our disposal. √ In 1951, Gill [6] proposed a new value b 3 = 13 + 62 , because this was more efficient in terms of memory than the classical choice proposed by Kutta 50 years earlier. But everyone seems to like the classical method of Kutta, with b 3 = 13 . It is simple, easy to remember and works well: 0 1 2 1 2 1

1 2 0 0 1 6

1 2 0 1 3

1 1 3

1 6

4.2 A method of order 4 or 5 When we were counting order conditions, we saw that, in the “First approach”, there were 16 conditions for order 5 but, in the “Second approach”, there

Numer Algor (2010) 53:153–170

165

were 17 conditions. The reason is quite simple: the two elementary differentials arising in the second approach

f  f, f  f  f

and



f  f  f, f  f

correspond to identical terms



f f yy + fxy f y f f y + fx ,

in the first approach formulation. The consequence is a confusion of two order conditions 

b i ci aija jk ck =

1 30

and

 1 b i aijc ja jk ck = , 40

(8)

which become the single condition   1 1 + . b i ci aija jk ck + b i aijc ja jk ck = 30 40

(9)

In a short 1995 paper [4], a method was constructed which satisfies all the equations for order 5 in the first approach but not for the second approach. A new method satisfying the same requirements is given by the tableau 0 1 10 1 5 3 10 3 4

1 10 11 260 3 − 170 13035 544 165709 1 − 918 5 54

41 260 63 9 340 68 75447 9009 − 1088 136 733747 638689 − 1326 1122 0

0

1287 64 32053 5320 − 162 11583 250 32 567 81

−

1 14

For this method, all order 5 conditions are satisfied except those appearing in (8). These are replaced by  1 − b i ci aijajk ck = 30  1 b i aijcjajk ck = + 40

127 , 20800 127 , 20800

166

Numer Algor (2010) 53:153–170

Fig. 3 Solutions to (10) and (11)

so that (9) is still satisfied. A suitable single differential equation to test the order of convergence of this method, together with a closely related autonomous system, are dy y−x = , dx y+x     1 d y y−x , = dt x x2 + y2 y + x

y(0) = 1, 

   y(1) 1 = . x(1) 0

(10) (11)

The solution of (10), written in parametric coordinates, is x = t sin(ln(t)), y = t cos(ln(t)), and this exactly the solution to (11). In numerical tests, we will approximate the solution when t = x = exp(π/2), y = 0. The solution on the interval t ∈ [0, exp(π/2)] is shown in Fig. 3. In two experiments, the same method was used to solve each of the problems and to find the absolute values of the computational errors in y(exp(π/2)), found by solving (10) for a range of values of n, the number of equally spaced steps. Denote this quantity by En(0) . Similarly, using the formulation (11), quantities En(1) and En(2) are computed as the absolute values of the errors for the two components y(exp(π/2)) and x(exp(π/2)) respectively. Since our purpose is to assess the rate of convergence as n increases, it is (i) appropriate to calculate En(i) /E2n for i = 0, 1, 2. If the orders are as claimed in this paper, the values of these ratios for i = 0 should become close to 32 as n increases and for i = 1 and i = 2, they should approximate 16. The actual

Table 5 Numerical calculations illustrating the different convergence rates for problems (10) and (11) (0) En(0) /E2n (1) (1) En /E2n (2) (2) En /E2n

n = 20

n = 40

n = 80

n = 160

n = 320

33.5068

33.0459

32.8983

32.5659

32.4374

13.4076

14.8227

15.4446

15.7309

15.8740

14.8548

15.5186

15.7833

15.8971

15.9678

Numer Algor (2010) 53:153–170

167

values are given in Table 5 for values of n from 20 up to 320, in powers of 2, and this claim is supported by these results.

5 Extensions of the theory 5.1 Algebraic interpretation We will introduce an algebraic system which represents individual RungeKutta methods and also compositions of methods. This centres on the meaning of order for Runge-Kutta methods and leads to the possible generalisation of “effective order”. Denote by G1 the group consisting of mappings of (rooted) trees to real numbers, where the group operation is defined according to the algebraic theory of Runge-Kutta methods [3, 5], or to the (equivalent) theory of B-series [7, 8]. We will illustrate this operation in Table 6, where we also introduce the special member E ∈ G, defined by E(t) = 1/γ (t). Hp will denote the normal subgroup defined by t → 0 for r(t) ≤ p. If α ∈ G1 then this maps canonically to α Hp ∈ G1 /Hp . If α is defined from the elementary weights for a Runge-Kutta method then order p can be written as α Hp = EHp . 5.2 Effective order of methods Effective order p is defined by the existence of β such that βα Hp = Eβ Hp . The computational interpretation of this idea is that we can carry out a starting step corresponding to β and a finishing step corresponding to β −1 , with Table 6 Group multiplication table

r(ti ) i

ti α(ti ) β(ti ) (αβ)(ti )

1

1

α1

β1

α1 + β1

2

2

α2

β2

α2 + α1 β1 + β2

3

3

α3

β3

α3 + α12 β1 + 2α1 β2 + β3

3

4

α4

β4

α4 + α2 β1 + α1 β2 + β4

4

5

α5

β5

α5 + α13 β1 + 3α12 β2 + 3α1 β3 + β5

4

6

α6

β6

α6 + α1 α2 β1 + (α12 + α2 )β2

E(ti ) 1 1 2 1 3 1 6 1 4 1 8

+ α1 (β3 + β4 ) + β6 4

7

α7

β7

α7 + α3 β1 + α12 β2 + 2α1 β4 + β7

4

8

α8

β8

α8 + α4 β1 + α2 β2 + α1 β4 + β8

1 12 1 24

168

Numer Algor (2010) 53:153–170

many steps in between corresponding to α. This is equivalent to many steps all corresponding to βαβ −1 , which has classical order p. Thus, the benefits of high order can be enjoyed by methods with high effective order. We analyse the conditions for effective order 4. Without loss of generality assume β(t1 ) = 0. i 1

(βα)(ti ) α1

2

β2 + α2

3

β3 + α3

4

β4 + β2 α1 + α4

5

β5 + α5

6

β6 + β2 α2 + α6

7

β7 + β3 α1 + α7

8

β8 + β4 α1 + β2 α2 + α8

(Eβ)(ti ) 1 1 + β2 2 1 + 2β2 + β3 3 1 + β2 + β4 6 1 + 3β2 + 3β3 + β5 4 1 3 + β2 + β3 + β4 + β6 8 2 1 + β2 + 2β4 + β7 12 1 1 + β2 + β4 + β8 24 2

Of these eight conditions, only five are conditions on α, while the remaining three are conditions on β, to be imposed once α is known. The five α order conditions, written in terms of the Runge-Kutta tableau, are  

b i = 1,

1 , 2  1 b i aijcj = , 6  1 b i aijajk ck = , 24   1 b i aijcj(2ci − cj) = . b i ci2 (1 − ci ) + 4 b i ci =

It is possible, in principle, to optimize in terms of metrics of the vector of fifth effective order error constants. Although it is not known if this can lead to more accurate methods, a similar question for the effective order three case leads to the interesting result that each of the error coefficients, except for tree number 8, can be made zero. When a similar analysis is carried out for effective order 5, it is found that this can be achieved with only 5 stages even though classical order 5 cannot be found with less than 6 stages.

Numer Algor (2010) 53:153–170

169

5.3 Application to multivalue methods The quantities involved in the computations carried out in the first step of a multivalue method have Taylor expansions about the initial value of the form β0 y(x0 ) + O(h), where the O(h) term can be further expanded in terms of elementary differentials. In the case β0 = 1, it is possible to represent the corresponding numerical quantity by a member of G1 . To allow for greater generality, the vector space G is introduced in which T is augmented by the addition of the empty tree which will be written as ∅. For β ∈ G, we can re-calculate Table 6, where now β0 = β(∅). It is found that β0 now appears as the coefficient of αi in the formula for (αβ)(ti ). It is now possible to define D ∈ G by the values D(∅) = 0, D(τ ) = 1, D(t) = 0 for r(t) > 1. The significance of D is that if α ∈ G1 represents y1 = y0 + O(h), then α D represents hf (y1 ). Write a multivalue, or general linear, method using a matrix quadruple (A, U, B, V), as in [5]. The vectors of stage values Y and stage derivatives F are related to the vectors of incoming approximations y[n−1] and outgoing approximations y[n] by the equations      Y AU hF = . y[n] BV y[n−1] In Fig. 4, order for this type of method is shown in terms of a starting method S which produces exactly the right input so that, to within order p, the result of carrying out a single step is the same as moving forward a unit stepsize and then applying the starting process from this point. If E, in this context, represents the exact solution through a standard step, and M represents one step of the method, then the order condition “M has order p relative to S” is represented by equating SM and ES to within O(h p+1 ). To carry out a Taylor series analysis of this condition, write ξ as an r dimensional vector of elements of G to represent S, and write η as an s dimensional vector of elements of G1 to represent the stage values. If Eξ and ηD are interpreted as

Fig. 4 Schema showing error of M relative to starting method S

170

Numer Algor (2010) 53:153–170

scalar multiplications of vectors, then up to trees of order p, the relative order can be written as      η AU ηD = . Eξ BV ξ This analysis can be applied generally to specify order conditions for a wide variety of general linear methods. In particular it can be used to analyze the order of two-step Runge–Kutta methods. Acknowledgement The author expresses his gratitude to an anonymous referee whose comments have led to improvements in the paper.

References 1. Butcher, J.C.: Coefficients for the study of Runge–Kutta integration processes. J. Aust. Math. Soc. 3, 185–201 (1963) 2. Butcher, J.C.: On the integration processes of A. Hutˇa. J. Aust. Math. Soc. 3, 202–206 (1963) 3. Butcher, J.C.: An algebraic theory of integration methods. Math. Comp. 26, 79–106 (1972) 4. Butcher, J.C.: On fifth order Runge–Kutta methods. BIT 35, 202–209 (1995) 5. Butcher, J.C.: Numerical Methods for Ordinary Differential Equations, 2nd edn. Wiley, Chichester (2008) 6. Gill, S.: A process for the step-by-step integration of differential equations in an automatic digital computing machine. Proc. Camb. Philos. Soc. 47, 96–108 (1951) 7. Hairer, E., Nørsett, S.P., Wanner, G.: Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd edn. Springer, Berlin (1993) 8. Hairer, E., Wanner, G.: On the Butcher group and general multi-value methods. Computing 13, 1–15 (1974) 9. Heun, K.: Neue methoden zur approximativen integration der differentialgleichungen einer unabhängigen Veränderlichen. Z. Math. Phys. 45, 23–38 (1900) 10. Hutˇa, A.: Une amélioration de la méthode de Runge–Kutta–Nyström pour la résolution numérique des équations différentielles du premier ordre. Acta Fac. Nat. Univ. Comen. Math. 1, 201–224 (1956) 11. Kutta, W.: Beitrag zur näherungsweisen integration totaler differentialgleichungen. Z. Math. Phys. 46, 435–453 (1901) 12. Nyström, E.J.: Über die numerische integration von differentialgleichungen. Acta Soc. Scient. Fenn. 50(13), 55pp. (1925) 13. Runge, C.: Über die numerische Auflösung von differentialgleichungen. Math. Ann. 46, 167–178 (1895)