Arieh Iserles
email: a.iserles@damtp.cam.ac.uk

Hans Z Munthe-Kaas

Syvert P Nørsett
email: norsett@math.ntnu.no

Antonella Zanna

Lie-group methods

Many differential equations of practical interest evolve on Lie groups or on manifolds acted upon by Lie groups. The retention of Lie-group structure under discretization is often vital in the recovery of qualitatively-correct geometry and dynamics and in the minimisation of numerical error. Having introduced requisite elements of differential geometry, this paper surveys the novel theory of numerical integrators that respect Lie-group structure, highlighting theory, algorithmic issues and a number of applications. 1 CONTENTS 1 Numerical analysts in Plato's temple 2 2 Theory and background 8 3 Runge-Kutta on manifolds and RK-MK 36 4 Magnus and Fer expansions 40 5 Quadrature and graded algebras 54 6 Alternative coordinates 70 7 Adjoint methods 88 8 Computation of exponentials 96 9 Stability and backward error analysis 103 10 Implementation, error control and DiffMan 107 11 Applications 112 References 132 A List of methods 137 B Fast computation of 3D rotations

1. Numerical analysts in Plato's temple "Ageometretos medeis eisito": let nobody enter who does not understand geometry. These were the words written at the entrance to Plato's Temple of the Muses. Are numerical analysts welcome in Plato's temple?

Historically, the answer is negative. Computational mathematics is all about rendering mathematical phenomena in an algorithmic form, amenable to sufficiently precise, affordable and robust number crunching. A mathematical phenomenon can be approached in one of two ways: either by exploring its qualitative features (which, to a large extent, are synonymous with geometry or, at the very least, can be formulated in geometric terminology) or by approximating its quantitative character. Although only purists reside completely at either end of the spectrum, it is fair to point out that numerical analysis, by its very 'rules of engagement', is what 'quantitative mathematics' is all about. Ask a numerical analyst "How good is the solution?" and the likely answer will address itself to a subtly different question, "How small is the magnitude of the error?".

In principle, the emphasis on quantitative aspects in mathematical computing served it well. It is hard to imagine modern technological civilisation without the multitude of silent computer programs in the background, flying the airplanes, predicting the weather, making sense of CAT scans, controlling robots, identifying fingerprints, keeping reactions from running away and predicting the behaviour of stock markets. This is the success story of numerical analysis, of this 'quantitative number crunching', and nothing should be allowed to obscure it. So, perhaps if we are doing so well everywhere else, we might cede Plato's temple to our more 'pure' brethren and sisters: let them engage in sterile intellectual discourse while we change the world!

The main contention of this review and of the emerging discipline of geometric integration is that this approach, although tempting, is at best incomplete, at worst badly misguided. The history is not just a heroic tale of numerical algorithms fleshing out mathematical concepts as numbers and graphs. Progress has always occurred along parallel, intertwined tracks: both better theoretical understanding of qualitative attributes of a mathematical construct and its better computation. The airplane-flying, weather-predicting and CAT-scanning programs can do their job only because they deliver an answer that explains in a satisfactory manner qualitative features, as well as producing the 'right' numbers! Indeed, an artificial dichotomy of quantitative and qualitative aspects of mathematical research is in our opinion misleading and it serves ill mathematical and applied communities alike.

On the one hand, computation tells pure mathematics what to prove. Phenomena are often initially identified when observed under discretization and subsequently subjected to the full rigour of mathematical analysis. A familiar case in point is the discovery of solitons in the solution of the Korteweg-de Vries equation by [START_REF] Zabusky | Interaction of solitons in a colisionless plasma and the recurrences of initial states[END_REF], an event which launched a whole new mathematical discipline, other examples abound. Indeed, we are so used to rely on the computer as a laboratory of pure mathematics that it is difficult to imagine the heroic work of Gaston [START_REF] Julia | Memoire sur l'iteration des fonctions rationalles[END_REF] on the geometry of fractals while bearing in mind that he has had no access to computers, never able to easily calculate a sequence of rational iterations or to visualise a fractal on a computer monitor! On the other hand, qualitative analysis tells computation, quite literally, what to compute. Every seasoned numerical analyst knows that the procedure of 'discretise everything in sight and throw it on a computer' works only with toy problems. The more we know about the qualitative behaviour of the underlying mathematical construct, the more we can identify the right computational approach, concentrate resources at the right place, focus on features that influence more the quantitative behaviour and, by the conclusion of the computation, have well-founded expectation that the graph on the computer monitor corresponds to a genuine solution of the problem in hand. Moreover, consumers of numerical calculations are not interested just in numbers, graphs and impressive visualisation. Very often it is the qualitative features, most conveniently phrased in the language of geometry, that draw genuine interest in applications: periodicity, chaoticity, conservation of energy or angular momentum, reduction to lower-dimensional manifolds, symmetry, reversibility, . . . . The contention of this review is not just that the contribution of geometry to computation in the special case of time-evolving systems of differential equations is absolutely crucial, but that the terminology of differential geometry, and in particular Lie groups, creates the right backdrop to this process. The following example will help to elucidate this point, while serving as a convenient introduction to the theme of this paper.

Let us denote by S N the set of all N × N real symmetric matrices and consider the solution of the isospectral flow

Y ′ = B(Y )Y -Y B(Y ), t ≥ 0, Y(0) = Y 0 ∈ S N , (1.1)
where the (sufficiently smooth) function B maps S N to N × N real skew-symmetric matrices. The solution itself remains in S N for all t ≥ 0. Such flows occur in a variety of applications. Perhaps the earliest (and the best known) is the Toda lattice of material points subjected to nearest-neighbour interaction. It has been demonstrated by Flaschka that, in the case of an exponential interaction potential, the underlying Hamiltonian system can be rendered in the form (1.1), where Y is tridiagonal and B maps ⎡ 1981). Another important application of (1.1) is to Lax pairs in fluid dynamics, whence S N need be replaced by a suitable function space, B is a differential operator and the outcome is a partial differential equation of a hyperbolic type [START_REF] Toda | Theory of Nonlinear Lattices[END_REF]). Before we mention another application of isospectral flows, we need to single out their most remarkable qualitative feature which, coincidentally, explains their name: as the time evolves, the eigenvalues of Y (t) stay put! Upon a moment's reflection, this renders such flows interesting in the context of numerical algebra. Indeed, the classical QR algorithm is intimately related to sampling the solution of (1.1) at unit intervals [START_REF] Deift | Ordinary differential equations and the symmetric eigenvalue problem[END_REF]. Many other iterative algorithms can be phrased in this terminology and, perhaps more importantly, many interesting algorithms rely on this construct in the first place. Pride of place belongs here to methods for the inverse eigenvalue problem: seeking a matrix of a given structure that possesses a specified set of eigenvalues (or singular values). Such problems are important in a wide range of applications, ranging from the theory of vibrations to control theory, tomography, system identification, geophysics, all the way to particle physics. Isospectral flows are a common denominator to perhaps the most powerful approach toward the design of practical algorithms for the inverse eigenvalue problem, which has been pioneered in the main by [START_REF] Chu | Inverse eigenvalue problems[END_REF]. Suppose that we are seeking a matrix in a class T ⊂ S N with the eigenvalues η ∈ R N . Often it is possible to design a matrix function B so that attractive fixed points of (1.1) lie in T . In that case, letting Y (0) = diag η results in a flow that converges to the solution of the inverse eigenvalue problem.

⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ α 1 β 1 0 • • • 0 β 1 α 2 β 2 . . . . . . 0 . . . . . . . . . 0 . . . . . . β N -2 α N -1 β N -1 0 • • • 0 β N -1 α N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ to ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 β 1 0 • • • 0 -β 1 0 β 2 . . . . . . 0 . . . . . . . . . 0 . . . . . . -β N -2 0 β N -1 0 • • • 0 -β N -1 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ( Toda
As an example of such a procedure we mention the inverse eigenvalue problem for Toeplitz matrices. Thus, T consists of symmetric N × N Toeplitz matrices:

X ∈ T ⇔ x k,l = t |k-l| , k,l = 1, 2, . . . , N,
where t 0 , t 1 , . . . , t N -1 are arbitrary real numbers. Such problems are important in the design of control systems but, remarkably, even the very existence of a solution has been until very recently an open problem, which has been answered by [START_REF] Landau | The inverse eigenvalue problem for real symmetric Toeplitz matrices[END_REF] in a beautiful, yet non-constructive, existence proof. Following [START_REF] Chu | On a differential equation dX dt = [X, k(X)] where k is a Toeplitz annihilator[END_REF][START_REF] Trench | Numerical solution of the inverse eigenvalue problem for real symmetric Toepliz matrices[END_REF], we let

b k,l (Y ) = ⎧ ⎨ ⎩ y k,l-1 -y k+1,l , 1 ≤ k < l ≤ N, 0, 1 ≤ k = l ≤ N, y k+1,l -y k,l-1 , 1 ≤ l < k ≤ N
be a Toeplitz annihilator. Note that B is indeed skew symmetric and that B(Y ) = O for Y ∈ T : thus, a solution is a fixed point. While remarking that many important questions with regard to the convergence of the above algorithm are still wide open, we should draw the reader's attention to a crucial observation. For numerical purposes, sooner or later we must replace (1.1) by a computational time-stepping scheme. Will such a scheme respect the eigenstructure of Y ? This is not simply an optional extra since the whole point of the exercise is to evaluate an answer in T . Yet, as proved in [START_REF] Calvo | Numerical solution of isospectral flows[END_REF]), the most popular numerical methods, multistep and Runge-Kutta schemes, do not respect isospectral structure and they fail to converge to the correct element of T : the error on the eigenvalues, of the same order of magnitude as the error in the numerical trajectory itself, is unacceptable.

An alternative, proposed by [START_REF] Calvo | Numerical solution of isospectral flows[END_REF], is to observe that all the elements of the isospectral manifold

I(η) = {X ∈ S N : σ(X) = η},
where σ(X) denotes the spectrum of X, can be written in the form X = QY 0 Q T , where Y 0 = diag η is our initial condition and Q ∈ SO(N ), the set of all N × N real orthogonal matrices with unit determinant. The main idea is to seek, in place of (1.1), a differential equation that is satisfied by Q(t) in the representation

Y (t) = Q(t)Y 0 Q(t) T , t ≥ 0.
(1.2)

It is easy to ascertain that this equation has the form

Q ′ = B(QY 0 Q T )Q, t ≥ 0, Q(0) = I (1.3)
and that, provided we can solve it while retaining orthogonality, we can easily recover the solution of the original isospectral flow.

To illustrate our point, let us consider a simple numerical experiment. We choose Fig. 1.1. Toeplitz error and error on the eigenvalues versus number of iterations for a 5 × 5 symmetric inverse eigenvalue problem when solved with the Forward Euler scheme (FE) and the Cayley-based Forward Euler (CFE). Although both methods converge asymptotically to a Toeplitz matrix, the error on the eigenvalues of CFE stays within machine accuracy while the error of FE is completely determined by the choice of integration step size.

η = [1, 2, 3, 4, 5] T , and solve the symmetric inverse eigenvalue problem (1.1), where the matrix function B is the Toeplitz annihilator introduced above. We consider first the standard Forward Euler (FE) scheme,

Y n+1 = Y n + h[B(Y n ), Y n ], n∈ Z + ,
with initial condition Y 0 = diag(η) and step size h = 1 10 . The Toeplitz error, ∥B(Y n )∥ 2 , is plotted in Figure 1.1 as a dotted line. Clearly, as n tends to infinity, the Toeplitz error becomes progressively smaller and Y n tends to a Toeplitz matrix, as the theory predicted. Do the eigenvalues stay put? The answer is negative, since the FE scheme is not isospectral, as proven in [START_REF] Calvo | Numerical solution of isospectral flows[END_REF]. Therefore the error on the eigenvalues after the first step is of the same order of the error of FE and is carried along the whole integration (dash-dotted line in Figure 1.1).

Consider next the iteration

Q n+1 = I -1 2 hB(Y n ) -1 I + 1 2 hB(Y n ) , Y n+1 = Q n+1 Y n Q T n+1 , n ∈ Z + ,
which is equivalent to solving (1.3) with a modified version of the Forward Euler scheme based on the Cayley expansion (note that Q n+1 is orthogonal) in tandem with a similarity transformation for the update Y n+1 . This scheme, to which we will refer as CFE (Cayley-type Forward Euler), is explicit, has the same order of accuracy and requires only slightly more computations than the more classical FE. However, unlike FE, it is isospectral by design and preserves the eigenvalues to machine accuracy. In conclusion, FE tends to a Toeplitz matrix with wrong eigenvalues, while CFE, a simple modification of FE that instead preserves the qualitative features of the flow, tends to a Toeplitz matrix with the right eigenvalues.

The Toeplitz and the eigenvalue error of CFE are displayed in Figure 1.1 and correspond to the dashed and 'plus' lines respectively. The set SO(N ) in which the matrix Q of the above example evolves is an instance of a Lie group, a concept that will be described and debated in great detail in Section 2, while (1.2) and (1.3) are special instances of a group action and a Liegroup equation.

Let us comment briefly on the contents of this review. In Section 2 we have assembled the common mathematical denominator underlying this paper: elements of differential geometry, Lie groups and algebras, homogeneous spaces and differential equations evolving on such objects. Section 3 is devoted to Runge-Kutta-Munthe-Kaas schemes, the most natural approach to Lie-group solvers in our setting. In Section 4 we describe expansions, originally due to Magnus and to Fer, which can be converted into interesting computational tools. Section 5 is concerned with a make-or-break issue for many Lie-group methods, multivariate quadrature of multilinear forms over polytopes. We demonstrate there that some very technical tools from Lie-algebra theory can be used to a great effect in reducing the numerical cost. Lie-group methods are typically based on local imposition of a convenient coordinate system in the group. In Section 6 we debate less conventional choices of the coordinate map, which are suitable for important Lie groups and equations of practical interest. The theme of Section 7 is time-symmetric methods that, by design, exhibit many favourable features, while the concern of Section 8 is the practical approximation of a matrix exponential from a Lie algebra, so that the result lies in the right Lie group. Stability issues are addressed in Section 9, while Section 10 reviews practical issues of implementation and error control and introduces the DiffMan package. A sample of the many applications of Lie-group methods is presented in Section 11. Finally, in Appendix A we list practical Lie-group methods, while Appendix B displays useful explicit formulae for integration in SO(3), perhaps the single Lie group with most important relevance to problems in science and engineering.

The purpose of this survey is not to cover the entire corpus of Lie-group methods but to present a unified introduction to a young discipline that is likely to undergo many exciting further developments. We have omitted many interesting methods and papers, with due apologies to their authors, to keep our narrative more focused and clear. Only the future can tell which methods and techniques will survive.

It is vital throughout the paper to distinguish what exactly is the type of objects under consideration. We will be often mixing in our formulae elements of Lie groups and Lie algebras, scalars, matrices and vectors. To assist the reader, we have adhered to a consistent naming convention:

• Elements in a Lie algebra are denoted by the Roman letters a, b, . . . , h and A, B, . . . , H and by the Greek letters α, β, . . . , ξ and ∆, Θ, Ξ. • Elements in a Lie group are denoted by the Roman letters p, q, . . . , z and P, Q, . . . , Z and by the remaining Greek letters: π, ρ, . . . , ω and Υ, Φ, Ψ, Ω. • Elements in an abstract construct (e.g., an abstract Lie group) are denoted mostly with lower-case letters. However, as soon as we are concerned with specific representation, we reserve lower-case letters for scalars, upper case for matrices and lower-case, boldfaced letters for vectors. • Special elements deserve special names. Thus, I is the identity in a Lie group, while O is the zero of a Lie algebra. A generic Lie group will be denoted by G and a generic Lie algebra by g. In general, we reserve Gothic font for Lie algebras.

• As in all naming conventions, we make obvious compromises with standard mathematical practice and common sense. Proper names remain unchanged: thus, sin t is the familiar sine function, not an element of a group, {B k } k∈Z + are Bernoulli numbers, Φ is the set of roots of a Lie algebra, h is the step size of a time-stepping algorithm and so on. We also employ a plethora of integration variables, summation indices, constants etc. Occasionally variables evolve in structures which are neither Lie groups nor Lie algebras, e.g. the isospectral manifold I(η). In all these cases, which should be obvious from the context, we use ad hoc notation.

Theory and background

Lie groups and Lie algebras are mathematical objects which have originated in the seminal work of Sophus Lie (1842-1899) on solving differential equations by quadrature, using symmetry methods. Originally these concepts were quite concrete, related to flows of differential equations on R N . Early in the twentieth century an abstract view of Lie group theory emerged, commencing from the work of Elie Cartan on the classification of Lie algebras. The advantage of abstract formulation is that it simplifies mathematical analysis, and hence this presentation has become dominant throughout mathematical literature. However, the abstract theory is concentrating on understanding mathematical structures rather than exposing applications in solving differential equations. Hence it is not at all clear to most applied mathematicians that Lie groups are really very useful objects also in applied and computational mathematics, and it might be difficult to find the motivation to learn an abstract theory.

We believe that the original idea of arriving at Lie algebras via continuous actions on a domain should be an excellent starting point for computationally-oriented mathematicians, and in fact for many applications it is important to keep this view in focus. In this presentation we will commence from this perspective and gradually move towards somewhat more abstract formulations. Eventually, in Section 2.5 we will return to concrete matrix formulation, concentrating on the numerical solution of matrix differential equations of the form

Y ′ = A(t, Y )Y, t ≥ 0, Y(0) = Y 0 ,
where Y and A(t, Y ) are N × N matrices. It turns out that all our solution techniques can be derived in this concrete matrix setting and without major modifications they can be applied to more general situations. All the algorithms of this paper will be derived within the matrix framework, and hence the theory from this point and up to Section 2.5 might be read in a relaxed manner, without the need to master all the details at first reading. Numerical integration of ordinary differential equations (ODEs) is traditionally concerned with solving initial value problems evolving on R N ,

y ′ = f (t, y), t ≥ 0, y(0) = y 0 , y(t) ∈ R N ,
where f is a vector field on R + × R N . Well-known numerical integrators, such as Runge-Kutta and multistep methods, advance a time-stepping procedure by adding vectors in R N , y n+1 = y n + ha n ,

where a n = a n (h, y n , . . .) is computed by the given numerical method and h is the time step. One might say that classical integrators are formulated using a set of 'basic motions' given by translations on R N to advance the numerical solution.

A major motivation for Lie group methods is the possibility of replacing the domain R N with more general configuration spaces and replacing translations on R N by more general families of 'basic motions' on the domain. For example, if y(t) is a vector known to evolve on a sphere, one might consider rotations y n+1 = Q n y n , where Q n is an orthogonal matrix, as basic motions. We have already encountered another example, isospectral flows (1.1). In that case Y (t) evolves on the isospectral manifold I(σ(Y 0 )) and this is the right configuration space. Moreover, the natural 'basic motions' are (1.2) and they rest upon the fact that any two elements in I(σ(Y 0 )) are similar via an orthogonal matrix. A generalisation of the last example does not require any more that Y 0 is symmetric and that the matrix function B is skew symmetric: any nonsingular initial value and sufficiently smooth matrix function will do. The flow remains isospectral but the elements of the configuration space are no longer orthogonally similar. The representation (1.2) is valid, however, if we allow Q(t) to range across all possible nonsingular real matrices and replace Q T with Q -1 .

An important reason why a manifold, rather than the entire R N , is a suitable configuration space is that it often expresses crucial geometric attributes of the underlying differential system, e.g. conservation laws, symmetries or symplectic structure. As will be seen in the sequel, an added bonus of this approach is that it frequently leads to interesting numerical advantages, in particular to slower error accumulation.

In seeking abstractions and generalisations of classical numerical methods it is important to bear in mind abstractions in pure-mathematical treatment of differential equations. However, the transition from pure to computational mathematics is not straightforward. Whenever a pure mathematician says "There exists an animal such that. . . ", an applied mathematician must add questions like "Can we compute this animal efficiently?" and "How can we represent it in software?".

Vector fields and flows on manifolds

A cornerstone of all abstract mathematical presentations of differential equations is the concept of differential manifolds as the definition of domains on which differential equations evolve. A good general introduction to manifolds and to differential geometry are [START_REF] Abraham | Foundations of Mechanics[END_REF] and [START_REF] Guillemin | Differential Topology[END_REF].

Intuitively one should think of a d-dimensional manifold as being a smooth domain which in a (small) neighbourhood of any point 'looks like' R d , but typically looks different globally. It is known that any d-dimensional manifold can be represented as a d-dimensional surface embedded in R N for some N ≥ d. This is a conceptually very useful view of manifolds. It is, however, also important to know that all geometric properties of manifolds exist independently of any particular representation1 .

It is possible to discuss these properties in a coordinate-free language, which is independent of any particular representation. Although this abstract presentation is mathematically very elegant and it provides vital clues toward writing good software in an object-oriented language [START_REF] Engø | DiffMan -an object oriented MATLAB toolbox for solving differential equations on manifolds[END_REF], it requires quite a bit of work to present and comprehend. The advantage of the abstract language is that it focuses on the essential structures. However, in most of the applications we will deal with, the manifolds exist naturally as surfaces embedded in R N , and furthermore it is fully possible to understand and use the numerical techniques we will discuss in this paper without knowing the abstract theory of manifolds. We have therefore decided to base our discussion on the following very concrete definition.

Definition 2.1 A d-dimensional manifold M is a d-dimensional smooth surface M ⊂ R N for some N ≥ d.
It should be made crystal clear that all the numerical techniques we are about to present rely solely on those properties of M that exist independently of any particular embedding in R N . We believe that a reader with knowledge of coordinate-free is a manifold for ρ ∈ (0, 1) but not when ρ = 1, because of a singularity at the origin. • An important representation of a manifold is as a smooth subset of solutions of a smooth algebraic equation, g(x) = 0. Thus, for example, g(x) = ∥x∥ 2 2 -1 defines a unit sphere.

• The algebraic-equation representation is of direct relevance to geometric integration, since conservation laws and integrals of differential systems are nothing else but algebraic equations constant along the solution trajectory. Thus, for example, a Hamiltonian system with Hamiltonian energy H(p, q) evolves on the manifold {H(p, q) = H(p(0), q(0)

) : p, q ∈ R N }. • The set O(N ) of all N ×N orthogonal matrices is a manifold, since X ∈ O(N ) is equivalent to g(X) = ∥X T X -I∥ 2 2 = 0.
So is also SL(N ), the set of all N × N matrices with unit determinant, since g(X) = det X -1 is a smooth function.

• A Stiefel manifold is the set of all real M ×N matrices X such that X T X = I: typically M > N, such matrices are 'long and skinny' and XX T ̸ = I. • A Grassmann manifold is a Stiefel manifold, equivalenced by O(M ). In other words, we identify

X 1 , X 2 such that X T 1 X 1 = X T 2 X 2 = I if there exists an orthogonal M × M matrix Q such that X 1 = QX 2 .
An alternative interpretation of a Grassmann manifold is as the set of all N -dimensional subspaces of R M .

The single most important property of a manifold is the existence of tangents to the manifold in any point p ∈ M. If we think of the manifold as a surface in R N , then a tangent at p can be defined as a vector a such that dist(p + εa, M) = O(ε 2 ). This construction of a tangent is relying on the embedding of M in R N , on the linear-space structure of R N and even on the metric structure of R N . An alternative way to define tangents is by differentiating a curve. This approach has the advantage of making no use of the embedding of M in R N and hence makes sense also on general manifolds. Definition 2.2 Let M be a d-dimensional manifold and suppose that ρ(t) ∈ M is a smooth curve such that ρ(0) = p. A tangent vector at p is defined as

a = dρ(t) dt t=0 .
The set of all tangents at p is called the tangent space at p and denoted by TM| p . It has the structure of a d-dimensional linear space: if a, b ∈ TM| p then a+b ∈ TM| p and α a ∈ TM| p for any real α. The collection of all tangent spaces at all points p ∈ M is called the tangent bundle of M and denoted by TM = p∈M TM| p .

Note that whereas it is fine to add tangents based at the same point , there is in general no rule for adding tangent vectors based at different points. Thus, to specify a tangent completely we need to provide both the basepoint p and the tangent itself. Hence TM is a 2d-dimensional space, with elements (p, a) consisting of all possible tangents a for all possible basepoints p. Definition 2.3 A (tangent) vector field on M is a smooth function F : M → TM such that F (p) ∈ TM| p for all p ∈ M. The collection of all vector fields on M is denoted by X(M).

Addition and scalar multiplication of vector fields are defined pointwise in a natural way as

(F + G)(p) = F (p) + G(p) and (αF )(p) = α(F (p)). If F, G ∈ X(M) then also F + G ∈ X(M) and αF ∈ X(M) for all real α.
Definition 2.4 Let F be a tangent vector field on M. By a differential equation (evolving) on M we mean a differential equation of the form

y ′ = F (y), t ≥ 0, y(0) ∈ M, (2.1)
where F ∈ X(M). Whenever convenient, we allow F in (2.1) to be a function of time, F = F (t, y). The flow of F is the solution operator Ψ t,F : M → M such that y(t) = Ψ t,F (y 0 ) solves (2.1).

Note that we can find the vector field F from Ψ t,F by differentiation,

F (y) = d dt Ψ t,F (y) t=0 .
F is often called the infinitesimal generator of the flow Ψ t,F . By reparametrizing time (or scaling the vector field) we can see that the flow operator satisfies the identity

Ψ α,F = Ψ 1,αF .
(2.

2)

The task of computing the flow of a given vector field is often called the exponentiation of the vector field . We will occasionally employ the notation

Ψ 1,F ≡ exp(F ) ⇔ Ψ t,F ≡ exp(tF ).
Computation of a flow is a particular example of an exponential map. We will return to a more general definition of exponential maps later. The notation Ψ will be used when we want to emphasise that we are discussing flows of vector fields, and exp whenever the map can be conveniently given the more general interpretation.

Recall our goal of integrating (2.1) numerically, using a chosen set of 'basic motions' to advance the numerical solution. If the analytical solution evolves on M it is natural to choose a set of basic motions which are everywhere tangent to M, which will produce a numerical solution also evolving on M. It is useful to consider also these basic motions as flows of a finite or infinite collection of vector fields, B 1 , B 2 , . . .. If we want an efficient solver, we must be able to compute the flow, i.e. exponentiate these B i s, efficiently.

From the standpoint of geometry, numerical integration of ODEs is concerned with the task of approximating the exponential of a general vector field F by exponentials originating in a family of simpler vector fields B 1 , B 2 , . . .. Example 2.2 Let M = R N and let T a stand for the constant vector field, T a (y) = a for some vector a ∈ R N . The flow of T a is translation along a, Ψ t,Ta (y 0 ) = y 0 + ta.

The set of all translations can be obviously used to advance the numerical solution in any desired direction on R N . Note that translations commute,

Ψ t1,Ta •Ψ t2,T b (y 0 ) = y 0 + t 1 a + t 2 b = Ψ t2,T b •Ψ t1,Ta (y 0 ).
Generally flows do not commute. We will later see that a major difference between traditional numerical integrators and Lie-group methods is that the former are

Φs,t(y0) • Ψs,F 2 2 y0 st[F,G](y0)+h.o.t c c r r r r r r r r r • Ψt,G K K 7 5 3 1 ) ' % Ψs,F , , • Ψt,G B B ! Ù × Fig. 2.2. A geometric interpretation of the commutator.
based on a set of commuting flows whereas the latter allow more general flows as basic movements to advance the solution.

The degree to which the flows of two vector fields fail to commute is measured by the commutator of the two vector fields. Consider two general vector fields F and G and their flows Ψ s,F and Ψ t,G . To investigate commutativity, we form the following composition of the flows, 2.2). Obviously the flows commute if Φ s,t (y) = y for all s, t ≥ 0 and all y ∈ R N . However, this is a nonlinear condition that in many cases may be difficult to compute or verify. One of the fundamental ideas due to Sophus Lie was the observation that such nonlinear conditions can be turned into equivalent linear infinitesimal conditions. To accomplish this we want to linearise Φ s,t for small s and t. Since Φ 0,t (y) = Φ s,0 (y) = y for all y, we must have

Φ s,t = Ψ s,F •Ψ t,G •Ψ -s,F •Ψ -t,G ≡ exp(sF )• exp(tG)• exp(-sF )• exp(-tG) (2.3) (see Figure
Φ s,t (y) = y + stH(y) + O(s 2 t) + O(st 2 ) ( 2 . 4 )
for some vector field H. This vector field H is called the commutator , or the Jacobi bracket of F and G, and it is written as

H = [F, G].
Lemma 2.1 Given two vector fields F, G on R N , the commutator H = [F, G] can be computed componentwise at a given point y ∈ R N as

H i (y) = N j=1 G j (y) ∂F i (y) ∂y j -F j (y) ∂G i (y) ∂y j .
(2.5) (Note that many authors define the commutator with an opposite sign. The commutator as given here is often called the (-)Jacobi bracket.)

Proof. From (2.3) and (2.4) we get

H = ∂ 2 ∂s∂t exp(sF )• exp(tG)• exp(-sF )• exp(-tG) s=t=0 .
Since ∂ ∂s exp(sF )• exp(tG)• exp(-sF ) s=0 = 0, this simplifies to

H = [F, G] = ∂ 2 ∂s∂t exp(sF )• exp(tG)• exp(-sF ) s=t=0 .
(2.6)

Neglecting higher order terms in s and t, Euler's integration scheme yields

Ψ s,F •Ψ t,G •Ψ -s,F (y) = y -sF (y) + tG(y 1 ) + sF (y 2 ) + h. o. t,
where y 1 = y -sF (y) , y 2 = y -sF (y) + tG(y 1 ).

Hence (2.6) implies that

H i (y) = ∂ ∂s G i (y 1 ) + ∂ ∂t F i (y 2 ) s=t=0 = N j=1 - ∂G i (y) ∂y j F j (y) + ∂F i (y) ∂y j G j (y) .

✷

We will now review the most salient properties of the commutator.

Lemma 2.2 If F, G ∈ X(M) then H = [F, G] ∈ X(M).
Proof. The function

ρ(t) = exp( √ tF )• exp( √ tG)• exp(- √ tF )• exp(- √ tG)(y 0 )
is a curve that evolves on M. The tangent defined by this curve is [F, G](y 0 ). ✷

Dividing the polygon of Figure 2.2 into infinitesimally small rectangles, we can verify that Lemma 2.3 Two flows Ψ s,F and Ψ t,G commute if and only if [F, G] = 0. From (2.5) one may prove the following important features of the commutator which should be familiar in the special case (which we will encounter soon again) of a commutator of two matrices.

Lemma 2.4 The commutator of vector fields satisfies the identities

[F, G] = -[G, F ] (skew symmetry), (2.7) [αF, G] = α[F, G] for α ∈ R, (2.8) [F + G, H] = [F, H] + [G, H] (bilinearity), (2.9) 0 = [F, [G, H]] + [G, [H, F ]] + [H, [F, G]] (Jacobi's identity). (2.10)
Example 2.3 Let L A denote the linear vector field on R N , given by some matrix A, that is L A (y) = Ay. The solution of the linear equation y ′ = Ay is given as

y(t) = ∞ j=0 (tA) j j! y 0 = expm(tA)y 0 ,
where expm denotes the classical matrix exponential. Hence

Ψ t,LA (y 0 ) ≡ exp(tL A )(y 0 ) = expm(tA)y 0 .
(This motivates the name 'exponentiation' for computing the flow.) Now let us compute the commutator of two linear vector fields from (2.5),

[L A , L B ] i (y) = j,k,l B j,k y k ∂A i,l y l ∂y j -A j,k y k ∂B i,l y l ∂y j = j,k (A i,j B j,k -B i,j A j,k ) y k . Thus [L A , L B ] = L C where C = AB -BA, (2.11)
the familiar definition of a commutator from linear algebra. Note that linear vector fields constitute a complete family of vector fields, closed under commutators and linear combinations, L A + L B = L A+B and αL A = L αA .

In applications of Lie group integrators to partial differential equations (PDEs) it is often useful to consider a more general version of (2.5). Let y be a point in a (finite or infinite-dimensional) linear space M. By a vector field F on M we mean some operator (linear or nonlinear) such that the (ordinary or partial) differential equation

∂y ∂t = F (y)
is well defined. An infinite-dimensional example is a parabolic PDE, where y belongs to some function space on a domain and F is a spatial differentiation operator, e.g.

F (y) = ∇ 2 y. If F and G are two vector fields on M then [F, G](y) = ∂ ∂s [F (y + sG(y)) -G(y + sF (y))] s=0 .
(2.12)

The proof is very similar to that of Lemma 2.1 and it is straightforward to verify that in the finite-dimensional case (2.12) reduces to (2.5). Note that if F and G are linear operators then (2.12) yields immediately [F, G](y) = F (G(y)) -G(F (y)), as we saw in Example 2.3.

Example 2.4 Let y ∈ R N and consider the set of all affine linear vector fields

F (A,a) (y) = Ay + a
where A is an N × N matrix and a ∈ R N . Let us compute the commutator of two such vector fields F (A,a) and F (B,b) . By inserting these vector fields in (2.12) we obtain

F (A,a) , F (B,b) (y) = ∂ ∂s {A(y + s(By + b)) + a -B(y + s(Ay + a)) -b} s=0 = ABy + Ab -BAy -Ba.
Thus

F (A,a) , F (B,b) = F (C,c) , where (C, c) = (AB -BA, Ab -Ba).
(2.13)

We note that the set of all affine linear vector fields is yet another example of a collection of vector fields closed under linear combination and commutation.

Lie algebras, Lie groups and Lie group actions

A problem, fundamental to numerical analysis of differential equations on manifolds, is to determine the set of all possible flows that can be obtained by composing a given set of basic flows. If we restrict the discussion to 'sufficiently small t', important information is provided by the so-called BCH formula.

Theorem 2.5. (Baker-Campbell-Hausdorff ) For sufficiently small t ≥ 0 we have exp(tF )• exp(tG) = exp(tH),

where H = bch(F, G) can be constructed from iterated commutators of F and G.

The first few terms are

H = F + G + 1 2 t[F, G] + 1 12 t 2 ([F, [F, G]] + [G, [G, F ]]) + O(t 3 ).
Higher order terms can be obtained by recursion [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF].

Definition 2.5 A Lie algebra of vector fields is a collection of vector fields which is closed under linear combination and commutation. In other words, letting g denote the Lie algebra,

B ∈ g ⇒ αB ∈ g for all α ∈ R. B 1 , B 2 ∈ g ⇒ B 1 + B 2 , [B 1 , B 2 ] ∈ g
Given a collection of vector fields B = {B 1 , B 2 , . . .}, the least Lie algebra of vector fields containing B is called the Lie algebra generated by B.

We arrive at the following conclusion. Let g be the Lie algebra generated by the set B = {B 1 , B 2 , . . .} of vector fields. For small t, the combination of flows of vector fields in B yields the flow of a vector field in g. Furthermore, the flow of any vector field in g can, provided t ≥ 0 is small enough, be approximated arbitrarily well by composing flows of vector fields in B. Thus the Lie algebra contains (for small t ≥ 0) all the information about composition of flows.

Until now we have discussed vector fields and flows on a manifold . Yet, whether we wish engage in mathematical analysis or produce software for solving ODEs, it is natural to ask are there important properties of vector fields and flows that can be specified (and possibly programmed) independently of which manifold they are acting upon? This question will lead us to the abstract definition of Lie algebras and Lie groups. The 'glue' that connects an abstract Lie algebra to concrete vector fields on a manifold is called a Lie algebra homomorphism, and abstract Lie groups are connected to flows on a manifold via a Lie group action. To elucidate this state of affairs, let us commence with an important example.

Example 2.5 In Example 2.3 we have shown that there exists natural correspondence between N ×N matrices and linear vector fields defined on R N . We can illustrate this as

A → L A αA → αL A A + B → L A + L B AB -BA → [L A , L B ]
A linear subspace of matrices closed under matrix commutation, [A, B] ≡ AB -BA, is called a matrix Lie algebra. The arrow is an example of a Lie algebra homomorphism, a linear map between two Lie algebras which preserves commutators.

Even the computation of flows of linear vector fields and compositions of such flows can be transformed into linear algebra operations. To achieve this, we must specify how a given N ×N matrix P corresponds to a motion on our domain R N . The simplest possible choice is motions by matrix-vector products. Define thus the map Λ(P, y) = P y.

Identifying a matrix P with the motion Λ(P, • ) leads to the correspondence

P → Λ(P, • ) expm(sA) → Λ( expm(sA), • ) = Ψ s,LA expm(sA) expm(tB) → Λ( expm(sA) expm(tB), • ) = Ψ s,LA •Ψ t,LB .
Note that the latter of these identifications rely on the associative property of the map Λ, Λ(P, Λ(R, y)) = Λ(P R, y).

Motivated by this example, we now proceed to precise mathematical definition of the underlying concepts in a more abstract setting. Definition 2.6 A Lie algebra is a linear space V equipped with a Lie bracket, a bilinear, skew-symmetric mapping

[ • , • ] : V ×V → V
that obeys identities (2.7-10) from Lemma 2.4. Definition 2.7 A Lie algebra homomorphism is a linear map between two Lie algebras, ϕ : g → h, satisfying the identity

ϕ([v, w] g ) = [ϕ(v), ϕ(w)] h , v,w∈ g.
An invertible homomorphism is called an isomorphism.

Definition 2.8 A Lie group is a differential manifold G equipped with a product

• : G ×G → G satisfying p•(q•r) = (p•q)•r ∀ p, q, r ∈ G (associativity), ∃I ∈ G such that I •p = p•I = p ∀ p ∈ G (identity element), ∀p ∈ G ∃ p -1 ∈ G such that p -1 •p = I (inverse),
The maps (p, r) → p•r and p → p -1 are smooth functions (smoothness).

Definition 2.9 An action of a Lie group G on a manifold M is a smooth map

Λ : G ×M → M satisfying Λ(I, y) = y ∀y ∈ M Λ(p, Λ(r, y)) = Λ(pr, y) ∀p, r ∈ G, y ∈ M. (2.14)
If this relation does hold only in a local sense, for all elements p and r sufficiently close to the identity I ∈ G, we say that Λ is local action.

From finite to infinitesimal and back

Given a set of flows on a domain we can find their vector fields by differentiation.

From the discussion above we know that if the flows are closed under composition then the vector fields are closed under Jacobi brackets and linear combinations, hence form a Lie algebra. On the other hand, provided we know the vector fields, we may recover the corresponding flows by integrating differential equations. Similar correspondence between the finite and the infinitesimal is fundamental in abstract Lie theory. We will review here a number of basic results, interpreting them in the terminology of group actions on a manifold. Let G be a Lie group, acting on a manifold M through Λ : G ×M → M and let ρ(t) ∈ G be a curve such that ρ(0) = I, the identity of G. This curve produces a flow Λ(ρ(t), • ) on M and by differentiation we find a vector field

F (y) = d dt Λ(ρ(t), y) t=0 .
The collection of all such vector fields forms a Lie algebra. Note that in order to produce F it is only necessary to know the tangent to ρ(t) at t = 0. Thus the set of all tangents at identity can be endowed with a structure of a Lie algebra.

Definition 2.10 The Lie algebra g of a Lie group G is defined as the linear space of all tangents to G at the identity I. The Lie bracket in g is defined as

[a, b] = ∂ 2 ∂s∂t ρ(s)σ(t)ρ(-s) s=t=0 (2.15)
where ρ(s) and σ(t) are two smooth curves on G such that ρ(0

) = σ(0) = I, ρ ′ (0) = a and σ ′ (0) = b.
Note that the bracket defined in (2.15) is essentially the same as (2.6). From this it is straightforward to verify that the correspondence between elements in g and vector fields on M is an algebra homomorphism.

Lemma 2.6 Let λ * : g → X(M) be defined as

λ * (a)(y) = d ds Λ(ρ(s), y) s=0 , (2.16) where ρ(s) is a curve in G such that ρ(0) = I and ρ ′ (0) = a. Then λ * is a linear map between Lie algebras such that [a, b] g = [λ * (a), λ * (b)] X(M) .
Thus we can go from the finite to the infinitesimal (from groups and group actions to algebras and algebra homomorphisms) by differentiation. To do the opposite and move from the infinitesimal to the finite, we must somehow compute flows of vector fields. A discussion of this process leads to the general definition of the exponential mapping.

Suppose G is a Lie group with Lie algebra g and let Λ : G×M → M be a group action. A given fixed element a ∈ g corresponds to a vector field λ * (a) ∈ X(M). We want to compute the flow of this field. Let us first assume for simplicity that G is a matrix group.

Lemma 2.7 For a fixed A ∈ g the flow of λ * (A), i.e. the solution of y ′ (t) = λ * (A)(y(t)) for y(0) = y 0 ∈ M can be expressed in the form

y(t) = Λ(S(t), y 0 )
where the curve S(t) ∈ G satisfies the matrix differential equation

S ′ (t) = AS, t ≥ 0, S(0) = I,
which has the explicit solution

S(t) = expm(tA), t≥ 0.
Proof. We assume that y(t) = Λ(S(t), y 0 ). Differentiation results in y ′ (t) = ∂ 1 Λ(S ′ (t), y 0 ), where ∂ 1 is the derivative with respect to the first argument. On the other hand, to compute λ * (A) we pick a curve R(s

) ∈ G such that R ′ (0) = A and R(0) = I. From (2.14) we get y ′ (t) = λ * (A)(y(t)) = ∂ ∂s Λ(R(s), Λ(S(t), y 0 )) s=0 = ∂ ∂s Λ(R(s)S(t), y 0 ) s=0 = ∂ 1 Λ(R ′ (0)S(t), y 0 ) = ∂ 1 Λ(AS(t), y 0 ). Thus S ′ (t) = AS(t).
Obviously S(0) = I. The explicit solution is easily verified.

✷

Lemma 2.7 holds unaltered for a general group G if we define the product of an element of an algebra a ∈ g with an element of a group σ ∈ G as

aσ ≡ d ds ρ(s)σ s=0 ,
(2.17) where ρ(s) ∈ G is a smooth curve such that ρ ′ (0) = a and ρ(0) = I. We also define the exponential mapping so that the flow of λ * (a) is of the form Λ(exp(ta), • ):

Definition 2.11 Let G be a Lie group and g its Lie algebra. The exponential mapping exp : g → G is defined as exp(a) = σ(1) where σ(t) ∈ G satisfies the differential equation σ ′ (t) = aσ(t), σ(0) = I.

These definitions lead to the following general form of Lemma 2.7.

Theorem 2.8 Let Λ : G ×M → M be a group action and λ * : g → X(M) the corresponding Lie algebra homomorphism (2.16). For any a ∈ g the flow of the vector field F = λ * (a), i.e. the solution of the equation

y ′ (t) = F (y(t)) = λ * (a)(y(t)), t ≥ 0, y(0) = y 0 ∈ M,
is given as y(t) = Λ(exp(ta), y 0 ).

Let us at this point make a small detour to introduce the adjoint representation which is fundamental in many contexts. By splitting (2.15) into two smaller steps, we obtain: Definition 2.12 Let p ∈ G and let σ(t) be a smooth curve on G such that σ(0) = I and σ ′ (0) = b ∈ g. The adjoint representation is defined as

Ad p (b) = d dt pσ(t)p -1 t=0 (2.18)
The derivative of Ad with respect to the first argument is denoted ad. Let ρ(s) be a smooth curve on G such that ρ(0) = I and ρ ′ (0) = a. Definition 2.10 now yields:

ad a (b) ≡ d ds Ad ρ(s) (b) s=0 = [a, b]. (2.19)
The following formulae show that Ad is both a linear group action (of G on g) and also that for a fixed argument p it is a Lie-algebra isomorphism of g onto itself:

Ad p (a) ∈ g for all p ∈ G, a ∈ g. ( 2 . 2 0 ) Ad p •Ad q = Ad pq (2.21) Ad p (a + b) = Ad p (a) + Ad p (b) ( 2 . 2 2 ) Ad p ([a, b]) = [Ad p (a), Ad p (b)] (2.23)
Note that according to (2.22) both Ad p and ad a are linear in their second argument, hence they may be regarded as matrices acting on the linear space g. This gives meaning to the following important formula relating Ad, ad and the exponential mapping: Ad exp(a) = expm(ad a ).

(2.24)

Differential equations on manifolds

We wish to return to general differential equations on manifolds, as given in (2.1).

In order to construct and implement numerical solvers for this equation, we require a concrete way of representing the vector field F (y). Herewith we describe a very general approach presented in (Munthe-Kaas and Zanna 1997).

Assumption 2.1 Given a differential equation y ′ (t) = F (t, y) on a manifold M, we assume the existence of a Lie algebra g, a Lie algebra homomorphism λ * : g → X(M) and a function a : t×M → g such that the equation can be written in the form y ′ (t) = λ * (a(t, y))(y) ( 2 . 2 5 )

If λ * is known from the context and no confusion is likely, we will usually write equation (2.25) in the shorthand form y ′ (t) = a(t, y)y.

In many important examples the function a depends only on t and not on y. These equations, y ′ (t) = a(t)y, are called equations of Lie type, or linear-type Liegroup equations. Some of the algorithms to be presented in the sequel are aimed at the general equation (2.25), while others are aimed at exploiting the special structure of linear equations.

Given an equation to be solved, an important challenge is to find a 'good' homomorphism λ * . It is not difficult to see that any differential equation can be written in the form (2.25). We might for example let g = X(M) and choose λ * as the identity map, which would trivially render any equation in this form. However, in order to construct practical solution algorithms we need to make some additional assumptions about g. We will usually assume that either all the elements of g, or at least a particular basis of g can be exponentiated efficiently. To achieve this, one might embed M in a linear space R n , and let g be the set of all translations on R n , since translations are trivial to exponentiate. This choise will, however, fail to capture much of the structure of the equations to be solved. In fact we will see that for this choice most of our numerical solution techniques will reduce to classical Runge-Kutta methods. The task of finding a 'good action' is in many respects similar to the task of finding a good preconditioner in the theory of iterative methods for solving linear algebraic equations Ax = b. In both cases we want to find some approximation to our original equation which is both simple to solve and which captures some important structural feature of the equation. The two extreme choices, on the one hand g = X(M) and on the other hand g as a set of all translations on R n , are similar to preconditioning a linear system with the matrix A itself or on the other hand choosing the identity matrix as preconditioner.

Let us now examine briefly a number of examples of equations presented in the form of Assumption 2.1. A useful approach to find a good action is the following:

1 Given a differential equation written in a familiar form, look at the terms and see if it is possible to find related equations which are simpler to integrate. 2 Check that the family of simpler equations forms a Lie algebra. Find a suitable representation g for the algebra and the corresponding homomorphism λ * . 3 Check that the original equation can be written in the form (2.25). This is not possible only in the case where there exist some points on M where the vector fields in the Lie algebra do not generate the direction of the original equation. In this case one must search for a larger Lie algebra.

Example 2.6 (Orthogonal matrix flows) Matrix differential equations of the form (2.26) where A : R + × O(N ) → so(N), are called orthogonal flows -we have already encountered such a flow in (1.3). It is well known that the exact solution Y (t) is an orthogonal matrix for all t ≥ 0 [START_REF] Dieci | Unitary integrators and applications to continuous orthonormalization techniques[END_REF]. A simpler family of equations is given by

Y ′ = A(t, Y )Y, t ≥ 0, Y(0) = Y 0 ∈ O(N ),
Y ′ = CY, t ≥ 0, Y(0) = Y 0 ∈ O(N ), (2.27)
where C is any constant matrix in so(N). The solution of (2.27) is given as

Y (t) = expm(tC)Y 0 , t ≥ 0.
From (2.11) we find that if

F (Y ) = CY and G(Y ) = DY then [F, G](Y ) = (CD - DC)Y .
Thus, the family of simple vector fields is a Lie algebra isomorphic to so(N).

We have

λ * (C)(Y ) = CY , hence (2.26) is in the form (2.25) if a(t, Y ) = A(t, Y ).
Note that also the flow of (2.27) is orthogonal, hence any numerical method based on the composition of such flows will yield a solution which retains orthogonality.

Example 2.7 (Isospectral matrix flows) We have already encountered such flows in Section 1, in (1.1). With slightly greater generality, we write them in the form (2.28) where B : R + × S N → so(N). The analytical solution Y (t) is a family of matrices with eigenvalues invariant under the flow. A simpler family of equations is given by

Y ′ = B(t, Y )Y -Y B(t, Y ), t ≥ 0, Y(0) = Y 0 ∈ S N ,
Y ′ (t) = CY -Y C, t ≥ 0, Y(0) = Y 0 ∈ S N , (2.29)
where again C ∈ so(N) is constant. The solution of this equation is given explicitly

in the form Y (t) = expm(tC)Y 0 expm(-tC), t ≥ 0.
(2.30)

We may now proceed, exactly like in the previous example, to find the brackets of such vector fields. An alternative route, already anticipated in (1.2), is to note that the basic flows in (2.30) are given by the orthogonal matrix group O(N ) acting on Y ∈ S N . We have the action Λ : O(N ) × S N → S N given by

Λ(Q, X) = QXQ T for any Q ∈ O(N ).
By differentiation, similarly to Lemma 2.6, we find that λ * (C)(Y ) = CY -Y C. Hence (2.28) is in the form (2.25) if a(t, Y ) = B(t, Y ). The basic flow (2.30) is isospectral since it is a similarity transformation, hence also any numerical method based on this flow is automatically isospectral.

Note that, essentially, these two examples are identical, except for the action involved. Hence in a computer implementation we might reuse the implementation of the Lie algebra and the Lie group, while simply changing the action: there is no need to develop separate computational approach to orthogonal and isospectral flows! This illustrates the importance of working with abstractly-defined groups and algebras rather than tying these concepts to flows and vector fields on particular manifolds.

Example 2.8 (ODEs on R N ) Consider an ordinary differential equation in the familiar form required by all classical numerical integrators,

y ′ (t) = g(t, y), t ≥ 0, y(0) = y 0 ∈ R N , (2.31)
where g : R × R N → R N . A simple family of basic flows is given by translations, e.g.

z ′ = a,
where a ∈ R N is constant. The algebra of these flows can be identified with R N and, since translations commute, we obtain the trivial bracket [a, b] = 0 for all a, b ∈ R N . The algebra homomorphism is given here simply as the identity map, hence the equation is in the form (2.25) if a(t, y) = g(t, y). This choice yields nothing new compared to classical integration schemes.

A more interesting choice is choosing linear vector fields as basic flows,

z ′ = Az.
We already know from Example (2.3) that the Lie algebra of these is gl(N), the set of all N ×N matrices with the usual matrix bracket. However, these flows cannot produce everywhere all possible tangent directions. Indeed, at y = 0 these flows cannot produce any nonzero tangent. Therefore, using this action, it is generally not possible to write (2.31) in the form (2.25).

We might finally consider a combination of translations and linear maps, i.e. the set of all affine linear maps given by equations of the form

z ′ = Az + b for A ∈ gl(N), b ∈ R N .
The Lie algebra can be identified with all pairs (A, b) where A is a matrix and b a vector and we have already seen in Example (2.4) that the bracket is given as

[(A, b), (C, d)] = (AC -CA, Ad -Cb). The algebra homomorphism is λ * (A, b)(y) = Ay + b.
In this situation there are many possible choices of a function a(t, y) such that (2.31) acquires the form (2.25).

No matter what we pick as the matrix part, it is always possible to adjust the vector so that λ * (a(t, y))(y) = g(t, y).

A natural choice is local linearization, namely letting the matrix part be the Jacobian of g at y, Jg(t, y) i,j = ∂g i (t, y) ∂y j .

The resulting a : R × M → g is a(t, y) = (Jg(t, y), g(t, y) -Jg(t, y)y) .

Possible advantages of using affine motions to advance the solution of stiff ODEs are a subject of ongoing research.

Example 2.9. (ODEs on a sphere) Many mechanical problems involve rotations in a 3-space. In Appendix B we list useful formulae for fast computations in SO(3). As a simple example we will consider a motion on the surface of a sphere.

y ′ (t) = a(y(t)) × y(t), (2.32)
where a, y ∈ R 3 . If ∥y(0)∥ 2 = 1 then y(t) evolves on the unit sphere S 2 ⊂ R 3 . A simpler system that also evolves on S 2 is given by

y ′ (t) = c × y(t), (2.33)
for any fixed vector c. To compute the commutator of two equations of the form (2.33), we employ the hat map (B.1) taking a 3-vector c to a 3 × 3 matrix c, such that c×y = cy. Using (B.2) we see that the commutator of the vector fields c×y and d×y is given by (c×d)×y. Thus the Lie algebra in this example may be identified with (R 3 , ×), the real 3-space with the Lie bracket

[c, d] = c×d
given by a vector product. The simplified equation (2.33) is a matrix equation y ′ = cy with solution y(t) = expm(t c)y 0 ,

where expm(t c) can be computed fast using the Rodrigues formula (B.10).

Example 2.10. (Parabolic PDEs) The final example in this section is chosen to illustrate the diversity of problems that may be tackled with the machinery of Lie-group methods. This example involves infinite-dimensional Lie algebras, a topic that is technically more demanding than the finite-dimensional case. There are several ways to circumvent the mathematical problems, we may either discuss the problem after it has been discretised in space (and has become finite dimensional), or we may plunge straight ahead using the available techniques, disregarding possible mathematical difficulties. We will henceforth follow the latter approach, being aware that resulting algorithms must be verified by other means.

Suppose that we wish to integrate a parabolic PDE with coefficients varying in space and time, e.g. the heat equation

∂u(t, x) ∂t = ∇ • (µ(x)∇u(t, x)) ,
where u is the temperature and µ the heat conductivity of the material. With greater generality, consider equations of the form

∂u(t, x) ∂t = L(u), (2.34)
where L is an elliptic operator. To simplify the discussion, we wish to neglect boundary conditions, so suppose that u is defined on the unit square with periodic boundary conditions. Thus, u should be thought of as a 'point' on the infinitedimensional manifold M = C ∞ (T), the collection of all smooth functions on a torus. It is well known that (classical) explicit integrators for parabolic PDEs are typically stiff and stability analysis leads to severe step size restriction: ∆t < c(∆x) 2 for explicit finite-difference methods. A family of simpler equations, which can be solved exactly, explicitly and very efficiently with the Fast Fourier Transform (FFT) is the set of all parabolic equations with constant coefficients of the form

∂u(t, x) ∂t = μ∇ 2 u(t, x),
where μ is constant. However, just like in Example 2.8, these equations cannot move an arbitrary point u ∈ M in an arbitrary direction. (In other words, they do not define a transitive action on M). Hence, we enlarge the family of simplified equations by adding an inhomogeneous term,

∂u(t, x) ∂t = μ∇ 2 u(t, x) + b(x), (2.35)
where b ∈ C ∞ (T). Also this equation is easy to solve using FFTs. Letting the flows of (2.35) define our group action on M, we see that the corresponding Lie algebra g can be identified with pairs (μ, b). The Lie-algebra action is given by

λ * ((μ, b))(u) = μ∇ 2 u + b.
Using this action we see that any equation of the form (2.34) can be cast into the form (2.25) by choosing the function a : R×M → g as

a(u) = (μ, L(u) -μ∇ 2 u) ( 2 . 3 6 )
for some choice of μ. For example, if L(u) = ∇ • (µ(x)∇u(t, x)), we would let μ be some averaged value of µ(x).

In order to define the whole structure of g, we need to determine the Lie bracket. Let F and G be two vector fields on M defined at a point u as

F (u) = μ∇ 2 u + f (x) G(u) = ν∇ 2 u + g(x).
Using (2.12), we obtain

[F, G](u) = ∂ ∂s μ∇ 2 (u + s(ν∇ 2 u + g)) + f -ν∇ 2 (u + s(μ∇ 2 u + f )) -g s=0 = 0 • ∇ 2 u + μ∇ 2 g -ν∇ 2 f.
Thus, the bracket on g = R×C ∞ (T ) is

[(μ, f ), (ν, g)] = (0, μ∇ 2 g -ν∇ 2 f ).
It is interesting to note that also the bracket can be computed efficiently using FFTs. This type of equations is considered in greater detail and various numerical examples are given in (Munthe-Kaas and Lodden n.d.). It turns out that, at least in some cases, it is possible to construct explicit integrators based on this action that are not subjected to any step size restriction involving the spatial discretization ∆x. Thus, the methods are stable regardless of how spatial discretization is chosen. This is a topic of ongoing research, and many aspects of these integrators are currently incompletely understood.

We will return to more examples and to numerical experiments in Section 11. Before discussing numerical algorithms, we need to study some important properties of the exponential map.

Much Ado about something

In this section we have emphasised a general view of differential equations on manifolds, based on Lie groups acting on manifolds. This outlook is important not just for the sake of mathematical beauty or abstraction but, as we hope to have persuaded the reader, also from the point of view of applications and computation. However, insofar as clarity of exposition is concerned, it is often better to restrict ourselves to the far simpler, familiar and more intuitive matrix theory.

In fact, it turns out that for all the algorithms that we present in this paper it is quite straightforward to translate results derived in matrix setting to the more general setting of local Lie-group actions on some domain. The following theorem of Ado underscores the importance of studying the matrix case [START_REF] Olver | Equivalence, Invariants, and Symmetry[END_REF][START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF].

Theorem 2.9. (Ado's theorem) Any finite-dimensional Lie algebra is isomorphic to a subalgebra of the matrix algebra gl(N ) for some N ≥ 1.

Although similar result does not hold for all finite-dimensional Lie groups, it is true that whenever we are given a finite-dimensional local Lie group action, we can always find an equivalent local action by a matrix Lie group. More information on these topics can be found in [START_REF] Olver | Equivalence, Invariants, and Symmetry[END_REF] and [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF].

Aware of the danger of rules-of-a-thumb being mathematically imprecise, it is nonetheless worthwhile to summarise these results as follows.

For practically any concept in general Lie theory there exists a corresponding concept within matrix Lie theory. Vice versa, practically any result which holds in the matrix case remains valid within the general Lie theory.

The above remark is even more important in numerical context, since computation always takes place in a finite-dimensional setting. Even if the original equation evolves on an infinite-dimensional manifold, its practical computation must sooneror-later involve a discretization to a finite-dimensional formulation.

At this point it is time to wake up the readers who have surfed through the general theory in a relaxed manner. We will restate the main definitions in the concrete form in which they appear within the matrix theory. It is worthwhile to compare these definitions to the corresponding general definitions above. Definition 2.13 A real matrix Lie group is a smooth subset G ⊆ R n×n , closed under matrix products and matrix inversion. We let I ∈ G denote the identity matrix.

Definition 2.14 The Lie algebra g of a matrix Lie group G is the linear subspace g ⊆ R n×n consisting of all matrices of the form

g = A ∈ R n×n : A = dρ(s) ds s=0 ,
where ρ(s) ∈ G is a smooth curve such that ρ(0) = I. The space g is closed under matrix additions, scalar multiplication and the matrix commutator

[A, B] = AB -BA.
(2.37)

Complex matrix Lie groups and algebras are defined similarly.

The time has come to introduce some of the main dramatis personae of our survey, concrete examples of Lie groups and algebras. In each case it is easy to verify that all the axioms of a group or an algebra, as the case might be, are fulfilled, and we leave this as an exercise to the reader.

• The set of all real N × N nonsingular matrices is a (multiplicative) Lie group, the general linear group GL(N ). The corresponding Lie algebra is the set R N ×N of all N × N real matrices which, in keeping with our terminology, we denote by gl(N).

The general linear group and algebra can be defined over other fields than R, in which case we communicate this in the second argument. For example, GL(N ; C) consists of all nonsingular N × N complex matrices.

• All members of GL(N ) with unit determinant form the special linear group SL(N ). Its Lie algebra, sl(N), consists of all matrices in gl(N) with zero trace.

• N × N real orthogonal matrices form the orthogonal group O(N ), whose Lie algebra so(N) consists of N × N skew-symmetric matrices.

The set SO(N ) = SL(N )∩O(N ), consisting of N ×N real orthogonal matrices with unit determinant, is the special orthogonal group. Its Lie algebra is so(N), whch we have just encountered. This is not contradictory: we never claimed that two different Lie groups must have different Lie algebras! As a matter of fact, more is true: If G is a Lie group and G Id its connected component such that I ∈ G Id (precisely the situation with O(N ) and SO(N ), respectively) then they produce the same Lie algebra.

• The set of all (2N ) × (2N ) real matrices X such that XJX T = J, where

J = O N I N -I N O N ,
is the symplectic group and denoted by Sp(N ). (The Jacobian of the flow of a Hamiltonian ODE system evolves in Sp(N ).) The corresponding Lie algebra, sp(N), consists of F ∈ gl(2N) such that F J + JF T = O.

• All the matrices X ∈ SL(4) such that XJX T = J, where

J = ⎡ ⎢ ⎢ ⎣ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 -1 ⎤ ⎥ ⎥ ⎦ ,
form the Lorenz group SO(3, 1). Its Lie algebra so(3, 1) is made out of all F ∈ gl(4) such that F J + JF T = O.

• As an example of complex Lie groups, we mention the unitary group U(N ; C) of all N × N complex unitary matrices: X ∈ U(N ; C) if and only if XX H = I.

The Lie algebra corresponding to U(N ; C) is the set u(N; C) of all skew-Hermitian matrices in gl(N; C).

The unitary group should not be confused with O(N ; C), the group of all N × N complex orthogonal matrices, whose Lie algebra is so(N; C).

• Similarly to the case of O(N ) and SO(N ), we obtain the special unitary group intersecting U(N ; C) with SL(N ; C). Its Lie algebra, su(N; C), is made out of N × N complex skew-Hermitian matrices with zero trace.

Definition 2.15 A differential equation on a matrix Lie group is an equation of the form (2.38) where A : R×G → g and AY is the usual matrix product between A ∈ g and Y ∈ G.

Y ′ = A(t, Y )Y, t ≥ 0, Y(0) ∈ G,
The reader may verify that this is the special case of the general form of a differential equation on a manifold, given in Assumption 2.1, where M = G is a matrix Lie group and the action Λ is taken to be the left (matrix) multiplication in G, Λ(R, Y ) = RY.

Using (2.16) we find

λ * (A)(Y ) = AY.
Since g is defined as the collection of all tangent directions at I ∈ G and matrix multiplication by Y is an invertible mapping, we see that any tangent at Y can be written in the form AY and all differential equations on G can be written in the form (2.38).

Definition 2.16 The exponential mapping expm : g → G is defined as

expm(A) = ∞ j=0 A j j! .
(2.39)

Note that expm(O) = I, and that for A sufficiently near O ∈ g the exponential has a smooth inverse given by the matrix logarithm logm : G → g.

Definition 2.17 The adjoint representation, Ad, and its derivative, ad, are defined as

Ad P (A) = P AP -1 (2.40) ad A (B) = AB -BA = [A, B].
(2.41)

AY , A ∈ g, Y ∈ G λ * (a)(y), a ∈ g, y ∈ M Y ′ = A(Y, t)Y y ′ = λ * (a(y, t))(y) RY , R, Y ∈ G Λ(r, y), r ∈ G, y ∈ M expm(A) = ∞ j=0 A j /j! Definition 2.11 [A, B] = AB -BA Definition 2.10 P AP -1 , P ∈ G, A ∈ g Definition 2.12
It is easy to verify that (2.20-23) hold in the matrix case, while (2.24) is far from being an obvious identity even in that case.

Table 2.1 summarises the correspondence between the matrix case and the general case.

The differential of the exponential map

We have introduced the exponential mapping exp : g → G as the fundamental solution of the equation y ′ = ay, or more explicitly in the matrix case as exp A = ∞ j=0 A j /j!. For development of numerical algorithms it is essential to discuss the derivative of the exponential map. This will first be used to deduce an infinitesimal version of the BCH formula and in the following chapters to derive a variety of different numerical algorithms for differential equations on Lie groups.

To simplify the exposition we will restrict the proofs to matrix theory. The results are, however, valid in an abstract setting. The reader is referred to [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF] for proofs in a general context.

Given a scalar function a(t) ∈ R, the derivative of the exponential is given as d exp(a(t))/dt = a ′ (t) exp(a(t)). One might have hoped for a similar result also when A(t) is a matrix. However, due to the fact that in general [A, A ′ ] ̸ = O, this is not the case and we must correct this formula. Note that d exp(A(t))/dt must be tangent to G in the point P (t) = exp(A(t)). We have seen in Section 2.5 that any such tangent can be written as C(t)P (t), where C(t) ∈ g. Furthermore, general properties of d/dt imply that B(t) must depend only on A(t) and A ′ (t), and that the dependence on A ′ is linear. This function is denoted by dexp.

Definition 2.18 The differential of the exponential mapping is defined as the 'right trivialised' tangent of the exponential map, i.e. as a function dexp :

g×g → g such that d dt exp(A(t)) = dexp A(t) (A ′ (t)) exp(A(t)).
(2.42)

Just like the functions Ad A and ad A defined in (2.40) and (2.41) respectively, also dexp A is linear in its second argument for a fixed A. Hence we may regard all these as being matrices acting on g. In fact dexp A is an analytic function of the matrix transformation ad A :

dexp A = expm(ad A ) -I ad A . (2.43)
This formula should be read as a power series in the following manner: Since

e x -1 x = 1 + 1 2! x + 1 3! x 2 + 1 4! x 3 • • • + 1 (j+1)! x j + • • • , we obtain dexp A (C) = C + 1 2! [A, C] + 1 3! [A, [A, C]] + 1 4! [A, [A, [A, C]]] + • • • = ∞ j=0 1 (j + 1)! ad j A C.
(2.44)

The fact that dexp A is an analytic function in ad a makes it easy to invert the matrix dexp A simply by inverting the analytic function,

dexp -1 A = ad A expm(ad A ) -I . (2.45) Recall that x e x -1 = 1 -1 2 x + 1 12 x 2 -1 720 x 4 + • • • = ∞ j=0 B j j! x j ,
where B j are Bernoulli numbers [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. Thus,

dexp -1 A (C) = C -1 2 [A, C] + 1 12 [A, [A, C]] + • • • = ∞ j=0 B j j! ad j A (C).
(2.46)

Note that, except for B 1 , all odd-indexed Bernoulli numbers vanish. Hence dexp -1 A + 1 2 ad A is an even function of ad A . We have based the formulae here on right trivialisations, i.e. tangents at a point P ∈ G being written as CP , C ∈ g. It is equally possible to derive formulae based on left trivialisations, tangents written in the form CP . If P C = CP , we observe that C = P -1 CP = Ad P -1 (C). Using (2.24), we compute the left-trivialised formulae as

Ad exp(-A) dexp A = exp(ad -A ) exp(ad A ) -I ad A = I -exp(ad -A ) ad A = dexp -A . Hence d dt exp(A(t)) = dexp A(t) (A ′ (t)) exp(A(t)) = exp(A(t)) dexp -A(t) (A ′ (t)). (2.47)
Thus we can arrive at the left versions by changing the sign of every commutator. Note that dexp -1 A (C) and dexp -1 -A (C) differ only in the sign of the term

± 1 2 [A, C].
The definition of dexp in Definition 2.18 can be generalised to any smooth function ψ : g → G:

Definition 2.19 Given a smooth function ψ : g → G, we define the right trivialized tangent of ψ as the function dψ : g×g → g defined such that

d dt ψ(A(t)) = dψ A(t) (A ′ (t))ψ(A(t)).
(2.48)

The function dψ is always linear in the second argument, A ′ .

Let us now apply dexp and dexp -1 to obtain a differential equation for the BCH formula. Going back to Theorem 2.5, we define a function bch : g×g → g such that expm( bch AB) = expm(A) expm(B).

We may compute C = bch(A, B) by integrating a differential equation. Let C(t) = bch(tA, B). Clearly C(0) = B, and we seek C(1). Writing expm(C(t)) = expm(tA) expm(B), we find by differentiation that

dexp Ct (C ′ (t)) exp(C(t)) = dexp tA (A) exp(tA) exp(B) = A exp(C(t)) ⇒ C ′ (t) = dexp -1 C(t) (A).
We have proved: the following result.

Lemma 2.10 The function C = bch(A, B) can be computed by integrating the differential equation

C ′ (t) = dexp -1 C(t) (A), 0 ≤ t ≤ 1, C(0) = B, (2.49) from t = 0 to t = 1.
In this light the dexp -1 function may be regarded as a kind of 'infinitesimal BCH generator'. A more concise presentation of this idea is given in [START_REF] Engø | On the construction of geometric integrators in the RKMK class[END_REF]. A symbolic algorithm to compute the BCH from this formula can be found in (Munthe-Kaas and Owren 1999).

Crouch-Grossman methods

The discipline of Lie-group methods owes a great deal to the pioneering work of Peter Crouch and his coworkers, who were the first to introduce in a systematic, mathematically-sophisticated manner ODE methods that evolve on manifolds. It is interesting to note that their work was primarily motivated by problems in robotics and control theory.

The main algorithm originating in this circle of ideas is the method of rigid frames of [START_REF] Crouch | Numerical integration of ordinary differential equations on manifolds[END_REF]. It has been originally stated in a more general formalism of differential equations evolving on manifolds. To fit the method into our narrative and to simplify its exposition we restrict our discussion to Lie-group equations.

In essence, the Crouch-Grossman approach is an attempt to apply a Runge-Kutta method to (2.1) (or, in a Lie-group context, to (2.38)) by repeatedly freezing and thawing coefficients and keeping the flow in the correct configuration space. The solution of a 'frozen' Lie-group equation,

Y ′ = A( t, Ỹ )Y t ≥ t * , Y(t * ) = Y * , is simply expm((t -t * )A( t, Ỹ ))Y * . Freezing (2.38) at t n and letting Y n+1 = dexp(hA(t n , Y n ))Y n , where Y m ≈ Y (t m ), t m = mh
, results in a first-order method of very little merit. This can be remedied in the following procedure. We choose constants c k , b l , a k,l , 1 ≤ l < k ≤ ν and let

X k = e ha k,k-1 F k-1 e ha k,k-2 F k-2 • • • e ha k,1 F1 Y n F k = A(t n + c k h, X k ), ⎫ ⎬ ⎭ k = 1, 2, . . . , ν, Y n+1 = e hbν Fν e hbν-1Fν-1 • • • e hb1F1 Y n .
(2.50)

In other words, we model the solution as a product of ν 'frozen' steps. Note that

X k ∈ G, F k ∈ g, hence, as required, Y n+1 ∈ G.
The method (2.50) might appear to the initiated as a 'Lie-group version' of a Runge-Kutta scheme, an analogy that we have deliberately reinforced by employing notation that will be reserved in the sequel to RK schemes. Yet, order conditions are considerably more challenging than in the classical RK case: they are nonlinear in the weights b 1 , b 2 , . . . , b ν and there are more of them! Moderate headway can be made by elementary means and a great deal of algebra. Thus, [START_REF] Crouch | Numerical integration of ordinary differential equations on manifolds[END_REF] have derived three-stage methods of order three, for example Inquiry into higher-order methods, though, requires more than algebra and elbow grease. The situation is further complicated by the fact that (2.50) is typically formulated in considerably more abstract manner, in a manifold setting: this does not make order analysis any simpler! The order of classical RK methods is nowadays determined by a method due to Butcher (1963), which identifies expansion terms of both the exact and the approximate solution with rooted trees. Remarkably, similar approach can be generalised to Crouch-Grossman methods and this has been accomplished by Owren and Marthinsen (1999b). Details of their work are outside the scope of this survey and we refer the readers to the primary source. Let us just mention that, unlike the classical RK case, there are no fourth-order methods of this kind with four stages and ν = 5 is required. Moreover, Owren and Marthinsen (1999b) extended (2.50) to implicit methods, whereby a k,l is given for all k, l = 1, 2, . . . , ν and

X 1 = Y n , F 1 = A(t n , X 1 ), X 2 = e 3 4 hF1 Y n , F 2 = A(t n + 3 4 h, X 2 ), X 3 = e
X k = e ha k,ν Fν e ha k,ν-1 Fν-1 • • • e ha k,1 F1 Y n , k= 1, 2, . . . , ν.
This allows for better order/stages ratio but the downside is the need to solve nonlinear equations in Lie groups. Classical methods, e.g. the Newton-Raphson technique, are of little use here since, unless iterated to convergence, they are unlikely to deliver a solution which resides in the Lie group G. Recently, however, [START_REF] Owren | The Newton iteration on Lie groups[END_REF] developed two variants of Newton's method that always evolve in G. This brings implicit versions of (2.50) within the realm of computation, although they are expensive.

The mainstream of research into Lie-group methods has moved in the last few years away from the Crouch-Grossman approach. The main reason is that the RK-MK methods, the theme of the next section, provide a considerably more convenient, intuitive and easy-to-analyse means of translating Runge-Kutta formalism to a Lie-group setting. Yet, it would be unfair to pronounce Crouch-Grossman methods as unviable or of purely historical interest. Firstly, in their more general setting, Crouch-Grossman methods can be made (with some effort!) to evolve on arbitrary smooth manifolds, while the scope of RK-MK is restricted to group actions. Secondly, in the present stage in the lifetime of geometric integration we are denied the comfort of discarding lines of inquiry simply because of our current, incomplete understanding.

Runge-Kutta on manifolds and RK-MK

In this section we describe a class of numerical integration schemes for computing (2.38), or more generally (2.25). The methods have become known under the name of RK-MK-type schemes. We will later see that they might just as well be called integration schemes based on canonical coordinates of the first kind. These methods were originally developed in a series of four papers, (Munthe-Kaas 1995, Munthe-Kaas and Zanna 1997, Munthe-Kaas 1998) and (Munthe-Kaas 1999).

A major motivation behind the first of these papers was an attempt to understand and specify the basic operations underlying classical Runge-Kutta methods for integration of differential equations. Abstract specifications of mathematical structures are fundamental in theoretical computer science as a tool for structuring software. An object-oriented program consists of a collection of program modules which interact in a well-specified manner. A module could e.g. represent some mathematical structure, like a linear space, a Lie algebra or a Lie group. The basic idea of object-oriented programming is that particular representations of the mathematical structure to be modelled should be hidden within the program module, and that interactions between different program modules should be independent of particular representations. Although this approach to programming has been very successful for discrete problems, considerably less has been done within areas of continuous mathematics, such as integration of differential equations. Much insight about the important underlying structures can be gained by studying 'pure' mathematical definitions, since these focus more on what the essential mathematical structures are, rather than on how they can be represented. Thus it is natural, for example, to specify that domains of differential equations should be differential manifolds and that the right hand side of the equation should be a vector field on the manifold.

Seen from this perspective, classical integration schemes such as Runge-Kutta methods contain 'type errors' in their formulation. As an example, consider the trapezoidal rule

y n+1 = y n + h 2 [f (y n ) + f (y n+1 )].
All the operations are valid if y and f are vectors. However, if y is a point on a manifold and f a vector field, then this expression involves addition of tangent vectors at different base points, and also the addition of a point on a manifold and a tangent to the manifold, both being invalid operations in the context of general manifolds.

In (Munthe-Kaas 1995) classical Runge-Kutta methods are reformulated using coordinate-independent operations on a Lie group. It is shown there that the Butcher theory for order conditions of Runge-Kutta methods (Butcher 1963, Hairer, Nørsett and[START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF] can be reformulated in a geometrical language, replacing the 'Butcher trees' with commutators in a Lie algebra. The outcome is a so-called Lie-Butcher theory. Although the resulting algorithms respect Lie-group structure, they can, in the simplest version, attain at most order two for a general non-commutative Lie group. In the sequel, (Munthe-Kaas 1998), the Lie-Butcher order theory is improved, order conditions derived to arbitrary order in general Lie groups and explicit methods of Runge-Kutta type presented up to order four. The paper (Munthe-Kaas and Zanna 1997) generalises the theory from equations on Lie groups to equations evolving on more general manifolds acted upon by a Lie group. In the final of these papers, (Munthe-Kaas 1999), it is shown that similar methods of arbitrarily high order can be constructed and analysed in a relatively simple manner, without employing the Lie-Butcher theory.

In this section our goal is to arrive at the main ideas of the algorithms while employing a minimal amount of formal theory. We have therefore decided not to discuss the general Lie-Butcher theory since this would require a significant amount of Lie theory, beyond what we have already introduced in Section 2. The interested reader is referred to (Munthe-Kaas 1998) for details on Lie-Butcher series. We will continue along the lines of (Munthe-Kaas 1999), but instead of proofs for the general equation (2.25) we restrict the discussion mainly to the simple matrix case (2.38).

We have seen that classical Runge-Kutta methods are valid and 'type-correct' only for differential equations evolving in a linear space V , since then the configuration space and the space of vector fields coincide. If the analytical solution of the differential equation evolves on some linear subspace W ⊂ V , then it is easy to show that a consistent numerical integrator will also evolve on W (up to a departure due to roundoff errors). On the other hand, if the differential equation evolves on some nonlinear manifold embedded in V , it is much more difficult to devise numerical algorithms that stay on the right submanifold. It is well known that traditional ν-stage Runge-Kutta methods preserve quadratic submanifolds if the coefficients satisfy the condition [START_REF] Cooper | Stability of Runge-Kutta methods for trajectory problems[END_REF], while [START_REF] Calvo | Numerical solution of isospectral flows[END_REF] show that this condition is also necessary. In the same paper it is also shown that it is in general impossible to devise classical Runge-Kutta methods that preserve arbitrary cubic submanifolds. Linear multistep methods or truncated Taylor expansions cannot in general preserve even quadratic submanifolds [START_REF] Iserles | Multistep methods on manifolds[END_REF]).

b k b l = b k a k,l + b l a l,k , k,l= 1, 2, . . . , ν (
In the case of equations on a Lie group G, recall that the local structure in a neighbourhood of any point can be described by the Lie algebra g, which is a linear space. Even better, if H is a Lie subgroup of G, then there exists a (linear!) subalgebra h of g describing the local structure of H. Given a differential equation evolving on H, it is in general impossible to devise classical integration scheme that will evolve on H. On the other hand, if an equation is evolving on h, so will almost any reasonable numerical integrator. It thus seems to be a good idea to try to solve a differential equation in the Lie algebra rather than in the Lie group! Given the equation

Y ′ = A(t, Y )Y , Y (0) = Y 0 , we call the map g ∋ A → expm(A)Y 0 ∈ G (3.1)
canonical coordinates of the first kind [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF]. This defines a smooth invertible map between a neighbourhood of O ∈ g and a neighbourhood of Y 0 ∈ G. We say that these coordinates are centred about Y 0 ∈ G. A crucial step is 'pulling back' the equation from G to g using this map.

Lemma 3.1 For small t ≥ 0 the solution of (2.38) is given by

Y (t) = expm(Θ(t))Y 0 .
where Θ ∈ g satisfies the differential equation

Θ ′ (t) = dexp -1 Θ(t) (A(t, Y )) Θ(0) = O (3.2)
and the dexp -1 operator has been defined in (2.46).

Proof. Differentiation of Y (t) = exp(Θ(t))Y 0 yields Y ′ (t) = dexp Θ(t) (Θ ′ (t)) exp(Θ(t))Y 0 = dexp Θ(t) (Θ ′ (t))Y (t).
The lemma follows from

Y ′ (t) = A(t, Y )Y (t). ✷ Equation (3.
2) is absolutely crucial to the entire business of Lie-group methods. It has been originally stated by Felix [START_REF] Hausdorff | Die symbolische Exponentialformel in der Gruppentheorie[END_REF], although some attribute it to John Edward Campbell, who might have published it few years earlier. (The names of both, together with Henry Frederick Baker, have been immortalised in the BCH formula, cf. Theorem 2.5.) The corresponding result for the general case (2.25) is given in [START_REF] Munthe-Kaas | High order Runge-Kutta methods on manifolds[END_REF].

Note that the proof uses no other property of the exponential mapping than the definition of dexp. Hence the argument can be easily generalised, replacing exp with a general coordinate map ψ and dexp with dψ as discussed in Definition 2.19. This will be used in Section 6.

The simplest and most natural solution strategy is to apply a classically-formulated Runge-Kutta method to (3.2), rather than to the original equation (2.38). At each step we choose a coordinate system of the form (3.1), centred at the last known point Y n . Let us consider briefly the details of this algorithm. Recall that a ν-stage Runge-Kutta method is defined by constants [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]. Applied to a standard vector equation y ′ = f (t, y), a single step of length h from y n to y n+1 is given by

{a k,l } ν k,l=1 , {b l } ν l=1 , {c k } ν k=1 , usually written as a Butcher tableau c 1 a 1,1 a 1,2 • • • a 1,ν c 2 a 2,1 a 2,2 • • • a 2,ν . . . . . . . . . . . . c ν a ν,1 a ν,2 • • • a ν,ν b 1 b 2 • • • b ν (
θ k = y n + ν l=1 a k,l f l , f k = hf (t n + c k h, θ k ), ⎫ ⎪ ⎬ ⎪ ⎭ k = 1, . . . , ν, y n+1 = y n + ν l=1 b l f l .
(3.3)

Applying this scheme to (3.2), we obtain the RK-MK algorithm. The following equations describes a single RK-MK step from

Y n ∈ G to Y n+1 ∈ G: Θ k = ν l=1 a k,l F l , A k = hA(t n + c k h, expm(Θ k )Y n ), F k = dexp -1 Θ k (A k ), ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ k = 1, . . . , ν, Θ = ν l=1 b l F l , Y n+1 = expm(Θ)Y n .
(3.4)

The same algorithm also integrates the general equation (2.28), provided we replace expm(Θ)Y n with its general form Λ(exp(Θ), Y n ).

In order to complete this algorithm, we need to provide practical means for computing dexp -1 Θ k (A k ). In some cases there exists fast direct algorithms for this, see Appendix B. Note that even if dexp -1 is approximated, the resulting algorithm will evolve on the correct manifold. In general one may use the expansion (2.46), truncated to the order of the underlying Runge-Kutta scheme, and the resulting algorithm will obtain the same order as the underlying Runge-Kutta scheme, while staying on the correct manifold. For high order methods, a significant number of commutators must be computed if dexp -1 is computed using (2.46). In Section 5.3 the structure of so called free Lie algebras is used to dramatically reduce the number of commutators.

Examples of specific methods based on (3.4) feature in Appendix A, where, in a more general setting, it is redesignated as (A.1). Here we just stress again the main difference between (3.3) and (3.4): the latter acts in the Lie algebra g, which is a linear space, thereby respecting Lie-group structure.

Magnus and Fer expansions

There are several possible points of departure for our description of Magnus and Fer expansions. Perhaps the simplest and the most intuitive is the scalar linear differential equation

y ′ = a(t)y, t ≥ 0, y(0) = y 0 .
Its solution, y(t) = exp t 0 a(ξ)dξ y 0 , is familiar to all well-trained mathematics undergraduates. Bearing in mind that the definition of the exponential can be easily extended from scalars to matrices, one might have perhaps hoped that its obvious generalisation,

expm t 0 A(ξ)dξ Y 0 , (4.1)
is the solution of the matrix linear system

Y ′ = A(t)Y, t ≥ 0, Y(0) = Y 0 . (4.2)
Unless A(t 1 ) and A(t 2 ) commute with each other for all t 1 , t 2 ≥ 0, this is, unfortunately, misplaced hope. Before offering possible remedies to this state of affairs, we mention that if Y 0 ∈ G, a Lie group, and A lies in its Lie algebra g (whence (4.2) is a Lie-group equation) then (4.1) evolves in G. Bearing in mind the advisability of respecting Lie-group structure, we need to 'correct' (4.1) without losing this feature. Two possible remedies suggest themselves. Firstly, we may seek a correction ∆(t)

evolving in the Lie algebra g so that

Y (t) = expm t 0 A(ξ)dξ + ∆(t) Y 0 .
Alternatively, we may correct with V (t) in the Lie group G,

Y (t) = expm t 0 A(ξ)dξ V (t).
This gives rise to Magnus and Fer expansions, respectively. Both [START_REF] Magnus | On the exponential solution of differential equations for a linear operator[END_REF] and [START_REF] Fer | Résolution del l'equation matricielle U = pU par produit infini d'exponentielles matricielles[END_REF] expansions originated within a non-numerical context and they have found extensive use, e.g. in mathematical physics, quantum chemistry, control theory and stochastic differential equations as a perturbative tool in the investigation of linear systems (4.2). Fashioning them into a numerical weapon is nontrivial and will occupy us in this and next sections.

Magnus expansions and rooted trees

Our point of departure is the dexpinv equation (2.46) which we recall for completeness of exposition: the solution of (4.2) can be written in the form Y (t) = expm(Θ(t))Y 0 , t ≥ 0, where

Θ ′ = dexp -1 Θ A = ∞ k=0 B k k! ad k Θ A, t ≥ 0, Θ(0) = O, (4.3) {B k } k∈Z + being Bernoulli numbers.
As a first step, we attempt to solve (4.3) by Picard iteration,

Θ [0] (t) ≡ O, Θ [m+1] (t) = t 0 dexp -1 Θ [m] (ξ) A(ξ)dξ = ∞ k=0 B k k! t 0 ad k Θ [m] (ξ) A(ξ)dξ, m = 0, 1, . . . .
Rearranging terms for simplicity, we obtain

Θ [1] (t) = t 0 A(ξ 1 )dξ 1 , Θ [2] (t) = t 0 A(ξ 1 )dξ 1 -1 2 t 0 ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 + 1 12 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 + • • • Θ [3] (t) = t 0 A(ξ 1 )dξ 1 -1 2 t 0 ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 + 1 12 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 + 1 4 t 0 ξ1 0 ξ2 0 A(ξ 3 )ξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 -1 24 t 0 ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 -1 24 t 0 ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 -1 24 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 + • • •
and so on. The Picard theorem implies that Θ(t) = lim m→∞ Θ [m] (t) exists in a suitably small neighbourhood of the origin and the above first few iterations indicate that it can be expanded as a linear combination of terms that are composed from integrals and commutators acting recursively on the matrix A. This is the Magnus expansion

Θ(t) = ∞ k=0 H k (t), (4.4)
where each H k is a linear combination of terms that include exactly k + 1 integrals (or -and later we will see that it boils down to the same thing -k commutators).

Thus,

H 0 (t) = t 0 A(ξ 1 )dξ 1 , H 1 (t) = -1 2 t 0 ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 , H 2 (t) = 1 12 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 + 1 4 t 0 ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 , H 3 (t) = -1 24 t 0 ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 -1 24 t 0 ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 -1 24 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 ξ2 0 A(ξ 3 )dξ 3 , A(ξ 2 ) dξ 2 , A(ξ 1 ) dξ 1 -1 8 t 0 ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , ξ1 0 A(ξ 2 )dξ 2 , A(ξ 1 ) dξ 1 .
It is possible to derive the next few sets H k with enough perseverance and perhaps a little help from a computer algebra package. Yet, it is evident that the terms are becoming increasingly complex. A considerably more transparent form of the Magnus expansion, amenable for easy recursive derivation and easier discussion of computational issues, can be obtained by associating each term in the expansion with a rooted binary tree, an approach that has been pioneered by Iserles and Nørsett (1999).

Let us recall briefly relevant terminology of graph theory [START_REF] Harary | Graph Theory[END_REF].

• Let V = {v 1 , v 2 , . . . , v r } be a finite set of distinct vertices and E ∈ V × V a set of edges. (The edges (v i , v j ) and (v j , v i ) are identified with each other.) We say that G = ⟨V , E ⟩ is a graph. • The ordered set {(v s l , v t l ) : l = 1, 2, . . . , r} of edges is a path from v i ∈ V to v j ∈ V , i ̸ = j, if s 1 = i, t l = s l+1 , l = 1, 2, .
. . , r -1 and t r = j. • The graph is said to be connected if there is a path between any two distinct vertices. It is a tree if exactly one path links every two vertices. • A rooted tree is the pair T = (G, w), where G is a tree and w ∈ V is its root.

There exists natural partial order on T : we say that v i ≺ v j if v i precedes v j in the unique path extending from the root w to v j . In that case v i is the ancestor of v j , while v j is the successor of v i . (Thus, the root is the ancestor of all vertices in V \ {w}.)

• If v i ≺ v j and there is no v k ∈ V such that v i ≺ v k ≺ v j ,
we say that v i is the parent of v j (most graph-theory texts adopt a more sexist definition) and v j the child of v i . Childless vertices are called leaves. • If each vertex in a rooted tree has at most two children, T is called a binary tree. If each vertex has either exactly two children or is a leaf, T is said to be a strictly binary tree.

It is always worthwhile to draw a pictorial representation of a graph, whereby edges are merely undirected lines extending between vertices, which are denoted by black discs. The graph-theoretical convention is that the root is always placed at the bottom, although computer scientists occasionally defy gravity by locating it at the top. We will follow the mathematical convention. We commence our investigation of the Magnus expansion by rewriting (4.4) in the form

Θ(t) = ∞ k=0 H k (t) = ∞ k=0 τ ∈T k α(τ ) t 0 C τ (ξ)dξ, t ≥ 0, (4.5)
Thus, each C τ for τ ∈ T k is made out of exactly k integrals and k commutators, while each α(τ ) is a scalar constant. Recalling how the expansion has been obtained from Picard's iteration, we observe that the terms C τ can be obtained by just two composition rules:

1 The index set T 0 is a simpleton, T 0 = {τ • }, say, and C τ• (t) = A(t); 2 If τ 1 ∈ T m1 and τ 2 ∈ T m2 then there exists τ ∈ T m1+m2+1 such that

C τ (t) = t 0 C τ1 (ξ)dξ, C τ2 (t) . (4.6)
Although this observation makes the expansion somewhat more transparent, much greater transparency is obtained by identifying the index sets T k with subsets of binary rooted trees, in a manner that makes the above composition rules stand out pictorially. We express the relationship between the index τ and the term C τ (t) by the map τ ❀ C τ (t).

1 We let T 0 consist of the single rooted tree with one vertex, r , and

r ❀ A(t). 2 Suppose that T m1 ∋ τ 1 ❀ C τ1 (t) and T m2 ∋ τ 2 ❀ C τ2 (t). Then T m1+m2+1 ∋ r r ❅ ❅ τ 1 τ 2 ❀ t 0 C τ1 (ξ)dξ, C τ2 (t) . (4.7)
Thus, (4.6) is 'coded' in graph terminology by denoting integration by adding a root to a tree, while commutation is denoted by joining two trees with a common root. It is possible to show that all the terms in the Magnus expansion can be obtained in this manner (Iserles and Nørsett 1999).

To derive T 1 we have just one option, m 1 = m 2 = 0, and the outcome is a single tree,

τ 1 = r , τ 2 = r ⇒ τ = r ❅ ❅ r r r .
There are two possibilities in T 2 , namely m 1 = 0, m 2 = 1 and m 1 = 1, m 2 = 0. They yield

τ 1 = r , τ 2 = r ❅ ❅ r r r ⇒ τ = r r r ❅ ❅ ❅ ❅ r r r r , τ 1 = r ❅ ❅ r r r , τ 2 = r ⇒ τ = r ❅ ❅ r r r ❅ ❅ r r r .
Next, we construct T 3 :

τ 1 = r , τ 2 = r r r ❅ ❅ ❅ ❅ r r r r ⇒ τ = r r r r ❅ ❅ ❅ ❅ ❅ ❅ r r r r r r , τ 1 = r , τ 2 = r ❅ ❅ r r r ❅ ❅ r r r ⇒ τ = r ❅ ❅ r r r r ❅ ❅ r r r ❅ ❅ r r r , τ 1 = r ❅ ❅ r r r , τ 2 = r ❅ ❅ r r r ⇒ τ = r ◗ ◗ ✑ ✑ r r ❅ ❅ r r r r ❅ ❅ r r r , τ 1 = r r r ❅ ❅ ❅ ❅ r r r r , τ 2 = r ⇒ τ = r ❅ ❅ r r r r r ❅ ❅ ❅ ❅ r r r r , τ 1 = r ❅ ❅ r r r ❅ ❅ r r r , τ 2 = r ⇒ τ = r ❅ ❅ r r ❅ ❅ r r r ❅ ❅ r r r r .
The principle should be quite clear by now, as should be the correspondence between trees and expansion terms which, indeed, can be 'read' directly from the graph. Thus, the last tree can be 'recited' as "the integral of A, commuted with A, integrated, commuted with A, integrated and commuted with A".

We now have a recursive algorithm to derive the C τ s but recall that, to complete our description of (4.5), we also require the constants α(τ ). Fortunately, also this can be deduced from the tree formalism. We thus commence by letting α( r ) = 1 and continue by observing that each tree in ∪ k∈N T k can be written in a unique form as

τ = r r r r ❅ ❅ ❅ ❅ ❅ ❅ r r r r ❅ ❅ r r τ 1 τ 2 τ 3 τ s . . . (4.8)
for some s ≥ 1. We can assume by induction that α(τ i ) are already known for i = 1, 2, . . . , s. In that case it has been proved in (Iserles and Nørsett 1999) that

α(τ ) = B s s! s i=1 α(τ i ).
(4.9)

Note that B 2m+1 = 0 for m ≥ 1. This implies that many (although, unfortunately, not too many. . . ) terms in (4.5) vanish.

We can now write the first terms of the Magnus expansion in a tree formalism, using the convention that linear combination of trees corresponds to a linear combination of expansion terms,

Θ(t) = r r -1 2 r r ❅ ❅ r r r + 1 12 r r r r ❅ ❅ ❅ ❅ r r r r + 1 4 r r ❅ ❅ r r r ❅ ❅ r r r -1 8 r r ❅ ❅ r r ❅ ❅ r r r ❅ ❅ r r r r -1 24 r r ❅ ❅ r r r r ❅ ❅ r r r ❅ ❅ r r r -1 24 r r ◗ ◗ ✑ ✑ r r ❅ ❅ r r r r ❅ ❅ r r r -1 24 r r ❅ ❅ r r r r r ❅ ❅ ❅ ❅ r r r r + • • • . (4.10)
Note that the coefficient corresponding to the last tree in T 3 vanishes, this being a consequence of (4.9) and of B 3 = 0.

It is relatively easy to continue the expansion to higher-order terms. Moreover, it is easy to prove that

τ ∈ T k ⇒ t 0 C τ (ξ)dξ = O t k+1 , k∈ Z + (4.11)
for every sufficiently-smooth matrix function A. On the face of it, this gives a handy device to truncate the Magnus expansion to obtain an approximant of given order -we will bother later about the calculation of multivariate integrals. This, however, is grossly misleading, since the naive estimate (4.11) can be improved a very great deal: far fewer terms are required!

Convergence of the Magnus expansion

Before we proceed to improve the estimate (4.11) and even consider the question of designing a realistic numerical algorithms based on the Magnus expansion, we need to examine the issue of convergence.

It has been proved in (Iserles and Nørsett 1999) that convergence takes place for sufficiently small t ≥ 0, but the result was unduly pessimistic. A considerably better (and in a well-defined sense optimal) estimate has been obtained by [START_REF] Blanes | Magnus and Fer expansions for matrix differential equations: The convergence problem[END_REF]. Herewith we present briefly a short and elegant proof due to [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF].

Theorem 4.1 Suppose that the Lie algebra g is equipped with the norm ∥ • ∥.

The Magnus expansion (4.4) absolutely converges in this norm for every t ≥ 0 such that

t 0 ∥A(ξ)∥dξ ≤ 2π 0 dξ 4 + ξ[1 -cot(ξ/2)]
≈ 1.086868702. (4.12)

Proof. Integrating (4.3) and taking norms, we have by the triangle inequality and the trivial bound

∥ad k B A∥ ≤ (2∥B∥) k ∥A∥ that ∥Θ(t)∥ = t 0 dexp -1 Θ(ξ) A(ξ)dξ ≤ t 0 ∥ dexp -1 Θ(ξ) A(ξ)∥dξ ≤ t 0 ∞ k=0 |B k | k! (2∥Θ(ξ)∥) k ∥A(ξ)∥dξ = t 0 g(2∥Θ(ξ)∥)∥A(ξ)∥dξ, where g(x) = 2 + x 2 1 -cot x 2 .
We now use a Bihari-type inequality from [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF]): suppose that h, g, v ∈ C(0, t * ) are positive and that g is nondecreasing. Then

h(t) ≤ t 0 g(h(ξ))v(ξ)dξ, t ∈ (0, t * ) implies that h(t) ≤ g-1 t 0 v(ξ)dξ , t ∈ (0, t * * ), where g(x) = x 0 dξ g(ξ)
and t * * ∈ (0, t * ] is such that g t 0 v(ξ)dξ is bounded in (0, t * * ). In our case h(t) = 2∥Θ(t)∥, v(t) = ∥A(t)∥ and g(t) are all positive and the latter is nondecreasing for t ∈ (0, 2π). Therefore,

∥Θ(t)∥ ≤ 1 2 g-1 t 0 ∥A(ξ)∥dξ
and ∥Θ(t)∥ is bounded, provided that g t 0 ∥A(ξ)∥dξ is bounded. The latter holds as long as condition (4.12) is satisfied.

✷

The condition of Theorem 4.1 has been recently improved for a more relaxed convergence framework. Moan (n.d.) proved that the Magnus expansion converges in norm for all t ∈ (0, t * ) with regard to the Euclidean norm,

lim m→∞ Θ(t) - m k=0 H m (t) 2 = 0, provided that t * 0 ∥A(ξ)∥ 2 dξ < π.
(4.13)

Magnus expansion cannot be expected to converge always. In its numerical implementation, this means that the expansion (like any other numerical method for ODEs) need be advanced in a time-stepping fashion, rather than being applied globally. Yet, the condition of Theorem 4.1 is not unduly restrictive and (4.13) is even less so. It might be problematic for stiff systems, a subject that has not received much attention in the study of Lie-group methods. An important exception is the use of Magnus expansions in the calculation of Sturm-Liouville spectra, where an elegant device allows to integrate the equations way beyond the formal upper bound (4.12) (cf. Section 11.2 and (Moan 1998)).

Power of a tree and time symmetry

Following Iserles, Nørsett and Rasmussen (1998), we say that a tree τ

∈ ∪ k∈Z + T k is of power m if m ≥ 0 is the least integer such that C τ (t) = O(t m )
for all sufficiently smooth matrix functions A. Letting F m be the set of all trees of power m and truncating the Magnus expansion,

Θ p (t) = p-1 m=0 τ ∈Fm α(τ ) t 0 C τ (ξ)dξ, (4.14)
it is trivial to verify that Θ p (t) = Θ(t) + O t p+1 and we have an order-p approximant. We already know from (4.11) that τ ∈ T k implies τ ∈ F m for some m ≥ k. How large can m get, though? The main mechanism that increases order is commutation. Thus, suppose that

C τi (t) = D i t mi + E i t mi+1 + F i t mi+2 + • • • , i= 1, 2. Then t 0 C τ1 (ξ)dξ, C τ2 (t) = 1 m 1 + 1 [D 1 , D 2 ]t m1+m2+1 + 1 m 1 + 1 [D 1 , E 2 ] + 1 m 1 + 2 [E 1 , D 2 ] t m1+m2+2 + 1 m 1 + 1 [D 1 , F 2 ] + 1 m 1 + 2 [E 1 , E 2 ] + 1 m 1 + 3 [F 1 , D 2 ] t m1+m2+3 + • • • .
In general, [D 1 , D 2 ] ̸ = O, hence the new term resides in F m1+m2+1 . However, cancellation takes place in the important special case τ 1 = τ 2 , whence

t 0 C τ1 (ξ)dξ, C τ1 (t) = 1 (m 1 + 1)(m 1 + 2) [D 1 , E 1 ]t 2m1+2 + 2 (m 1 + 1)(m 1 + 3) [D 1 , F 1 ]t 2m1+3 + • • •
and the new term is in F 2m1+2 , a gain of one unit. While F 0 = { r }, we observe that F 1 = ∅ and

F 2 = r ❅ ❅ r r r F 3 = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ r r r ❅ ❅ ❅ ❅ r r r r , r ❅ ❅ r r r ❅ ❅ r r r ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ , F 4 = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ r r r r ❅ ❅ ❅ ❅ ❅ ❅ r r r r r r , r ❅ ❅ r r r r ❅ ❅ r r r ❅ ❅ r r r , r ❅ ❅ r r r r r ❅ ❅ ❅ ❅ r r r r , r ❅ ❅ r r ❅ ❅ r r r ❅ ❅ r r r r ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ .
In general, the number of terms counted according to power is significantly smaller than similar enumeration by the number of commutators. It is possible to prove that lim sup and Nørsett 1999, Iserles et al. 1998). As an example of the reduction in cardinality, compare #T 6 = 132 with #F 6 = 21. Using the 'truncation by power' (4.14) is thus aptly justified. However, before we rush to pronounce this as the 'correct' truncated Magnus expansion, we need to pay attention to yet another device that reduces the number of terms in the expansion.

k→∞ (#T k ) 1/k = 4, lim sup m→∞ (#F m ) 1/m ≈ 3.11674 (4.15) (Iserles
Let Φ t be the flow corresponding to the differential equation (4.2), Φ t (Y 0 ) = Y (t). It is obvious that the flow is time symmetric, Φ -t • Φ t = Id, since integrating from 0 to t and back to 0 returns us to the original initial value. What is far less obvious, yet has been proved in [START_REF] Iserles | Time-symmetry and highorder Magnus methods[END_REF], is that the truncation by power (4.14) respects time symmetry. In other words, let Φt (Y 0 ) = e Θp(t) Y 0 , t≥ 0.

Then Φ-t • Φt = Id. This is remarkable, since any analytic time-symmetric map S t can be represented in the form S t = e Ft where the map F t is expandable in odd powers of t only [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]). This fits our framework perfectly.

Theorem 4.2 The function Θ p can be expanded in odd powers of t and

Θ 2q-1 (t) = Θ(t) + O t 2q+1 , q∈ N.
Therefore, truncating by power with odd p leads to a gain of an extra unit of order! This is a critical observation which leads to substantial savings in high-order Magnus expansions.

It is important to realise that it is not true that individual elements in H 2q-1 , q ≥ 2, are O t 2q+1 : it is their linear combination that knocks out the O t 2q term!

We may now re-examine the expansion (4.10), truncating by power and identifying the order. Reverting from trees to standard notation, we have

Θ(t) =
+ 1 12 t 0 ξ1 0 ξ1 0 [A(ξ 2 ), [A(ξ 3 ), A(ξ 1 )]]dξ + 1 4 t 0 ξ1 0 ξ2 0 [[A(ξ 3 ), A(ξ 2 )], A(ξ 1 )]dξ -1 24 t 0 ξ1 0 ξ1 0 ξ3 0 [A(ξ 2 ), [[A(ξ 4 ), A(ξ 3 )], A(ξ 1 )]]dξ -1 24 t 0 ξ1 0 ξ2 0 ξ2 0 [[A(ξ 3 ), [A(ξ 4 ), A(ξ 2 )]], A(ξ 1 )]dξ -1 8 t 0 ξ1 0 ξ2 0 ξ3 0 [[[A(ξ 4 ), A(ξ 3 )], A(ξ 2 )], A(ξ 1 )]dξ . . . . order 6 (4.18) + • • •
Very often entries of A and its commutators can be integrated explicitly, e.g. when they are polynomials or trigonometric functions. In that case the truncated Magnus expansions (4.16-18), say, can be computed explicitly. However, a more comprehensive numerical approach to Magnus expansions requires the computation of multivariate integrals. Although at the first glance this may appear as a very formidable task, it turns out that the special structure of 'Magnus integrals' renders them amenable to very effective and affordable numerical treatment. We defer the discussion of this issue to Section 5.

Fer expansions

At the first instance, we wish to represent the solution of (4.2) in the form

Y (t) = e R t 0 A(ξ)dξ V (t), t≥ 0. (4.19)
Direct differentiation yields

V ′ = d dt e - R t 0 A(ξ)dξ Y + e - R t 0 A(ξ)dξ Y ′ = d dt e - R t 0 A(ξ)dξ e R t 0 A(ξ)dξ V + e - R t 0 A(ξ)dξ Ae R t 0 A(ξ)dξ V = Ad expm(- R t 0 A(ξ)dξ) A(t) -dexp expm(- R t 0 A(ξ)dξ) A(t) V.
Recalling from (2.24) and (2.43) that Ad expm(E) D = e adE D,

dexp E D = e adE -I ad E D,
we deduce that the correction term V itself obeys a linear differential equation,

V ′ = (I + ad E ) e -adE -I ad E D V, t ≥ 0, V(0) = Y 0 ,
where

D = A(t), E = t 0 A(ξ)dξ.
This is indeed the main step in constructing the Fer expansion: the correction V in (4.19) satisfies the equation

V ′ = ∞ k=1 (-1) k k (k + 1)! ad k R t 0 A(ξ)dξ A(t) V, t ≥ 0, V(0) = Y 0 . (4.20)
The idea is now to iterate (4.20). Thus, we let B 0 = A and generate the sequence {B m } m∈Z + recursively, where

B m (t) = ∞ k=1 (-1) k k (k + 1)! ad k R t 0 Bm-1(ξ)dξ B m-1 (t), t≥ 0, m ∈ N. (4.21)
The Fer expansion of the solution of (4.2) is

Y (t) = e R t 0 B0(ξ)dξ e R t 0 B1(ξ)dξ e R t 0 B2(ξ)dξ • • • Y 0 , t≥ 0. (4.22)
This expansion was introduced by [START_REF] Fer | Résolution del l'equation matricielle U = pU par produit infini d'exponentielles matricielles[END_REF] who, remarkably, did not recognise that it respects Lie-group structure. It was rediscovered by [START_REF] Iserles | Solving linear ordinary differential equations by exponentials of iterated commutators[END_REF] in a numerical context but, again, Lie groups were not mentioned. Finally, [START_REF] Zanna | The method of iterated commutators for ordinary differential equations on Lie groups[END_REF] recognised (4.22) as a Lie-group solver. In [START_REF] Zanna | Iterated commutators, Lie's reduction method and ordinary differential equations on matrix Lie groups[END_REF] it is shown that this expansion can be understood as a version of so-called Lie reduction.

From this it follows that if g is solvable (cf. Section 6.5) then the expansion in (4.2) is exact for a finite product of exponentials. The first step in a numerical implementation of the Fer expansion (4.22) is truncation of the infinite product. To this end it is vital to recognise the rate of decay of the matrices B m . Assuming that B m (t) = O(t pm ), m ∈ Z + , it is easy to verify from (4.21) that p m = 2p m-1 + 2. This, in tandem with p 0 = 0, yields p m = 2 m+1 -2 and we deduce that e

R t 0 B0(ξ)dξ e R t 0 B1(ξ)dξ • • • e R t 0 Bs-1(ξ)dξ Y 0 = Y (t) + O t 2 s+1 -1 . (4.23)
Thus, we obtain an approximant of order 2 s+1 -2: the order grows exponentially with s (Iserles 1984)! Of course, if we are interested in an order-p Fer approximant to the solution of (4.2), there is no need to carry out the summation in (4.21) ad infinitum. Systematic analysis of order conditions necessary for the formation of Fer approximants of various orders has been carried out by [START_REF] Zanna | On the Numerical Solution of Isospectral Flows[END_REF] using the same binary rooted trees that we have already encountered in our analysis of the Magnus expansion. Identifying A with the single-vertex tree, we have the following explicit form of a Fer approximant of any given order:

Order 2: s = 1, B 0 (t) : r . Order 3: s = 2, B 0 (t) : r , B 1 (t) : 1 2 r ❅ ❅ r r r . Order 4: s = 2, B 0 (t) : r , B 1 (t) : 1 2 r ❅ ❅ r r r + 1 3 r r r ❅ ❅ ❅ ❅ r r r r . Order 5: s = 2, B 0 (t) : r , B 1 (t) : 1 2 r ❅ ❅ r r r + 1 3 r r r ❅ ❅ ❅ ❅ r r r r + 1 8 r r r r ❅ ❅ ❅ ❅ ❅ ❅ r r r r r r . Order 6: s = 2, B 0 (t) : r , B 1 (t) : 1 2 r ❅ ❅ r r r + 1 3 r r r ❅ ❅ ❅ ❅ r r r r + 1 8 r r r r ❅ ❅ ❅ ❅ ❅ ❅
1 8 r ◗ ◗ ✑ ✑ r r ❅ ❅ r r r r ❅ ❅ r r r .
and so on. Note that just a subset of 'Magnus trees' occurs in the expansion. This representation of expansion terms as linear combination of trees is central to the application of the multivariate quadrature algorithms of the next section and derivation of practical numerical methods for the Fer expansion.

5. Quadrature and graded algebras

Multivariate quadrature over polytopes

Casting our eyes again over the Magnus method (4.16-18) and considering its numerical implementation, let us discuss the computation of the first four integrals, noting that each need be carried out over a different polytope:

I 1 (t) = t 0 A(ξ)dξ over the line segment 0 t I 2 (t) = t 0 ξ1 0 [A(ξ 2 ), A(ξ 1 )]dξ over the triangle ❅ ❅ ❅ 0 t 0 t I 3 (t) = t 0 ξ1 0 ξ1 0 [A(ξ 2 ), [A(ξ 3 ), A(ξ 1 )]]dξ over the prism ❅ ❅ ❅ ✏ ✏ ✏ ✏ ❅ ❅ ❅ 0 t 0 t 0 t I 4 (t) = t 0 ξ1 0 ξ2 0 [[A(ξ 3 ), A(ξ 2 )], A(ξ 1 )]dξ over the pyramid ❅ ❅ ❅ ☞ ☞ ❍ ❍ ❍ 0 t 0 t 0 t
Unless we can replace all these integrals by affordable and accurate quadrature, the Magnus method (and by the same token the Fer method) of non-trivial order is of little but theoretical value. Yet, multivariate quadrature is notoriously expensive in terms of function evaluations [START_REF] Cools | Constructing cubature formulas: the science behind the art[END_REF]. Fortunately, the special nature of integrals occurring within the context of Magnus and Fer expansions renders them particularly suitable for numerical quadrature with a remarkably small number of function evaluations (Iserles and Nørsett 1999).

We commence by observing that, time-stepping with step h > 0, each Magnus or Fer expansion term is of the form

I(h) = S L(A(ξ 1 ), A(ξ 2 ), . . . , A(ξ s ))dξ, (5.1)
where L is a multilinear form, while S is a polytope of a special form,

S = {ξ ∈ R s : ξ 1 ∈ [0, h], ξ l ∈ [0, ξ m l ], l = 2, 3, . . . , s},
where m l ∈ {1, 2, . . . , l -1}, l = 2, 3, . . . , s. Thus, for example,

t 0 ξ1 0 ξ1 0 [A(ξ 2 ), [A(ξ 3 ), A(ξ 1 )]]dξ and t 0 ξ1 0 ξ2 0 [[A(ξ 3 ), A(ξ 2 )], A(ξ 1 )]dξ yield s = 3 and L(E 1 , E 2 , E 3 ) = [E 2 , [E 3 , E 1 ]], m 2 = 1, m 3 = 1 and L(E 1 , E 2 , E 3 ) = [[E 3 , E 2 ], E 1 ], m 2 = 1, m 3 = 2
respectively. Following Iserles and Nørsett (1999), we propose to discretise I(h) as follows: Choose ν distinct quadrature points, c 1 , c 2 , . . . , c ν ∈ [0, 1], evaluate A k = hA(c k h), k = 1, 2, . . . , ν and form the quadrature

K(h) = k∈C ν s b k L(A k1 , A k2 , . . . , A ks ), (5.2)
where C ν s is the set of all combinations of length s from the set {1, 2, . . . , ν}. The weights are

b k = S s i=1 ℓ ki (ξ i )dξ, (5.3) where S = {ξ ∈ R s : ξ 1 ∈ [0, 1], ξ l ∈ [0, ξ m l ], l = 2, 3, . . . , s}
is the polytope S scaled to the unit cube and

ℓ j (x) = ν i=1 i̸ =j x -c j c i -c j , j = 1, 2, . . . , ν,
are the familiar cardinal polynomials of Lagrange's interpolation. Note that (5.3) are interpolatory weights: they follow naturally by substituting the interpolation polynomial at the quadrature points,

Ã(t) = h -1 ν k=1 ℓ k t h A k (5.4)
in place of A(t) in (5.1) and carrying out the integration explicitly. Conceptually, (5.2) recycles ν function values at all possible combinations at the s 'slots' of the multilinear function L.

How well does the quadrature (5.2) approximate the integral (5.1)? The answer is straightforward in the case s = 1, L(E) = E, since we recover standard univariate interpolatory quadrature, which is of order ν +m, where m ≥ 0 is the largest integer so that

1 0 ξ i-1 c(ξ)dξ = 0, i = 1, 2, . . . , m where c(t) = ν k=1 (t -c k ) (5.5)
is the collocation polynomial.

Theorem 5.1 The orthogonality condition (5.5) implies that the quadrature rule (5.2) is of order ν + m for all polytopes S and all multilinear forms L. In particular, if c 1 , c 2 , . . . , c ν are the roots of the Legendre polynomial P ν , shifted to the interval [0, 1] (Gauss-Legendre points), then the quadrature is of order 2ν.

We do not propose to present the proof from (Iserles and Nørsett 1999), which is long, technical and not particularly illuminating. In the sequel, in the context of Magnus methods for nonlinear Lie-group equations, we describe a much clearer argument due to [START_REF] Zanna | Collocation and relaxed collocation for the Fer and the Magnus expansions[END_REF] that explains why a suitable linear combination of quadratures (5.2) approximates the truncated Magnus expansion to the order reported in Theorem 5.1. Instead, we present an example, the quadrature of the four integrals that have opened this subsection using (5.2) at three Gauss-Legendre points.

Letting ν = 3, hence order six, we have

c 1 = 1 2 - √ 15 10 , c 2 = 1 2 , c 3 = 1 2 + √ 15 10 , whence ℓ 1 (x) = 5 + √ 15 6 - 10 + √ 15 3 x + 10 3 x 2 , ℓ 2 (x) = - 2 3 + 20 3 x - 20 3 x 2 , ℓ 3 (x) = 5 - √ 15 6 - 10 - √ 15 3 x + 10 3 x 2 .
Insofar as the univariate integral I 1 is concerned, we have the familiar Gauss-Legendre quadrature, I 1 (h) ≈ 1 18 (5A 1 + 8A 2 + 5A 3 ), while, after taking into account skew-symmetry of commutators, the quadrature of the planar integral is

I 2 (h) ≈ √ 15 54 (2[A 1 , A 2 ] + [A 1 , A 3 ] + 2[A 2 , A 3 ]].
Insofar as the two cubic integrals are concerned, we are interested (cf. (4.18)) in their linear combination which, after a great deal of simplification, yields

1 12 I 3 (h) + 1 4 I 4 (h) ≈ 1 27216 (94[A 1 , [A 1 , A 2 ]] + 45[A 1 , [A 1 , A 3 ]] + 194[A 1 , [A 2 , A 3 ]] -152[A 2 , [A 1 , A 2 ]] + 152[A 2 , [A 2 , A 3 ]] (5.6) -194[A 3 , [A 1 , A 2 ]] -45[A 3 , [A 1 , A 3 ]] -94[A 3 , [A 2 , A 3 ]]).
Before any attempts are made either to anoint (5.2) as a new wonder-quadrature or to expand efforts to optimise its calculation, we hasten to say that it is suboptimal! Indeed, it is an immediate consequence of Theorem 4.2 that 1 12 I 3 (h)+ 1 4 I 4 (h) = O h 4 , rather than O h 3 . Moreover, some commutators can be expressed in terms of other commutators. Although all this can be done on an ad hoc basis, it is significantly better and more efficient to understand this phenomenon mathematically and apply the fruits of our understanding not just to truncated Magnus method but also to other Lie-group solvers.

Quadrature (5.2) scores exceedingly well in terms of function evaluations: the number of function evaluations required to compute all the integrals in a Magnus expansion to requisite order is the same as the cost of the corresponding univariate Gauss-Legendre quadrature! The trade-off, though, is that this approach requires a very large volume of linear-algebra calculations, since the number of all combinations is inordinately large. Fortunately, the cost of linear algebra can be reduced a very great deal by exploiting the theory of free Lie-algebras. This salutary example of a mathematical theory purer than a driven snow finding a very practical application in the design of numerical algorithms is told in the remainder of this section.

The self-adjoint basis

The first step in the effort to improve the multivariate quadrature formula (5.2) is a change of basis. As suggested first by Munthe-Kaas and Owren (1999), we choose c 1 < c 2 < • • • c ν symmetric with respect to 1 2 and replace the function values A 1 , A 2 , . . . , A ν with the solution of the Vandermonde system

ν l=1 (c k -1 2 ) l-1 B l = A k , k= 1, 2, . . . , ν.
(5.7)

We say that {B 1 , B 2 , . . . , B ν } is a self-adjoint basis and note for future reference that Gauss-Legendre points in [0, 1], being symmetric with respect to 1 2 , lead to such basis.

Proposition 5.2 Given a sufficiently smooth matrix function A, it is true that

B l = c l h l A (l-1) ( 1 2 h) + O h l+1
, where c l ̸ = 0 is a scalar constant, l = 1, 2, . . . , ν. Moreover, each h -l B l can be expanded in even powers of h.

It follows from the proposition that the interpolating polynomial can be written in the form

Ã(t) = h -1 ν-1 l=0 B l ( t h -1 2 ) l .
(5.8) Substituting this into (5.1) allows us to rephrase the quadrature formula (5.2) in a considerably more convenient form, (ξ i -1 2 ) li-1 dξ.

K(h) =
(5.10)

Proposition 5.3. (Munthe-Kaas and Owren (1999)) Suppose that a linear combination of integrals I(h) can be expanded in odd powers of h. Then so can be the linear combination of quadratures (5.9).

Recall from Section 4.3 that Magnus series, truncated by power, conform with the assumptions of Proposition 5.3: at a stroke, roughly half commutators go away. Yet, this is but the first of three important savings that are a consequence of the change of a basis.

An element F ∈ g constructed from the basis terms B 1 , B 2 , . . . , B ν by the standard Lie-algebra operations of linear combination and commutation is said to be of grade m if F = O(h m ) for all sufficiently smooth matrix functions A. This is denoted by ω(F ) = m. We note from Proposition 5.2 that ω(B l ) = l, l = 1, 2, . . . , ν. Moreover, the grade is inherited under commutation, e.g. ω([B k , B l ]) = ω(B k ) + ω(B l ) and, with greater generality,

ω(L(B l1 , B l2 , . . . , B ls )) = |l| = ν i=1 l i .
By the definition of the grade, this is equivalent to

L(B l1 , B l2 , . . . , B ls ) = O h |l| .
Thus, as long as we are interested in an order-p quadrature, we can discard higherorder terms in (5.10). The outcome is

K(h) = l∈ Ĉν,p s bl L(B l1 , B l2 , . . . , B ls ),
(5.11)

where Ĉν,p

s ⊆ C ν s such that l ∈ Ĉν,p s ⇔ |l| ≤ p.
Let us recall the four integrals from Section 5.1. We presently obtain using (5.11) the following order-6 quadrature formulae using Gauss-Legendre points with ν = 3. In line with (5.7), we let

B 1 = A 2 , B 2 = √ 15 10 (A 3 -A 1 ) and B 3 = 20 3 (A 3 -2A 2 +A 1 ). I 1 (h) ≈ B 1 + 1 12 B 3 , I 2 (h) ≈ 1 6 [B 2 , B 1 ] -1 120 [B 3 , B 2 ], I 3 (h) ≈ -1 8 [B 1 , [B 2 , B 1 ]] + 1 80 [B 2 , [B 2 , B 1 ]] -1 120 [B 1 , [B 3 , B 1 ]] + 1 480 [B 1 , [B 3 , B 2 ]] + 1 240 [B 2 , [B 3 , B 1 ]] -7 480 [B 3 , [B 2 , B 1 ]] -1 160 [B 1 , [B 3 , B 1 ]], I 4 (h) ≈ 1 24 [B 1 , [B 2 , B 1 ]] + 1 80 [B 2 , [B 2 , B 1 ]] -1 120 [B 1 , [B 3 , B 1 ]] -1 1440 [B 1 , [B 3 , B 2 ]] -1 720 [B 2 , [B 3 , B 1 ]] + 7 1440 [B 3 , [B 2 , B 1 ]] + 1 480 [B 1 , [B 3 , B 1 ]].
Moreover, consistently with Proposition 5.3, we have

1 12 I 3 (h) + 1 4 I 4 (h) ≈ 1 240 [B 2 , [B 2 , B 1 ]] -1 360 [B 1 , [B 3 , B 1 ]].
Just two terms survive! Taken together, throwing away terms of high enough grade and the removal of terms with even |l| after summation removes a high proportion of commutators. Having said this, the most important feature of the self-adjoint basis that allows us to reduce the computational cost has not been mentioned mentioned yet!

We have already exploited skew symmetry of the commutator in the derivation and 'beautification' of our integration formulae. This is a fairly transparent procedure. However, let us recall that the commutator is also subject to the Jacobi identity (2.10). This allows for a very powerful mechanism to express commutators as linear combinations of other commutators and leads to results that are of importance not just to Magnus expansions, but also to RK-MK methods and the evaluation of the BCH formula. This is the theme of the next subsection.

Free Lie algebras

In the previous subsections we have seen how to construct an approximation to the solution Y of the Lie-group differential equation Y ′ = A(t)Y by means of linear combinations of matrices G 1 , G 2 , . . . , G ν ∈ g, whereby the G i are either 'samples' of the matrix function A(t), in which case G i = hA(c i h), or terms of the self-adjoint basis, in which case the G i s coincide with the matrices B i of Section 5.2.

It is clear that, if we want to make the best out of the properties of the commutator (skew-symmetry and Jacobi identity) it is useful to depart from specific representations (A i s and B i s) and treat the G i s as abstract objects in an abstract algebra g that embodies the structure that is common to all Lie algebras but nothing more. This is the main idea behind free Lie algebras, a formalisation of a Lie algebra whose terms can be generated by means of brackets of pairwise elements and such that there are no reducing mechanisms other than skew-symmetry and the Jacobi identity of the commutator.

More precisely, the following definition formalises the concept of free Lie algebras presented above (Munthe-Kaas and Owren 1999) .

Definition 5.1 Let I be a set of indices, either finite or countable. A Lie algebra g is free over the set I if i) For every index i ∈ I there exists G i ∈ g; ii) For any Lie algebra h and any function i ∈ I → H i ∈ h there exists a unique Lie algebra homomorphism π : g → h such that π(G i ) = H i for all i ∈ I.

Moreover, S = {G i } i∈I is called the set of generators of the free Lie algebra g.

In our exposition it is useful to think of a free Lie algebra (FLA) as a linear space and to describe it in terms of a basis. The Hall basis H can be endowed with a total ordering defined recursively on the length ℓ of its elements. In general, we say that

G ≺ H if ℓ(G) < ℓ(H). If ℓ(G) = ℓ(H) then G ≺ H if G
precedes H in lexicographic order. Moreover, to take into account skew-symmetry and the Jacobi identity, we require that

• elements of length two [G i , G j ] are included in H if G i ≺ G j ;
• elements of length greater or equal to three are included in H if they are of the form

[H i , [H j , H k ]], with H i , H j , H k , [H j , H k ] ∈ H and moreover H j ≼ H i ≺ [H j , H k ].
An example of the first terms of a Hall basis generated by three elements

G 1 , G 2 , G 3 is given by G 1 , G 2 , G 3 , [G 1 , G 2 ], [G 1 , G 3 ], [G 2 , G 3 ], [G 1 , [G 1 , G 2 ]], [G 1 , [G 1 , G 3 ]], [G 2 , [G 1 , G 2 ]], [G 2 , [G 1 , G 3 ]], [G 2 , [G 2 , G 3 ]], [G 3 , [G 1 , G 2 ]], [G 3 , [G 1 , G 3 ]], [G 3 , [G 2 , G 3 ]], ....
We shall not go into details of algorithmic construction of the Hall basis, which can be found in (Bourbaki 1975, Munthe-Kaas andOwren 1999). It is interesting, however, to mention how fast the number of elements of the Hall basis grows: assuming that I is finite and consists of ν indices (equivalently, S consists of ν generators), the linear subspace of terms of length exactly equal to m has dimension

ρ m = 1 m d|m µ(d)ν m/d ,
the sum being carried over all integers d dividing m, a result known as Witt's formula. The function µ is the Möbius function, defined as follows: assume that d can be factorized as

d = d n1 1 d n2 2 • • • d nq q
, with each d i a prime number and n i ≥ 1. Then

µ(d) = ⎧ ⎨ ⎩ 1,
d= 1, (-1) q , n i = 1 for all i = 1, 2, . . . , q, 0.

otherwise.

The number ρ m grows quite fast, as illustrated in Table 5.1 for ν = 3. Why is all this relevant to our discussion? Let us represent

G 1 = A 1 , G 2 = A 2 , . . . , G ν = A ν , where A i = hA(c i h), for i = 1, . . . , ν. Since A i = O(h), a
term containing exactly m commutators corresponds to a combination that is at least of order O(h m ). Thus, in order to have a numerical approximation of order m=1 ρ m if the c i are symmetrically distributed with respect to 1 2 . For instance, for a sixth-order Gauss-Legendre scheme (ν = 3), this leads to 80 linearly independent terms of the Hall basis, not counting the number of commutators involved!

The growth of ρ m can be reduced introducing a grading ω of the FLA g, as suggested by Munthe-Kaas and Owren (1999). Assume that

ω(G i ) = ω i ∈ N, G i ∈ S, i = 1, 2, . . . , ν.
The grading propagates in a natural manner in the Hall basis

H: for all H ∈ H of the form H = [H 1 , H 2 ] we let ω(H) = ω(H 1 ) + ω(H 2 ).
A consequence of the grading is that the Hall basis H splits into a disjoint union of sets of grade m,

H = ∞ m=1 H m , H m = {H ∈ H : ω(H) = m}, consequently g becomes a direct sum of subspaces, g = ∞ m=1 g m , g m = span H m ,
and we say that g is a graded FLA algebra.

The fundamental result on the dimension of g m is due to Munthe-Kaas and Owren (1999).

Theorem 5.4 Let g be the graded FLA generated by S = {G 1 , G 2 , . . . , G ν }, with grades ω 1 , ω 2 , . . . , ω ν respectively. Denote by λ 1 , λ 2 , . . . , λ r the roots of the rth degree polynomial

p(z) = 1 - ν i=1 z ωi , r= max 1≤i≤ν ω i . Then dim g m = ρm = 1 m d|m r i=1 λ m/d i µ(d).
(5.12)

To illustrate the benefits of the grading, Table 5.2 displays dim g m , for m = 1, 2, . . . , 10 for a graded algebra generated by G 1 , G 2 , G 3 with weights ω 1 = 1, ω 2 = 2, ω 3 = 3. Comparison with Table 5.1 reveals a significant reduction in the number of linearly independent terms. The results have been obtained using the Matlab package DiffMan, which will be further discussed in Section 10.2 (Engø et al. 1999). Linking again with the theory of Section 5.1-2, the case of a graded algebra with generators G 1 , . . . , G ν and weights 1, 2, . . . , ν, corresponds to the realization

G i = h i (i -1)! A (i-1) (ξ i ), i= 1, 2, . . . , ν,
where A (i-1) (ξ i ) is the (i -1)th derivative of the function A for some ξ i ∈ (0, h).

The weight ω i merely indicates that G i is an O h i term. Moreover,

Ã(t) = h -1 ν i=1 G i ( t h ) i ,
is exactly the collocation polynomial (5.4), the information about the nodes c i being hidden in the G i s. The equivalence is revealed by means of the Vandermonde transformation

A i = ν j=1 c j-1 j G j , i= 1, . . . , ν.
This procedure applies to all kinds of collocation, whether with symmetric nodes or otherwise. In particular we deduce from Table 5.2 that with three collocation points is is possible to obtain order six with at most 15 terms! A substantial saving, compared with 80 using an ungraded algebra. . . . To reduce further the dimension of the graded FLA, we exploit the argument of Theorem 5.2, assuming that the collocation points are symmetric in [0,1] with respect to 1 2 . We construct

Ã(t) = h -1 ν i=1 G i ( t h -1 2 ) i ,
where G i ≡ B i are the self-adjoint basis introduced in Section 5.2. Theorem 5.2 implies that only terms with odd grades need be considered! As an example, if S = {G 1 , G 2 , G 3 }, with weights 1, 2 and 3 respectively, then the growth of the dimension of the graded FLA is given by m 1 2 3 4 5 6 7 8 9 10

m i=1 i odd
ρi 1 1 3 3 7 7 11 11 37 37 and comparison with Table 5.2 is remarkable. In particular, we deduce that for methods based on three symmetric collocation points (for instance a sixth-order Gauss-Legendre), we need at most seven terms of the graded Hall basis H. Letting G i = B i , the seven terms are

B 1 , B 3 , [B 1 , B 2 ], [B 2 , B 3 ], [B 1 , [B 1 , B 3 ]], [B 2 , [B 1 , B 2 ]], [B 1 , [B 1 , [B 1 , B 2 ]]].
Bounds on the number of independent terms for different orders are given below for methods based on Gauss-Legendre quadrature. The sharpest bound corresponds to the case of a graded FLA including only odd terms.

ν (stages) 1 2 3 4 5 2ν (order of the method) 2 4 6 8 10

2ν-1 i=1 ρ i 1 5 80 3304 > 10000 2ν-1 i=1 ρi 1 3 10 33 111 2ν-1 i=1 i odd ρi 1 2 7 22 73
It is important to note that p-1 i=1 ρi odd p i=1 ρ i respectively in the case of the self-adjoint (non self-adjoint resp.) basis is an upper bound on the number of linearly-independent commutators required for a method of order p.

We have seen the advantage of changing the basis in the case of linear equations Y ′ = A(t)Y . Similar savings can be also achieved in the case of explicit RK-MK methods for the general equation Y ′ = A(t, Y )Y . By combining the stage values A i computed by the algorithm, we seek linear combinations B i of highest-possible grade. However, in order to obtain an explicit method, we must require that B i are related to A i by a triangular matrix. Optimal combinations can be found using linear algebra and the theory of B-series. The details can be found in (Munthe-Kaas and Owren 1999). It turns out that in this case it is in general not possible to change basis in such a way that ω(B i ) = i: the weights grow slower. The method (A.7) originates in the classical four-stage RK4 scheme. The resulting basis B i has the grading 1, 2, 3, 3, and the scheme has just two commutators. Similar savings have been realised for the seven-stage DOPRI5(4) method in (Munthe-Kaas and Owren 1999).

Reducing further the number of commutators

The theory of free Lie algebras allows us to derive an upper bound on the number of linearly independent terms required to obtain numerical methods of given order, taking into account skew-symmetry and the Jacobi identity of the commutator. It also provides algorithms to compute the requisite pattern of dependency, e.g. in terms of the Hall basis.

Although the theory of free Lie algebras estimates the numbers of commutators for a method of order p, this by no means indicates the least number of commutators required for a method of a given order.

For general N × N matrices, the computation of the commutator is an O N 3 operation, a cost that quickly adds up when we consider methods of order three and higher.

At present there is no theory that systematically reduces the number of commutators to minimum. However, we shall present a technique, due to Blanes, Casas and Ros (1999) that gives, true for today, the least number of commutators for methods based on Gauss-Legendre and Newton-Cotes quadratures up to order ten. Before proceeding further, it is important to remark that the content of this section applies to linear Lie group problems Y ′ = A(t)Y solved with Magnus-type/RK-MK methods based on symmetric nodes in [0, 1]. It is convenient to illustrate the procedure with an example. Let us thus construct a sixth-order Gauss-Legendre method based on the collocation nodes

c 1 = 1 2 - √ 15 10 , c 2 = 1 2 , c 3 = 1 2 + √ 15 10 .
Using the toolbox of graded algebras, we obtain

Θ = B 1 + 1 12 B 3 -1 12 [B 1 , B 2 ] + 1 240 [B 2 , B 3 ] + 1 360 [B 1 , [B 1 , B 3 ]] -1 240 [B 2 , [B 1 , B 2 ]] + 1 720 [B 1 , [B 1 , [B 1 , B 2 ]]].
Since in the self-adjoint basis B i = O h i , note that Θ includes terms with odd powers of h only and, in this form, can be evaluated by computing just seven commutators.

Let us focus on the portion of Θ consisting of single commutators,

C 2 = -1 12 [B 1 , B 2 ] + 1 240 [B 2 , B 3 ].
The first fundamental observation is that terms of the form [B i , B j ] can only be obtained if one of the indices is even and the other is odd. Therefore, we look for a linear combination

[b 1 B 1 + b 3 B 3 , b 2 B 2 ] ( 5 . 1 3 )
for some real coefficients b i , that equals C 2 . This can be achieved by choosing for instance b 1 = -1 12 , b 3 = -1 240 and b 2 = 1. Note that computing (5.13) requires one commutator only instead of two. In general, given B 1 , B 2 , . . . , B ν , terms of the form

ν-1 i=1 ν j=i+1 k i,j [B i , B j ] may be replaced by a single commutator, ⎡ ⎣ ν/2 i=1 b 2i-1 B 2i-1 , ν/2 j=1 b 2j B 2j ⎤ ⎦
provided that the coefficients {b i } 1≤i≤ν and {k i,j } 1≤i<j≤ν are compatible up to the order of the method.

Next, let us consider terms with double commutators. Let us focus on

-1 240 [B 2 , [B 1 , B 2 ]]. Since C 2 = -1 12 [B 1 , B 2 ] + O h 5 , evaluating the term -1 20 [B 2 , C 2 ] in place of -1 240 [B 2 , [B 1 , B 2 ]] amounts to the calculation of just one commutator. Note that -1 20 [B 2 , C 2 ] = -1 240 [B 2 , [B 1 , B 2 ]] + O h 7
, hence the approximation retains the odd-power of h expansion and the error introduced is subsumed in the local truncation error of the method.

Finally, let us consider the combination

1 360 [B 1 , [B 1 , B 3 ]] + 1 720 [B 1 , [B 1 , [B 1 , B 2 ]]].
(5.14)

Clearly, [B 1 , B 3 ] is not obtained from (5.13) using the linearity of the bracket, and needs to be taken into account. However, [B 1 , B 2 ] conforms with (5.13) and we can replace it with C 2 , introducing only odd-powered error in h which is also subsumed in the local truncation error. In summary, (5.14) can be replaced by the term

C 3 = [B 1 , [B 1 , 1 360 B 3 -1 60 C 2 ]]
which requires the computation of two commutators. Therefore, Θ can be computed in the form

Θ = B 1 + 1 12 B 3 + C 2 + C 3 , requiring four commutators only.
It is difficult to formalise the last two steps, involving two or more commutators. The reduction in the number of commutators has, in present state of knowledge, to be done on a case-by-case analysis.

Some examples of methods for linear problems obtained by means of graded algebras, made more economic with the technique of [START_REF] Blanes | Improved high order integrators based on Magnus expansion[END_REF], are described in Appendix A.2.

Nonlinear problems: collocation methods

Magnus and Fer expansions presented in Section 4 have been designed for linear problems, when the matrix function A depends on time only. When A = A(t, Y ), though, multivariate integrals appearing in (5.1) depend on A k = hA(c k h, Y (c k h)), namely also on the value of the unknown variable Y at quadrature points. Assume that the quadrature points c 1 , c 2 , . . . c ν obey the orthogonality conditions (5.5) for some m ≤ ν, hence the corresponding univariate interpolatory quadrature has order p = ν + m. Since each A k is a multiple of h, it is clear that the values A k can be replaced by hA(c k h, X k ), where X k is an approximation to Y (c k h) of order at least p -1. Following this point of view, it is possible to derive X k ≈ Y (c k h) with a numerical method of order p -1. For instance, using an RK-MK method of order p -1 and a Magnus (or Fer) expansion of order p is a Lie-group equivalent of a predictor-corrector method in classical numerical ODE theory.

A more elegant approach, due to [START_REF] Zanna | Collocation and relaxed collocation for the Fer and the Magnus expansions[END_REF], is to construct suitable approximants to X k s using directly the underlying principle of collocation methods. Proceeding as in Section 5.1, we replace the function A(t, Y ) by its Lagrangian interpolating polynomial

Ã(t, Y ) = h -1 ν k=1 ℓ k t h A k where A k = hA(c k h, X k ), k = 1, 2, . . . , ν.
The corresponding dexpinv equation is integrated by means of Picard iterations á la Section 4, obtaining the usual order trees for the Magnus and the Fer expansion. The main difference with the linear case is that now the same order trees are also employed to evaluate the integration coefficients for the internal stages X k . At each internal stage we need to calculate quadratures of the form (5.2),

K l (h) = j∈C ν s a l;j L(A j1 , A j2 , . . . , A js ),
where C ν s is the set of all combinations of length s from the set {1, 2, . . . , ν}. The integration weights are different from those in (5.2) and are defined as

a k;j = Sk s i=1 ℓ ji (ξ i )dξ,
where

Sk = {ξ ∈ R s : ξ 1 ∈ [0, c k ], ξ l ∈ [0, ξ m l ], l = 2, 3, . . . , s}
is the polytope S scaled to the [0, c k ]-cube instead of the unit cube. The weights b j are recovered by substituting c k = 1.

Theorem 5.5 Let c 1 , c 2 , . . . , c ν be ν collocation nodes and let p = ν + m, where m is the largest index such that (5.5) is satisfied. Assume that the Magnus or Fer expansion is truncated to include all trees of power q ≤ p for the evaluation of Y n+1 , and of power q ≤ p -1 for the intermediate stages. Then the resulting scheme has order p.

We sketch the main idea and refer the reader to [START_REF] Zanna | Collocation and relaxed collocation for the Fer and the Magnus expansions[END_REF] for details. The starting point is the Alekseev-Gröbner lemma, a nonlinear version of the variation of constants formula whose proof can be found in [START_REF] Nørsett | Perturbed collocation and Runge-Kutta methods[END_REF], stating that, if y is the solution of the differential equation y ′ = f (t, y) with initial condition y 0 , and if w(t) is a C 1 approximation to y such that w(t 0 ) = y 0 , then

y(t) -w(t) = t t 0 Φ(t, ξ, w(ξ))[f (ξ, w(ξ)) -w ′ (ξ)]dξ,
where Φ(t, ξ, w(ξ)) is the partial derivative of the the solution passing through (ξ, w(ξ)) with respect to the initial condition w(ξ). In the usual classical collocation setting for RK methods, w corresponds to the case when the function f is replaced by a collocation polynomial: at the nodes it is true that f (t n + c k h, w(t n + c k h))w ′ (t n + c k h)) = 0, hence the error reduces solely to quadrature error.

In our Lie-group setting the main difference consists in the fact that only the function A(t, Y ) is collocated, and not the whole right-hand side A(t, Y )Y . However, this can be viewed as a collocation method in the algebra g, where classical analysis remains valid.

A typical example of such collocation methods is the fourth-order scheme

X 1 = Y n , A 1 = hA(t n , X 1 ), X 2 = expm{ 5 24 A 1 + 1 3 A 2 -1 24 A 3 -1 2 ( 11 240 [A 1 , A 2 ] + 5 576 [A 1 , A 3 ] + 1 72 [A 2 , A 3 ])}Y n , A 2 = hA(t n + 1 2 h, X 2 ), X 3 = expm{ 1 6 A 1 + 2 3 A 2 + 1 6 A 3 -1 2 ( 2 15 [A 1 , A 2 ] + 1 30 [A 1 , A 3 ] + 2 15 [A 2 , A 3 ])}Y n , A 3 = hA(t n + h, X 3 ), Y n+1 = expm{ 1 6 A 1 + 2 3 A 2 + 1 6 A 3 -1 2 ( 2 15 [A 1 , A 2 ] + 1 30 [A 1 , A 3 ] + 2 15 [A 2 , A 3 ])}Y n
for n ∈ Z + , with the Gauss-Lobatto quadrature points c 1 = 0, c 2 = 1 2 and c 3 = 1 [START_REF] Zanna | On the Numerical Solution of Isospectral Flows[END_REF]. Note that the coefficients of the A i 's are the classical Runge-Kutta coefficients of Lobatto collocation scheme with the same quadrature points [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF], while the coefficients of the commutator terms [A i , A j ] are evaluated by integrating

a k;i,j -a k;j,i = c k 0 ξ1 0 ℓ i (ξ 1 )ℓ j (ξ 2 )dξ 2 dξ 1 - c k 0 ξ1 0 ℓ j (ξ 1 )ℓ i (ξ 2 )dξ 2 dξ 1 ,
a term corresponding to the power-three tree in the Magnus expansion (cf. Section 4).

A useful formula for evaluating the a k;i,j s corresponding to the power-three tree is given by {a

k;i,j } ν i,j=1 = a(c k ), k= 1, 2, . . . , ν, (5.15)
where a(θ) is the ν × ν matrix function of the scalar argument θ defined as

a(θ) = V -T T (θ)HJT (θ)V -1 , where J = diag(1, 1 2 , 1 3 , . . . 1 ν ), T (θ) = diag(θ, θ 2 , . . . , θ ν ), V is the Vandermonde matrix V = ⎡ ⎢ ⎢ ⎢ ⎣ 1 c 1 • • • c ν-1 1 1 c 2 • • • c ν-1 2 . . . . . . . . . 1 c ν • • • c ν-1 ν ⎤ ⎥ ⎥ ⎥ ⎦ ,
and finally H is the Hilbert matrix with entries H i,j = 1 i+j , i, j = 1, 2, . . . , ν.

This formula is reminiscent of the matrix representation of the standard Runge-Kutta matrix of a collocation method [START_REF] Nørsett | Perturbed collocation and Runge-Kutta methods[END_REF]. This is not a coincidence: the methods of [START_REF] Zanna | On the Numerical Solution of Isospectral Flows[END_REF] generalise the concept of collocation to the special multivariate integrals that occur in Magnus or Fer expansions.

Alternative coordinates

All the methods so far, whether applied to Lie groups or in a homogeneous-space setting, have been based on the exponential map. In other words, we represented the solution as an exponential (RK-MK and Magnus methods) or as a product of exponentials (Crouch-Grossman and Fer methods). It is entirely legitimate to query to which extent this renders such methods unduly expensive and non-competitive.

Sometimes the exact computation of the exponential is easy: a case in point is the application of Magnus expansions to the computation of Sturm-Liouville spectra in Section 11.2, since the exponential of an element in sl(2) can be evaluated exactly with great ease. In other cases the exponential can be replaced by a suitable approximation δ : g → G, δ(z) ≈ e z . This is the case with quadratic Lie groups: A

Lie group G is quadratic if G = {X ∈ GL(N ) : XP X T = P }, (6.1)
where P ∈ GL(N ) is a given matrix. Many Lie groups that appear in applications are of this kind, for example O(N ), Sp(N ) and O(N, M ). Moreover, some complex groups can be brought into this framework by replacing the transpose T by the Hermitian (i.e., conjugate) transpose H . The Lie algebra of the quadratic Lie group (6.1) is

g = {B ∈ gl(N) : BP + P B T = O} (6.2)
and it is easy to prove that δ : g → G whenever δ(z) = e γ(z) and γ is an odd function, analytic in the neighbourhood of the origin [START_REF] Celledoni | Approximating the exponential from a Lie algebra to a Lie group[END_REF]).

An important case occurs when γ(z) = log q(z)-log q(-z), where q is a polynomial, and it leads to rational functions δ(z) = q(z)/q(-z). In particular, diagonal Padé approximants to the exponential are of this form [START_REF] Baker | Essentials of Padé Approximants[END_REF]) and they can be applied very effectively in place of the exponential.

Yet, for some algebras there exists no analytic nonconstant function δ : g → G except for the (scaled) exponential. This, in particular, is the case with SL(N ) [START_REF] Kang | Volume-preserving algorithmms for source-free dynamical systems[END_REF]. In yet other cases, although we may replace the exponential with, say, a Padé approximant with impunity, the sheer size of the system renders this impractical when the number of variables is large. In that case there exist two possibilities. Firstly, we may endeavour to approximate the exponential of a matrix by some nonstandard means while keeping the outcome in a Lie group. This is the theme of Section 8. In the present section we consider another approach, which disposes of the dexpinv equation (2.46) altogether.

The research into 'alternative coordinates' is in a fairly preliminary stage and just two surrogates to the dexpinv equation have been identified so far, the Cayley transform and canonical coordinates of the second kind. We consider them in detail in the remainder of this section.

The Cayley transform and RK-Cayley methods

Let G be a quadratic Lie group. The main idea is to replace the exponential with the Cayley transform. Thus, given the Lie-group equation Y ′ = A(t, Y )Y , where A : R + × G → g, we seek a solution in the form

Y (t) = cay[∆(t)]Y 0 = [I -1 2 ∆(t)] -1 [I + 1 2 ∆(t)]Y 0 , t≥ 0, (6.3)
and at the first instance seek a differential equation for ∆. More generally, we may solve the homogeneous-space equation (2.26) replacing the exponential with the Cayley action but, for the sake of simplicity, the discussion is restricted to the 'straight' Lie-group case.

It is important to realise that our approach has no connection whatsoever with approximating the exponential. True, δ(z) = (1+ 1 2 z)/(1 -1 2 z) is a special case of a Padé approximant to the exponential, δ(z) = expm(z) + O z 3 , but this is entirely coincidental: as a matter of fact, we could have replaced, at the cost of slightly more complicated coefficients, the number 1 2 with an arbitrary nonzero constant. So far, everything is exact and no approximation has taken place.

It is an easy exercise, left to the reader, to ascertain that the function ∆ in (6.3) obeys the differential equation

∆ ′ = dcay -1 ∆ A(cay(∆)Y 0 , t) = A-1 2 [∆, A]-1 4 ∆A∆, t ≥ 0, ∆(0) = O. (6.4)
Note the presence of the term ∆A∆ in the above equation. In general, we cannot expect such a term to reside in g, but quadratic Lie algebras (6.2) are an exception. For future reference we note that the more general symmetric triple product [[D, E, F ]] 3 = DEF +F ED resides in g for all D, E, F ∈ g and quadratic Lie groups g. (The reason for this notation will be clear in Section 6.3.) Applying (6.2) thrice,

P (DEF ) T = P F T E T D T = -F P E T D T = F EP D T = -F EDP.
Therefore (DEF + F ED)P + P (DEF + F ED) T = (DEF + F ED)P -(F ED + DEF )P = O and indeed [[D, E, F ]] 3 ∈ g. This confirms that no illicit terms have crept into (6.4) and ∆ ′ evolves in g. (To be more precise, it evolves in Tg, except that the latter can be identified with g.)

The simplest implementation of (6.3) to quadratic Lie-group solvers is by employing Runge-Kutta methods in the Lie algebra, similarly to RK-MK, except that expm and dexp -1 in (3.4) need be replaced by cay and dcay -1 respectively. This has been accomplished systematically by [START_REF] Engø | On the construction of geometric integrators in the RKMK class[END_REF].

Cayley expansions

Proceeding like in Section 4.1 and subjecting (6.4) to Picard iteration, we observe that the solution ∆ can be expanded quite similarly to the Magnus expansion of the dexpinv equation (4.3),

∆(t) = t 0 A(ξ)dξ -1 2 t 0 ξ1 0 [A(ξ 2 ), A(ξ 1 )]dξ + 1 4 t 0 ξ1 0 ξ2 0 [[A(ξ 3 ), A(ξ 2 )], A(ξ 1 )]dξ 3 dξ 2 dξ 1 -1 4 t 0 ξ1 0 ξ1 0 A(ξ 2 )A(ξ 1 )A(ξ 3 )dξ 3 dξ 2 dξ 1 + • • • .
We seek to expand ∆, in greater generality, in a manner similar to (4.5),

∆(t) = ∞ k=0 τ ∈S k α(τ )D τ (t),
t≥ 0, (6.5)

where each D τ for τ ∈ S k is made out of exactly k + 1 integrals. (Note that, unlike (4.5), the 'exterior' integral is already included in the expansion term -intuitively speaking, D τ (t) = t 0 C τ (ξ)dξ. This makes the notation somewhat simpler and more transparent.) Following upon the construction of [START_REF] Iserles | On Cayley-transform methods for the discretization of Lie-group equations[END_REF], we identify three composition rules that are needed to assemble the terms in (6.5).

1 S 0 = {τ • }, and

D τ• (t) = t 0 A(ξ)dξ; 2 If τ 1 ∈ S k-1 , k ≥ 1, then there exists τ ∈ S k such that D τ (t) = t 0 t 0 D τ1 (ξ), A(ξ) dξ; ( 6 . 6 ) 3 If k ≥ 2 and τ 1 ∈ S k-j , τ 2 ∈ S j-1 for some 1 ≤ j ≤ k then there exists τ ∈ S k such that D τ (t) = t 0 D τ1 (ξ)A(ξ)D τ2 (ξ)dξ. (6.7)
Note that the outcome resides in the Lie algebra, as long as α(τ ) = α(τ ) where τ has been given in (6.7) and

D τ (t) = t 0 D τ2 (ξ)A(ξ)D τ1 (ξ)dξ
is the conjugate term of D τ (t). (The existence of such a term is assured by the third composition rule.)

Similarly to the association between rooted binary trees and terms in the Magnus expansion, we wish to render the structure of the above composition rules clearer by using graph theory. The presence of three, rather then two, composition rules makes this goal different and 'plain' rooted binary trees are no longer adequate for the task in hand. Instead, following [START_REF] Iserles | On Cayley-transform methods for the discretization of Lie-group equations[END_REF], we employ rooted bicolour binary trees: each vertex can be of one of two colours, black or white. The composition rules are interpreted in the following manner, borrowing as much as possible from the construction in Section 4.

1 S 0 = { r r } and r r ❀ t 0 A(ξ)dξ; 2 If S k-1 ∋ τ 1 ❀ D τ1 (t) then (6.6) corresponds to S k ∋ r r ❅ r τ 1 ❀ t 0 [D τ1 (ξ), A(ξ)]dξ; 3 Letting S k-j ∋ τ 1 ❀ D τ1 (t) and S j ∋ τ 1 ❀ D τ2 (t), (6.7) corresponds to S k ∋ r ❜ ❅ τ 1 τ 2 ❀ D τ (t) = t 0 D τ1 (ξ)A(ξ)D τ2 (ξ)dξ.
Unlike the case of the Magnus expansion, the derivation of the coefficients α(τ ) is straightforward and does not require any recursion. Given τ ∈ ∪ S k , we denote by γ(τ ) the number of white nodes therein. It is possible to prove that α(τ ) = (-1) k+γ(τ ) 2 -k and the outcome is the Cayley expansion

∆(t) = ∞ k=0 (-1) k 2 k τ ∈S k (-1) γ(τ ) D τ (t) ( 6 . 8 ) = r r -1 2 r r ❅ ❅ r r r + 1 4 r r ❅ ❅ r r r ❅ ❅ r r r -1 4 r ❜ ❅ ❅ r r r r -1 8 r r ❅ ❅ r r ❅ ❅ r r r ❅ ❅ r r r r r r ❅ ❅ r r ❜ ❅ ❅ r r r r + 1 8 r ❜ ❅ ❅ r r r ❅ ❅ r r r r + 1 8 r ❜ ❅ ❅ r r r r ❅ ❅ r r r + • • • .
Note that the last two trees above are conjugate and that they have the same weights. This is true in general, since if τ and τ are conjugate then they have the same number of white vertices, γ(τ ) = γ(τ ). Therefore, conjugate trees translate into (scaled) symmetric triple products and we stay safely within the Lie algebra. Absolute convergence of the Cayley expansion (6.8) was proved in (Iserles 1999b) for t ∈ (0, t * ), provided that t 0 ∥A(ξ)∥dξ 2, a result that can be somewhat improved for some norms and Lie algebras. However, inasmuch as convergence in norm with respect to ∥ • ∥ 2 is concerned, the Magnus-expansion condition (4.13) of Moan (n.d.) remains valid in the present setting.

Similarly to the case of Magnus expansions, it makes sense to truncate the series (6.8) by power,

∆(t) ≈ p-1 m=0 τ ∈Gm (-1) β(τ )+γ(τ ) 2 β(τ ) D τ (t), (6.9)
where β(τ ) + 1 is the number of integrals in D τ (in other words, τ ∈ S β(τ ) ), while G m stands for the set of trees of power m,

τ ∈ G m ⇐⇒ D τ (t) = O t m+1
for all sufficiently-smooth matrix functions A. The mechanism that allows m to exceed β(τ ) is subtly different from that of Magnus expansions (4.14) since, except for the second tree in the expansion (6.9), we never encounter an instance of

D τ (t) = [D τ1 , D ′ τ1 ]
. Instead, we say that a tree is basic if it has no black nodes with two children (equivalently, if the corresponding expansion term contains no commutators). The first few basic trees are

r r , r ❜ ❅ ❅ r r r r , r ❜ ❅ ❅ r r r ❜ ❅ ❅ r r r r , r ❜ ❅ ❅ r r ❜ ❅ ❅ r r r r r , ....
Each basic tree has an even number of vertices and it is easy to verify that if τ ∈ S 2m and τ is basic then

A(t) = A 0 + O(t) ⇒ D τ (t) = cA 2m+1 0 t 2m+1 = O t 2m+2
where c ̸ = 0 is scalar. In other words, τ ∈ G 2m and nothing is gained. However, as soon as we form the tree

τ = r r ❅ ❅ r τ 1
, where τ 1 is basic, it is trivial to notice that τ ∈ G 2m+2 , a 'gain' of one unit in power. Needless to say, this gain is inherited each time τ features as a component of a larger tree.

Truncating by power economises on the number of components: it has been proved in [START_REF] Iserles | On Cayley-transform methods for the discretization of Lie-group equations[END_REF]) that lim sup

k→∞ (#S k ) 1/k = 3 and lim sup m→∞ (#G m ) 1/m = 2.69805 . . .
(in either case there is substantial saving in comparison with the Magnus expansion, cf. (4.15)). However, a very important feature of Magnus expansions is unfortunately lost: The Cayley expansion (6.9), truncated by power, is no longer time symmetric! (We should perhaps emphasise that the role of time symmetry survives when the exponential is replaced with the Cayley transform: it still implies even order.) In other words, if we truncate the Cayley expansion (with all the integrals evaluated exactly) so that p = 3, say, in (6.9), the order will be just three.

Similarly to (4.16-18), we conclude by presenting Cayley expansions in standard notation, rather than in a tree terminology.

∆(t) =

[A(ξ 3 )A(ξ 2 )A(ξ 4 ), A(ξ 1 )]dξ + 1 8 t 0 ξ1 0 ξ2 0 ξ1 0 [A(ξ 3 ), A(ξ 2 )]A(ξ 1 )A(ξ 4 )dξ + 1 8 t 0 ξ1 0 ξ1 0 ξ3 0 A(ξ 2 )A(ξ 1 )[A(ξ 4 ), A(ξ 3 )]dξ + 1 16 t 0 ξ1 0 ξ2 0 ξ3 0 ξ4 0 [[[[A(ξ 5 ), A(ξ 4 )], A(ξ 3 )], A(ξ 2 )], A(ξ 1 )]dξ -1 16 t 0 ξ1 0 ξ2 0 ξ3 0 ξ1 0 [[A(ξ 4 ), A(ξ 3 )], A(ξ 2 )]A(ξ 1 )A(ξ 5 )dξ -1 16 t 0 ξ1 0 ξ1 0 ξ3 0 ξ4 0 A(ξ 2 )A(ξ 1 )[[A(ξ 5 ), A(ξ 4 )], A(ξ 3 )]dξ . . . . order 6 + • • • .
The above expansion underscores the importance of time symmetry in reducing the number of terms: compare the order-6 truncation with (4.18). We hasten to reassure the disappointed reader that not all is lost: time symmetry and even order will be regained in the next subsection.

Quadrature of the Cayley expansion and hierarchical algebras

In principle, the terms in the Cayley expansion (6.9) can be approximated exactly like 'Magnus integrals', since they are all consistent with (5.1): integrals of a multilinear form L over a polytope S. The theory of Section 5.1-2 is robust enough to cater for symmetric triple products, not just commutators.

Thus, to obtain a third-order method we truncate by power,

∆(t) ≈ r r -1 2 r r ❅ ❅ r r r -1 4 r ❜ ❅ ❅ r r r r , (6.10)
and replace integrals by quadrature using a self-adjoint basis. More specifically, we evaluate hA at the Gauss-Legendre points ( 1 2 ± quadratures are

r r : B 1 , r r ❅ ❅ r r r : -1 6 [B 2 , B 1 ], r ❜ ❅ ❅ r r r r
:

1 3 B 3 1 + 1 12 B 1 B 2 B 1 -1 24 B 2 B 2 1 -1 24 B 2 1 B 2
but we can throw away the last three terms with complete impunity: after all, we want a third-order method! The outcome is (6.11) Just to be on the safe side, we expand the solution, only to find that, lo and behold, the order of (6.11) is four. Not the order of (6.10), we hasten to say: the miracle has occurred just as the integrals have been replaced by quadrature! This is not a serendipitous coincidence. Quadrature recovers time symmetry, thereby boosting the order of an odd-order truncation (6.9) [START_REF] Iserles | On Cayley-transform methods for the discretization of Lie-group equations[END_REF]). Thus, herewith for example a sixth-order method, where B 0 , B 1 , B 2 have been obtained from order-six Gauss-Legendre quadrature:

B 1 + 1 12 [B 2 , B 1 ] -1 12 B 3 1 .
B 1 + 1 12 B 3 + 1 12 [B 2 , B 1 ] -1 12 B 3 1 -1 240 [B 3 , B 2 ] -1 240 [[B 3 , B 1 ], B 1 ] -1 240 [[B 2 , B 1 ], B 2 ] -1 48 B 1 B 3 B 1 -3 320 [B 2 , B 3 1 ] + 7 960 B 1 [B 2 , B 1 ]B 1 + 1 960 [[[B 2 , B 1 ], B 1 ], B 1 ] + 1 120 B 5 1 .
(6.12)

Having hopefully learnt something from our analysis of discretised Magnus expansions in Section 5, our next question is whether all the terms in (6.12) are necessary or can we perhaps replace some with linear combinations of other terms. In other words, we wish to repeat here the discussion from Section 5.3, except that in the present situation we should reckon with two operations: commutation and the symmetric triple product, the latter characteristic of quadratic Lie algebras. Wishing to derive the dimension of linear spaces of graded free algebras, along the lines of Section 5.3, a natural temptation is to express symmetric triple products in terms of commutators, but this soon leads to mushrooming complexity. More effective approach is described in [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF].

Let (g, +) be an Abelian group over a field of zero characteristic and introduce a countable family of m-nary operations

[[ • , . . . , • ]] m : m times g × g × • • • × g → g, m∈ N,
which is subject to the following three axioms:

1 Alternate symmetry: For all F 1 , F 2 , . . . , F m ∈ g

[[F 1 , F 2 , . . . , F m ]] m + (-1) m [[F m , F m-1 , . . . , F 1 ]] m = O; 2 Multilinearity: [[F 1 , F 2 , . . . , F m ]
] m is linear in each of its m arguments; 3 Hierarchy condition: For all F 1 , . . . , F m , E 1 , . . . , E s ∈ g and 1 ≤ l ≤ m it is true that [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF] a hierarchical algebra. It is easy to see that each hierarchical algebra is a Lie algebra (with the commutator defined as [

[[F 1 , . . . , F l-1 , [[E 1 , . . . , E s ]] s , F l+1 , . . . , F m ]] m = [[F 1 , . . . , F l-1 , E 1 , . . . , E s , F l+1 , . . . , F m ]] m+s-1 -(-1) s [[F 1 , . . . , F l-1 , E s , . . . , E 1 , F l+1 , . . . , F m ]] m+s-1 . The triple (g, +, {[[• • •]] m } m∈N ) has called by
• , • ] = [[ • , • ]] 2 , while each quadratic Lie algebra is hierarchical with [[F 1 , . . . , F m ]] m = F 1 F 2 • • • F m -(-1) m F m F m-1 • • • F 1 , F 1 , . . . , F m ∈ g.
We now proceed along the same path as Section 5.3, choose a set G 1 , G 2 , . . . , G ν of generators and define a free hierarchical algebra (FHA) similarly to Definition 5.1. We endow the generators with grading ω and extend it to FHA in a natural manner,

ω([[H i1 , H i2 , . . . , H ir ]] r = r k=1 ω(H k ).
One should not take the analogy with FLA too far, since FHA require subtly different approach. At the heart of the discussion of Section 5.3 it the fact that, using for example the Hall basis, we can express every element of a FLA as a linear combination of primitive terms of the form

[G i1 , [G i1 , [. . . , [G ir-1 , G ir ] • • • ]]].
In the case of FHA [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF] prove that the primitive 'building blocks' can be chosen of the form

[[G i1 , G i2 , . . . , G ir ]] r .
The method of proof is constructive, repeatedly using the three axioms. Anticipating future discussion, we exemplify it with one of the terms in the sixth-order Cayley expansion (6.12), assuming that ν ≥ 2:

[[[[[[G 2 , G 1 ]] 2 , G 1 ]] 2 , G 1 ]] 2 Axiom 3 = [[[[G 2 , G 1 ]] 2 , G 1 , G 1 ]] 3 -[[G 1 , [[G 2 , G 1 ]] 2 , G 1 ]] 3 Axiom 3 = ([[G 2 , G 1 , G 1 , G 1 ]] 4 -[[G 1 , G 2 , G 1 , G 1 ]] 4 ) -([[G 1 , G 2 , G 1 , G 1 ]] 4 -[[G 1 , G 1 , G 2 , G 1 ]] 4 ) Axiom 1 = [[G 2 , G 1 , G 1 , G 1 ]] 4 -3[[G 1 , G 2 , G 1 , G 1 ]] 4 .
Let g m be the set of all the grade-m elements in the FHA g. Similarly to FLA, we can express g as a direct sum of g m for m ∈ N. The dimension of each g m has been characterised in [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF]. Theorem 6.1 Let g be the graded FHA generated by S = {G 1 , G 2 , . . . , G ν }, with grades ω 1 , ω 2 , . . . , ω ν respectively. Denote by λ 1 , λ 2 , . . . , λ r the roots of the rth degree polynomial

p(z) = 1 - ν i=1 z ωi , r= max 1≤i≤ν ω i ,
and assume that they are all distinct. Then

σ2m = 1 2 r l=1 λ -m-1 l p ′ (λ l ) 2 -λ -m l -1 2 [p(λ 1/2 l ) + p(-λ 1/2 l )] σ2m+1 = 1 2 r l=1 λ -m-3/2 l p ′ (λ l ) -λ -m-1/2 l + 1 2 [p(λ 1/2 l ) -p(-λ 1/2 l )] ,
(6.13)

where σm = dim g m , m ∈ N.

The proof is long and technical. Its main step is in demonstrating that

∞ m=1 t m dim g m = 1 2 1 p(t) + p(t) p(t 2
) . Comparison with Tables 5.1-2 demonstrates that the dimension of graded subspaces of FHA grows larger than in the FLA case. This is hardly a surprise. A Lie algebra is closed with respect to just one binary operation, commutation, in addition to the usual linear-space operations. A hierarchical algebra, however, is closed with respect to a countable number of operations! On the face of it, there are infinitely more ways of forming terms in FHA. Fortunately, the hierarchy conditions mean that the operations are interconnected and the growth in dimension is not as bad as we might have expected.

The ratio of the dimensions ρm and σm from (5.12) and (6.13) respectively can be determined asymptotically. The dominant zero of polynomial p, λ 1 , say, is in (1, ∞) and simple. [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF] proved that

σm ρm = - λ 1 m 2p ′ (λ 1 ) [1 + o(1)], m≫ 1.
The method of proof of Theorem 6.1 is constructive and it naturally leads to a basis and to algorithmic means of its construction. We refer the reader to [START_REF] Iserles | On the dimension of certain graded Lie algebras arising in geometric integration of differential equations[END_REF] for details, here just presenting the results for ν = 2 and the grades ω(G i ) = i. The basis of g m is denoted by B m .

B 1 : {[[G 1 ]] 1 }, B 2 : {[[G 2 ]] 1 }, B 3 : {[[G 3 ]] 1 , [[G 1 , G 2 ]] 2 , [[G 1 , G 1 , G 1 ]] 3 }, B 4 : {[[G 1 , G 3 ]] 2 , [[G 1 , G 1 , G 2 ]] 3 , [[G 1 , G 2 , G 1 ]]} B 5 : {[[G 2 , G 3 ]] 2 , [[G 1 , G 2 , G 2 ]] 3 , [[G 1 , G 1 , G 3 ]] 3 , [[G 2 , G 1 , G 2 ]] 3 , [[G 1 , G 3 , G 1 ]] 3 , [[G 1 , G 1 , G 1 , G 2 ]] 4 , ]]G 1 , G 1 , G 2 , G 1 ]] 4 , [[G 1 , G 1 , G 1 , G 1 , G 1 [[ 5 }.
We conclude by going back to the sixth-order Cayley expansion (6.12) and representing it in the FHA basis. As before, we let G i = B i and ω i = i. The outcome, after long but straightforward algebra, is Note that, thanks to time symmetry, only odd-grade elements enter the expansion.

+ 1 12 [[B 3 ]] 1 -1 12 [[B 1 , B 2 ]] 2 -1 24 [[B 1 , B 1 , B 1 ]] 3 . . . . . . . . . . . . . . . . . . . . . . . . . . grade 3 + 1 240 [[B 2 , B 3 ]] 2 + 1 240 [[B 1 , B 2 , B 2 ]] 3 -1 240 ]]B 1 , B 1 , B 3 ]] 3 -1 240 [[B 2 , B 1 , B 2 ]] 3 -1 160 [[B 1 , B 3 , B 1 ]] 3 + 1 120 [[B 1 , B 1 , B 1 , B 2 ]] 4 -1 240 [[B 1 , B 1 , B 2 , B 1 ]] 4 + 1 240 [[B 1 , B 1 , B 1 , B 1 ,

Canonical coordinates of the second kind I: A naive approach

The main idea of the present section is to consider alternatives to the standard exponential map expm : g → G, which we can write in the form (6.14) where d = dim g and C = {C 1 , C 2 , . . . , C d } is a basis of the Lie algebra. The map (6.14) induces (at least locally, near the identity) a coordinate system in the Lie group which has been termed by [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF] the canonical coordinates of the first kind. An alternative to (6.14) (which, incidentally, explains why we have insisted to write it in such a strange form) are the canonical coordinates of the second kind (CCSK),

g ∋ d k=1 θ k C k ∼ (θ 1 , θ 2 , . . . , θ d ) → expm d k=1 θ k C k ∈ G,
g ∋ d k=1 θ k C k ∼ (θ 1 , θ 2 , . . . , θ d ) → e θ1C1 e θ2C2 • • • e θ d C d ∈ G. (6.15)
Why should we consider (6.15)? On the face of it, we have replaced a single exponential with d exponentials (and the whole exercise becomes really interesting when d ≫ 1!), hardly a sensible point of departure. However, as long as C is appropriately chosen, the computation of each expm(θ k C k ) can be exceedingly cheap and, moreover, the approach lends itself naturally to the exploitation of sparsity: as long as we can expect that there should be no component in the C k direction, say (or that it is suitably small), we can drop the relevant exponential from the product.2 A useful analogy is the distinction between Householder reflections and Givens rotations in numerical algebra.

It is possible to approach the issue of CCSK within the context of this survey from two distinctive points of view. Although ultimately they are closely related, they follow different philosophy, the first 'naive' and the other more mathematically sophisticated. This subsection is devoted to the more 'naive' approach, which associates CCSK with splittings.

Splitting methods have rich history throughout numerical analysis of differential equations and they are exceedingly useful in geometric integration, e.g. in the computational of Hamiltonian systems [START_REF] Serna | Numerical Hamiltonian Problems, Chapman & Hall. C. Størmer[END_REF]Calvo 1994, Yoshida 1990) and in the recovery of integrals and conservation laws [START_REF] Mclachlan | A unified approach to Hamiltonian systems, Poisson systems, gradient systems, and systems with Lyapunov functions and/or first integrals[END_REF]). Yet, the splitting of a flow into components corresponding to elements of a basis allows a significant enhancement of the technique. For simplicity, let us assume that we are solving the linear Lie-group equation (4.2), namely Y ′ = A(t)Y , t ≥ 0. We express the solution in the form (6.16) where θ 1 , θ 2 , . . . , θ d are scalar functions. It has been proved by [START_REF] Wei | On global representations of the solutions of linear differential equations as a product of exponentials[END_REF] that such functions always exist locally (and, in the case of solvable Lie algebras or for 2 × 2 real matrices, globally). The exact derivation of θ 1 , θ 2 , . . . , θ d is, needless to say, impossible in general, otherwise we could have written down the solution of (4.2) explicitly! Instead, we replace the θ k s with polynomials, which are chosen so as to match suitable order conditions at t = 0. Letting t = 0 in (6.16), we note that θ k (0) = 0, k = 1, 2, . . . , d. To obtain more useful order conditions we differentiate Y . Brief algebra confirms that

Y (t) = e θ1(t)C1 e θ2(t)C2 • • • e θ d (t)C d Y 0 , t≥ 0,
A(t) = Y ′ (t)Y -1 (t) = d k=1 θ ′ i (t)e θ1(t)C1 • • • e θ k-1 (t)C k-1 C k e -θ k-1 (t)C k-1 • • • e -θ1(t)C1 = d k=1 θ ′ k (t)Ad expm[θ1(t)C1] • • • Ad expm[θ k-1 (t)C k-1 ] C k .
(6.17)

Letting t = 0 in (6.17) we obtain the first-order condition

A(0) = d k=1 θ ′ k (0)C k .
Recalling that A(0) ∈ g, we can expand it in the elements of C and this yields θ ′ k (0) explicitly. Higher-order conditions can be obtained by differentiating (6.17) and massaging the formulae with a great deal of (fairly unpleasant) algebra. Thus, for example,

A ′ = d k=1 θ ′′ k Ad e θ 1 C 1 • • • Ad e θ k-1 C k-1 C k + d k=1 k-1 l=1 θ ′ k θ ′ l Ad e θ 1 C 1 • • • Ad e θ l C l [C l , Ad e θ l+1 C l+1 • • • Ad e θ k-1 C k-1 C k ]
and, letting t = 0, we have

d k=1 θ ′′ (0)C k = A ′ (0) - d k=1 k-1 l=1 θ ′ k (0)θ ′ l (0)[C l , C k ].
Recall that C is a basis of g, hence there exist scalars c j k,l such that

[C k , C l ] = d j=1 c j k,l C j , k,l= 1, 2, . . . , d.
They are called the structure constants of g and play an important role in the theory of Lie algebras [START_REF] Olver | Equivalence, Invariants, and Symmetry[END_REF][START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF]. Using structure constants and observing that A ′ (0) ∈ g can be expanded in elements of C, we obtain (6.18) hence second-order conditions.

d k=1 θ ′′ k (0)C k = A ′ (0) + d j=1 d k=1 k-1 l=1 θ ′ k (0)c j k,l θ ′ l (0)C j ,
Typically, the dimension d is quite large, for example dim so(N) = 1 2 (N -1)N and dim sl(N) = N 2 -1. Thus, in principle it might be costly to evaluate θ ′′ k (0), k = 1, 2, . . . , d in (6.18). Higher-order conditions are substantially costlier still. Yet, the cost can be reduced a very great deal by the right choice of the basis C.

The most suitable basis C is provided by a root-space decomposition of the (nonnilpotent) Lie algebra g. Deferring our discussion of this construct to the next subsection, we describe in a more nontechnical setting the special case of so(N).

To this end we choose the basis

C = {C k,l = e k e T l -e l e T k : 1 ≤ k < l ≤ N },
where e j ∈ R N is the jth unit vector. Note that

B ∈ so(N) ⇒ B = N -1 k=1 N l=k+1 b k,l C k,l
and that U = expm(tC k,l ) is a rigid rotation in the (k, l) plane: it coincides with the identity matrix, except for

u k,k u k,l u l,k u l,l = cos t sin t -sin t cos t .
Therefore, multiplying a matrix with expm(tC k,l ) is cheap. Moreover, it is easy to verify that

[C r,s , C k,l ] = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ C k,s , r= l, s ̸ = k, C l,r , r̸ = l, s = k, C r,k , r̸ = k, s = l, C s,l , r= k, s ̸ = l, O, otherwise, r < s, k < l,
where we identify C i,j with -C j,i for i > j. The condition (6.18) simplifies to

θ ′′ k,l (0) = a ′ k,l (0) - N i=1 a k,i (0)a i,l (0), 1 ≤ k < l ≤ N,
where

A(t) = N -1 k=1 N l=k+1 α k,l (t)C k,l
. This requires O N 3 flops altogether, in comparison with O N 6 if sparsity of structure constants is disregarded. Similar construction can be applied in other Lie algebras.

In principle, we can go on to derive θ (i) k (0) for i = 0, 1, . . . , p, but this procedure, even while utilising root-space decomposition, becomes progressively more expensive for larger orders p. Herewith we present a device which, to our knowledge, is new and which allows one to obtain order p while computing one less derivative. Observing that (6.16) is sensitive to the ordering of the basis, the main idea is to alternate the order of elements of C while time-stepping the numerical method. Thus, suppose that t m = mh and

Y 2n+1 = e θ2n,1(t2n+1)C1 e θ2n,2(t2n+1)C2 • • • e θ 2n,d (t2n+1)C d Y 2n , Y 2n+2 = e θ 2n+1,d (t2n+2)C d e θ 2n+1,d-1 (t2n+2)C d-1 • • • e θ2n+1,1(t2n+2)C1 Y 2n+1 ,
where θ m,k are p-degree polynomials. Without loss of generality, we assume that θ m,k are consistent with order p for m = 0, 1, . . . , 2n. Let

X(t) = e θ 2n+1,d (t)C d e θ 2n+1,d-1 (t)C d-1 • • • e θ2n+1,1(t)C1 Y 2n+1 , hence Y 2n+1 = X(t 2n+1 ) and Y 2n+2 = X(t 2n+2
). Repeatedly multiplying by inverted exponentials, we obtain

Y 2n+1 = e -θ2n+1,1(t)C1 e -θ2n+1,2(t)C2 • • • e -θ 2n+1,d (t)C d X(t)
Assuming that the θ 2n+1,k are chosen consistently with order p and letting t = t 2n we deduce that

Y 2n+1 = e -θ2n+1,1(t2n)C1 e -θ2n+1,2(t2n)C1 • • • e -θ 2n+1,d (t2n)C d Y 2n + O h p+1 .
In other words, we may take θ 2n+1,k (t 2n ) = -θ 2n,k (t 2n+1 ), k = 1, 2, . . . , d. This, together with the values of θ (i) 2n+1,k (t 2n+1 ), i = 0, 1, . . . , p-1, is just right to determine the θ 2n+1,k s consistently with order p.

Canonical coordinates of the second kind II: Admissible bases

The technique of canonical coordinates of the second kind can be enhanced a great deal at the cost of increased mathematical sophistication. The point of departure for our discussion, presently based on the important paper of Owren and Marthinsen (1999a), is the equation (6.17). We rewrite it in the form

d k=1 θ ′ k Ad expm(θ1C1) Ad expm(θ2C2) • • • Ad expm(θ k-1 C k-1 ) C k = A(t). (6.19)
Here θ 1 , θ 2 , . . . , θ d are known and we seek the scalars θ ′ 1 , θ ′ 2 , . . . , θ ′ d . Suppose that we have the means to solve (6.19) for arbitrary inputs θ 1 , θ 2 , . . . , θ d (in (6.17) we needed just θ = 0). This yields a differential equation

θ ′ = g(θ, A(t)), t ≥ 0, θ(0) = 0. (6.20)
Once the solution of (6.20) is known (or adequately approximated), we can use it to advance the solution of the Lie-group equation (4.2) through the CCSK representation (6.16). Moreover, the argument extends at once to nonlinear Lie-group equations Y ′ = A(t, Y )Y , where A : G × R + → g. Again, we represent the solution in the CCSK form (6.16), except that now

g = dccsk -1 θ A = g(θ, A(t, e θ1C1 e θ2C2 • • • e θ d C d Y 0 )) in the dccskinv equation (6.20).
Applying a Runge-Kutta method, say, to (6.20) results in time-stepping scheme that is guaranteed to respect Lie group structure. In effect, the only difference between the RK-MK methods of Section 3 and such methods is that, in place of canonical coordinates of the first kind and the dexpinv equation, they utilise canonical coordinates of the second kind and the dccskinv equation. All this motivates a thorough discussion of the problem how to invert an equation of the form (6.19).

Letting

F k = Ad expm(θ k C k ) , v k = θ ′ k , k = 1, 2, .
. . , d, we commence by writing (6.20) as

d k=1 v k F 1 F 2 • • • F k-1 C k = A.
Let P l be a projection on the trailing dl coordinates,

P l d k=1 c k C k = d k=l+1 c k C k ,
and set Fl = I -P l + P l A l , l = 1, 2, . . . , d. We can easily verify that

Fl C k = F l C k , l<k, C k , l≥ k.
Owren and Marthinsen (1999a) say that C is an admissible ordered basis (AOB) if for every θ 1 , θ 2 , . . . , θ d it is true that

F 1 F 2 • • • F k P k = F1 F2 • • • Fk P k , k= 1, 2, . . . , d -1. (6.21)
Provided that C is an AOB, it is simple to prove that our equation can be rewritten in the form (6.22) where

d k=1 v k F1 F2 • • • Fd-1 C k = F1 F2 • • • Fd-1 F = A,
F = d k=1 v k C k .
In the sequel we will see that in a number of important cases AOB implies that each Fk can be inverted very cheaply.

At a first glance, the AOB condition (6.21) is exceedingly demanding. Surprisingly, it is often achievable but we need to introduce a little bit more Lie-algebra theory before being in position to describe exactly how. The following brief extract should be ideally supplemented by perusing a Lie-algebra monograph: The book of [START_REF] Varadarajan | Lie Groups, Lie Algebras, and Their Representations[END_REF] is a good place to start.

• A subalgebra h of a Lie algebra g is an ideal if [h, g] ⊆ h.

• The Lie algebra g is solvable if there exists m ∈ Z + such that g (m) = {0}, where g (0) = g and g (i+1) = [g (i) , g (i) ] ⊆ g (i) . • The radical of g, denoted by Rad g, is the maximal solvable ideal in g. We say that the Lie algebra is semisimple if Rad g = {0}.

If definitions have become hazy by now, let us just point out that all specific Lie algebras in this survey (and in known applications within its framework) are semisimple, inclusive of sl(N), so(N) and sp(N). (All these three Lie algebras are, as a matter of fact, simple: their only ideals are {0} and the algebra itself.)

• An element in g is semisimple if all the roots of its minimal polynomial are distinct: in a matrix representation it means that the element can be diagonalised. • A subalgebra is toral if all its elements are semisimple. It is easy to see that every toral algebra must be abelian: in a matrix representation we can restate this by saying that the elements of the subalgebra share all eigenvectors, hence they commute.

• Unless g is nilpotent, it possesses a nonzero maximal toral subalgebra. Such subalgebra, which we denote by h, is unique up to an isomorphism. If g is a simple algebra, h is also known (subject to an equivalent definition) as a Cartan subalgebra.

• Suppose that g is a linear space over C. We denote by h * the dual space of a maximal toral subalgebra. It consists of all linear functionals h → C. The nonzero functional α ∈ h * is a root if there exists f ∈ g \ {0} such that

[h, f ] = α(h)f, h ∈ h.
In a matrix representation, α(H) is an eigenvalue of the commutator operator generated by H ∈ h, while F can be 'translated' into its eigenvector. • Denote the set of all roots of g by Φ. It is possible to prove that g can be subjected to the root-space decomposition

g = h ⊕ α∈Φ g α , (6.23) where g α = {f ∈ g : [h, f ] = α(h)f, h ∈ h} ̸ = {0}.
• The decomposition (6.23) motivates the choice of a Chevalley basis of the Lie algebra g: we choose one basis vector for each one-dimensional subspace g α , α ∈ Φ, and combine it with an arbitrary basis of h.

• There exists an integer k * ≥ 1 such that ad k * +1 h = 0 for every h ∈ g α , α ∈ Φ.

Many of the above concepts can be illustrated briefly with an example, and we choose sl(N, C). It is semisimple (as a matter of fact, we have already mentioned that it is a simple Lie algebra). Using the standard representation of sl(N, C) as matrices of zero trace, we can easily identify a maximal toral subalgebra h with diagonal zero-trace matrices. Setting E k,l = e k e T l , k, l = 1, 2, . . . , N, we may choose the basis {E k,k -E k+1,k+1 : k = 1, 2, . . . , N -1} for h. Moreover, given

h ∋ H = N k=1 h k E k,k , N k=1 h k = 0,
and letting h 0 = 0, we verify easily that

[H, E r,s ] = (h r -h r-1 -h s + h s-1 )E r,s , r,s= 1, 2, . . . , N, r ̸ = s.
Hence we identify the root α(H) = h rh r-1h s + h s-1 and construct our basis by placing there E r,s for every r ̸ = s and appending to this the above basis of h.

Note that this results in N 2 -1 terms, matching exactly the dimension of sl(N, C).

To determine k * we compute

ad 2 E r,s E k,l = E r,s , k= s, l = r, O, otherwise, r ̸ = s, k ̸ = l ⇒ ad 3 E r,s E k,l = O and ad 2 E r,s (E k,k -E k+1,k+1 ) = O. Therefore ad 3 E r,s = O, r ̸ = s,
and we deduce that k * = 2.3 Theorem 6.2 Let {ρ 1 , ρ 2 , . . . , ρ d * }, where d * = ddim h, be the set of roots Φ of a semisimple non-nilpotent Lie algebra g. Suppose that the Chevalley basis C is ordered so that the basis of h comes last. Then this basis is AOB if

kρ i + ρ j = ρ m , m < i < j ≤ d * , 1 ≤ k ≤ k * ⇒ ρ m + ρ n ̸ ∈ Φ ∪ {0} (6.24)
for all n = m + 1, m + 2, . . . , i -1. Moreover, in that case (6.25) Returning to so(N, C), it has been proved in (Owren and Marthinsen 1999a) that conditions of Theorem 6.2 are satisfied as long as super-diagonal elements are ordered lexicographically by rows in front of the elements underneath the diagonal, which are ordered lexicographically by columns: thus, the ordered basis is

F -1 k = I + k * l=1 (-1) l θ l k l! ad l C k P k , k= 1, 2, . . . , d -1.
{E k,l : 1 ≤ k < l ≤ N } ∪ {E k,l : 1 ≤ l < k ≤ N } ∪ {E k,k -E k+1,k+1 : 1 ≤ k ≤ N -1}.
We omit the largely-technical proof. Likewise, it is possible, using Theorem 6.2, to identify AOB of sp(N, C) and so(N, C). This, however, is probably of less importance than in the case of so(N, C), since the increase in dimension due to the replacement of R with C leads to a significant increase in the volume of computations. In the present stage of the development of Lie-group methods it is fair to say, we believe, that the Cayley-transform-based techniques from Section 6.1-2 are the method of choice for quadratic Lie groups, while CCSK should be used with the special linear group.

Adjoint methods

In the previous sections we have encountered a number of numerical integrators for Lie groups. Although such methods produce solutions that stay on a given Lie group G by design, it is of interest to study how well such schemes respect other qualitative features of the underlying equations: the retention of a symplectic form, Lie-Poisson structure, conservation of energy, time symmetry, time reversibility et cetera. Given the novelty of the proposed schemes, many of the above features and their implications on the 'quality' of the solution are still under investigation (Engø and Faltinsen 1999, Faltinsen n.d.). For this reason, we shall focus here just on time symmetry for Lie-group methods, which is at present one of few features that are better understood, deriving adjoint and self-adjoint Lie-group methods. Before proceeding further, let us recall that the flow Φ of the differential equation in R N y ′ = f (t, y), t ≥ t 0 , y(t 0 ) = y 0 , defined as Φ(t, t 0 , y 0 ) = y(t) obeys the following conditions:

i. Φ(t 0 , t 0 , y 0 ) = y 0 , ii. Φ(t + τ, t 0 , y 0 ) = Φ(τ, t, Φ(t, t 0 , y 0 )), provided that the above function f is Lipschitz with respect to y [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]. In particular, the second condition implies that

Φ(-τ, t + τ, Φ(τ, t, y(t))) = y(t),
a condition that in literature is mostly known as time symmetry or self adjointness of the exact flow Φ.

Adjoint methods in classical setting

Numerical methods usually approximate the flow Φ by a discrete flow, say Ψ, so that

y n+1 = Ψ(t n + h, t n , y n ), n∈ Z + ,
approximates the exact solution y(t n + h) to given order p. Although numerical integrators for ODEs always obey the condition Ψ(t n , t n , y n ) = y n , they usually fail to satisfy condition ii. However, its weaker variant, time symmetry, is easier to impose and numerical methods such that

Ψ(-h, t n + h, Ψ(h, t n , y n )) = y n , n∈ Z + ,
are usually called self-adjoint or time-symmetric methods. If a method is not self adjoint, its adjoint Ψ * is defined as the map

Ψ * (-h, t n + h, Ψ(h, t n , y n )) = y n .
In shorthand notation, we write Ψ h for Ψ(h, t n , • ), so that Ψ * h denotes the adjoint method of Ψ h and moreover Ψ * h • Ψ h = Id, the identity map, or equivalently

Ψ -h = Ψ -1
h , if and only if the method is self adjoint. The theory of adjoint and self-adjoint numerical methods for ODEs in R N is well established and we refer the reader to [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]) for further reading.

One might question why self adjointness of a numerical integration scheme is desirable. It turns out that for numerical integration schemes that are self adjoint it is possible to develop a theory analogous to the KAM theory for Hamiltonian and time-reversible problems [START_REF] Moser | Stable and Random Motion in Dynamical Systems[END_REF], which usually implies better approximation of the solution and slower accumulation of error over long integration intervals [START_REF] Estep | The Rate of Error Growth in Hamiltonian-Conserving Integrators[END_REF]Stuart 1995, Reich 1996).

With regards to Lie-group methods and time symmetry, we have already shown in Section 4 that the Magnus expansion truncated by power is time symmetric when applied to linear Lie-group differential equations Y ′ = A(t)Y , provided that the underlying quadrature is based on quadrature nodes in [0, 1] which are symmetric with respect to 1 2 . A similar result applies also to the RK-MK methods, provided that the underlying Runge-Kutta scheme is self adjoint.

However, it can be easily verified by means of numerical experiments that the above-mentioned Lie-group methods are not self adjoint for nonlinear problems. Hence, using self-adjoint methods to solve a Lie-algebra differential equation is not sufficient to derive self-adjoint Lie-group methods!

In this section we shall be discussing a more general procedure, derived by [START_REF] Zanna | Adjoint and selfadjoint Liegroup methods[END_REF], that allows us to construct self-adjoint Lie-group methods for linear and nonlinear problems alike and for all type of coordinate maps φ : g → G that we may use to represent the solution. In lifting the Lie-group equation from G to g, we make an implicit choice of a coordinate map, which has to be taken into account in the construction of self-adjoint methods.

Coordinate maps centred at arbitrary points

Schematically, the first step in the development of Lie-group schemes introduced in Sections 3-6 is the choice of a smooth map, say φ : g → G, such that φ(O) = I, the identity of the Lie group G, and φ ′ (O) = I (more precisely Tφ(O, B) = B, where B ∈ g). In other words, φ is a diffeomorphism mapping a neighbourhood of O ∈ g into a neighbourhood of I in G. Thus, φ(B) = expm(B), B ∈ g, is an example of such a map, which in (6.14) we have termed canonical coordinates of the first kind. Similarly, we might consider canonical coordinates of the second kind (6.15), namely

φ(B) = ccsk(B) = expm(β 1 B 1 ) expm(β 2 B 2 ) • • • expm(β d B d ), B ∈ g, where B = d i=1 α i B i , α i s
being real coefficients, the B i s are basis elements of the algebra g, and the β i s real functions of α 1 , α 2 , . . . , α d . Yet another example of this kind of map is the Cayley transform (6.3), which in the current formalism reads

φ(B) = (I -1 2 B) -1 (I + 1 2 B), B ∈ g,
that maps g into G whenever G is a quadratic group and g is its corresponding quadratic algebra. Secondly, the Lie-group differential equation

Y ′ = A(t, Y )Y, t ≥ 0, Y(0) = Y 0 ∈ G, (7.1)
is lifted by means of the inverse of the map dφ to an ordinary differential equation in g. Thus, for coordinates of the first kind, one has dφ -1 = dexp -1 , the dexpinv equation (3.2) that we have already encountered time and again in the course of the present article. Similarly, we have derived the expressions dcay -1 and dccsk -1 in Section 6.

Finally, the dφ -1 equation is solved in g with either a Runge-Kutta method or with a Magnus or a Cayley-type expansion. Assuming that an approximation Y n ∈ G has been already derived, we typically solve

Θ ′ = dφ -1 Θ (A(t, φ(Θ)Y n )), Θ(t n ) = O, t ∈ [t n , t n+1 ],
where

t n+1 = t n + h. The choice of the initial condition Θ(t n ) = O is equivalent to 'centring' the coordinate map φ at Y n . Instead, coordinates centred at any point X ∈ G can be obtained inverting the map B ∈ g → φ(B)X ∈ G.
In the general case we write

X = φ(C) -1 Y n ,
for some C ∈ g to be specified in the sequel. Note that φ(C) -1 is the inverse (in G) of the group element φ(C). Before proceeding further, we observe that both canonical coordinates of the first kind and the Cayley transform obey the relation

φ(B) -1 = φ(-B),
for all B ∈ g. For canonical coordinates of the second kind one has instead

ccsk(B) -1 = expm(-β d B d ) expm(-β d-1 B d-1 ) • • • expm(-β 1 B 1 ).
We seek a solution of (7.1) of the form

Y (t) = φ(Θ(t))φ(C) -1 Y n (7.2) for t ∈ [t n , t n+1 ].
Differentiating in the usual fashion we obtain a differential equation for Θ(t)

Θ ′ = dφ -1 Θ (A(t, Y )), t∈ [t n , t n+1 ], (7.3)
where Y is as in (7.2), in tandem with the initial condition

Θ(t n ) = C. (7.4)
Thus, changing the centre of coordinate map does not affect the differential equation obeyed by Θ, just its initial condition.

The adjoint of Lie-group methods

Assume that the Lie-algebra differential equation ( 7.3), with the initial condition (7.4), is computed with a numerical method Ψ h , and denote by Ψ * h its adjoint in the classical sense of Section 7.1, namely

(Ψ * -h • Ψ h )B = B
for all B ∈ g. The corresponding Lie-group method is such that

Y n+1 = Ψh Y n = φ(Υ h,n+1 )φ(C h,n ) -1 Y n Υ h,n+1 = Ψ(h, t n , C h,n ),
where C h,n is the initial condition of Θ in the interval [t n , t n+1 ], and we allow it to depend on the interval of integration and on the step size h. In order to obtain the adjoint of the Lie-group method Ψh , we need not just to use Ψ * h in g, but also to make sure that the coordinate map employed while stepping forward with the method Ψh is the same as when stepping backward with the adjoint method Ψ * h . Define a pair of methods Ψ and Ψ * on G as

Ψ(t n + h, t n , Y n ) = φ(Υ n,n+1 )φ(C h,n ) -1 Y n , Ψ * (t n + h, t n , Y n ) = φ(Θ * h,n+1 )φ(C * h,n ) -1 Y n , C * -h,n+1 = Υ h,n+1 , (7.5)
where

Υ h,n+1 = Ψ(h, t n , C h,n ) and Θ * h,n+1 = Ψ * (h, t n , C * h,n ).
C h,n = ϑ(F 1 , F 2 , . . . , F ν ) = -1 2 ν i=1 b i F i .
Flow-symmetric coordinates are instead more naturally defined for methods based on collocation. If ℓ i (x), i = 1, . . . ν, denote the familiar cardinal polynomials of Lagrangian interpolation, already introduced in Section 5.1, we set

w i = 1 2 0 ℓ i (τ ) dτ, i = 1, 2, . . . , ν, (7.7)
and the choice

C h,n = ϑ(F 1 , F 2 , . . . , F ν ) = - ν i=1 w i F i .
corresponds to the flow midpoint.

There exist also other choices of functions ϑ that obey (7.6) and it is possible to show that the set of such functions is convex. See [START_REF] Zanna | Adjoint and selfadjoint Liegroup methods[END_REF] for further examples of coordinate maps that yield self-adjoint schemes.

To conclude this section, we illustrate with a numerical experiment the benefits of using the geodesic and flow-symmetric coordinates instead of classical coordinates centred at Y n . We consider the Euler equations for the rigid body,

y ′ = y × M y, t ≥ 0, y(0) = y 0 , (7.8)
where y ∈ R 3 (we assume that ∥y 0 ∥ 2 = 1), the symbol '×' denotes the classical vector product os R 3 and M is a diagonal matrix, M = diag (m 1 , m 2 , m 3 ). This system has the Hamiltonian function

H(y) = 1 2 (m 1 y 2 1 + m 2 y 2 2 + m 3 y 2 3
) and obeys ∥y∥ 2 = 1. It can be represented by means of Lie-group action of SO(3) on R 3 by representing the solution y(t) as Ω(t)y n , n = 0, 1, 2, . . ., with Ω ∈ SO(3), t ∈ [t n , t n+1 ]. Hence, in each interval [t n , t n+1 ] we solve the differential equation

Ω ′ (t) = A(y(t))Ω(t), t ≥ t n , Ω(t n ) = I, (7.9)
where

A(y(t)) = - ⎡ ⎣ 0 -m 3 y 3 m 2 y 2 m 3 y 3 0 -m 1 y 1 -m 2 y 2 m 1 y 1 0 ⎤ ⎦ . 93 3 5
In this numerical experiment, m 1 =1,m 2 = 1 and m 3 = 1 and h = t n

+1

t n = 1 10 , while the initial condition is a random 3-vector with unit norm. We remark that such action automatically obeys the homogeneous-space condition ∥y∥ 2 = 1 whenever a Lie-group method is applied to (7.9). Although generally the two latter choices would correspond to different error, in this case the errors are very similar. The bottom plot corresponds to the difference between the errors in geodesic and flow-symmetric coordinates.

We compare an RK-MK method of order four based on Gauss-Legendre quadrature, using coordinates centred at y n with geodesic and flow-symmetric coordinates. The error in the Hamiltonian function, evaluated as

err H = H(y n ) -H(y 0 ), is displayed in Figure 7.2.
Similar comparison is displayed in Figure 7.3 for a method based on Magnus expansion of order four.

Precise details of numerical schemes used in the above example can be found in Appendix A.

Computation of exponentials

8.1. Six dubious ways to compute the exponential of a matrix in a Lie algebra Most (but by no means all) Lie-group methods require repeated calculation of exponentials or, in a more realistic setting, an approximation of exponentials. In principle, this is a well-tried and familiar task in numerical analysis and can be accomplished in one of several ways: rational approximants (Baker 1975, Iserles and[START_REF] Iserles | Order Stars[END_REF], Krylov-subspace methods [START_REF] Hochbruck | On Krylov subspace approximations to the matrix exponential operator[END_REF], Schur decomposition [START_REF] Golub | Matrix Computations[END_REF] and so on. Although such methods have occasionally attracted healthy scepticism [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix[END_REF], it is fair to say that they have a distinguished track record across numerical analysis. However, as we have already commented in Section 6, our present task is subject to a crucial restriction: Our approximant must map the Lie algebra g to the Lie group G! Low-dimensional algebras are easy and often we can evaluate the exponential explicitly. In particular, the following two cases are of practical importance.

• Firstly, given A = a b c -a ∈ sl(2), we can easily establish that e A = cosh ω I + sinh ω ω A, where ω = a 2 + bc. (8.1)
This will be of use in Section 11, in our discussion of the application of Magnus expansions to the calculation of Sturm-Liouville spectra. • Secondly, the exponential of

A = ⎡ ⎣ 0 a b -a 0 c -b -c 0 ⎤ ⎦ ∈ so(3) is I + sin σ σ A + 1 -cos σ σ 2 A 2 , where σ = √ a 2 + b 2 + c 2 .
Given the number of spatial dimensions in our universe, it will come as little surprise that many useful equations, e.g. (7.8), can be formulated in so(3). We will return to the above expression, known as the Rodrigues formula, in Appendix B. These, however, are the exceptions.

An exact formula being unavailable, an appealing alternative is to compute the exponential to machine accuracy. This, however, is neither affordable not always reliable. The Matlab function expm computes the exponential by scaling-andsquaring a diagonal Padé approximant: the procedure is very expensive for large dimensions and the outcome often falls short of machine accuracy and is subject to fast error accumulation.

When neither an explicit formula nor computation to machine accuracy are feasible, we must resort to approximation. Any such procedure must conform with two conditions: The outcome lies in the correct Lie group G and it departs from the exact exponential only to an extent consistent with the order of the Lie-algebra method. A standard means of approximation is to replace e z with a function r(z), analytic in the neighbourhood of the origin. The action of such function can be extended from C to gl(N), hence to any matrix Lie algebra, by elementary means.

Our two desiderata can be now reformulated by requiring that r(z) = e z +O z p+1 , where p ≥ 1 is the order of the Lie-algebra time-stepping procedure, and r(g) ⊆ G.

As we have already mentioned in Section 6, above conditions might be much too restrictive. As has been proved by [START_REF] Kang | Volume-preserving algorithmms for source-free dynamical systems[END_REF], the only analytic function that maps sl(N) to SL(N ) for every N ∈ N and takes zero to identity is r(z) = e αz for α ∈ R. Requiring consistency means that we must choose the exact exponential! On the other hand, in a quadratic Lie algebra we are faced with an abundance of riches: given an arbitrary odd function f , analytic about the origin, it is true that e f (g) ⊆ G [START_REF] Celledoni | Approximating the exponential from a Lie algebra to a Lie group[END_REF]. In particular, this is the case with all diagonal Padé approximants,

r(z) = p m (z) p m (-z) , where p m (z) = m k=0 m k (2m -k)! (2m)! z k , m ∈ N.
Yet, all this is of lesser utility since, arguably, the method of choice for quadratic Lie algebras rests upon the use of the Cayley transform (i.e., the diagonal Padé approximant with m = 1) as an alternative action, thereby avoiding altogether the need to approximate the exponential function! Yet another option is to evaluate the exponential with a Krylov-subspace method.

Assuming for simplicity that we wish to approximate e A v, where A ∈ gl(N)

and v ∈ R N , such techniques choose the approximant from the space K N,M = span {v, Av, . . . A M -1 v}. Surprisingly small values of M produce remarkably good and affordable approximants [START_REF] Hochbruck | On Krylov subspace approximations to the matrix exponential operator[END_REF]). Yet, there is absolutely nothing in this approach to guarantee that the outcome resides in the correct Lie group.

The last (and perhaps the most obvious) candidate for our list of alternatives to the matrix exponential is projection. For example, to travel from sl(N) to SL(n), we may employ a diagonal Padé approximant. The outcome, V = r(A), say, cannot be expected to reside in the special linear group. However, replacing V with V /(det V ) 1/N produces an element in SL(n). Unfortunately, experience tells that this procedure is prone to instability (cf. Section 11.4).

Splitting methods

Let A ∈ g. In the spirit of [START_REF] Celledoni | Approximating the exponential from a Lie algebra to a Lie group[END_REF], we wish to approximate

e tA ≈ R(tA) = e tB1 e tB2 • • • e tBs , (8.2)
where the matrices B = {B 1 , B 2 . . . , B s } are subject to the following requirements.

1 Each B l resides in g;

2 It is cheap to evaluate expm(tB l ) ∈ G exactly for each l;

3 It is cheap to multiply exponentials in (8.2); 4 The error is suitably small, R(tA) = e tA + O t p+1 .

We choose the splitting B so that it consists of low-rank matrices.

Let K and L be two N × r matrices, where r ≥ 1 is small, and assume that C = B l = KL T . Then e tC = I + tKD -1 (e tD -I)L T , where D = L T K [START_REF] Celledoni | Approximating the exponential from a Lie algebra to a Lie group[END_REF]. Note that D is just r × r and the cost of evaluating e tC exactly is modest for small values of r.

As an example, let us consider g = so(N). We let r = 2, s = N -1, set

B [0] = A = [b [0] 1 , b [0] 2 , . . . , b [0] N ], and choose B 1 = b [0] 1 e T 1 -e 1 b [0] 1 T ∈ so(N), where e k ∈ R N is the kth unit vector. Letting B [1] = B [0]
-B 1 , we observe that its first row and column vanish. We continue in a manner similar to LU factorisation, letting

B [i] = B [i-1] -B i and
B i = b [i-1] i e T i -e i b [i-1] i T ∈ so(N), i= 1, 2, . . . , N -1.
We refer to [START_REF] Celledoni | Approximating the exponential from a Lie algebra to a Lie group[END_REF] for precise estimation of cost, implementation details and a similar example for sl(N), as well as for an example of a splitting, again with r = 2, that eliminates two rows and columns of B ∈ so(N) at a time.

Although the number of exponentials in (8.2) is generally quite large for low-rank splittings, the underlying linear algebra carries a reasonable price tag.

The main disadvantage of low-rank splitting methods is the quality of approximation: in general, we can expect order p = 1. The second-order condition is

s-1 k=1 s l=k+1 [B k , B l ] = O
and it is easy to verify that it is satisfied when (8.2) is the Strang splitting: s = 2s+1 and B s+i = B s-i , i = 1, 2, . . . , s. In principle, it is easy to convert any low-rank matrix so that it becomes a Strang splitting by first approximating 1 2 A with a firstorder splitting (8.2), next approximating the same matrix with the same splitting but with the matrices B l arranged in reverse order and finally 'aggregating' the middle two terms. The outcome, S(tA) = e tB1/2 e tB2/2 • • • e tBs-1/2 e tBs e tBs-1/2 • • • e tB2/2 e tB1/2 (8.3) costs twice as much as (8.2) but it has a crucial advantage: it is not just secondorder but also time symmetric. This renders (8.3) amenable to the application of the Yoshida device (Sanz [START_REF] Serna | Numerical Hamiltonian Problems, Chapman & Hall. C. Størmer[END_REF]Calvo 1994, Yoshida 1990). Thus, the function S(αtA)S((1 -2α)tA)S(αtA), where

α = 2 3 + 3 √ 2 3 + 3 √ 4 6
approximates e tA to order four. Similar procedure can we used to increase the order further in increments of two.

A most welcome feature of low-rank splittings is that they can be implemented to take advantage of sparsity. Provided that A is banded, say, all the computations can be confined to the relevant band and sparsity is inherited as we are 'mopping up' rows and columns similarly to the so(N) algorithm above.

Canonical coordinates of the second kind

Our point of departure is similar to the reasoning behind the CCSK representation (6.16). Again, C = {C 1 , C 2 , . . . , C d } is a basis of the Lie algebra g and we seek polynomials θ 1 , θ 2 , . . . , θ d so that e θ1(t)C1 e θ2(t)C2 • • • e θ d (t)C d = e tA + O t p+1 .

(8.4)

C being a basis, there exist scalars a 1 , a 2 , . . . , a d so that A = d k=1 a k C k . Letting θ k (t) = a k t, k = 1, 2, . . . , d, gives us a first-order splitting (8.2). With greater generality, we set θ k (0) = 0, θ ′ k (0) = a k , k = 1, 2, . . . , d, to guarantee p ≥ 1. To obtain higher-order conditions in (8.4) we proceed like in Section 6.4. Differentiation and further algebra produce, equivalently to (6.17), the equation

A = d k=1 a k C k = d k=1 θ ′ k (t)Ad expm[θ1(t)C1] • • • Ad expm[θ1(t)C1] C k + O(t p ) . (8.5)
We go on differentiating (8.5) and letting t = 0. This yields order conditions, which need be unscrambled by further algebra, exploiting the structure constants (cf. Section 6.4 for the definition) of C. The outcome is [START_REF] Celledoni | Methods for the approximation of the matrix exponential in a Lie-algebraic setting[END_REF]. The cost of this procedure is, at the first glance, prohibitive. We might just about get away with computing order-two conditions at the cost of O d 3 operations, but a price tag of O d 5 operations for order three, to say nothing of higher orders, is out of the question. This naive impression is misleading, since we are absolutely free to exploit, along the lines of Sections 6.4-5, sparsity in structure constants in a serendipitously-chosen basis C. 4 Again, Chevalley bases present the best choice, as well as leading to very easy computation of exponentials in (8.4). For example, choosing the basis {E k,l -E l,k : 1 ≤ k < l ≤ N } of so(N), an order-two approximation is attained by letting

p ≥ 2 : θ ′′ k (0) = d l=1 l-1 j=1 a l c k l,j a j , k= 1, 2, . . . , d p ≥ 3 : θ ′′′ k (0) = 2 d l=1 l-1 j=1 c k l,j [θ ′′ l (0)a j + a l θ ′′ j (0)] + 2 d l=1 l-1 j=1 j-1 i=1 d m=1 c m l.j c k i,m a l a j a i + d l=1 d j=1 l-1 i=1 c j l,i c k i,j β l β 2 i , k= 1, 2, . . . , d (
θ k,l (t) = a k,l t + 1 2 N i=1 a k,i a i,l t 2 , 1 ≤ k < l ≤ N,
an O N 3 = O d 3/2 procedure. In the case of sl(N) the approach advocated in [START_REF] Celledoni | Methods for the approximation of the matrix exponential in a Lie-algebraic setting[END_REF] is to choose the ordered basis C = C 1 ∪ C 2 , where

C 1 = {E k,l : k ̸ = l}, C 2 = {E k,k -E k+1,k+1 : k = 1, 2, . . . , N -1},
where each set is ordered lexicographically. Let

A = N k,l=1 k̸ =l a k,l E k,l + N -1 k=1 b k (E k,k -E k+1,k+1 )
and denote the coefficients corresponding to terms in C 1 and C 2 by θ k,l and η k respectively. The second-order conditions are

θ k,l (t) = a k,l t + 1 2 k-1 i=1 a k,i a i,l + a k,l (b k-1 -b k + b l -b l-1 ) - N i=k+1 a k,i a i,l t 2 , k̸ = l, η k (t) = b k t -1 2 k i=1 N j=k+1 a i,j a j,i t 2 , k= 1, 2, . . . , N -1, again just O N 3 = O d 3/2 operations.
An intriguing aspect of approximants based on CCSK is that they might well be ideally suitable to handle sparsity in the matrix A. Although this issue is by no means fully understood, there are enough encouraging pointers to justify a brief discussion. One mechanism that exploits sparsity is that the latter can be taken into account in the evaluation of the θ k s because of the association between terms in a Chevalley basis and the entries of A. For example, the second-order conditions for a tridiagonal matrix A ∈ sl(N) reduce to

θ k,l (t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 2 a k,k-1 a k-1,k-2 t 2 , l = k -2, a k,k-1 t -1 2 a k,k-1 (b k -2b k-1 + b k-2 )t 2 , l= k -1, a k,k+1 t + 1 2 a k,k+1 (b k+1 -2b k + b k-1 )t 2 , l= k + 1, -1 2 a k,k+1 a k+1,k+2 t 2 , l = k + 2, 0, |k -l| ̸ = 1, η k (t) = b k t -1 2 a k,k+1 a k+1,k t 2 .
This entails just O(N ) operations and, equally importantly, just O(N ) terms survive in the product (8.4). The cost scales with the number of nonzero elements in the matrix, rather than with dimension, a hallmark of a good method for sparse matrices.

Another mechanism is of a more tentative value, yet we believe that it deserves mentioning. Even if the matrix A is sparse, its exponential is dense. However, it is possible to prove, at least in the case of banded matrices, that most of the entries are exceedingly small and that the loci of large elements are predictable (Iserles 1999a). As an example, we have averaged 1000 exponentials of 100 × 100 matrices with the cruciform sparsity pattern displayed in Figure 8.1a and with random elements uniformly distributed in (-1, 1) and normalised so that max k,l |a k,l | = 1. Let W be the average of all the exponentials. The matrix is dense, yet most of the elements of W are tiny! Thus, Figure 8.1b displays the sparsity pattern of W without all its entries that are smaller than 10 -6 in magnitude. Although the cruciform shape 'swells', most of the matrix consists of zero entries. This observation is affirmed in Figure 8.1c, where we have plotted the matrix log 10 |W |. The vertical axis tells the magnitude of the entries in terms of significant (decimal) digits. The decay outside the original cruciform shape is evident. Using upper bounds from (Iserles 1999a), it is possible to say how fast entries decay for banded A but computer experiments (and, indeed, Figure 8.1) indicate that similar behaviour takes place for more exotic sparsity patterns.

Choosing a Chevalley basis, members of C mostly correspond to elements of A. Given a tolerance ε > 0 and knowing which elements of e A are bound to be smaller than ε in magnitude, we are free to remove them altogether from the product (8.4). The outcome departs entry-by-entry from the exact exponential by at most ε and, by design, it resides in the Lie group G.

The sparsity pattern of a matrix exponential of a sparse matrix after the excision of small entries is at present unknown. (Banded matrices are an exception.) Moreover, full implications of this phenomenon to the subject matter of this sec- Having said this, our analysis emphasises a very important and welcome feature of methods based upon canonical coordinates of the second kind, applied in conjunction with Chevalley bases: the connection between basis elements and entries of the matrix is a powerful tool, which we can bring to bear with pointwise precision upon the exploitation of sparsity.

Stability and backward error analysis

How stable are Lie-group methods? A combined wisdom of half-a-century of computational analysis of ODEs is that numerical algorithms for initial-value problems are to all intents and purposes useless unless they exhibit favourable stability properties. Indeed, much of the narrative of modern numerical analysis of ODEs is the tale of stability, linear and nonlinear alike, culminating in a profound understanding of the subject. We refer the reader to the monograph of [START_REF] Stuart | Dynamical Systems and Numerical Analysis[END_REF] for a comprehensive review of this important subject area. The word 'stability' has been so far conspicuously absent from our exposition. A partial reason is ignorance: much remains to be done in the realm of stability investigations in a Lie-group setting. Interesting results abound which cannot yet be fitted into general theory. Thus, we can learn from computation that merely projecting a solution into the right manifold often leads to instabilities, while intrinsic Lie-group methods exhibit much more favourable behaviour. Much of the advance in stability theory for classical numerical ODEs was concerned with identifying appropriate stability models: broad enough to provide insight about many differential systems of interest, yet sufficiently focused and narrowly defined to be amenable to rigorous analysis. Thus, the linear model, the monotone model and the many more advanced models from [START_REF] Stuart | Dynamical Systems and Numerical Analysis[END_REF], originating in the theory of nonlinear dynamical systems. This chapter in the narrative of Lie-group methods cannot yet be written, except for the observation that some Lie-group solvers, e.g. RK-MK, Magnus and Fer expansions with exactly-calculated exponentials compute the solution of linear equations with constant coefficients exactly: in that case there is no need for stability analysis!

The last few years have seen the emergence of an alternative stability theory, mainly within the context of symplectic integration of Hamiltonian ODEs and discretization of dynamical systems. In addition to asking "how near is the numerical solution to the exact one and how influenced is it by small perturbations?", the new breed of stability researchers also poses a different query: "what is it that our numerical method solves exactly? And how far apart is it from the equation that we wish to solve?". This is precisely the question of backward error analysis that James Hardy Wilkinson made into the centrepiece of modern numerical linear algebra. Arguably, it is just as relevant in the ODE setting and it has already led to impressive new insights [START_REF] Benettin | On the Hamiltonian interpolation of near-tothe-identity symplectic mappings with application to symplectic integration algorithm[END_REF][START_REF] Hairer | Backward analysis of numerical integrators and symplectic methods[END_REF][START_REF] Hairer | The life-span of backward error analysis for numerical integrator[END_REF][START_REF] Neishtadt | The separation of motions in systems with rapidly rotating phase[END_REF][START_REF] Reich | Backward error analysis for numerical integrators[END_REF].

In this section we report briefly on recent work of [START_REF] Faltinsen | Backward error analysis for Lie-group methods[END_REF], who has generalised backward error analysis to a Lie-group setting.

Discussion of stability is meaningless without the concept of distance. Thus, given a Lie group G, we seek d : G ×G → R + which is consistent with the standard axioms of a metric in an Euclidean space and, in addition, compatible with the topology of G. A good way to ensure compatibility is to require that d is left invariant, i. An example is the geodesic metric in O(N ) (which, as a matter of fact, is biinvariant): d(X, Y ) = ∥η∥ 2 , where the eigenvalues of X T Y are σ 1 , σ 2 , . . . , σ N and η k = -iσ k , k = 1, 2, . . . , N. (Note that the spectrum of elements in O(N ) lives on the complex unit circle, hence the η k s are real.) Other examples are more complicated to derive explicitly, but this is of little consequence since we do not require d in a closed form.

We are concerned with the solution of the Lie-group equation

Y ′ = A(t, Y )Y, t ≥ t 0 , Y(t 0 ) = Y 0 ∈ G, (9.2)
where A : [t 0 , ∞) × G → g, by a Lie-group method. The flow corresponding to (9.2) is denoted by Ψ t,t0,A , therefore Y (t) = Ψ t,t0,A (Y 0 ).

Assumption 9.2 The matrix function A is real analytic and there exist constants α, β, t * > 0 such that

∥A(t, Ŷ ) Ŷ ∥ F ≤ α, t ∈ [0, t * ], Ŷ ∈ B β (Y 0 ), where ∥ • ∥ F is the Forbenius norm and B β (Y 0 ) = { Ŷ ∈ G : d(Y 0 , Ŷ ) ≤ β}.
We are interested in Lie-group solvers that lift the solution to the corresponding Lie algebra g, whether once or repeatedely in the course of every time step. All the methods that we have described in this survey: Crouch-Grossman and Runge-Kutta-Munthe-Kaas schemes, Magnus, Fer and Cayley expansions and methods based on canonical coordinates of the second kind fit this framework. An example of a method that is outside the scope of the theory of this section is projection. For example, in the case G = O(N ) we might time-step from Y n to Y n+1 , say with an arbitrary ODE method which produces a new value Ȳn+1 . We subject Ȳn+1 to a polar decomposition and retain the orthogonal part as our new Y n+1 [START_REF] Higham | Time-stepping and preserving orthonormality[END_REF]. We hasten to acknowledge that this is a perfectly valid procedure, except that it is outside the scope of our present discussion. The map induced by the numerical method will be denoted by Φ h,tn,A , hence

Y n+1 = Φ h,tn,A (Y n ), n ∈ Z + .
Assumption 9.3 The Lie-group method is accurate to order p ≥ 1, i.e.

d(Ψ h,t0,A (Y 0 ), Φ h,t0,A (Y 0 )) = O h p+1 .
What is it that Φ h,tn,A solves exactly? In linear algebra this is precisely the question of backward error analysis. Insofar as ODEs are concerned, however, the situation is slightly more complicated.

Theorem 9.1. [START_REF] Faltinsen | Backward error analysis for Lie-group methods[END_REF]) Subject to Assumptions 9.1-3, there exists (9.4) where γ > 0 is a constant.

a matrix function A h : [t 0 , ∞) × G → g such that ∥A h (t, X) -A(t, X)∥ F = O(h p ) ( 9 . 3 ) and d(Ψ h,t0,A h (Y 0 ), Φ h,t0,A (Y 0 )) = O e -γ/h ,
Note that Ψ h,t0,A h evolves on the very same Lie group G and, because of (9.3), the modified equation Y ′ h = A(t, Y h )Y h 'approximates' the ODE (9.2). The above theorem argues that a single step of the numerical method departs to an exponentially small extent from the exact solution of a nearby equation! It is reminiscent of similar results in symplectic integration and its method of proof generalises the work of [START_REF] Reich | Backward error analysis for numerical integrators[END_REF] to Lie-group setting.

A good approximation across a single step does not tell much. Ideally, we wish to extend the scope of Theorem 9.1 to [t 0 , ∞) or, at the very least, to a large number of steps. This, however, requires further conditions. Similarly to backward error analysis for symplectic integration, the verification of such conditions in a nonlinear case might be difficult and require bespoke analysis for different methods. Matters simplify a great deal, though, for linear equations Y ′ = A(t)Y . Suppose that the exact solution is Y (t) = expm[Θ(t)]Y 0 , t ≥ t 0 . Let µ ∈ R be the least real number such that ∥Ad expmΘ B∥ F ≤ c e µ(t-t0) ∥B∥ F , t≥ t 0 for some c > 0 which may depend on Θ. Then (9.4) can be extended to a longer interval. Specifically, it is true that

d(Ψ tm,t0,A h (Y 0 ), Φ mh,t0,A (Y 0 )) = O e -γ * /h ,
where γ * > 0 and m ≤ M (h), where It is possible to prove that g = so(N) implies that the operator Ad expmΘ is itself skew symmetric and µ = 0. This implies that the scope of backward error analysis is quite significant at least in this case.

M (h) = ⎧ ⎨ ⎩ O(1) , µ>0, O(h -p ) , µ= 0
Another interesting phenomenon originally identified in the analysis of Hamiltonian problems is linear error growth [START_REF] Estep | The Rate of Error Growth in Hamiltonian-Conserving Integrators[END_REF]. Provided that the exact solution of a Hamiltonian differential system is periodic with period T > 0 and a pth-order numerical method satisfies convenient requirements (e.g., reversibility or conservation of Hamiltonian energy), it is possible to prove that ∥y nT /hy(nT )∥ ≤ cnh p , (9.5)

where c > 0 (for simplicity we assume that T /h is integer). This is a very important feature of successful long-term integration, since lesser methods typically produce quadratic error growth and the solution is unlikely to remain periodic for long. Does (9.5) remain valid in a Lie-group setting? Unfortunately, not always. In a recent paper, though, [START_REF] Engø | Numerical integration of Lie-Poisson systems while preserving coadjoint orbits and energy[END_REF] proved that solving a Lie-Poisson system with a method that is both a Lie-group action and conserves Hamiltonian energy results in linear error growth. As an example, let us recall the Euler equations for a rigid body (7.8). This is a Lie-Poisson system which conserves the Hamiltonian energy H(y) = 1 2 (m 1 y 2 1 + m 2 y 2 2 + m 3 y 2 3 ), as well as evolving on the unit sphere in R 3 . Therefore y(t) lives on a circle (an intersection of a sphere and an ellipsoid), a feature shared by energy-conserving methods based on group actions, and is periodic there. We have already seen in Figures 7.2 and 7.3 that self-adjoint methods do well in recovering a periodic solution.

In Figure 9.1 we have integrated the rigid body equations (7.8) with three RK-MK methods: Forward Euler, Runge-Kutta-Heun and the trapezoidal rule, all applied in the Lie algebra. Only the latter method conserves energy [START_REF] Engø | Numerical integration of Lie-Poisson systems while preserving coadjoint orbits and energy[END_REF]). The results have been displayed as point-plots of y n in three dimensions. First note that all the solution trajectories evolve on the unit sphere: unsurprising, since we are using Lie-group methods, yet beyond the reach of most classical algorithms. Secondly, forward Euler is no respecter of periodicity and its trajectory spirals to a fixed point. Runge-Kutta-Heun is much better, yet more careful examination demonstrates how the error accumulates and periodicity is lost. The solution is qualitatively correct for a while, but long integration leads to false dynamics also in this case. The energy-conserving trapezoidal rule, though, produces a trajectory which to all intents and purposes is periodic.

10. Implementation, error control and DiffMan 10.1. Implementation and error control of Lie-group solvers Practical implementation of Lie-group solvers requires much more than merely to program a numerical method. We must address ourselves to issues like error control and variable-step implementation. As is perhaps natural in a new subject, implementation details have received less attention so far than theoretical issues, a situation that is likely to be remedied in the next few years.

In this section we survey the little that is presently known about implementational issues, commencing with the welcome observation that variable-step procedure does not interfere with the retention of Lie-group structure. This is important, since it is known that another important geometric-integration technique, symplectic solution of Hamiltonian systems, loses many of its most favourable features unless implemented with (essentially) constant step size [START_REF] Serna | Numerical Hamiltonian Problems, Chapman & Hall. C. Størmer[END_REF].

There are two levels of discretization in Lie-group solvers and each should be monitored in variable-step implementation and contribute to the estimate of local error:

1 The error committed in the evaluation of the coordinate map; and 2 The error incurred in the solution of the Lie-algebraic equation.

Insofar as the coordinate map is concerned, the situation is simple. The Cayley map (6.3) requires an inversion of a matrix. Unless its size is large, this can be accomplished by direct methods, otherwise it requires iteration. In the first case the only source of error is roundoff and this issue is well understood by classical numerical algebra. In the second case the error is determined by the termination criteria and the issue is, again, transparent. The use of techniques based on coordinates of second kind (6.15) generates roundoff error and nothing else, since each individual exponential can be evaluated easily in exact arithmetic. Unfortunately, the situation is different with the most important coordinate map, expm. To date, there are no efficient means to monitor the error of methods from Section 8. The last statement refers not just to approximation methods but also to 'exact' calculation of the exponential: it will be seen in Section 11.5 that the Matlab function expm is a stumbling block to high-precision long-term integration of some systems. 5It is entirely conceivable that good techniques to monitor the error in expm should depend on the Lie group in question: intuitively, a compact object like O(N ) is easier to handle than SL(N ), say.

Estimation of the error in Lie algebra presents different challenge, more akin to classical error-control theory for numerical ODEs.

A fair share of quality software for 'classical' ODEs is based on multistep methods. There exist multistep Lie-group solvers [START_REF] Faltinsen | Multistep methods integrating ordinary differential equations on manifolds[END_REF], but their implementation requires travelling forwards and backwards between the group and the algebra in a manner which is, arguably, too cumbersome for the use of 'automatic differentiation' techniques of [START_REF] Gear | Numerical Initial Value Problems in Ordinary Differential Equations[END_REF], the cornerstone of most multistep ODE software.

In principle, error control of explicit RK-MK can be accomplished in a straightforward manner, using embedded Runge-Kutta schemes [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]. The sole difference with the classical framework is that, instead of estimating the error in the original configuration space, we do so in the algebra. An embedded RK scheme has the Butcher tableau

c 1 a 1,1 a 1,2 • • • a 1,ν c 2 a 2,1 a 2,2 • • • a 2,ν . . . . . . . . . . . . c ν a ν,1 a ν,2 • • • a ν,ν b 1 b 2 • • • b ν b1 b2 • • • bν
and the approximants

y n+1 = y n + h ν l=1 b l f l ȳn+1 = y n + h ν l=1
bl f l (compare with (3.3)!) are of order p and p ≥ p + 1 respectively. The higher-order approximant is used for error control, y n+1y(t n+1 ) ≈ y n+1 -ȳn+1 . Applying similar trick in the Lie algebra readily yields an estimate of the local error at relatively little extra cost for low-order methods. If the order is high, typically the embedded RK scheme requires increasingly large number of additional stages to evaluate ȳn+1 and the cost mounts. This is a problem common to high-order RK methods, also in classical setting.

We have neglected in our discussion of RK-MK one important source of error:

in practical applications the dexp -1 operator (2.46) is truncated consistently with the order of the method, typically replaced with

dexp -1 A (C, p) = p-1 j=0 B j j! ad j A C.
This carries error which we are forbidden to neglect, yet can estimate easily from the leading term of dexp

-1 A C -dexp -1 A (C, p), |B q | q! ∥ad q A C∥, where q = 1, p = 1, 2⌊(p + 1)/2⌋, p≥ 2.
Much more challenging is error control of Magnus expansions, the subject of [START_REF] Iserles | On the implementation of the method of Magnus series for linear differential equations[END_REF]. Again, there are two discretization steps: truncation by power (4.14) of the Magnus expansion and the replacement of integrals by quadrature. In the sequel we restrict the discussion to the linear case Y ′ = A(t)Y .

To estimate the error in the leading truncation term, we commence by assuming that the leading terms in the expansion of the matrix A are known,

A(t) = C 0 + C 1 t + C 2 t 2 + • • • .
The error term in t 0 C τ (ξ)dξ can be evaluated from the tree τ in the following manner. We label each leaf of the tree by few leading terms of the expansion and prune the tree according to the original composition rules from Section 4.1:

1 If two leaves share a parent, they are excised and the parent is labelled by the commutator of their labels.

2 If a leaf is the only child of a parent then it is eliminated and its parent is labelled with the integral of its label.

By the end of this procedure we throw away all the terms except for the leading one. An example will clarify this procedure,

r r ◗ ◗ ✑ ✑ r r ◗ ◗ ✑ ✑ r r ◗ ◗ ✑ ✑ r C0 C0+C1t C0+C1t C0 ⇒ r r ◗ ◗ ✑ ✑ r ◗ ◗ ✑ ✑ r r ❜ ❜ ✧ ✧ C0t C 0 t+ 1 2 C1t 2 C0+C1t C0 ⇒ r r ◗ ◗ ✑ ✑ r ◗ ◗ ✑ ✑ r C0t 1 2 [C0,C1]t 2 C0 ⇒ r r ◗ ◗ ✑ ✑ r ◗ ◗ ✑ ✑ C0t 1 6 [C0,C1]t 3 C0 ⇒ r r ❍ ❍ ❍ ✟ ✟ ✟ C0t - 1 6 [C 0,[C0,C1]]t 3 ⇒ r - 1 6 [C 0,[C0,[C0,C1]]]t 4 ⇒ -1 30 [C0,[C0,[C0,C1]]]t 5 .
This procedure has been automated in (Iserles et al. 1999), using a Maple program.

Often it is easy to find derivatives of A explicitly but in general-purpose software one needs to approximate them by finite differences. The scaled function values A 1 , A 2 , . . . , A ν are a perfectly good starting point for this elementary calculation. Note that very few derivatives are required, e.g. just C 0 and C 1 for the power-5 tree in the above example, and the approximants need not be very precise. If we are using the adjoint basis, we might just as well use (appropriately scaled) h -l B l in place of C l-1 : the difference between expansions at 0 and 1 2 h is subsumed into higher-order terms.

The computation of quadrature error is more challenging and the results of (Iserles et al. 1999) are of a more tentative nature. This is hardly surprising, since even the estimation of the error of univariate Gauss-Legendre quadrature is difficult [START_REF] Davis | Methods of Numerical Integration[END_REF]. Indeed, perhaps paradoxically, the univariate quadrature is the most problematic also in our setting. Of course, we may evaluate each integral by two quadratures, e.g. Gauss-Legendre and Gauss-Lobatto, or perhaps consecutive Gauss-Legendre rules, but this doubles the cost. The remedy, proposed in (Iserles et al. 1999), is to use again the derivatives of A: recall that in many instances they are easy to evaluate explicitly.

Straightforward, yet messy, calculation, expanding everything in sight into Taylor series in h, demonstrates that the leading error terms in the first four integrals in (4.10), using two-point Gauss-Legendre quadrature, are

E 1 = -1 180 C 4 h 5 E 2 = -( 1 240 [C 0 , C 3 ] + 1 1080 [C 1 , C 2 ])h 5 , E 3 = 1 1080 ([C 0 , [C 0 , C 2 ]] + [C 1 , [C 0 , C 1 ]])h 5 , E 4 = -( 1 270 [C 0 , [C 0 , C 2 ]] + 1 720 [C 1 , [C 0 , C 1 ]])h 5 .
respectively. Note that only E 1 depends on the fourth derivative. A Matlab program incorporating variable time-stepping strategy with above error estimators for the fourth-order Magnus method has been applied in [START_REF] Iserles | How large is the exponential of a bounded matrix?[END_REF]) to a large number of test problems, in each case behaving predictably and producing error consistent with the specified tolerance. Yet, much work remains to be done in this subject area, not just to estimate quadrature error without the need for higher derivatives, but also to develop the right strategies to the methods of [START_REF] Blanes | Improved high order integrators based on Magnus expansion[END_REF] from Section 5.4.

The DiffMan package

DiffMan (Engø et al. 1999) is an object-oriented Matlab toolbox for solving differential equations on manifolds. The package embodies most of the algorithms discussed in this survey. This software is in the public domain and it can be downloaded from the DiffMan home page at http://www.ii.uib.no/diffman. Some of the numerical examples of Section 11 are distributed with DiffMan. Since the package is still in the stage of intensive development, we do not wish to elaborate excessively upon the details of DiffMan and of its usage. Instead, we refer the reader to the latest version of the DiffMan user manual, which can be found at the above-mentioned homepage.

The basic philosophy behind DiffMan is the idea of Coordinate-free numerics (Munthe-Kaas and Haveraaen 1996), a research programme devoted to the study of the role of abstract formulations, independent of particular representations, in computational and applied mathematics. We have touched briefly upon this topic in Section 3, emphasising the importance of using abstractions as a tool for organising object-oriented software.

The DiffMan package is built upon classes which are modelling continuous mathematical structures. There are currently three main types of classes: domains, fields and flows. Domains consist of Lie algebras, Lie groups and homogeneous manifolds. Fields are structures built over manifolds. Currently the only fields are vector fields, but more general tensor fields (and possibly the even more general fibre bundles) are likely to be included in the future. Numerical algorithms are collected in the class of flows on manifolds.

The DiffMan package displays the importance of integrating symbolic and numerical techniques in the same software. For example, a particular Lie algebra implemented is the free Lie algebra class. It is capable of symbolic computations which remain valid in any particular finite-dimensional Lie algebra. Hence, numerical algorithms may be developed using this package, and the resulting expressions may be evaluated later, substituting data from a concrete Lie algebra. Abstract concepts play a central role, characteristic not just of pure mathematics but also of modern non-numerical software.

Another important insight that has emerged from this research is the observation that some numerical algorithms require detailed knowledge of representations of objects, while others can be completely formulated using general, coordinateindependent operations. This distinction is most easily seen in the area of solving linear algebraic equations. For example, Gaussian elimination requires detailed knowledge of matrices in terms of components (there is a major difference between sparse and dense Gaussian elimination!), while other algorithms such as conjugate gradients just require access to matrices via their properties as linear operators. Algorithms which can be formulated in terms of general, coordinate-independent operations are much more flexible in their use than algorithms tied to particular representations. The fact that Runge-Kutta methods can indeed be phrased as coordinate-free algorithms has many implications that we are only now beginning to understand.

Applications

The heavy top

Many systems evolving in the physical space can be modelled by motions consisting of rotations, like the rigid body that we have already encountered in (7.8), or translations, or combination thereof, in which case we say that the motion is described by the special Euclidean group SE(3).

The heavy top equations are an important instance in which we make use of SE(3) action. It differs from the rigid body (7.8) because of the presence of gravity.

The equations of motion can be described in either space or body coordinates. In this section we briefly discuss both, considering space coordinates first. We denote by Π s the angular momentum of the heavy top, by Ω s its angular velocity, by I s the inertia tensor and by λ s the unit vector pointing toward the centre of mass of the heavy top. The equations of motion are

Π ′ s = Mglk × λ s , I ′ s = [I s , Ω s ], λ ′ s = Ω s × λ s ,
where M is the mass of the top, g = 9.81 is the gravitational constant, l is the distance from the origin of the space-coordinate system to the centre of mass of the heavy top and k is the unit vector along the z-axis. The angular velocity Ω s is related to the angular momentum Π s by the identity

Ω s = I -1 s Π s .
We refer the reader to [START_REF] Arnold | [END_REF], Goldstein 1980) for background, theory and notation.

It is well known that the axis of the top displays a very special motion composed of a precession about the vertical axis z and nutations, nodding up and down between two bounding angles θ 1 and θ 2 as displayed in Figure 11.1 for a top with initial angular velocity Ω b = [0, 0, 100] T about its figure axis, l = √ 3/2, and mass M = 20/(gl). The inertia tensor in body coordinates is I b = diag(1, 1, 1/5) and the initial position of the axis of the top is rotated at an angle θ = -π/10 with respect to the z-axis. The same motions can be described in body coordinates, in terms of two vectors Π b and Γ b in R 3 , the angular momentum and the gravity vector in body coordinates respectively [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF]. The equations simplify to

Π ′ b = Π b × Ω b + MglΓ b × χ, Γ ′ b = Γ b × Ω b , (11.1)
where χ is the unit vector on the figure axis of the top. Note that in body coordinates the inertia tensor is constant. Also in this case, Ω b = I -1 b Π b . Although the vector Γ b is parallel to the z axis of the space-coordinate system, its motion does not correspond to the motion of the figure axis λ s of the heavy top. To reconstruct the motion of the axis we integrate the configuration matrix R ∈ SO(3), so that Ω s = RΩ b and in general the same type of transformation of body coordinates to space coordinates holds for all relevant vectors. The configuration matrix R obeys the differential equation

R ′ = -R Ω b . (11.2)
Once the configuration matrix R is known, the position of the axis of the top can be recovered by means of the transformation

λ s = R ⎡ ⎣ 0 0 1 ⎤ ⎦ .
In our first experiment, the equations (11.1)-(11.2) have been solved with the Matlab routine ode45 with variable step size and error control. Since in body coordinates the gravity vector Γ b rotates very fast, the Matlab routine employs a very small step size in the integration interval [0, 10], ranging from h min ≈ 0.5 × 10 -6 to h max ≈ 0.00214. The average value is h mean = 0.0018 and the standard deviation σ = 1.972 × 10 -4 . The motion of the axis is quite similar to the one observed in Figure 11.1. However, when the integration is performed over longer time, the numerical approximation of λ s is not longer on the unit sphere and the amplitude of precessions shrinks (see Figure 11.2). The routine ode45 does not preserve much of the geometry of the underlying problem, and many of the conserved quantities of the problem (Casimirs and first integrals, see [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF]) are not conserved (see Figure 11.3).

To illustrate the advantages of the methods discussed in this article, we solve numerically the equations of the heavy top (11.1)-(11.2) using the coadjoint action of SE(3) on the dual of the algebra se(3) * (Marsden and Ratiu 1994),

Λ (R, Ω b ), (Π b , Γ b ) = (R Π b + Ω b × R Γ b , R Γ b ).
Using numerically the coadjoint action implies that all the Casimirs, in this case the projection of the angular momentum on the gravity axis, Π b • Γ b , and the norm of the gravity vector ∥Γ b ∥, are automatically preserved to machine accuracy [START_REF] Engø | On the construction of geometric integrators in the RKMK class[END_REF]Faltinsen 1999, Zanna et al. 1999). In Figure 11 fourth-order RK-MK method based on the Butcher tableau 0

1 2 1 2 1 2 0 1 2 1 0 0 1 1 6 1 3 1 3 1 6
115 (cf. Appendix A), with step size h = 1 100 . We refer to this numerical scheme as RKMK4. Comparing Figure 11.4 with Figure 11.2, we observe that the tip of the axis remains on the unit sphere (with five times as large step size as the average one used by ode45), as we expected from theory. There is only a very slight variation in the angles θ 1 and θ 2 , bounding the nutations. This comes as no surprise: Lie-group type methods perform very well for systems with oscillatory behaviour, as will be further discussed in Section 11.5. Note that this scheme is neither time-reversible nor energy preserving. Using self-adjoint schemes, as those discussed in Section 7, or energy preserving schemes, as those of [START_REF] Engø | Numerical integration of Lie-Poisson systems while preserving coadjoint orbits and energy[END_REF], it is possible to further improve other geometrical features of the numerical solution.

We remark that for this type of actions, solved by Lie-group schemes based on canonical coordinates of the first kind, the integrator requires numerous evaluations of exponentials and dexp -1 of 3 × 3 skew-symmetric matrices, which can be computed exactly and very fast using the formulae described in Appendix B.

Sturm-Liouville problems

Numerous problems in applied mathematics require the solution of Sturm-Liouville problems, finding λ so that Ly = -y ′′ + q(t)y = λy, t ∈ (0, α), (11.3) with the sufficiently smooth potential q ∈ C m (0, α) and the boundary conditions (11.4) where rank [b 0 , b α ] = 2. Examples range from fluid flow to Schrödinger spectra to geophysics to NMR imaging and beyond. It is well known from classical functional analysis that the spectrum σ(L), i.e. the set of all λ that solve (11.3), is real, countable and accumulates at +∞. Its computation has attracted significant effort since the very dawn of numerical analysis and led to an impressive array of methods and software. We refer the reader to [START_REF] Pryce | Numerical Solution of Sturm-Liouville Problems[END_REF] and to references therein.

y(0)b 0 + y(α)b α = 0,
The problem (11.3) can be always formulated in SL(2), (11.5) and this makes it a 'natural' for Lie-group methods. The main idea is to approximate the fundamental solution of (11.5) by a product of exponentials and impose the boundary conditions (11.4). This technique has been introduced and thoroughly analysed by [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF] and it constitutes a far-fetched generalisation of Preuss-Fulton methods [START_REF] Preuss | Solving linear boundary value problems by approximating the coefficients[END_REF]). The outcome is a scalar nonlinear equation for the spectral parameter λ. For example, suppose that (11.4) simplifies to y(0) = y(α) = 0 and our approximation to the fundamental solution of (11. Suppose for example (and a very trivial example it is!) that Y is the second-order truncation of the Magnus expansion,

y ′ = 0 1 q(t) -λ 0 y, t ∈ (0, α), where y = y y ′ ,
Y (t) = expm 0 t t 0 q(ξ)dξ 0 . Therefore (11.6) becomes α sin α 2 λ -α α 0 q(ξ)dξ α 2 λ -α α 0 q(ξ)dξ = 0
and the approximate eigenvalues are

λ m ≈ m 2 π 2 + α α 0 q(ξ)dξ α 2 , m∈ N.
Adding the next term to the Magnus expansion yields a somehow better approximation,

λ m ≈ m 2 π 2 + α α 0 q(ξ)dξ -1 4 α 0 (2ξ -α)q(ξ)dξ 2 α 2 , m∈ N.
We can continue in this vein, possibly replacing integrals by quadrature, except that the outcome of this naive approach is very poor. Magnus expansions are an excellent means to approximate the fundamental solution locally and we can hardly expect a single exponential to approximate the solution well in the entire interval (0, α). The situation is further exacerbated when α is large or when (as is the case in many instances of practical interest) we desire to approximate a large number of eigenvalues. It is possible to show that ω is a polynomial in λ, of a degree that grows with the order of the Magnus expansion. Bearing in mind Theorem 4.1, we can thus hardly expect convergence for large λ. A superior alternative is to partition the interval into small subintervals, where convergence and adequate precision can be assured. The details of this procedure are reported in [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF].

Let us restrict the discussion for the sake of simplicity to a fourth-order method. Letting q n,1 = q((n

+ 1 2 - √ 3 6 )h), q n,2 = q((n + 1 2 + √ 3 6 )h), where h = α/n * , we set Θ n (λ) = - √ 3 12 h 2 (q n,1 -q n,2 ) h 1 2 h(q n,1 + q n,2 ) -hλ √ 3 12 h 2 (q n,1 -q n,2 ) , n= 0, 1, . . . , n * -1.
The fourth-order approximant to the solution of (11.5) at α is Y λ y(0), where λ) .

Y λ = e Θ n * -1 (λ) • • • e Θ1(λ) e Θ0 (
To force the boundary conditions y(0) = y(α) = 0, say, we seek λ so that Y λ 1,2 = 0. This nonlinear equation can be solved by Newton-Raphson iteration which, eigenvalues being simple, converges quadratically near the solution. [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF] recommends a procedure based on the Newton-Raphson method being applied to an augmented equation, thereby simplifying the computation of the derivative.

The outcome is a very efficient method for the computation of Sturm-Liouville problems and it can be easily generalised to more elaborate eigenvalue problems and boundary conditions. However, if the interval in question is long or a large number of eigenvalues is desired, efficiency suffers. This is in particular the case if high-order Magnus methods are used, since the degree of elements of Θ n as polynomials in λ grows. As a matter of fact, the error in a pth-order method grows as O h p+1 λ p/2-1 . To overcome this phenomenon, which occurs also in other methods for the computation of Sturm-Liouville problems, [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF] introduced an interesting device which might be relevant to other applications of Lie-group methods. We observe first that, given matrix values in a self-adjoint basis B 0 , B 1 , . . . , B ν-1 (cf. Section 5), it follows at once from (5.7) that only B 0 depends upon λ.

Rather than presenting the construction from [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF], we introduce a conceptually similar idea of Iserles (n.d.). Our point of departure is a decomposition of the Magnus expansion (4.5) into streamers: partial sums of the form

H τ [0] (t) = ∞ k=0 α(τ [k] ) t 0 C τ [k] (ξ)dξ, where τ [k+1] = r ❅ ❅ r r τ [k] , k∈ Z + . (11.7) Let τ = τ [0] be in the form (4.8). It is possible to show that α(τ [k] ) = B s+k (s + k)! α(τ ), where α(τ ) = s i=1 α(τ i )
(cf. (4.9)). Therefore we can write the streamer as

H τ [0] (t) = α(τ [0] ) t 0 ∞ k=0 B s+k (s + k)! ad k R ξ 0 A(η)dη C τ [0] (ξ)dξ.
(11.8)

We say that a tree is primitive if it cannot be written in the form (11.7). The set of primitive trees in F m is denoted by F p m , whereby we might replace the truncated Magnus expansion (4.14) with the Magnus streamer expansion

Ξ p (t) = p-1 m=0 τ ∈F p H τ (t).
(11.9) Like (4.14), this is a pth-order approximant, except that each tree therein is accompanied by an infinitely-long streamer and there is one less level of truncation in the Magnus expansion (4.5). Using (11.9) makes sense only if there exists a good method to evaluate the streamer (11.8) without truncating the expansion. This is the case if the Lie algebra is of sufficiently small dimension.

(H τ (t)) = α(τ ) t 0 ∞ k=s B k k! C k-s D(ξ) π(C τ (ξ))dξ = α(τ ) t 0 C -s D(ξ) C D(ξ) e C D(ξ) -I - s-1 k=0 B k k! C k D(ξ) π(C τ (ξ))dξ.
where D(t) = t 0 A(η)dη. Provided that d is small, typically we can obtain the commutator matrix and its exponential explicitly. This, together with quadrature, leads to an explicit formula for the calculation of the streamer that does not require any truncation of the series.

Similar idea, reported in [START_REF] Moan | Efficient approximation of Sturm-Liouville problems using Liegroup methods[END_REF], ameliorates the deterioration in accuracy for large λ. For example, the error of sixth-order Magnus streamer is O h 7 λ , hence it behaves like a fourth-order method when h 2 λ ≈ 1, when the standard sixth-order Magnus, which carries an error of O h 7 λ 2 , reduces to an order-two method.

Charged particles in a magnetic field

The motion of charged particles moving in a magnetic field was first studied numerically by Carl Størmer (1907) as a part of his work on explaining the origin of northern light (aurora borealis).

A particle (of unit mass and unit charge) is moving in a magnetic field b according to the equations

y ′ (t) = v v ′ (t) = b(t, y) × v
where v is the velocity, y is the position and b is the magnetic field. We attempt to obtain a simpler system by assuming that b = const, neglecting the dependence of b on space and time. Written in a matrix form, this yields summarises the results. We see that the Lie-group method is much more accurate than the classical integrator for the same step size. Figure 11.6 displays the time step selection as a function of time for the two methods for the simulations given in the first line of Table 11.1. The dip in the middle corresponds to the part of the trajectory where the particle is bouncing back. Both methods must reduce the step size in this region, but the classical one does so relatively more than the Lie group method. Note that, whereas the Lie group integrator preserves ∥v∥ exactly, the ode45 integrator does so only up to the the order of the method.

We conclude that for this problem, where both commutators and exponentials can be computed quickly, the total cost of a Lie-group integrator is significantly less than that of a classical integrator. The time interval of integration is relatively small in this example: the importance of preserving geometric properties becomes increasingly more crucial for very long integration intervals.

Toda-lattice equations

Imagine a regular lattice of N particles, all of unit mass, and assume that each particle interacts with just its nearest neighbours, subject to an exponential interaction potential. The outcome is the Toda lattice, governed by the Hamiltonian equations of motion dq k dt = p k , dp k dt = e q k-1 -q ke q k -q k+1 , ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ k = 1, 2, . . . , N, (11.12)

123
where q k and p k are the generalised coordinates and momenta, respectively [START_REF] Toda | Theory of Nonlinear Lattices[END_REF], where we have let q 0 = -∞, q N +1 = ∞ in the linear case and q N +k = q k when the particles are arranged in a ring (cf. Figure 11.7). Equation (11.12) corresponds to the Hamiltonian potential

H(p, q) = 1 2 N k=1 p 2 k + N k=1
(e q k -q k+1 -1).

It is a special case of a Fermi-Pasta-Ulam (FPU) flow. Letting α k = 1 2 e (q k -q k+1 )/2 , β k = 1 2 p k , [START_REF] Flaschka | The Toda lattice i[END_REF] showed that (11.12) can be recast in the Lax form (11.13) where

Y ′ = [B(Y ), Y ], t ≥ 0, Y(0) = Y 0 ,
Y = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ β 1 α 1 0 • • • α N α 1 β
α N • • • 0 α N -1 β N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
, with α N = 0 for linear lattice configuration, and B(Y ) = Y + -Y -∈ so(N), where M + and M -denote the upper-triangular and the lower-triangular portions of the matrix M , respectively. We immediately identify (11.13) with the isospectral flow (1.1).

A special case of (11.12) that has elicited much interest is three particles in a ring. In that case, the linear transformation p = T p, q = T q where T = ⎡ ⎢ ⎣ in the new variables p, q ∈ R 2 (Zanna 1998). 7 The solution now evolves on a compact submanifold of R 4 , more specifically on a 2-torus. Therefore, a Poincaré section consists of two closed curves and a good numerical scheme should retain this important property.

The following calculation has been performed in [START_REF] Zanna | On the Numerical Solution of Isospectral Flows[END_REF]) for a variety of methods which follow a set pattern: First we transform the three-particle lattice with the initial conditions p = [1, 1, 0] T , q = 0, to the Lax form (11.13). We then solve the isospectral form with the constant step size h = 1 10 for 10 4 steps, transform the variables to (p, q) and sketch the Poincaré section. The following methods have been used:

1 The three-stage, third-order explicit Runge-Kutta method with the Butcher tableau 0

1 2 1 2 1 -1 2 1 6 2 3 1 6
, (11.15) applied directly to (11.13). Note that this method is not isospectral. The Poincaré section is displayed in Figure 11.8a and we note at once that, instead of periodic trajectories, the motion collapses rapidly to a spurious fixed point. 2 The Runge-Kutta method (11.15) applied to the orthogonal flow (11.16) translating back to (11.13) with Y n+1 = Q n+1 Y n Q T n+1 . Note that the numerical approximation Q n+1 ≈ Q((n + 1)h), produced by the method (11.15), does not evolve on SO(3). Figure 11.8b demonstrates that the solution again spirals to a fixed point, although less rapidly than in the former case. 3 The former method can be 'orthogonalised' by projection, subjecting Q n+1 to polar decomposition and discarding the nonorthogonal part. This can be further enhanced by subjecting the intermediate stages to polar decompositions. Relevant Poincaré sections are displayed in Figure 11.9. Evidently, the behaviour improves, yet the solution goes on spiralling to a spurious fixed point. 4 We solve (11.16) with the RK-MK method, using the Runge-Kutta scheme (11.15). The Poincaré section for this computation, the first with a 'proper' Lie-group method, is displayed in Figure 11.10a. We have not managed to get rid of the unwelcome spurious convergence to a fixed point, yet the trajectory spirals significantly slower. 5 Finally, (11.16) is solved with the fourth-order collocated Magnus method from Section A.3, applied to the orthogonal flow. The outcome is displayed in Figures 11.8-10 display consistent gradual improvement. This is not evident at all from the retention of the Lie-group per se, since the qualitative feature under examination, invariant tori, is of a different flavour. Yet, it is clear that conservation of Lie-group structure is advantageous and that purpose-designed Lie-group solvers hold the edge over projection methods.

Q ′ = B(QY n Q T )Q, t ≥ nh, Q(nh) = I,
Note that invariant tori would have been retained by a symplectic method, applied directly to the Hamiltonian equations induced by the potential (11.14). This is hardly a surprise: we have, after all, subjected methods applied in a Lie-group formalism to a 'Hamiltonian test'. It is possible to devise a test according to 'isospectral' ground rules, e.g. by monitoring global error in Lax-equation coordinates. This has been done in [START_REF] Calvo | Conservative methods for the Toda lattice equations[END_REF], demonstrating a marked advantage of isospectral methods (originating in Lie-group solvers) over symplectic algorithms.

Highly-oscillatory equations

Highly-oscillatory systems of ODEs feature in many applications and their computation currently absorbs much of the overall effort devoted to numerical analysis of ODEs. It is well known that classical solvers performs poorly, and this has motivated many novel and ingeneous techniques [START_REF] Petzold | Numerical solution of highly oscillatory ordinary differential equations[END_REF].

Past experience indicates that Lie-group methods might be a suitable means to solve highly-oscillatory systems (Iserles andNørsett 1999, Iserles et al. 1999). We commence from an example that has already featured in (Iserles and Nørsett 1999), the solution of the Airy equation y ′′ + ty = 0, t ≥ 0, y(0) = 1, y ′ (0) = 1.

(11.17)

The exact solution of (11.17) can be represented as a linear combination of Airy functions y(t) = 1 2 Γ( 2 3 )[3 2/3 Ai(-t) + 3 1/6 Bi(-t)], t≥ 0.

It is easy to prove, e.g. by the WKB technique, that the trajectory is a bounded function that oscillates like sin t 3/2 : the frequency increases with time. (Cf. Fig. 11.11.) This indicates that long-time integration is difficult and this can be confirmed by endeavouring to solve (11.17) with any popular, general-purpose solver. We have solved the Airy equation with two Matlab routines, ode113 and ode15s, the first for nonstiff and the latter for stiff problems. Both routines employ variable-order methods in tandem with sophisticated error control. Yet, although we have attempted a wide range of possible error tolerances, including the least that the software would accept (relative tolerance of 2.2 × 10 -14 ), the pointwise error at t = 100 consistently exceeded 9.42 × 10 -6 and 5.43 × 10 -4 respectively: very poor performance, in particular in the case of the stiff solver ode15s.8

We have solved the Airy equation (11.17) with Magnus methods of orders 2, 4 and 6 respectively, taken from (A.9), with the constant step-size sequence h = 1 10 , 1 20 , 1 40 , 1 80 . No attempt has been made to monitor the error or to optimize the solution. In particular, no advantage has been taken of the fact that (11.17) has been rendered as an SL(2) equation and, rather than using (8.1), we have computed the exponential with the Matlab function expm. Fig. 11.12. The solution y(t) of the Airy equation (11.17) at t = 100.

of our calculations. The plot on the left compares the absolute error at t = 100 for different values of h on a doubly-logarithmic scale. For comparison purposes we have included the optimal results for ode113 and ode15s, assigning to them h equal to the average step size. The plot on the right compares the logarithm of the absolute error with the number of flops expanded on the solution: a very imprecise measure, yet a fair reflection of the computational effort.

The main observation is that, unlike the Matlab routines, Magnus methods perform very well indeed. The deterioration in the accuracy of the sixth-order method for small step size can be explained by the accumulation of small errors in the evaluation of the expm function. In principle, we could have replaced it by the exact formula (8.1), except that, while nominally reducing the number of flops, this does not enhance the solution because of inexact computation of hyperbolic functions. In all likelihood, the behaviour in Figure 11.12 is as good as one can expect with IEEE computer arithmetic.

Another example of a highly-oscillatory ODE is t 2 y ′′ty ′ + (1 + t 2 )y = 0, t ≥ 1, y(1) = J 0 (1), y ′ (1) = J 0 (1) -J 1 (1), (11.18) whose solution is y(t) = tJ 0 (t). Here J ν is a Bessel function of a first kind, whose oscillatory behaviour is well known. Note that we commence integration at t = 1, to avoid a regular singularity at the origin. The top plot of Figure 11.13 depicts the exact solution of (11.18). The oscillatory behaviour is evident: as a matter of fact, it is elementary to prove that the solution behaves like 2t/π sin t for t ≫ 1. The bottom plot displays in a logarithtmic scale the maximal error in the interval [0, 250] in second-order, fourth-order and sixth-order Magnus methods for a range of different step sizes. For comparison purposes, we have also included the same data for the Matlab routine ode113 for different tolerances, down to the least permitted by the software. As in the case of Fig. 11.12, we have assigned to each such computation an average step size. Note that the step sizes are significantly larger for Matlab routines, but so is the error! Even though it uses error control, ode113 appears incapable of driving pointwise accuracy below quite a large treshold. We note in passing that the constant-step-size Magnus methods were (for the largest step size, h = 1 10 ) 244%, 317% and 476% more expensive than ode113, but no attempt has been made to optimize the calculations. As in the case of the Airy equation (11.17), the sixth-order method 'hits the buffers' for small step sizes and it appears that there exists a limit on attainable accuracy, probably caused by accumulation of roundoff error. Yet, this limit is twice the number of significant digints attainable by the general-purpose ode113.

The above comparison might appear as unfair to ode113 which, after all, has not been constructed with highly-oscillatory equations in mind. Such criticism misses the point altogether. The whole purpose of geometric integration is to take advantage of structure, rather than employing general-purpose methods! Moreover, our Magnus expansions have been implemented with constant step sizes. There is little doubt (and much evidence in (Iserles et al. 1999)) that variable-step implementation of Magnus methods, e.g. with the technique from Section 10.1, would have improved the odds drastically in their favour. Of course, there are other bespoke methods for highly-oscillatory systems [START_REF] Petzold | Numerical solution of highly oscillatory ordinary differential equations[END_REF]. It is not our claim that Lie-group methods are superior, since this has never been investigated comprehensively. Our aim is more modest, namely to argue that such an investigation might be very interesting.

What is the explanation for the remarkable performance of Magnus methods (and other Lie-group methods that we have not implemented in this section) for a highlyoscillatory problem? Although this issue constitutes an interesting open research problem, our suspicion is that the conservation of SL(2) structure has absolutely nothing to do with it. Classical numerical methods invariably employ the ansatz that locally the solution of an ODE behaves like a polynomial in t. Unfortunately, polynomials are a very poor means to approximate rapidly-varying and highlyoscillating functions. We either require minute step size or a very high degree of a polynomial -and the latter might lead to ill conditioning. All the Lie-group methods of this paper, however, are based on an entirely different representation of the solution, as an exponential of a matrix with polynomial entries. Unlike polynomials, exponentials of matrices with polynomial entries can easily model exponential change, high oscillation and changes in amplitude and frequency. This might well be the mechanism explaining the superior performance of Lie-group methods.

• Order-four Gauss-Legendre method based on geodesic coordinates

F 1 = 1 4 A 1 + 1 4 - √ 3 6 A 2 + 5 144 - √ 3 24 [A 1 , A 2 ], A 1 = hA t n + c 1 h, expm(F 1 -C h,n ) expm(C h,n )Y n , F 2 = 1 4 + √ 3 6 A 1 + 1 4 A 2 -5 144 + √ 3 24 [A 1 , A 2 ], A 2 = hA t n + c 2 h, expm(F 2 -C h,n ) expm(C h,n )Y n , C h,n = 1 4 (A 1 + A 2 ) - √ 3 24 [A 1 , A 2 ], Θ = 2C h,n , Y n+1 = expm(Θ)Y n .
• Order-four Gauss-Legendre methods with flow coordinates

F 1 = 1 4 A 1 + 1 4 - √ 3 6 A 2 + 1 288 - √ 3 96 [A 1 , A 2 ], A 1 = hA t n + c 1 h, expm(F 1 -C h,n ) expm(C h,n )Y n , F 2 = 1 4 + √ 3 6 A 1 + 1 4 A 2 -1 288 + √ 3 96 [A 1 , A 2 ], A 2 = hA t n + c 2 h, expm(F 2 -C h,n ) expm(C h,n )Y n , C h,n = 1 4 + √ 3 8 A 1 + 1 4 - √ 3 8 A 2 - √ 3 96 [A 1 , A 2 ], Θ = 1 2 (A 1 + A 2 ) - √ 3 48 [A 1 , A 2 ], Y n+1 = expm(Θ -C h,n ) expm(C h,n )Y n .
Although Magnus-type methods with geodesic/flow coordinates exist for every order p, devising such schemes for nonlinear problems and order greater than four is hard. For this reason we shall restrict our attention to order-four methods. We let c 1 , c 2 , . . . , c ν be the collocation nodes.

• Collocation order-four Magnus method with geodesic coordinates:

F k = ν j=1 a i,j A l + 1 2 ν l,j=1
(a k;l,j + 1 2 b l a k,l )[A l , A j ],

A k = hA(t n + c k h, expm(F k -C h,n ) expm(C h,n )Y n ), ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ k = 1, . . . , ν,

B. Fast computation of 3D rotations

Rotations in three dimensions are ubiquitous in computational mechanics, hence it is important to have fast algorithms for their computation. The Lie algebra so(3) can be realized either as the set of all skew-symmetric 3×3 matrices with the matrix commutator as the bracket, or as the Euclidean space R 3 with the vector product as the bracket. As we will see, some formulae are most easily expressed by representing so(3) as 3-vectors, while other formulae appear more naturally in terms of skewsymmetric 3×3 matrices. It is convenient to switch back and forth between these forms. The Lie algebra isomorphism between these two representations is given by the hat map, (B.9)

x = ⎡ ⎣ 0 -x 3 x 2 x 3 0 -x 1 -x 2 x 1 0 ⎤ ⎦ (B.
The first of those, (B.6), is well known in literature as the Rodrigues formula for the exponential mapping [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF].

In many instances we need to compute repeatedly expressions in the general form f (ad x )(y) for some function f (x). These expressions in so(3) can be computed fast in a very similar manner. Since ad x (y) = xy, it is true that f (ad x )(y) = f ( x)y, and we may use the same technique as above to simplify f ( x). where α 1 =tan -1 x 3 θ -1 sin θ + x 1 x 2 θ -2 (1cos θ) cos θ + x 2 1 θ -2 (1cos θ)

, α 2 =sin -1 [x 2 θ -1 sin θ + x 1 x 3 θ -2 (1cos θ)],

α 3 =tan -1 x 1 θ -1 sin θ + x 2 x 3 θ -2 (1cos θ) cos θ + x 2 3 θ -2 (1cos θ)

.

(B.16)

Fig. 2

 2 Fig. 2.1. Examples of manifolds embedded in R 3 : a sphere, a doubly-twisted Möbius-strip-like torus and a twisted ribbon.

  )dξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . order 2 (4.16) ξ 2 ), A(ξ 1 )]dξ . . . . . . . . . . . . . . . . . . . . . . . . . . order 4 (4.17)

  dξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .order 2 (ξ 3 ), A(ξ 2 )], A(ξ 1 )]dξ . . . . . . . . . . . . . . . . . . . . . . . order 3 2 )A(ξ 1 )A(ξ 3 )dξ . . . . . . . . . . . . . . . . . . . . . . . . . . . order 4 A(ξ 4 ), A(ξ 3 )], A(ξ 2 )], A(ξ 1 )]dξ . . . . . . . . . . .

[

  [B 1 ]] 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . grade 1

  B 1 ]] 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . grade 5

  Fig. 7.1. Representation of the geodesic and flow-symmetric midpoint

  Fig. 7.2. Error in the Hamiltonian versus time for an order-four RK-MK method based on Gaussian nodes. The top plot corresponds to coordinates centred at y n while the second plot corresponds to geodesic and flow-symmetric coordinates.Although generally the two latter choices would correspond to different error, in this case the errors are very similar. The bottom plot corresponds to the difference between the errors in geodesic and flow-symmetric coordinates.

  Fig. 7.3. Error in the Hamiltonian versus time for an order-four Magnus method based on Gaussian nodes. The top plot corresponds to coordinates centred at y n while the second plot corresponds to geodesic and flow-symmetric coordinates. The bottom plot corresponds to the difference between the errors in geodesic and flow-symmetric coordinates.

  Fig. 8.1. How large is the exponential of a sparse matrix?

  would have liked d to be bi-invariant: satisfying both (9.1) and the condition d(XZ, Y Z) = d(X, Y ) for all X, Y, Z ∈ G. This, however, is not always possible. On the other hand, according to the Birkhoff-Kakutani theorem, every Lie group G admits a left-invariant, almost right-invariant metric which, in addition to (9.1), obeys d(XZ, Y Z) ≤ ρ(Z) d(X, Y ), where the function ρ is finite. Assumption 9.1 The metric d is left-invariant and almost right-invariant.

  Fig. 9.1. Rigid body equations (7.8), as solved by three different methods with h = 1 10 , integrating for 100000 steps.

Fig. 11

 11 Fig. 11.1. Precession and nutations for t ∈ [0, 10] of the axis of a heavy top.

  Fig. 11.2. Numerical precession and nutations projected on the (x, y)-plane for t ∈ [0, 100] of the axis of a heavy top. The equations are solved in body coordinates with the Matlab routine ode45.

Fig. 11

 11 Fig. 11.4. Numerical precession and nutations projected on the (x, y)-plane for t ∈ [0, 100] of the axis of a heavy top. The equations are solved in body coordinates (employing the coadjoint action) with the Lie-group scheme RKMK4.

  2 + bc.Note that a, b, c, ω are functions of both t and λ. We impose the boundary conditions and the outcome is the equation

  Let dim g = d and let C = {C 1 , C 2 , . . . , C d } be a basis of the algebra. Then there exists a natural map π : g → R d defined by π(B) = θ where B = d k=1 θ k C k . The action of the ad operator is a linear transformation, thus π(ad D B) = C D π(B), where C D ∈ gl(d) is a commutator matrix. 6 It follows from (11.8) and the definition of Bernoulli numbers that

  π

  Fig. 11.6. Time-step selection for RK-MK (upper) and classical RK (lower).

  Fig. 11.7. Different Toda-lattice configurations with N = 5.

  tandem with rescaling and elimination of one of the variables by employing the linear conservation law p 1 + p 2 + p 3 ≡ const, results in a simplified two-

  Fig. 11.8. Three-particle lattice solved by the Runge-Kutta method (11.15), applied a. directly to the isospectral flow (11.13) (on the left); and b. to the orthogonal flow (11.16) (on the right).

  Fig. 11.10. Three-particle lattice solved by Lie-group methods: a. The RK-MK method with the tableau (11.15); and b. fourth-order Magnus collocation from Section A.3.

Figure 11 .

 11 Figure 11.10b: the Poincaré section consists of the qualitatively-correct two closed curves, a section across a 2-torus!

Figure

  Fig. 11.11. The Airy equation (11.17).

  Fig. 11.13. The solution of the Bessel-like equation (11.18).

  b l,j [A l , A j ], Θ = 2C h,n , Y n+1 = expm(Θ)Y n .

2 .

 2 1) mapping x ∈ R 3 into a 3×3 skew-symmetric matrix x such that xy = x×y. (Note that, for clarity sake, we have abandoned our convention of denoting matrices with upper-case letters.) In particular, the identities[ x, y] = x yy x = x×y, (B.2) x y x = -(x T y) x, (B.3) x 2 y + y x 2 = -(x T x) y -(x T y) x, (B.4) x 2 y 2y 2 x 2 = -(x T y) x×y (B.5)are obeyed. For future convenience, we letθ = ∥x∥ = (x T x) 1/2 , ϑ= θNote that since x 3 = -θ 2 x, we deduce that for any real analytic function f (z) we can easily obtain real functions c 0 (z), c 1 (z) and c 2 (z) such that f ( x) = c 0 (θ)I + c 1 (θ) x + c 2 (θ) x 2 .

  x = c I -1 2 x + 1 4 xx T ,(B.13) with c as in (B.9). These formulae are based on the representation of so(3) as 3-vectors, so e.g. (B.12) should read dcay x (y) = c(y + 1 2 xy) = c(y + 1 2 x×y). If U = expm( x) is an orthogonal matrix, the matrix x can be obtained by means of the matrix logarithm as x = logm(basis, we have expm( x) = expm(α 1 C 1 ) expm(α 2 C 2 ) expm(α 3 C 3 ),

Table 2 .

 2 1. Correspondence between the matrix case and general Lie theory

	Matrix case	General case

  One of the most popular is a Hall basis H that contains the generators, S ⊆ H[START_REF] Bourbaki | Coefficients for the study of Runge-Kutta integration processes[END_REF]. All elements of H are produced by recursive commutation of generators. We can associate a length function to each element in the following fashion: ℓ(G i ) = 1 for G i ∈ S and ℓ(H) = ℓ(H 1 ) + ℓ(H 2 ) for all H ∈ H \ S, where H = [H 1 , H 2 ]. Intuitively, we may say that the length of H corresponds to how many commutators of generators are needed to construct H. In other words, the length function merely counts commutators.

Table 5 .

 5 1. Dimension of linear spaces of words of length equal to m in the Hall basis generated by G 1 , G 2 and G 3 .

	m 1 2 3	4	5	6	7	8	9	10
	ρ m 3 3 8 18 48 116 312 810 2184 5580
	p, the number of linearly-independent terms that we need to take into account is bounded by p m=1 ρ m or, because of Theorem 4.2, by p-1

Table 5 .

 5 2. Dimension of linear spaces of terms of weight exactly equal to m (top) and weight at most m (bottom) in the Hall basis generated by G 1 , G 2 and G 3 with weights 1, 2 and 3 respectively.

	m	1 2 3 4	5	6	7	8	9	10
	ρm	1 1 2 2	4	5 10 15 26	42
	m
	ρi 1 2 4 6 10 15 25 40 66 108
	i=1

Table 6

 6

	m 1 2	3	4	5	6	7	8	9	10
	ω i ≡ 1 σm 3 3 18 36 135 351 1134 3240 9963 29403 ω i = i σm 1 1 3 3 8 11 25 39 80 134

.1. Dimensions (6.13) of graded FHA for ν = 3 in two cases: ω i ≡ 1 and ω i = i.

Existence independently of representations should in fact be taken as the very definition of a geometric property.

)h, denote these function values by A 1 , A 2 and letB 1 = 1 2 (A 1 + A 2 ), B 2 = √ 3(A 2 -A 1). The relevant fourth-order

We return to this point in far greater detail in Section 8.

There are easier ways to determine k * but they require more Lie-algebra theory.

Sparsity in structure constants has no connection whatsoever with sparsity (or otherwise) of the matrix A.

In fairness to Matlab, in our experience expm performed better than alternatives.

Note that the usual definition of a commutator matrix is as an object in gl(N 2 ) for g ⊆ gl(N).

One reason why this case is interesting is that truncation of H to cubic terms results in the famous Hénon-Heiles Hamiltonian, which is known to be nonintegrable[START_REF] Berry | Regular and irregular motion[END_REF][START_REF] Toda | Theory of Nonlinear Lattices[END_REF]).

Highly-oscillating ODEs should be never confused with stiff ODEs. Our calculation merely confirms this.

Acknowledgements

The authors wish to thank many of their colleagues and fellow dexperts who have read portions of this survey as it was taking shape, in particular Krister Åhlander, Sergio Blanes, Elena Celledoni, Kenth Engø, Stig Faltinsen, Arne Marthinsen, Per Christian Moan, Ranjiva Munasinghe, Brynjulf Owren and Reinout Quispel. Their comments and perceptive remarks have contributed to diminish the number of mistakes, typos and infelicities and to improve the cohesion and the presentation of this paper.

HMK, SPN and AZ wish to acknowledge the support of Norwegian Research

Proof. With the same notation as above, we have

Because of (7.5), one has φ(C h,n ) -1 Y n = φ(C * -h,n+1 )Y n+1 , from which we deduce that φ(C * -h,n+1 ) -1 φ(Υ h,n+1 ) = I. Furthermore Υ h,n+1 and Υ * -h,n+2 are solution of the same differential equation (7.3) whereby the initial condition of Θ * is the endpoint of Θ and Θ * is obtained by means of Ψ * -h , where Ψ * h is the adjoint of Ψ h in the classical sense. Thus ( Ψ * -h • Ψh )Y n = Y n and the assertion follows.

✷

We have not specified yet what is C h,n . In general we let it be a function of the step size of integration h and of the current stage values F i of a Runge-Kutta method Ψ, C h,n = ϑ(F 1 , F 2 , . . . , F ν ),

ν being the number of the stages of the scheme (cf. Appendix A for notation). Thus, C * -h,n+1 = ϑ * (F * 1 , F * 2 , . . . , F * ν ), therefore the functions ϑ and ϑ * obey the fundamental adjointness condition

for the centre of the coordinate map.

The following result, which can be found in [START_REF] Zanna | Adjoint and selfadjoint Liegroup methods[END_REF], characterises self-adjoint Lie-group methods.

Theorem 7.2 Assume that Ψ is self adjoint. Then Ψ is self adjoint provided that

(7.6)

Geodesic and flow-symmetric coordinate maps

Basically, there are two distinct ways to generate coordinate maps such that (7.6) is obeyed. One way to achieve this goal is to choose ϑ so that the value

Y n is a 'midpoint in space' between Y n and Y n+ , which will generate what we call the geodesic midpoint method. An alternative is to choose ϑ so that Y n+1/2 = φ(C h,n ) -1 Y n is instead a 'midpoint in time', thus generating a flow midpoint. The situation is schematically represented in Figure 7.1.

In the first instance we say that the coordinate map is geodesic symmetric, while in the second case we say that the coordinate map is flow symmetric.

where b is a 3×3 skew symmetric matrix given by the hat map (B.1) in Appendix B. We compute the bracket of the vector fields in (11.10). Since they are linear, we find that 0

Hence, the vector fields do not form a Lie algebra. The remedy is to enlarge the family of 'simple' equations into (11.11) where B is a general 3 × 3 matrix. This results in the commutator

It is worth remarking on the similarity between this bracket and the brackets in Example 2.8, Example 2.10 and in Section 11.1. These are all special instances of so-called semidirect product brackets. Semidirect products are one of the most fundamental ways of constructing Lie algebras and Lie groups from simpler algebras and groups.

The solution of the simplified equation (11.11) is given in terms of initial conditions y 0 and v 0 as

a formula that can be computed exactly and efficiently using (B.6) and (B.10).

In a numerical experiment we have taken b(t, y) = b(y) as the magnetic dipole field b(y) = 3e T y me ym ∥y∥ 3 , where m = [0, 0, 1] T is the magnetic dipole vector. Figure 11.5 displays the orbit of a single particle with initial conditions y 0 = [0, -5 2 , 0] T , v 0 = 12 • 10 -3 [0, 0, 1] T . The time interval is [0, 10000]. The particle is moving back and forth in the van Allen belt around the Earth. This system is Hamiltonian, the Hamiltonian energy being H = ∥v∥ 2 . The preservation of H follows also trivially from the observation that the acceleration is always orthogonal to v. We wish to compare step sizes in classical integrators and RK-MK-type methods for this problem. The classical integrator used in the study is the ode45 code of Matlab, based on DOPRI5(4), an embedded (4,5) pair of Runge-Kutta methods due to Dormand and Prince. The code uses variable step size control. The RK-MK method we use is also based on the DOPRI5(4) method and we have employed the step-size controller of DiffMan. The interval of integration ranged from 0 to 500, corresponding to a motion of the particle from the equator up towards the north pole, bouncing back once towards the equator. We varied the tolerance of the step-size controllers and measured the relative error of the answer at the endpoint. For each simulation we have counted the number of steps taken.

A. List of methods

In this appendix we list few of the methods that have been described in the survey.

No attempt has been expanded to explain the methods, beyond references to the relevant earlier material.

A.1. RK-MK methods

All classical RK methods can be translated into Lie-group methods. Assume that the Butcher tableau

defines a Runge-Kutta method of order p in the classical sense [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF] and that φ is a map from g to G, for instance φ = expm, the exponential mapping, or φ = cay, the Cayley mapping for quadratic Lie groups. Then the corresponding order-p RK-MK algorithm for the Lie-group equation

for n ∈ N, and it is explicit provided that the underlying RK scheme is explicit.

The function dφ -1 (B, C, p) is a truncation of dφ -1 B (C) to order p -1, which is usually sufficient for a method of order p, given that the error is subsumed in the O h p+1 term. In some instances, like when g = so(3), the function dφ -1 B (C) can be evaluated exactly (see Appendix B).

Some popular schemes of the type (A.1) are Lie-group versions of

• forward Euler,

• the implicit midpoint rule,

• the trapezoidal rule,

The scheme (A.1) employs coordinates centred at Y n . Schemes with coordinates centred at an arbitrary point can be obtained in a similar manner,

for n = 0, 1, 2, . . .. In particular, for geodesic symmetric coordinates we have

while flow-symmetric coordinates yield

where the weights w l s are obtained integrating the Lagrangian cardinal polynomials, as in (7.7).

Herewith a number of important examples, originating in familiar classical RK methods.

• Order-four Gauss-Legendre scheme,

• Order-six Gauss-Legendre scheme, -1 24 (A.5)

• Order-four Gauss-Lobatto scheme, 0 0 0 0

Note that both in the case of the Gauss-Legendre of order six and Gauss-Lobatto of order four it is true that C h,n = -F 2 , hence C h,n need not be explicitly evaluated in (A.3). Algorithms (A.1-3) are practical for methods of order p ≤ 4, since for higherorder methods the computation of dφ -1 (A, B, p) requires the evaluation of a large number of commutators. For higher-order schemes it is recommended instead to use schemes based on graded Lie algebras, along the lines of Section 5.3.

Graded Lie algebras can be also effectively used to optimise existing RK-MK methods. The scheme

requires only two commutators instead of six if implemented using (A.1).

A.2. Higher-order methods for linear equations using graded Lie algebras

In this subsection we present methods for linear equations Y ′ = A(t)Y , whereby the number of commutators is reduced using graded algebras and a further technique due to [START_REF] Blanes | Improved high order integrators based on Magnus expansion[END_REF] and described briefly in Section 5.4. We do not distinguish between RK-MK and Magnus-type methods because, subject to these reductions, the two methods produce very similar results.

We consider collocation-type schemes. Denote by

the Vandermonde matrix generated by the vector d. In particular,

• set V = VDM(c), c T = (c 1 , c 2 , . . . c ν ) for non-symmetric collocation schemes;

• for symmetric collocation schemes, we take full advantage of symmetry, setting

where dϕ -1 is an order-p truncation to the dφ -1 -equation in the graded basis B 1 , . . . , B ν . For symmetric collocation schemes, the function dϕ -1 in (A.8) depends on terms of odd grade only. Specifically, for φ = expm,

• For order-six Gauss-Legendre scheme

(for order two and order four it suffices to consider the first one and two terms only, respectively, and the nodes of the corresponding quadrature), which can be evaluated with just four commutators using the technique of [START_REF] Blanes | Improved high order integrators based on Magnus expansion[END_REF],

producing an order-six time-symmetric truncation of dϕ -1 (B 1 , B 2 , B 3 ). • For order-four Gauss-Lobatto scheme (Ehle III A ),

• For order-six Gauss-Lobatto scheme, with nodes

which can be evaluated to correct order as

with just four commutators.

It may strike the reader that, after the transformation in the self-adjoint basis {B i }, the expression of Θ has always the same type of expansion, regardless of the choice of collocation nodes. This is no surprise, the information about the nodes being hidden in the basis elements B i . Changing into the self-adjoint basis implies that we integrate a Taylor-type expansion, in combination with a truncation of dexp -1 , which is independent of the nodes, depending only on the order of the method.

A.3. Methods using the Magnus expansion

Magnus-type methods are well suited to collocation-type techniques. This results in implicit schemes for nonlinear Lie-group equations.

• Order-four Gauss-Legendre:

6 , i = 1, 2. • Order-four Gauss-Lobatto:

Relaxing the collocation conditions, it is possible to obtain explicit methods.

• An explicit order-three scheme:

A.4. Magnus-type methods with geodesic/flow coordinates Magnus-type methods introduced in Appendix A.3 employ coordinates at Y n , hence they cease to be self adjoint for nonlinear problems. Below we describe their selfadjoint modification.

• Collocation order-four Magnus method with flow coordinates:

In both cases the a k;l,j s are evaluated according to (5.15) and the b l,j s and w l,j s are evaluated from (5.15) for θ = 1 and θ = 1 2 respectively. The weights w l are given by (7.7).

To conclude, it should be noted that all Lie-group methods for nonlinear problems require a number of exponential evaluations in the internal stages. This is a consequence of the fact that the funtion A(t, Z) may fail to be an element of g for arguments Z ̸ ∈ G. If A(t, Z) ∈ g for all matrices Z then exponentiations (or, with greater generality, evaluations of the map φ) in the internal stages may be disregarded, at the cost of a minor increase of local truncation error. If A(t, Z) ∈ g only when Z ∈ G, however, disregarding the evaluation of the map φ in the internal stages would compromise the assurance that the numerical approximation Y n+1 resides in G. However, as observed by [START_REF] Liu | Projected Runge-Kutta methods for differential equations on matrix Lie groups[END_REF], in some cases it is possible to devise simplified versions of the methods, whereby A(t, Z) is projected on g according to need. Specifically, setting A = A(t, Z), we note that

is a projector onto so(N), the algebra of skew-symmetric matrices,

is a projector onto sl(N), the algebra of matrices with zero trace, and

is a projector onto sp(N), the algebra of symplectic matrices. Using these projectors in the internal stages of Lie-group methods may significantly reduce the number of evaluations of the map φ.

Different ordering of C 1 , C 2 , C 3 leads to similar formulae. Note that, to avoid loss of significant digits for small θ, it is convenient to implement 1cos θ, a recurring factor in (B.16), using an angle-doubling formulae. The outcome, 1cos θ = 2 sin 2 ϑ, is more stable when numerical methods are implemented with small step size.