Applied Mathematics
Runge{Kutta Methods, Trees, and Maple
On a Simple Proof of Butcher's Theorem and the
Automatic Generation of Order Conditions
Folkmar Bornemann
Center of Mathematical Sciences, Munich University of Technology, 80290 Munich, Germany e-mail:bornemann@ma.tum.de
February 9, 2001
Summary.
This paper presents a simple and elementary proof of Butcher's theorem on the order conditions of Runge-Kutta methods. It is based on a recursive denition of rooted trees and avoids com-binatorial tools such as labelings and Faa di Bruno's formula. This strictly recursive approach can easily and elegantly be implemented using modern computer algebra systems like Maple for automatically generating the order conditions. The full, but short source code is presented and applied to some instructive examples.Key words:
numerical solution of ODEs, Runge-Kutta methods, recursive representation of rooted trees, Butcher's theorem, automat-ic generation of order conditions, computer algebra systems.Mathematics Subject Classication (1991): 65-01, 65L06, 65Y99
Vol. 2, No. 1, pp. 3{15, 2001
1. Introduction
A rst step towards the construction of Runge-Kutta methods is the calculation of the order conditions that the coe cients have to obey. In the old days they were obtained by expanding the error term in a Taylor series by hand, a procedure which for higher orders sooner or later runs into di culties because of the largely increasing com-binatorial complexity. It was a major break-through when Butch-er 1] published his result of systematically describing ordButch-er condi-tions by rooted trees. The proof of this result has evolved very much in meantime, mainly under the inuence of Butcher's later work 3]
and the contributions of Hairer and Wanner 7,6]. In this paper we will present a simple and elementary proof of Butcher's theorem by using very consequently the recursive structure of rooted trees. This way we avoid lengthy calculations of combinatorial coe cients, the use of tree-labelings, or Faa di Bruno's formula as in 3,6]. Our proof is very similar in spirit to the presentation of B-series by Hairer in
Chapter 2 of his lecture notes 5].
As early as 1976 Jenks 8] posed the problem of automatically generating order conditions for Runge-Kutta methods using comput-er algebra systems, but no replies wcomput-ere received. In 1988 Keipcomput-er of Wolfram Research concluded that the method of automatically cal-culating Taylor expansions by brute force was bound to be very ine -cient. Naturallyhe turned to the elegant results of Butcher's. Utilizing them, he wrote the Mathematica packageButcher.m, which has been
available as part of the standard distribution of Mathematica since then. This package was later considerably improved by Sofroniou 9] and oers a lot of sophisticated tools.
While teaching the simple proof of Butcher's result in a rst course on numerical ODEs, the author realized that the underlying recursive structure could also be exploited for a simple and elegant computer implementation. This approach diers from the work of Sofroniou in various respects. We will present the full source code in Maple and some applications.
2. Runge-Kutta methods
The Runge-Kutta methods are one-step discretizations of initial-value problems for systems of dordinary dierential equations,
x 0= f(tx) x(t 0) = x 0
where the right-hand sidef : t 0
T] R Rd!Rd is assumed
to be su ciently smooth. The continuous evolution x(t) = tt 0
x 0
of the initial-value problem is approximated in steps of length by t
+t
x t +t
x. This discrete evolution t +t
x is dened as an
approximation of the integral-equation representation
t+t x=x+ Z t + t f( t x)d
by appropriate quadrature formulas: (1) t +t x=x+ s X i=1 biki ki=f 0 @ t+cix+ s X j=1 aijkj 1 A :
The vectors ki 2 Rd, i = 1:::s, are called stages, s is the stage
number. Following the standard notation, we collect the coe cients of the method into a matrix and two vectors
A= (aij)ij 2R ss b= (b 1 :::bs) T 2R s c= (c 1 :::cs) T 2R s :
The method is explicit, ifAis strictly lower triangular. The method
hasorder p2N, if the error term expands to t+t x; t+t x=O( p+1) :
In terms of the Taylor expansion of the error;, the vanishing of
all lower order terms in just denes the conditions which have to
be satised by the coe cients A,band cof a Runge-Kutta method.
If we choose ci = s X j=1 aij
it can be shown 6], that there is no loss of generality in considering
autonomous systems only, i.e., those with no dependence of f on t.
Doing so, the expressionst +t
xandt +t
xare likewise independent
of t. We will writex and x for short, calling them theowand
thediscrete owof the continuous resp. the discrete system.
3. Elementary dierentials and rooted trees
The Taylor expansions of both the phase ow x and the discrete
owxare linear combinations of elementarydierentials like f 000( f 0 ff 0 ff) = X ijklm @ 3 f @xi@xj@xk @fi @xl fl @fj @xm fmfk:
We will use the short multilinear notation of the left hand side for the rest of the paper.
An elementary dierential can be expressed uniquely by the struc-ture of how the subterms enter the multilinear maps. For instance, looking at the expressionf
000( f 0 ff 0 ff) we observe thatf 000must be
a third derivative, since three arguments make it three-linear. This structure can be expressed in general byrooted trees, e.g.,
; ; @ @ e r r r r expressing f 0 f 00( ff) ; ; @ @ e r r r r expressing f 00( f 0 ff):
Every node with n children denotes a nth derivative of f, which is
applied as a multilinear map to further elementary dierentials, ac-cording to the structure of the tree. We start reading o this structure by looking at the root. This denes a recursive procedure, if we ob-serve the following: Having removed the root and its edges, a rooted tree decomposes into rooted subtrees
1
:::n with strictly less
nodes. The roots of the subtrees 1
:::nare exactly thenchildren
of's root. This way a rooted tree can be dened as theunordered
list of its successors
(2) =
1
:::n] #= 1 + # 1+
:::+ #n:
Here, we denote by # theorder of a rooted tree, i.e., the number
of its nodes. The root itself can be identied with the empty list,
= ].
An application of this procedure shows for the examples above thatf 0( f 00( ff)) is expressed by ]] andf 00( f 0( f)f) is expressed
by ]]. The reader will observe the perfect matching of
paren-theses and commas. In general the relation between a rooted tree
= 1
:::n] and its corresponding elementary dierentialf ()( x) is recursively dened by f ()( x) =f (n)( x) f ( 1 )( x):::f (n)( x) :
The dot of multiplication denotes the multilinear application of the derivative to the ngiven arguments. Due to the symmetry of the n
-linear mapf
(n), the order of the subtrees
1
:::ndoes not matter,
which means, that f
() depends in a well-dened way on
as an
unordered list only.
From = ] we deduce f () =
f. Analogously, each of the
re-cursive denitions in the following will have a well-dened meaning if applied to the single root = ], mostly by using the reasonable
convention that empty products evaluate to one and empty sums to zero|a convention that is also observed by most computer algebra systems.
4. A simple proof of Butcher's theorem
We are now in a position to calculate and denote the Taylor expansion of the continuous ow in a clear and compact fashion.
Lemma 1.
Given f 2Cp(Rd) the ow x expands to x=x+ X #p # ! f ()( x) +O(p +1) :The coe cients!and belonging to a rooted tree = 1
:::n]
are recursively dened by
(3) ! = (#) 1! :::n! = n! 1 ::: n :
By we denote the number of dierent ordered tuples ( 1
:::n)
which correspond to the same unordered list = 1
:::n].
Proof. The assertion is obviously true for p = 0. We proceed by
induction on p. Using the assertion for p, the multivariate Taylor
formula and the multilinearity of the derivatives we obtain
f(x)=f 0 @ x+ X #p # ! f ()+ O(p +1) 1 A =Xp n=0 1 n! f (n) 0 @ X # 1 p # 1 1! 1 f ( 1 ) ::: X #np # n n! n f (n) 1 A +O(p +1) =Xp n=0 1 n! X #1+:::+#np #1+:::+# n 1! :::n! 1 ::: n f (n) f ( 1 ) :::f (n) +O(p +1) =Xp n=0 X = 1::: n] #p+1 # #;1 ! n! 1 ::: n | {z } = f ()+ O(p +1) = X #p+1 # #;1 ! f ()+ O(p +1) :
Plugging this into the integral form of the initial value problem we obtain x=x+ Z 0 f( x)d=x+ X #p+1 # ! f ()+ O( p+2)
which proves the assertion forp+ 1. ut
A likewise clear and compact expression can be calculated for the Taylor expansion of the discrete ow.
Lemma 2.
Given f 2Cp(Rd) the discrete ow x expands to x=x+ X #p # bTA () f ()( x) +O(p +1) : The vector A () 2Rs, = 1 :::n], is recursively dened by (4) A () i = AA (1) i::: AA (n) i i= 1:::s:Proof. Because of the denition (1)of the discrete ow we have to prove that the stageski expand to
ki= X #p #;1 A () i f ()+ O(p):
This is obviously the case for p= 0. We proceed by induction on p.
Using the assertion for p, the denition of the stages ki, the
multi-variate Taylor formula and the multilinearity of the derivatives we obtain ki=f 0 @ x+ 0 @ X #p #;1 AA () i f ()+ O(p) 1 A 1 A =Xp n=0 1 n! f (n) 0 @ X #1p # 1 1 AA ( 1 ) i f ( 1 ) ::: ::: X # n p # n n AA (n) i f (n) 1 A + O( p+1) =Xp n=0 1 n! X # 1 +:::+# n p # 1 +:::+# n 1 ::: n AA ( 1 ) i ::: AA ( n ) i f (n) f ( 1 ) :::f ( n ) +O(p +1) =Xp n=0 X = 1::: n] #p+1 #;1 n! 1 ::: n | {z } = A () i f ()+ O(p +1) = X #p+1 #;1 A () i f ()+ O( p+1)
which proves the assertion forp+ 1. ut
Comparing the coe cients of the elementary dierentials in the expansions of both the phase ow and the discrete ow, we immedi-ately obtain Butcher's theorem 1].
Theorem 1 (Butcher 1963).
A Runge-Kutta-method (bA) is oforder p2N for all f 2Cp(Rd), if the order conditions bTA
()= 1 =!
are satised for all rooted trees of order # p.
5. Generating order conditions with Maple
The recursive constructions underlying the proof of Butcher's theo-rem can easily be realized using modern computer algebra systems like Maple. We assume that the reader is familiar with this particular package.
The recursive data-structure of a rooted tree itself can be realized by self-reference as a list of rooted trees,
> `type/tree`:=list(tree): f:= ]:
Conveniently, we have given a name to the root = ]. Because
Maple's lists areordered, we have to sort the trees recursively before testing for equality:
TreeSort:=proc(beta::tree) sort(map(TreeSort,beta)) end: # TreeSort
Now, here is an example for the tree representing the elementary dif-ferentialf 00( f 00( ff 0( f))f) =f 00( ff 00( f 0( f)f)). The corresponding
tree is obtained by erasing all the derivativesf
(n)from the expression
and replacing parentheses (:::) by :::]:
> beta 1]:= f, f]],f]: beta 2]:= f, f],f]]:
> evalb(TreeSort(beta 1])=TreeSort(beta 2])) true
The denition (2)of the order # is simply expressed by the recursive
procedure
TreeOrder:=proc(beta::tree) option remember
1+`+`(op(map(TreeOrder,beta))) end: # TreeOrder
As an example, we take the tree that represents the elementary dif-ferentialf 00( f 000( f 0( f)f 0( f)f)f): > beta:= f], f],f],f]: TreeOrder(beta) 8
TreeFactorial:=proc(beta::tree) option remember
TreeOrder(beta)*`*`(op(map(TreeFactorial,beta))) end: # TreeFactorial
The above used tree of order 8 gives: > TreeFactorial(beta)
192
For the sake of completeness we also express the recursive deni-tion (3)of using Maple:
TreeAlpha:=proc(beta::tree) option remember local n n:=nops(beta) nops(combinat permute](beta,n))/n!* `*`(op(map(TreeAlpha,beta))) end: # TreeAlpha
An application to the above example yields:
> TreeAlpha(beta)
1 2
We are now ready to map the recursive denition (4)of A
() and of
the order conditionbTA ()= 1 =! to Maple: TreeA:=proc(beta::tree,no::posint) option remember local vars,val,son vars:= i,j,k,l,m,p,q,r,u,v,w] val:=1
for son in beta do
val:= val*Sum(a vars no],vars no+1]]* TreeA(son,no+1),vars no+1]=1..s) od val end: # TreeA TreeOrderCondition:=proc(beta::tree) Sum(b i]*TreeA(beta,1),i=1..s)= 1/TreeFactorial(beta) end: # TreeOrderCondition
The coordinate index for the vector A
() can be chosen from the
sec-ond argument toTreeA. The order condition belonging to the above
example is obtained as follows:
> TreeOrderCondition(beta) s X i=1 bi 0 @ s X j=1 aij s X k=1 ajk( s X l=1 akl) ! 2 (Xs k=1 ajk) 1 A ( s X j=1 aij) = 1 192 Even the typesetting of this formula was done completely automati-cally, using Maple's ability to generate TEX-sources.To generate all the order conditions for a given order
p, we need a
device that constructs the set of all trees with # p. There are,
in principle, two dierent recursive approaches:
{
root-oriented: generate all trees of order # =pby rst, listingall integer partitions p;1 = p 1 +
:::+pn, n = 1:::p;1,
and next, setting = 1
:::n] for all trees 1 :::n of order # 1 = p 1
<p:::# =pn < p. These trees have already been
generated by the recursion.
{
leaf-oriented: Add a leaf to each node of the trees ^ = 1:::n]
of order #^ =p;1, increasing thereby the order exactly by one.
This can be done recursively by adding a leaf to every node of the subtrees
1
:::n.
The root-oriented approach was chosen by Sofroniou 9] in his Math-ematica package Butcher.m. It requires an e cient integer partition
package and the handling of cartesian products. The leaf-oriented approach is as least as e cient as the other one, but much easier to code: Trees:=proc(order::posint) option remember local Replace,AddLeaf,all,trees Replace:=proc(new::tree,old::posint,beta::tree) sort(subsop(old=new,beta)) # order-independent... end: # Replace AddLeaf:=proc(beta::tree) option remember local val,child,new val:=fsort( ],op(beta)])g
for child from 1 to nops(beta) do new:=AddLeaf(beta child])
val:=val union map(Replace,new,child,beta) od
val end: # AddLeaf trees:=f ]g all:= trees] to order-1 do trees:=`union`(op(map(AddLeaf,trees))) all:= op(all),trees] od: all end: # Trees
Given an order p this procedure generates a list of the sets of trees
for each order q p, e.g., > Trees(4)
f]gf]]gf]]]]]]gf]]]]]]]]]]]]]]]]g]
For instance, the number of trees for each order p 10 is given by
the entries of the following list:
> map(nops,Trees(10))
112492048115286719]
The number of order conditions for the order p = 10 can thus be
obtained by:
> `+`(op(%))
1205
Finally, for concrete calculations one has to specify the number s
of stages. The following procedure then generates the specic set of equations forexplicitRunge-Kutta methods:
OrderConditions:=proc(order::posint,stages::posint) option remember local eqs,vars,auto,explicit explicit:=seq(seq(a i,j]=0,j=i..stages), i=1..stages) vars:=eval(seq(b i],i=1..stages), seq(seq(a i,j],j=1..stages),i=1..stages), seq(c i],i=1..stages),explicit) minus f0g auto:=eval(seq(sum(a i,j],j=1..stages)=c i], i=1..stages),explicit) eqs:=value(Eval(map(TreeOrderCondition, `union`(op(Trees(order)))),s=stages)) eqs:=eval(eval(eqs,explicit),auto) eqs,auto,vars end: # OrderConditions
This way, we can automatically generate and typeset the order con-ditions for the classical explicit 4-stage Runge-Kutta methods of or-der 4: > OrderConditions(4,4): % 1] fb 1+ b 2+ b 3+ b 4 = 1 b 4 a 43 a 32 c 2= 124 b 2 c 2+ b 3 c 3+ b 4 c 4 = 12 b 2 c 2 3+ b 3 c 3 3+ b 4 c 4 3 = 14 b 3 a 32 c 2+ b 4( a 42 c 2+ a 43 c 3) = 16 b 3 c 3 a 32 c 2+ b 4 c 4( a 42 c 2+ a 43 c 3) = 18 b 3 a 32 c 2 2+ b 4( a 42 c 2 2+ a 43 c 3 2 ) = 112 b 2 c 2 2+ b 3 c 3 2+ b 4 c 4 2 = 13g
Remark 1.Even in the more recent literature one can nd examples like 4], where order conditions for Runge-Kutta methods are gen-erated by using a computer algebra system to calculate the Taylor expansions of the ow and the discrete ow directly. This approach is typically bound to scalar non-autonomous equations, i.e., d = 1.
Besides being ine cient for higher orders, it is well-known 2] that for p 5 additional order conditions for general systems make an
appearance, which do not show up in the scalar case.
6. Examples of usage
The following simple procedure tempts to solve the order conditions for a given order p and stage number s by using brute force, i.e.,
Maple'ssolve-command. To simplify the task, the user is allowed to
supply a set pre of a priori chosen additional equations and
assign-ments that he thinks to be helpful.
RungeKuttaMethod:=proc(p::posint,s::posint,pre::set) # explicit methods only
local conds,auto,vars,eqs,sols,sol,val conds,auto,vars:=OrderConditions(p,s) eqs:=conds union auto union pre
sols:=solve(eqs,vars) val:=NULL
val:=val, a=matrix( seq( seq(eval(a i,j],sol),j=1..i-1), seq(0,j=i..s)],i=1..s)]), b=vector( seq(eval(b i],sol),i=1..s)]), c=vector( seq(eval(c i],sol),i=1..s)])] od: val end: # RungeKuttaMethod
This way we can simply generate the general explicit 2-stage Runge-Kutta method of orderp= 2:
> RungeKuttaMethod(2,2,fb 2]=thetag) " a= " 0 0 1 21 0 # b= ;+ 1]c= 0 1 21 #
The next example is more demanding. In his book 3, p. 199] Butcher describes an algorithm for the construction of explicit 6-stage meth-ods of order p= 5. The choicesc
6= 1 and b
2 = 0 together with the
free parametersc 2 c 3 c 4 c 5 and a
43 yield a unique method. Butcher
provides a two-parameter example by choosingc 2= uc 3= 1 =4c 4= 1=2c 5= 3 =4a 43=
v. By just passing this additional information to
Maple'ssolve-command we obtain the following solution
> pre:=fc 2]=u,c 3]=1/4,c 4]=1/2,c 5]=3/4,c 6]=1, > b 2]=0,a 4,3]=vg: > RungeKuttaMethod(5,6,pre): % 1] > %% 2],%% 3] a= 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 0 0 0 0 0 0 u 0 0 0 0 0 1 32 ;1 + 8u u 1 32 1u 0 0 0 0 ; 1 8 ;2v+ 1 + 8vu;4u u ; 1 8 2v;1 u v 0 0 0 3 16 ;v+ 1;3u+ 4vu u 3 16 v;1 u 3 4; 3 4v 9 16 0 0 ; 1 14 ;6v+ 7 + 24vu;22u u ; 1 14 6v;7 u 12 7 v ;12 7 87 0 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 b= 7 900 16 45 2 15 16 45 7 90 c= 0u 1 4 1 2 3 41
This result shows that the coe cients a 51 and
a
52 of Butcher's
so-lution 3, p. 199] are in error, a fact that was already observed by Sofroniou 9] using the Mathematica packageButcher.m.
References
1. Butcher, J. C. (1963): Coecients for the study of Runge-Kutta integration processes, J. Austral. Math. Soc.3, 185{201.
2. Butcher, J. C. (1963): On the integration processes of A. Huta, J. Austral. Math. Soc.3, 202{206.
3. Butcher, J. C. (1987):The Numerical Analysis of Ordinary Dierential Equa-tions. Runge-Kutta Methods and General Linear Methods, John Wiley & Sons, Chichester.
4. Gander, W. and Gruntz, D. (1999): Derivation of numerical methods using computer algebra, SIAM Review41, 577{593.
5. Hairer, E. (1999):Numerical Geometric Integration, Lecture notes of a course given in 1998/99,http://www.unige.ch/math/folks/hairer/polycop.html. 6. Hairer, E., Nrsett, S. P., and Wanner, G. (1993): Solving Ordinary Dier-ential Equations I. Nonsti Problems, 2nd Edition, Springer-Verlag, Berlin, Heidelberg, New York.
7. Hairer, E. and Wanner, G. (1974): On the Butcher group and general multi-value methods, Computing 13, 1{15.
8. Jenks, R. J., (1976): Problem # 11: Generation of Runge-Kutta equations, SIGSAM Bulletin 10, 6.
9. Sofroniou, M. (1994): Symbolic derivation of Runge-Kutta methods, J. Symb. Comp.18, 265{296.
