BIT 0006-3835/98/3901-0116 $15.001999, Vol. 39, No. 1, pp. 116–142 c© Swets & Zeitlinger
RUNGE–KUTTA METHODS ADAPTED TOMANIFOLDS AND BASED ON RIGID FRAMES ∗
B. OWREN and A. MARTHINSEN †
Department of Mathematical Sciences, NTNU, N-7034 Trondheim, Norway
email: [emailprotected], [emailprotected]
Abstract.
We consider numerical integration methods for differentiable manifolds as proposedby Crouch and Grossman. The differential system is phrased by means of a systemof frame vector fields E1, . . . , En on the manifold. The numerical approximation isobtained by composing flows of certain vector fields in the linear span of E1, . . . , En
that are tangent to the differential system at various points. The methods reduce totraditional Runge–Kutta methods if the frame vector fields are chosen as the standardbasis of euclidean R
n. A complete theory for the order conditions involving orderedrooted trees is developed. Examples of explicit and diagonal implicit methods arepresented, along with some numerical results.
AMS subject classification: 65L06, 34A50.
Key words: Geometric integration, numerical integration of ordinary differentialequations on manifolds, Runge–Kutta methods.
1 Introduction.
For many years openings like the following:“Consider the initial value problem
y = f(t, y), y(0) = y0, f : R × Rn → R
n(1.1)
. . . ” have been fairly common in papers within the field of numerical ODEs.For readers who normally solve differential equations in the vector space R
n,(1.1) is of course an excellent point of departure. Those whose problems arephrased on a differentiable manifold M may still be happy with the formulation(1.1), since a manifold by definition has a local representation in R
n for some n.Certainly, there may be difficulties when the numerical method needs to switchfrom one coordinate chart to another, but experience will perhaps reveal thatthis is not a serious obstacle in most practical problems. It is possibly a moreserious difficulty that for many manifolds, the local coordinate charts may becomplicated and that the derivation of a representation of the form (1.1) may callfor lengthy calculations by hand or by a computer algebra system. In many such
∗Received January 1997. Revised July 1998. Communicated by Stig Skelboe.†This work was sponsored in part by The Norwegian Research Council under contract no.
111038/410, through the SYNODE project (http://www.math.ntnu.no/num/synode/).
INTEGRATION METHODS BASED ON RIGID FRAMES 117
cases it is preferable to embed the manifold in RN for some N > n. This may
lead to simpler expressions, and the representation can be made global. But adifficulty that arises here is that the computed solution may drift off the manifoldunless information about it is encoded in the algorithm. Recently, some authorslike Crouch and Grossman [3], and Munthe-Kaas with co-authors [10, 11, 12]have adopted a different view on differential equations on manifolds, taking onestep up the ladder of abstraction. Rather than introducing coordinates alreadyin the definition of the differential equation, like in (1.1), the solution to the ODEis considered as a smooth curve on the manifold and the differential equationis interpreted as a vector field for which this smooth curve is an integral curve.The integration method is then formulated in this abstract setting. It shouldbe emphasized that our motivation for this approach is not to develop a moregeneral type of method, but rather the opposite; to tailor integration methods byallowing the user to take into consideration the particular manifold on which thesolution is known to evolve. The method is designed such that the numericalapproximation may be forced to evolve on the manifold by means of flows ofcertain vector fields. The methods proposed by Munthe-Kaas are phrased in themost general case for hom*ogeneous manifolds; see [12], and he also shows that themethods by Crouch and Grossman [3] can be interpreted in the setting of suchmanifolds. Perhaps the most important examples of hom*ogeneous manifolds arethe Lie groups. The configuration of a suspended rigid body can be thought ofas an element of the Lie group SO(3), the group of rotations in three dimensions.Hence, many applications can be found for instance within the field of structuraldynamics.
An important class of differential equations on manifolds are the differential-algebraic problems. One may for instance have a set of algebraic equationsthat define the manifold as an embedding in a higher dimensional vector space.In such a case, it may appear natural to look for a system of vector fields ininvolution that spans the tangent bundle of this manifold. A numerical methodfor such a problem may then typically use as primitives flows of these vectorfields.
In this paper we shall study the approach proposed by Crouch and Grossmanin [3] and applied to a problem in mechanics in [4]. They develop methodsfor manifolds both as generalizations of multistep methods and Runge–Kuttamethods. Here we will focus on the latter type, and give a complete descriptionof the conditions that must be imposed by the method coefficients to ensure thatthe order of consistency is equal to a given integer. We begin in Section 2 byintroducing some necessary notation and present the format of the methods. InSection 3, we introduce the concept of ordered rooted trees and we consider theirconnection to the Lie series of the exact solution. By developing the numericalsolution in a similar series, we obtain in Section 4 a characterization of the orderconditions. It turns out that the conditions are not all independent, and wepresent a way of reducing them to a minimal set in Section 5. In Section 6, wepresent a few examples of Crouch–Grossman methods, to our knowledge the firstexplicit one of order 4. Finally, we show how the methods work for Lie groups,
118 B. OWREN AND A. MARTHINSEN
and pay special attention to the case of matrix groups. We give numericalexamples to demonstrate that the order theory is correct for the spinning topproblem.
2 Background and notation.
We introduce only the concepts from differential geometry that are strictlynecessary for the discussion that follows. The reader may find that the text byWarner [18] is quite compatible with the notation used here.
Let M be a differentiable manifold. A tangent vector v at the point p ∈ M is alinear derivation, acting on functions on M , such that for any f, g ∈ C∞(M, R)and λ ∈ R
v[f + λg] = v[f ] + λv[g]v[f · g] = f(p)v[g] + g(p)v[f ] (Leibniz’ rule).
The set of all tangent vectors to M at p is called the tangent space to M atp, and is denoted TMp. It can be imposed the structure of a vector space inthe natural way, in fact, with a choice of local coordinates, say (x1, . . . , xn), atangent vector v can be expressed as
v = v1∂
∂x1+ v2
∂
∂x2+ · · · + vn
∂
∂xn.
Thus, a tangent vector v can be associated with a direction (v1, . . . , vn) ∈ Rn
and we can think of its action on f at p as the rate of change in f from p alongthis direction.
The tangent bundle of M is defined as
TM =⋃
p∈M
TMp.
A differentiable structure on TM is induced in a natural way from that of M .A vector field F on M is a section of TM , i.e.
F : M → TM
p → v(p), v(p) ∈ TMp.
If φ : M → N is a C∞ map between the two manifolds M and N , we define thedifferential of φ at the point p to be the linear map dφ : TMp → TNφ(p)
dφ(v)[g] = v[g φ], v ∈ TMp and g ∈ C∞(N, R).
A smooth curve on M is a C∞ map σ : (a, b) → M where (a, b) is an interval.The tangent vector of σ at t = t0 is an element, σ(t0) ∈ TMσ(t0), given as
σ(t0) = dσ
(d
dt
∣∣∣∣t0
).
INTEGRATION METHODS BASED ON RIGID FRAMES 119
We may interpret this in coordinates (x1, . . . , xn) as
σ(t0) =n∑
i=1
d(xi σ)dt
∣∣∣∣t0
∂
∂xi.
A smooth curve σ(t) is an integral curve of a vector field F if
σ(t) = F (σ(t)), t ∈ (a, b).
If F is a smooth vector field on M , then there exists, for each p ∈ M , an interval(a(p), b(p)) and a smooth curve, σp : (a(p), b(p)) → M with σ(0) = p such thatσp(t) is an integral curve of F . Now, fixing t, we can obtain a diffeomorphismof M by setting
φF,t : Dt → M, where Dt = p ∈ M : t ∈ (a(p), b(p))p → σp(t).
This map is called the flow of the vector field F .Let E1, . . . , En be smooth vector fields on M . In the obvious way, they gen-
erate a vector space, V over R, and a module, Ω over the ring C∞(M, R). Inwhat follows, we shall be interested in approximating flows of vector fields inΩ with flows of vector fields in V . Notice that for any F ∈ Ω and p ∈ M wecan associate a vector field Fp ∈ V with frozen coefficients relative to the framevector fields E1, . . . , En where
Fp =n∑
i=1
fi(p)Ei.
For any vector field F in either Ω or V we shall use the notation
φF,t(p) =: etF p = exp(tF )p
for its flow. By a system of differential equations relative to the vector fieldsE1, . . . , En, we shall mean an equation of the form
y = F (y) =n∑
i=1
fi(y)Ei, y(0) = y0, F ∈ Ω,(2.1)
whose solution is an integral curve y : [0, T ] → M of F . It is well known that fora function ψ ∈ C∞(M, R) and h ∈ R, the composition ψ ehF p has a Lie series
ψ(p) + hF [ψ](p) +12h2F 2[ψ](p) + · · · + 1
r!hrF r[ψ](p) + · · ·(2.2)
where we use the convention F r+1[ψ] = F [F r[ψ]] , r = 1, 2, . . . . Observe thatthe operator F r, r ≥ 2, will not normally be a vector field, actually it can beregarded as a map p → vr(p) with vr(p) ∈ TM r
p , the vector space of rth order
120 B. OWREN AND A. MARTHINSEN
tangent vectors at p. With a choice x1, . . . , xn of local coordinates, a basis forTM r
p is∂α
∂xα
∣∣∣∣p
: 1 ≤ |α| ≤ r
, where α = (α1, . . . , αn) and |α| =
n∑i=1
αi.
Clearly one has TMp = TM1p ⊂ TM2
p ⊂ · · · . The union⋃
p TM rp is called the
rth tangent bundle of M .Crouch and Grossman propose to approximate the solution of (2.1) by a gen-
eralization of explicit Runge–Kutta methods, under the assumption that flows ofvector fields in V can be computed. In this paper, we also allow for the methodsto be implicit, hence we propose a procedure as follows
p = yk
Yr = ehar,sF sp · · · ehar1F 1
p p r = 1, . . . , s
F rp = FYr =
n∑i=1
fi(Yr)Ei r = 1, . . . , s
yk+1 = ehbsF sp ehbs−1F s−1
p · · · ehb1F 1p p.
(2.3)
The scheme is well defined for sufficiently small values of the step size h. Asit is implicit, both in the approximations Yr ∈ M and the vector fields F r
p ∈ V ,one needs a procedure for the solution of such equations on manifolds. As of yet,it is not clear to us how this should be accomplished for general differentiablemanifolds. However, in a recent paper by Owren and Welfert [14] a version ofthe Newton iteration is proposed in the case that M is a Lie group. We willsubsequently refer to methods of the format (2.3) as Crouch-Grossman methods.
Note that the notation for the coefficients of the method (arj) and (br) isconsistent with the usual Runge–Kutta methods, and hence they can be rep-resented in a Butcher tableau. We will subsequently refer to the coefficientscr, r = 1, . . . , s, where cr =
∑sj=1 arj .
Also, as pointed out in [3], the method reduces to a traditional Runge–Kuttamethod if the manifold is the vector space R
n and the vector fields Ei are chosenas the standard euclidean basis. To see this, let (x1, . . . , xn) be coordinates on R
n
and let the frame vector fields be defined as Ei : y → ∂/∂xi|y , 1 ≤ i ≤ n. For q ∈R
n we let Fq : y →∑
fi(q) ∂/∂xi|y. We then find that exp(aFq) p = p + af(q),where a ∈ R and f(q) = (f1(q), . . . , fn(q))T . By applying this repeatedly to theformulas in (2.3) we recover the usual Runge–Kutta method.
The method defined by (2.3) is a one-step method, hence we can write yk+1 =Ψ(h, yk) = Ψ(h, p). The order of such a method can be defined as the greatestinteger q for which
dr
dhr
(ψ ehF p
)∣∣∣∣h=0
=dr
dhr(ψ Ψ(h, p))
∣∣∣∣h=0
, r = 0, . . . , q, ∀ ψ ∈ C∞(M, R),
(2.4)
INTEGRATION METHODS BASED ON RIGID FRAMES 121
and note again that this definition of order is consistent with the usual definitionwhen the frame vector fields are chosen as the standard euclidean basis.
3 Trees and their relation to the Lie series of the exact solution.
We begin by quoting a definition from [6] for labelled rooted trees.Definition 3.1. Let A be an ordered chain of indices A = r < i1 < i2 < · · ·
and denote by Aq the subset consisting of the first q indices. A (rooted) labelledtree of order q is a mapping
t : Aq\r → Aq such that t(z) < z ∀ z ∈ Aq\r.
The set of labelled trees of order q is denoted LT q.While this definition is the same as in the classical case, we need to modify
the concept of equivalence classes of labelled trees as defined in [6].Definition 3.2. Two labelled trees t and u are equivalent if they have the
same order, say q, and if there exists a permutation σ : Aq → Aq such that
1. σ(r) = r;
2. tσ = σu on Aq\r;
3. for all z1, z2 ∈ Aq\r such that t(z1) = t(z2),
σ(z1) < σ(z2) ⇔ z1 < z2.
The equivalence classes induced by this relation are called ordered (rooted) trees(T q
O). The number of labelled rooted trees belonging to an equivalence class t ∈ T qO
is denoted α(t). The set of all rooted trees is TO = T 1O ∪ T 2
O ∪ · · · . The functionρ : TO → Z
+ is defined as ρ(t) = q if t ∈ T qO.
Definition 3.3. The set of special ordered (rooted) trees SqO with q nodes is
defined as follows.
1. SqO ⊆ T q
O.
2. A tree s ∈ LT q belongs to an equivalence class of s ∈ SqO if and only if the
set s−1(z) contains at most one element for all z ∈ Aq\r. In otherwords, Sq
O is the set of all ordered trees, with exactly q nodes, that have noramifications except at the root.
There is a more convenient notation for ordered rooted trees. We denote thetrivial one-node tree τ and write any ordered tree with two or more nodes on theform t = [t1, . . . , tµ]. Such a tree is obtained by adjoining each root of the treest1, . . . , tµ to a new node that becomes the root of t. Of course, the ordering ofthe trees in the bracket is now significant, the chosen convention being that forall labelled trees which correspond to the ordered tree t = [t1, . . . , tµ] we havei < j ⇒ ri > rj , where rk is the index of the root of tk.
122 B. OWREN AND A. MARTHINSEN
It is well known (see e.g. [2]) that the number of ordered trees with exactly rnodes is equal to the Catalan number
Nr =1r
(2r − 2r − 1
).
As for the number of labelled trees, α(t), belonging to an equivalence class t ∈ T ρO
we have the following recursion formula.Proposition 3.1. α(τ) = 1. Let t ∈ T ρ
O be of the form t = [t1, t2, . . . , tµ]where each subtree ti has ρi nodes. Then
α(t) =µ∏
=1
(∑i=1 ρi − 1ρ − 1
)α(t).
Proof. α(τ) = 1 is obvious. We must consider the number of ways todistribute ρ = 1 +
∑µi=1 ρi indices that are eligible according to Definition 3.2.
Clearly, the two smallest indices must be assigned to the roots of t and tµrespectively. There are now
(ρ−2
ρµ−1
)ways to choose the remaining ρµ−1 nodes of
tµ and these can be distributed internally on tµ in α(tµ) ways. If µ > 1 we mustproceed by assigning to the root of tµ−1 the smallest of the remaining indices.There are then ρµ−1 indices to choose from the remaining ρ − 2 − ρµ indices.This argument is used iteratively along with the relation ρ = 1 +
∑µi=1 ρi to
yield the required result.We now proceed to establish the correspondence between the set of ordered
trees and a counterpart to the elementary differentials known from the classicaltheory of Runge–Kutta methods. These constitute the terms obtained from thepowers of the vector field F in (2.2). We compute the first few powers of F usingonly Leibniz’ rule and linearity of the vector fields Ei.
F 1 = F =∑i1
fi1Ei1 ,
F 2 =∑i1,i2
(fi2fi1Ei2Ei1 + fi2Ei2 [fi1 ]Ei1 ) ,
F 3 =∑
i1,i2,i3
( (1)
fi3fi2fi1Ei3Ei2Ei1 +(2)
fi3fi2Ei3 [fi1 ]Ei2Ei1 +(3)
fi3Ei3 [fi2 ]fi1Ei2Ei1
+(4)
fi3fi2Ei2 [fi1 ]Ei3Ei1 +(5)
fi3Ei3 [fi2 ]Ei2 [fi1 ]Ei1 +(6)
fi3fi2Ei3Ei2 [fi1 ]Ei1
).
It is natural to identify the expression for F 1 with the 2-node tree in Figure 3.1.Similarly, we identify the two terms of F 2 with the two 3-node trees. In F 3 thereare six terms, each corresponding to a labelled rooted tree with 4 nodes in thefigure below. Terms 2 and 4 correspond to the labelled trees depicted in theframe. They only differ in the sense that the indices i2 and i3 are interchanged.But term 3 is in general different since Ei1Ei2 = Ei2Ei1 whenever i1 = i2. Thus,
INTEGRATION METHODS BASED ON RIGID FRAMES 123
when identifying a labelled tree with a term in F q one must make sure that theorder in the composition of frame vector fields Eij Eik
· · · is retained, above weuse the convention that the labels should be decreasing from left to right.
i1
i2i3
r
i1
i2
i3
r
i1
i2
r
i1i2
r
i1
rr
i2i1
i3
r
i3
i1
i2
r
i1 i1
i2 i3
i2i3
r r
Figure 3.1: Labelled rooted trees with q ≤ 4.
We see that the result of applying F (differentiation) to an elementary differ-ential is that a new branch is added at each node, introducing a new greaterlabel. In this way we find that there is a one-to-one correspondence betweenthe set of labelled trees with ρ(t) nodes and the set of elementary differentialsconstituting the differential operator F ρ(t)−1. To see this correspondence whendrawing the labelled trees in a diagram, it may be helpful to always add the newbranch to the left in each ramification (although no such distinction is made inthe definition of a labelled tree).
Definition 3.4. We define the elementary differentials, relative to the framevector fields E1, . . . , En in the following way.
F(τ) = I,
F([t1, . . . , tµ]) =∑
i1,...,iµ
F(t1)[fi1 ] · F(t2)[fi2 ] · · ·F(tµ)[fiµ ] Ei1Ei2 · · ·Eiµ ,
where I denotes the identity operator.We find that
ψ ehF p =∑t∈TO
hρ(t)−1
(ρ(t) − 1)!α(t) F(t)[ψ](p),(3.1)
which is similar to the classical theory.Notice that due to the structure of the elementary differentials, we have that
for any t ∈ TO of the form t = [t1, . . . , tµ],
F(t)[ψ](p) =∑
i1,...,iµ
F(t1)[fi1 ](p) · · ·F(tµ)[fiµ ](p)Ei1 · · ·Eiµ [ψ](p),
124 B. OWREN AND A. MARTHINSEN
and for this reason we define the frozen elementary differentials Fp(t) by
Fp(τ) = I,
Fp([t1, . . . , tµ]) =∑
i1,...,iµ
F(t1)[fi1 ](p) · F(t2)[fi2 ](p) · · ·F(tµ)[fiµ ](p)Ei1Ei2 · · ·Eiµ .
In particular,
Fp([τ ]) =∑i1
F(τ)[fi1 ](p)Ei1 =∑i1
fi1(p)Ei1 .
Clearly, from (3.1) we must have
dq
dhq
(ψ ehF p
)∣∣∣∣h=0
=∑
t∈T q+1O
α(t)Fp(t)[ψ](p), q = 0, 1, . . . .(3.2)
Remark 3.1. In (3.2) the qth derivative of the exact solution is related toordered trees with q + 1 nodes, contrary to the classical theory where the qthderivative corresponds to trees with only q nodes. The restriction of our approachto the euclidean case amounts to choosing coordinate functions φi : y → yi, andframe vector fields Ei : y → ∂/∂yi. By applying the elementary differential(operator) to ψi, we find that all Fp(t) vanish when t is of the form [t1, . . . , tµ]with µ ≥ 2. We are left simply with the Fp(t) corresponding to trees of the formt = [u], hence by identifying each such t with u we obtain a representation ofF ρ(t) by means of trees with ρ(t) nodes.
4 The derivatives of the approximate solution.
We need to find expressions for the rth derivative of ψ Ψ(h, p) as it occursin (2.4). This can be done in different ways, one of them is proposed in [3].Another technique is based on the use of the Baker–Campbell–Hausdorff formula[17], i.e., the composition of flows occurring in (2.3) can be replaced by the flowof one single vector field, say G(h). Then exp(G(h)) can be differentiated withrespect to h through the introduction of the differential of the exponential map.Although this approach is straightforward, it seems to generate extremely longand complicated expressions involving iterated Lie brackets of vector fields, andwe have not been able to discover any tractable way of dealing with them.Instead, we have chosen to base the expansion of the numerical approximationsolely on ordered rooted trees arising from the use of Lie series, not very unlikethe expansion of the exact solution as given in the previous section. This resultsin a nice and clean theory, and it is consistent with the theory for classicalRunge–Kutta methods. A disadvantage is that one obtains several redundantconditions, such that only a subset of them needs to be considered. We presentin the next section some results that can be used to select those conditions thatcan be discarded.
INTEGRATION METHODS BASED ON RIGID FRAMES 125
Let j = (j1, . . . , jr) be a multi-index of integer elements with 1 ≤ jr ≤ s. Wedefine
j! := q1!q2! · · · qs!, where qi is the number of occurrences of i in (j1, . . . , jr).(4.1)
For example, with j = (1, 1, 2, 3, 3) we have j! = 2! · 1! · 2! = 4.We have the following result.Lemma 4.1. For any real s-plet
v = (v1, . . . , vs), ψ ∈ C∞(M, R),
and smooth vector fields F 1p , . . . , F s
p ∈ V depending on h, the composition of theirflows is defined for sufficiently small values of h and can be written as follows:
ψ(ehvsF s
p · · · ehv1F 1p p
)= Wv(h)[ψ](p)
where Wv(h) is the formal operator
Wv(h) = I +∞∑
r=1
hrs∑
j1=1
s∑j2=j1
· · ·s∑
jr=jr−1
1j!
vj1 · · · vjr F j1p · · ·F jr
p .(4.2)
Furthermore, there is a Faa di Bruno type of formula for the qth derivative ofthe operator Wv(h),
W (q)v (h)
∣∣∣h=0
=∑
u∈Sq+1O
q!µ∏
ν=1
(δν − 1)!
×∑
1≤j1≤···≤jµ≤s
1j!
vj1 · · · vjµ
(F j1
p
)(δ1−1) · · ·(F jµ
p
)(δµ−1)∣∣∣h=0
,
for q = 1, 2, . . . . Here, u ∈ Sq+1O is of the form [u1, . . . , uµ] and δi = ρ(ui)
(compare [6, Lemma 2.8 (p. 150)]).Proof. We compute
ψ(ehvsF sp · · · ehv1F 1
p p)
by using formula (2.2) repeatedly, beginning with F = vsFsp and p replaced by
evs−1F s−1p · · · ehv1F 1
p p.This leads to the expression
∞∑r1=0
· · ·∞∑
rs=0
hr
r1! · · · rs!vr11 · · · vrs
s
(F 1
p
)r1 · · ·(F s
p
)rs [ψ](p),
where r = r1 + · · · + rs, which is equal to Wv(h)[ψ](p).
126 B. OWREN AND A. MARTHINSEN
For sufficiently small h we can do a term-wise differentiation of Wv(h), andwe recall the formula
(hr g(h))(q)∣∣∣h=0
=
0 if q < r
q!(q − r)!
g(h)(q−r)
∣∣∣∣h=0
if q ≥ r
for any smooth function g(h). We shall apply this formula to (4.2). Repeateduse of Leibniz’ rule yields(
F j1p · · ·F jr
p
)(q−r) =∑
δ1+···+δr=q
(q − r
δ1 − 1, . . . , δr − 1
)(F j1
p
)(δ1−1) · · ·(F jr
p
)(δr−1).
Any ordered r-plet (δ1, . . . , δr) identifies uniquely a term in the (q− r)th deriva-tive of ∑
1≤j1≤···≤jr≤s
1j!
(F j1
p · · ·F jrp
)(q−r).
On the other hand, the same r-plet identifies uniquely the tree u ∈ Sq+1O with
subtrees ui such that ρ(ui) = δi. Thus, by summing over all trees in Sq+1O we
obtain the last part of the lemma.We have the following main result for the derivatives of the approximate so-
lution.Theorem 4.2. For a given Crouch–Grossman method (2.3), let Φr : TO → R
be defined recursively for each r = 1, . . . , s as
Φr(τ) = 1, and for t = [t1, . . . , tµ]
Φr(t) =s∑
j1=1
s∑j2=j1
· · ·s∑
jµ=jµ−1
1j!
ar,j1 · · ·ar,jµΦj1(t1) · · ·Φjµ(tµ),(4.3)
let Φ : TO → R be defined as
Φ(τ) = 1, and for t = [t1, . . . , tµ]
Φ(t) =s∑
j1=1
s∑j2=j1
· · ·s∑
jµ=jµ−1
1j!
bj1 · · · bjµΦj1(t1) · · ·Φjµ(tµ),(4.4)
and let γ : TO → Z+ be defined recursively as
γ(τ) = 1, γ(t) =µ∏
r=1
γ(tr)r∑
i=1
ρ(ti) when t = [t1, . . . , tµ].(4.5)
Then
dq
dhqψ Ψ(h, p)
∣∣∣∣h=0
=∑
t∈T q+1O
α(t)γ(t)Φ(t)Fp(t)[ψ](p), q = 1, 2, . . . ,(4.6)
with α(t) as given in Proposition 3.1.
INTEGRATION METHODS BASED ON RIGID FRAMES 127
It now follows directly from (2.4) and (3.2) thatCorollary 4.3. A Crouch–Grossman method (2.3) is of order q if
Φ(t) =1
γ(t)for all t ∈ T 2
O ∪ T 3O ∪ · · · ∪ T q+1
O .(4.7)
Remark 4.1. Note that Corollary 4.3 provides sufficient conditions for aCrouch–Grossman method to have order q. The necessity of these conditionsdepends on which assumptions that are imposed with respect to the linear in-dependence amongst the elementary differentials
Fp(t) : t ∈ T 2
O ∪ · · · ∪ T q+1O
.
Such assumptions were made with respect to the frame vector fields E1, . . . , En
in [3] in order to obtain also necessary order conditions. We shall not make anysuch assumptions here, but rather leave the choice open to accommodate theapplication at hand.
Remark 4.2. The reader who is familiar to the classical Butcher theorymay be puzzled to observe that the conditions for order q are characterized interms of trees of orders 2, . . . , q + 1 rather then 1, . . . , q as in the classical case.Certainly, the coefficients of a Crouch–Grossman method of order q must satisfythe classical order q conditions since the theory is also valid in the case that theframe vector fields coincide with a basis for euclidean R
n, i.e. Ei : y → ∂∂xi
∣∣y. To
establish this connection, let t = [t1, . . . , tµ]no be a non-ordered rooted tree, suchthat ρ(t) = q. It is possible to show that the classical condition correspondingto t will hold if the order conditions arising from the set of ordered rooted trees
[[tσ(1), . . . , tσ(µ)]] ∈ T q+1O : σ ∈ Sµ
are satisfied, where Sµ is the group of permutations on 1, . . . , µ.
Prior to proving Theorem 4.2 we shall prove the following lemma.Lemma 4.4.
α(t)γ(t) = (ρ(t) − 1)! for all t ∈ TO.
Proof. α(τ)γ(τ) = 1 by definition. We proceed by induction and assume thatthe result is true for any t ∈ T 1
O ∪· · ·∪T q−1O with q > 1. Let t = [t1, . . . , tµ] ∈ T q
O
and set ρi := ρ(ti). We get from Proposition 3.1 and (4.5)
α(t)γ(t) =µ∏
=1
((∑i=1 ρi − 1ρ − 1
) ∑i=1
ρi
)α(t1) · · ·α(tµ) · γ(t1) · · · γ(tµ)
=µ∏
=1
(∑
i=1 ρi)!
(ρ − 1)!(∑−1
i=1 ρi)!(ρ − 1)! = (ρ(t) − 1)!
where as usual sums over empty index sets are taken as zero.
128 B. OWREN AND A. MARTHINSEN
Proof of Theorem 4.2. The theorem is proved by induction. We willobtain the auxiliary formula
dq
dhq(ψ Yr)
∣∣∣∣h=0
=∑
t∈T q+1O
α(t)γ(t)Φr(t)Fp(t), 1 ≤ r ≤ s, q = 1, 2, . . . ,
(4.8)
as a part of the induction proof. In fact, the formula (4.6) can be written in thisform if we define as+1,j := bj , 1 ≤ j ≤ s, Ys+1 := Ψ(h, p), and Φs+1(t) := Φ(t).We shall adopt this notation in the rest of the proof. Note also that (4.8)immediately implies
(F r
p
)(q)∣∣∣h=0
=dq
dhq
∣∣∣∣h=0
∑i
(fi Yr)Ei =∑
i
∑t∈T q+1
O
α(t)γ(t)Φr(t)Fp(t)[fi]Ei
=∑
t∈T q+1O
α(t)γ(t)Φr(t)Fp([t]).(4.9)
Let ar = (ar1, . . . , ars), 1 ≤ r ≤ s + 1. From (2.3) and Lemma 4.1 we findthat
ψ Yr = War (h)[ψ], 1 ≤ r ≤ s + 1,
for small h and that
dq
dhq(ψ Yr)
∣∣∣∣h=0
=∑
u∈Sq+1O
q!∏µν=1(δν − 1)!
s∑j1=1
· · ·s∑
jµ=jµ−1
1j!
ar,j1 · · · ar,jµ
×(F j1
p
)(δ1−1) · · ·(F jµ
p
)(δµ−1)∣∣∣h=0
[ψ](p)(4.10)
for q = 1, 2, . . . . First we set q = 1. The only tree in T 2O (and S2
O) is [τ ], so weget
d
dh(ψ Yr)
∣∣∣∣h=0
=s∑
j=1
arj (F jp )(0)
∣∣∣h=0
[ψ](p) =s∑
j=1
arjFp([τ ])[ψ](p)
for 1 ≤ r ≤ s + 1. Moreover Φr([τ ]) =∑s
j=1 ar,j, 1 ≤ r ≤ s + 1, and α([τ ]) =γ([τ ]) = 1, thus formulas (4.8) and (4.6) hold for q = 1.
We proceed by induction: Substitute q by k in (4.8) and (4.9) and assumethat the formulas hold for 1 ≤ k ≤ q − 1 where q > 1 is an arbitrary integer.Note that q + 1 = 1 + δ1 + · · · + δµ in (4.10), thus each δν − 1 ≤ q − 1. We usethe induction hypothesis and substitute (4.9) in (4.10) to obtain
INTEGRATION METHODS BASED ON RIGID FRAMES 129
dq
dhq(ψ Yr)
∣∣∣∣h=0
=∑
u∈Sq+1O
q!∏µν=1(δν − 1)!
s∑j1=1
· · ·s∑
jµ=jµ−1
1j!
ar,j1 · · · ar,jµ
×∑
t1∈Tδ1O
α(t1)γ(t1)Φj1(t1)Fp([t1]) · · ·∑
tµ∈TδµO
α(tµ)γ(tµ)Φjµ(tµ)Fp([tµ])[ψ](p)(4.11)
where 1 ≤ r ≤ s + 1.We now make two observations. First, to each u ∈ Sq+1
O characterized by theordered µ-plet (δ1, . . . , δµ), there is a bijection between T δ1
O × · · · × Tδµ
O andthe subset of T q+1
O consisting of trees with exactly µ ramifications at the root,namely the one that maps (t1, . . . , tµ) to t = [t1, . . . , tµ]. Consequently, we canreplace ∑
u∈Sq+1O
∑t1∈T
δ1O
· · ·∑
tµ∈TδµO
with∑
t∈T q+1O
.
Secondly, all the operators (vector fields) Fp([ti]) in the above expression havefrozen coefficients, hence we find
Fp([t1]) · · ·Fp([tµ])[ψ](p) =∑k1
Fp(t1)[fk1 ]Ek1 · · ·∑kµ
Fp(tµ)[fkµ ]Ekµ [ψ](p)
=∑
k1,...,kµ
Fp(t1)[fk1 ] · · ·Fp(tµ)[fkµ ]Ek1 · · ·Ekµ [ψ](p)
= Fp([t1, . . . , tµ]).(4.12)
We can now employ (4.12) and Lemma 4.4 to (4.11), keeping in mind thatq = ρ(t) − 1 where t = [t1, . . . , tµ] to obtain the stated formulas (4.8) and (4.6).
5 Reduction of order conditions.
It turns out that several of the conditions in Corollary 4.3 depend on eachother, hence only a subset of them needs to be considered. The key observation isthe following: Let Sµ be the symmetric group of order µ and for each i = 1, . . . , µ,let σi = (12 · · · i) ∈ Sµ written in cyclic notation, i.e.
σi : (p1, . . . , pi, pi+1, . . . , pµ) → (p2, . . . , pi, p1, pi+1, . . . , pµ).
Let ti ∈ TO, i = 1, . . . , µ. By a straightforward induction proof, we can establishthe identity
Φ([t1])Φ([t2, . . . , tµ]) =µ∑
i=1
Φ([tσi(1), . . . , tσi(µ)]), µ ≥ 2.(5.1)
130 B. OWREN AND A. MARTHINSEN
It is easily verified (again by induction) that (5.1) is consistent with the condi-tions of Corollary 4.3, i.e.
(γ([t1])γ([t2, . . . , tµ]))−1 =µ∑
i=1
(γ([tσi(1), . . . , tσi(µ)])
)−1.
Actually, Tracogna and Welfert [15] have shown that (5.1) is a special case ofthe more general relation∑
σ=σim···σi1
m≤i1≤···≤im≤µ
Φ([tσ−1(1), . . . , tσ−1(µ)]) = Φ([t1, . . . , tm])Φ([tm+1, . . . , tµ]),
1 ≤ m < µ. For our purpose it is sufficient to consider (5.1), but the interestedreader should consult [15] for the complete solution of the reduction problem.
To make use of the identity (5.1), let t1, . . . , tµ ⊂ TO be any collection oftrees and assume that the coefficients of the Crouch–Grossman method (2.3) arechosen such that
Φ([ti]) =1
γ([ti]), i = 1, . . . , µ,
andΦ([tσ(1), . . . , tσ(µ−1)]) =
1γ([tσ(1), . . . , tσ(µ−1)])
, ∀σ ∈ Sµ.
We now consider the order conditions corresponding to the set of trees
J =t = [tσ(1), . . . , tσ(µ)] : σ ∈ Sµ
.
For each permutation σ ∈ Sµ we obtain a linear equation in the quantitiesΦ(t) : t ∈ J. More precisely, from (5.1) we get
µ∑i=1
Φ([tσiσ(1), . . . , tσiσ(µ)]) =1
γ([tσ(1), . . . , tσ(µ−1)])γ([tσ(µ)]), ∀σ ∈ Sµ.(5.2)
Clearly, of the µ! linear equations above, some will be identical unless all oft1, . . . , tµ are distinct. Specifically, if there are ν distinct trees t1, . . . , tν and thetree ti appears i times there will be a total of
n :=µ!
1! · 2! · · · ν !
distinct linear equations in (5.2). In any case, given the set t1, . . . , tµ, (5.2)defines a linear transformation L : R
n → Rn . It is then possible to partition
the set J into subsets J1 and J2 and express each Φ(t), t ∈ J1, in terms of theelements of the set Φ(t) : t ∈ J2. Then J1 has precisely rank L elements andJ2 has dim(kerL) elements. In conclusion, it is at most necessary to considerNL := dim(kerL) order conditions from the set J . Note that if there are νdistinct trees among t1, . . . , tµ, and the ith of these distinct trees occurs i times,
INTEGRATION METHODS BASED ON RIGID FRAMES 131
NL depends only on 1, . . . , ν , i.e. NL = NL(1, . . . , ν). The precise expressionfor NL can be found in another paper by Tracogna and Welfert [16].
For the single tree generated by µ ≥ 2 copies of the tree t1, i.e. t = [t1, . . . , t1],we clearly have NL(µ) = 0, thus the order condition corresponding to t is satisfiedif Φ([t1]) = 1/γ([t1]).
Next we consider the case with µ− 1 copies of the tree t1 and one occurrenceof a different tree t2. The corresponding set of trees J has µ elements. By anatural choice of basis for the linear transformation above, we obtain the matrixrepresentation L of L:
L =
1 1 1 · · · 1
µ − 1 1 0 · · · 00 µ − 2 2 · · · 0...
. . . . . . . . ....
0 · · · 0 1 µ − 1
.
Clearly rank L = µ − 1, the first row multiplied by µ − 1 equals the sum of theremaining rows. Any choice of µ−1 columns of L results in a linearly independentset. In conclusion: It is sufficient to include one of the µ order conditionscorresponding to the set J , and the choice of order condition is arbitrary.
These observations can be used to find the minimal number of order conditionsthat must be satisfied by methods with q ≤ 5. For methods of higher order onemust use the theory of [15] to determine the minimal number of conditions thatmust be considered.
We illustrate the above discussion by applying the results to the order con-ditions up to order 4, conditions corresponding to trees with q + 1 nodes,q = 1, . . . , 4. Referring to Table 5.1 we recognize the trees t31, t41, t501, t509as being of the form [t, t, . . . , t] with at least two occurrences of t in the bracket.
Furthermore from each of the sets t42, t43, t501, t502, t503, t505, t506, andt507, t508 we need only consider the order condition corresponding to one (ar-bitrary) member.
In conclusion, a sufficient selection of order conditions to be satisfied are thosecorresponding to the following trees.
Order 1.
t21 :s∑
r=1
br = 1.
Order 2.
t32 :s∑
r=1
br cr =12.
Note that the conditions for orders 1 and 2 are exactly the same as forstandard Runge–Kutta methods, thus, by choosing the coefficients of anysuch method with (classical) order ≥ 2 we obtain at least order 2 also inthe Crouch–Grossman case.
132 B. OWREN AND A. MARTHINSEN
Order 3.
t42 :s∑
r=1
br cr
(12br +
s∑j=r+1
bj
)=
16,
t44 :s∑
r=1
s∑j=1
brarjcj =16,
t45 :s∑
r=1
br c2r =
13.
Note that the only extra condition compared to the standard case stemsfrom t42. This condition is nonlinear in the br-coefficients and it can alsobe found in [3] for explicit methods with three stages.
Order 4. A sufficient set is t502, t505, t507, t510, t511, t512, t513, t514. We writeout only those of the conditions that differ from the standard case.
t502 :s∑
r=1
br cr
(16b2r +
s∑j=r+1
bj
(12br +
12bj +
s∑k=j+1
bk
))=
124
,
t505 :12
s∑r=1
br c2r
(12br +
s∑j=r+1
bj
)=
124
,
t507 :s∑
r=1
s∑j=1
brarjcj
(12br +
s∑k=r+1
bk
)=
124
,
t511 :s∑
r=1
s∑j=1
brarjcj
(12arj +
s∑k=j+1
ark
)=
124
,
t512 :s∑
r=1
s∑j=1
brarj
(12arjcj +
s∑k=j+1
arkck
)=
18.
Note that if one adds up the conditions corresponding to t511 and t512 oneobtains the classical order condition corresponding to the (non-ordered)tree [τ, [τ ]]; see also Remark 4.2.
INTEGRATION METHODS BASED ON RIGID FRAMES 133
Table 5.1: Trees of order ≤ 5 and some associated functions.
lbl t γ(t) α(t) Φ(t)
τ 1 1 1
t21 1 1∑s
r=1br
t31 2 1∑s
r1=1
∑sr2=r1
1r!br1br2 = 1
2 (∑s
r=1br)2
t32 2 1∑s
r=1brcr
t41 6 1∑s
r1=1
∑sr2=r1
∑sr3=r2
1r!br1br2br3 = 1
6 (∑s
r=1br)3
t42 6 1∑s
r1=1
∑sr2=r1
1r!br1br2cr1
t43 3 2∑s
r1=1
∑sr2=r1
1r!br1br2cr2
t44 6 1∑s
r=1
∑sj=1brarjcj
t45 6 1∑s
r=1br12c2
r
t501 24 1∑s
r1=1
∑sr2=r1
∑sr3=r2
∑sr4=r3
1r!br1br2br3br4 = 1
24 (∑s
r=1br)4
t502 24 1∑s
r1=1
∑sr2=r1
∑sr3=r2
1r!br1br2br3cr1
t503 12 2∑s
r1=1
∑sr2=r1
∑sr3=r2
1r!br1br2br3cr2
t504 8 3∑s
r1=1
∑sr2=r1
∑sr3=r2
1r!br1br2br3cr3
t505 24 1∑s
r1=1
∑sr2=r1
1r!br1br2
12c2
r1
t506 8 3∑s
r1=1
∑sr2=r1
1r!br1br2
12c2
r2
t507 24 1∑s
r1=1
∑sr2=r1
∑sj=1
1r!br1br2ar1jcj
t508 8 3∑s
r1=1
∑sr2=r1
∑sj=1
1r!br1br2ar2jcj
t509 8 3∑s
r1=1
∑sr2=r1
1r!br1br2cr1cr2 = 1
2 (∑s
r=1brcr)2
t510 24 1∑s
r=1br16c3
r
t511 24 1∑s
r=1
∑sj1=1
∑sj2=j1
1j! brarj1arj2cj1
t512 12 2∑s
r=1
∑sj1=1
∑sj2=j1
1j! brarj1arj2cj2
t513 24 1∑s
r=1
∑sj=1
∑sk=1brarjajkck
t514 24 1∑s
r=1
∑sj=1brarj
12c2
j
134 B. OWREN AND A. MARTHINSEN
6 Examples of methods.
6.1 Explicit methods.
Explicit methods of order three have already been developed by Crouch andGrossman [3]. We give one example method with three stages in the Butchertableau below.
0 034
34
1724
119216
17108
1351 − 2
32417
(6.1)
For methods of order 4 there are a total of 13 order conditions versus 8 inthe standard case. The search for solutions of the order conditions is furthercomplicated by the nonlinearity in the br-coefficients. By careful investigationof the order conditions, we found that it is not possible to obtain explicit methodsof order 4 with only 4 stages. The proof is technical and is not included here.We therefore tried to find methods with 5 stages. To reduce the complexity inthe conditions, we introduced simplifying assumptions, requiring that
s∑j=1
arjcj =12c2r, ∀ r = 2.
Under these assumptions we found certain solutions; an example is:
c2 = 32 , c3 = 1
3κ + 16κ2 + 2
3 ,
c4 = − 13κ − 1
6κ2 + 13 , c5 = 1,
b1 = b5 =12· 1 + κ + κ2
κ + κ2, b2 = 0,
b3 = −16· 1 + 2κ + κ2
2 + κ + κ2, b4 = −1
2· 1κ + κ2
,
a32 = 118 (4 + 3κ + 2κ2), a42 = (1 + κ + κ2)θ − 1
18 (4 + 3κ + 2κ2),
a52 = θ, a43 =−9(1 + κ + κ2)θ + 3 + κ + κ2
4 + 2κ + κ2,
a54 = − κ + κ2
4 + 2κ + κ2, a53 =
−9(1 + κ + κ2)θ + 3 + 2κ + 2κ2
10 + 8κ + 7κ2,
(6.2)
where κ = 21/3 and θ is a root of the quadratic equation
81θ2 − 9(1 + κ + κ2)θ − (25 + 21κ + 17κ2) = 0.
INTEGRATION METHODS BASED ON RIGID FRAMES 135
6.2 SDIRK and ESDIRK methods.
An SDIRK method with three stages can be represented in the followingButcher tableau (see e.g. [6]):
γ γ
c2 a21 γ
c3 a31 a32 γ
b1 b2 b3
1 158 − 3
8 12348 − 125
432 − 25108 1
3175 − 4
34825
(6.3)
It is not clear to us what the desired properties of a Crouch–Grossman type ofSDIRK method should be, for an example, we derived a method of order 3 whichis A-stable (but not L-stable) in the standard case given to the right above.
Recently, Kværnø [8] has derived improved candidates among a type of meth-ods, denoted ESDIRK methods. These seem to have better properties than theSDIRK methods for classical problems. The methods are stiffly accurate anduse an explicit first stage, thereby acquiring the FSAL property, hence they ef-fectively cost one less stage per step. With 4 (effectively 3) stages, the methodsare of the form
0 0
c2 a21 γ
c3 a31 a32 γ
1 b1 b2 b3 γ
b1 b2 b3 γ
0 034 − 1
4 12524
59216 − 25
108 1
1 2275
23 − 24
25 12275
23 − 24
25 1
(6.4)
An A-stable (but not L-stable) example of such a Crouch–Grossman methodof order 3 with rational coefficients is given to the right. It turns out that theextra order condition corresponding to t42 of Table 5.1 can be fulfilled withoutmaking any major restrictions on the parameter γ. By choosing γ ≈ 0.43 onecan obtain a method which is L-stable in the classical sense.
7 Application to Lie groups.
A Lie group G is a differentiable manifold equipped with a group structuresuch that the map
G × G → G,
(y, z) → y · z−1, y, z ∈ G,
is smooth. By using this group operation, we define a diffeomorphism of G bymeans of right multiplication, i.e. for any y ∈ G, let Ry : G → G be the map
Ry(z) = z · y, ∀z ∈ G.
136 B. OWREN AND A. MARTHINSEN
There is an identity element in G, henceforth denoted e. The tangent space ate, g := TGe, equipped with a bilinear, skew-symmetric operation denoted theLie bracket, [·, ·] : g × g → g, is called the Lie algebra of G. This operation canbe defined as follows: For any v ∈ g let Xv be the right invariant vector field
Xv : y → dRy(v).
The set of right invariant vector fields is isomorphic to g. Let Z be the operator
Z = Xv Xw − Xw Xv.
It can then be shown (see e.g. [18]) that Z is a smooth right invariant vectorfield, hence there is an element u ∈ g such that Z = Xu. We define the bracket[v, w] := u.
The map dRy : g → TGy is an isomorphism of the tangent spaces TGe andTGy. Therefore, any smooth vector field on G can be written in the form
F : y → dRy(f(y)), f : G → g.(7.1)
Let v1, . . . , vn be a basis for g. We may now choose the frame vector fields Ei
as right invariant vector fields
Ei : y → dRy(vi), i = 1, . . . , n,
and it follows that f(y) =∑n
i=1 fi(y)vi with fi(y) defined by (2.1).In the case that G is a matrix Lie group (see e.g. [5]) the Crouch–Grossman
methods take a particularly simple form through the use of the matrix exponen-tial. The group elements are now invertible m × m-matrices, and we associateeach element in the Lie algebra with an m × m-matrix v such that the ma-trix exponential ev ∈ G. The Lie bracket is simply the matrix commutator[v, w] = v ·w−w · v and the flow of a vector field Xv ∈ V can be computed withthe formula
etXvy = etv · y,
where the matrix exponential is used on the right hand side. A consequenceof this is that the flows occurring in the Crouch–Grossman methods (2.3) areobtained by computing the matrix exponential of the function f(y) occurring in(7.1) for various values of y.
Examples of matrix Lie groups and their Lie algebras are:
1. the general linear group GL(m), consisting of all invertible m×m matrices.The corresponding Lie algebra is gl(m) the vector space of all m × mmatrices;
2. the group of m × m orthogonal matrices with determinant 1, denotedSO(m). This is a connected subgroup of GL(m). The Lie algebra is so(m),the vector space of skew-symmetric m × m matrices and its dimension ism(m − 1)/2;
INTEGRATION METHODS BASED ON RIGID FRAMES 137
Any vector space V (not necessarily of matrices) has the structure of an abelianLie group where the binary operation is addition, v−1 = (−v) for each v ∈ V ,and the identity element e is the zero vector. The Lie algebra is isomorphic toV itself and the Lie bracket [x, y] = 0 for all x, y ∈ V .
The cartesian product G1 ×G2 between two Lie groups G1 and G2 can againbe made into a Lie group through the multiplication rule
(σ1, σ2) · (τ1, τ2) = (σ1 · τ1, σ2 · τ2), σ1, τ1 ∈ G1, σ2, τ2 ∈ G2.
If g1 and g2 are the Lie algebras of G1 and G2 respectively, then g1 × g2 is theLie algebra of G1 × G2.
g
GlobalLocal B(t)ω
Figure 7.1: A schematic picture of the spinning top in local and global coordinates..
Example. We shall apply our theory to the spinning top in Figure 7.1, see e.g.[1, 9]. The top is modeled as a rigid body, suspended at one point. Its positioncan be represented as a rotation which transforms a point in local coordinates toglobal coordinates. Such a rotation B = B(t) can be represented as an elementof the Lie group SO(3), the group of orthogonal 3×3-matrices with determinantone. The tangent to B(t) is B(t) ∈ TSO(3)B(t) and can therefore be writtenin the form B(t) = ω(t) · B(t), where ω(t) ∈ so(3). We call ω(t) the angularmomentum of the spinning top. From a physical perspective it is natural tothink of ω(t) as a 3-vector ω(t) defined by the identification
ω =
ω1
ω2
ω3
↔
0 −ω3 ω2
ω3 0 −ω1
−ω2 ω1 0
= ω.
The configuration of the spinning top can be described by the pair(B(t), ω(t)
)∈ SO(3) × so(3) =: G.
Thus, G is a cartesian product of two Lie groups. The Lie algebra can be iden-tified with g := so(3)× so(3). We list some properties of this Lie group/algebra:
Product in G: (a, v) · (b, w) = (a · b, v + w).Addition in g: (u, v) + (u, v) = (u + u, v + v).Multiplication by scalar in g: α(u, v) = (αu, αv).Lie bracket in g: [(u, v), (u, v)] = ([u, u], 0).Exponential map from g to G: exp(u, v) = (exp(u), v).
138 B. OWREN AND A. MARTHINSEN
10−5
10−4
10−3
10−2
10−1
100
10−15
10−10
10−5
100
105
stepsize
glob
al e
rror
− CG4
− CG3
− RK4
Figure 7.2: Global error versus stepsize for the schemes CG3, CG4 and RK4 .
The dynamics of the spinning top is discussed in many text books; see forinstance [1]. Denote by g the external forces, e.g. gravitation, let A be theinertia tensor, and C the center of mass of the spinning top. Define
f : G → g,
(r, ω) → (ω, ω′),
where
ω′ = −BA−1BT ω BABT ω + BA−1BT (B C × g).
We obtain the differential equation
d
dt(B(t), ω(t)) = dR(B(t),ω(t)) (f(B(t), ω(t))) .
If we write out the right hand side in terms of matrix products we get
dR(B,ω) (f(B, ω)) = (ω · B, ω′).
The flows of vector fields appearing in (2.3) are now computed using the expo-nential map defined above for suitable values of the function f : G → g.
We have implemented the schemes CG3 (6.1), CG4 (6.2), SDIRK3 (6.3) andESDIRK3 (6.4) in MATLAB and applied them to the spinning top test problem,the main purpose being to verify the order theory presented in this paper. Theinitial data of the simulations were
B0 =
1.0 0.0 0.00.0 cos(φ) sin(φ)0.0 − sin(φ) cos(φ)
and ω0 =
0.0 −1.0 0.01.0 0.0 0.00.0 0.0 0.0
INTEGRATION METHODS BASED ON RIGID FRAMES 139
10−5
10−4
10−3
10−2
10−1
100
10−12
10−10
10−8
10−6
10−4
10−2
100
102
stepsize
glob
al e
rror
− SDIRK3
− ESDIRK3
Figure 7.3: Global error versus stepsize for the schemes SDIRK3 and ESDIRK3 .
with φ = π/16. The inertia tensor, the center of mass and the external forceswere given as
A =18
7 0 00 7 00 0 2
, C =
00√3/2
, and g =
00
−9.81
.
Figure 7.2 shows the global error versus stepsize of CG3 and CG4. For com-parison, we implemented “the Runge–Kutta method”, RK4 [6, p. 138], using theformat (2.3). The lines in the plot are included for reference, they have slopes2, 3 and 4. Note that even though RK4 has order 4 in the classical sense, itexhibits only order 2 in the Crouch–Grossman setting. This is consistent withthe fact that additional order conditions are introduced only from order 3.
The same type of data is plotted in Figure 7.3 for the SDIRK3 and ESDIRK3methods. The straight line has slope 3. The non-linear equations are solved ineach step by means of the Newton iteration on Lie groups as presented in [14],iterating to machine accuracy for each stage. The experiments show that for thespinning top problem, implicit methods are very inefficient, but the tests stillserve their purpose here.
In all the experiments, the global error was computed by first applying theMATLAB library routine ODE45 with tolerance 10−14. The error was measuredusing the sum of the 2-norm of each solution component.
We remark that the above explicit methods applied to the spinning top prob-lem are examples of explicit orthogonal integrators, see [7] for a discussion of or-thogonal integrators. It is known that classical methods cannot both be explicitand preserve orthogonality. However, it is known that for instance the classical
140 B. OWREN AND A. MARTHINSEN
0 0.5 110
−16
10−14
10−12
10−10
10−8
10−6
10−4
integration time
dist
ance
from
the
man
ifold
RK4
CG3
Figure 7.4: ‖B(t)T B(t)−B(0)T B(0)‖2 for CG3 and RK4 implemented in the classicalway.
Gauss methods do preserve orthogonality. This fact may also indicate that theimplicit Crouch–Grossman method is a poor choice for this type of problems. Todemonstrate the preservation of orthogonality in the explicit Crouch–Grossmanmethod, we have plotted in Figure 7.4 the quantity ‖B(t)T B(t) − B(0)T B(0)‖2
as a function of time for the method CG3 and for RK4 implemented in theclassical way. For the approximations obtained with CG3, the drift from themanifold SO(3) is of a magnitude consistent with the machine accuracy.
8 Concluding remarks.
We have presented a type of method which is a generalization of Runge–Kuttamethods. Rather than being able to move only along straight lines in a vectorspace, the primitives of the Crouch–Grossman methods consist of flows belongingto the linear span of certain frame vector fields that can be chosen at liberty bythe user to accommodate any particular application. It is thus assumed that theuser has the means to evaluate or approximate these flows, and the theory doesnot consider the effect of any inaccuracies in the computation of these flows. Itis not clear how the frame vector fields should be chosen in general. For certaintypes of manifolds, there will be a natural choice of frame vector fields, forinstance, for Lie groups one will typically choose a basis for the correspondingLie algebra. If the Lie group is a matrix group, the evaluation of a flow willamount to the computation of a matrix exponential. In certain applications thiscan be done cheaply, for instance if the group consists of rotations in 3D, onecan make use of the Rodrigues formula [5] for this purpose. In the general case,the computation of matrix exponentials may be time consuming.
As for classical Runge–Kutta methods, implementation issues include errorcontrol and stepsize variation strategies. It it natural to look for embeddedpairs of method coefficients. This, however, requires a metric on the manifold,typically realized by using a specific coordinate system.
It seems that the type of method presented here is only one of many possible
INTEGRATION METHODS BASED ON RIGID FRAMES 141
generalizations of Runge–Kutta methods. A few other approaches have beenproposed by Munthe–Kaas [10, 11]. They differ from the Crouch–Grossmanmethods in the sense that the coupling between internal stages is performedon vector fields rather than on their flows. Generally the two approaches arenot equivalent. We believe that there may be other feasible approaches, notyet discovered, but perhaps the analysis of similar techniques can be done in asimilar way to the one presented here. The order conditions derived here aregeneral in the sense that they do not depend on the choice of frame vector fields.We believe that such a criterion may be too restrictive, since our intention is notto design more general schemes, but rather to accommodate integration methodsfor particular manifolds found in applications. Therefore, it would be importantto answer the question of whether the set of order conditions can be reduced forparticular choices of manifolds and frame vector fields and thereby obtain lessexpensive schemes tailored for a particular type of application. For instance, inmany problems in structural dynamics, the configuration manifold consists ofsemidirect products [17, 13] between the rotation group and translations in R
3.One of the benefits of these new approaches for numerical integration on man-
ifolds is their suitability for object-oriented implementations. For instance, inthe definition of the methods, no explicit reference is made to coordinates or tothe particular representation of objects like flows or vector fields.
Acknowledgements.
We would like to thank Finn Faye Knudsen, Anne Kværnø, Hans Munthe-Kaas, and Bruno D. Welfert for useful discussions and comments during thework with this paper.
REFERENCES
1. V. I. Arnold, Mathematical Methods of Classical Mechanics, 2nd ed., Springer-Verlag, Berlin, 1989.
2. D. I. A. Cohen, Basic Techniques of Combinatorial Theory, Wiley, New York,1978.
3. P. E. Crouch and R. Grossman, Numerical integration of ordinary differentialequations on manifolds, J. Nonlinear Sci., 3 (1993), pp. 1–33.
4. P. E. Crouch, Y. Yan, and R. Grossman, On the numerical integration of the rollingball equations using geometrically exact algorithms, J. Mech. Struct. and Mach.,23 (1995), pp. 257–272.
5. M. L. Curtis, Matrix Groups, 2nd ed., Springer-Verlag, Berlin, 1984.
6. E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential EquationsI, Springer Verlag, 1993.
7. D. J. Higham, Runge–Kutta type methods for orthogonal integration, in AppliedNumerical Mathematics (special issue celebrating the Runge–Kutta centenary),J. C. Butcher, ed., vol. 22 (1996), pp. 217–223.
8. A. Kværnø, ESDIRK methods, manuscript, 1996.
9. A. Marthinsen, H. Munthe-Kaas, and B. Owren, Simulation of ordinary differentialequations on manifolds—some numerical experiments and verifications, Modeling,Identification and Control, 18 (1997), pp. 75–88.
142 B. OWREN AND A. MARTHINSEN
10. H. Munthe-Kaas, Lie–Butcher theory for Runge–Kutta methods, BIT, 35:4 (1995),pp. 572–587.
11. H. Munthe-Kaas, Runge–Kutta methods on Lie groups, BIT, 38:1 (1998), pp. 92–111.
12. H. Munthe-Kaas and A. Zanna, Numerical integration of differential equations onhom*ogeneous manifolds, in Foundations of Computational Mathematics, F. Cuckerand M. Shub, eds., Springer-Verlag, Berlin, 1997, pp. 305–315.
13. H. Munthe-Kaas and A. Zanna, Semidirect products, manuscript, 1997.
14. B. Owren and B. Welfert, The Newton iteration on Lie groups, Technical ReportNumerics No. 3/1996, Norwegian University of Science and Technology, Trondheim,Norway, 1996.
15. S. Tracogna and B. Welfert, Notes on riffle shuffles, to appear.
16. S. Tracogna and B. Welfert, Spectral analysis of generalized top to random shuffles,Technical Report TW-97-04, Department of Mathematics and Computer Science,Leiden University, The Netherlands, 1997.
17. V. S. Varadarajan, Lie Groups, Lie Algebras, and their Representations, GTM102, Springer-Verlag, Berlin, 1984.
18. F. W. Warner, Foundations of Differentiable Manifolds and Lie Groups, GTM 94,Springer-Verlag, Berlin, 1983.
