Variable stepsize commutator free Lie group integrators
Charles Curry Brynjulf Owren
math.NA

[
    Henry Kvinge, Can Ozan Oğuz, and Michael Reeks
    December 30, 2023
==================================================




We introduce variable stepsize commutator free Lie group integrators, where the error control is achieved using embedded Runge-Kutta pairs.





§ INTRODUCTION

Off the shelf computer software for the numerical solution of ordinary differential equations usually comes with built-in error control and variable step size. Typically, such codes compute an estimate of the local truncation error in each step. This estimate is then compared to a user specified tolerance, . If its norm is smaller than , the step will be accepted and a step size to be used in the succeeding step is computed. If the estimate is larger than , the step will be rejected and attempted again with a new reduced step size. Several adjustments to this procedure can be made in order to improve the behaviour of the integrator. For instance one can seek to avoid an excessive number of rejected steps, or an oscillating behaviour of the sequence of step sizes.

For Runge-Kutta schemes, the most popular method for obtaining an error estimate is by using the method of embedded pairs, see for instance the monograph <cit.> for details.
For a given initial point y_n  and a proposed step size h_n,  y_n+1 is computed with a method of convergence order p, together with an auxiliary second approximation ŷ_n+1 of convergence order p̂≠ p. The resulting pair is usually labelled p(p̂) in the literature.
Historically, it was customary to use p̂>p and take y_n+1-ŷ_n+1 as an estimate of the error in the lower order approximation y_n+1, a popular scheme of this form is the Runge–Kutta-Fehlberg 4(5) method. Yet, it seems unfortunate to not make use of the more accurate approximation, and it looks like the preferable choice these days is  to use p̂<p but adjust the step size changing formula in such a way that the global error behaves proportionally to the user given tolerance . This is done in several codes, such as the Dormand-Prince 5(4) method <cit.> implemented in the MATLAB solver  <cit.>.
An important feature of embedded Runge–Kutta pairs is that they share the same 
 internal stages. In this way, the cost of computing with the pair of method is not much higher than computing only the principal approximation y_n+1 in each step.
In the last few decades a subfield of numerical analysis called geometric integration has been established, see <cit.> for an exhaustive account. The main purpose is to develop numerical methods that preserve certain underlying geometric structures of the differential equation. Examples of such structures are symplecticity, volume, reversibility, and first integrals. In this paper we consider schemes which were designed for the situation where the differential equation has a natural formulation by means of a Lie group action on a smooth manifold, these schemes are called Lie group integrators, <cit.>, and the particular type considered here are called commutator-free Lie group integrators <cit.>.
Many of the aforementioned methods have primarily been implemented with constant step size.
In fact, for certain geometric properties of a scheme, such as symplecticity, it is not straightforward to vary the step size in such a way that the property is preserved.
Turning  now to Lie group integrators, it is for our purpose useful to divide them into two categories. One  consists of schemes which can be interpreted as applying a standard Runge–Kutta scheme to a local coordinate representation of the original differential equation, this category includes the Runge–Kutta–Munthe–Kaas methods. For such methods, error control and variable step size can be implemented simply by applying the technique described above to the coordinate representation. The second category, however, does not have such a natural coordinate representation, the schemes are typically constructed by composing flows of simple vector fields. In this case, embedded pairs can only be derived by solving order conditions for Lie group integrators. For a general and modern exposition of order conditions of such integrators, see for example <cit.>. A detailed account of the order conditions for commutator-free Lie group integrators was given in <cit.>.
The problem of deriving embedded pairs of efficient commutator-free schemes requires careful consideration. Not only is it important that the stages of the pair can be shared, in order for such schemes to be competitive to other Lie group integrators, they also need to reuse flow calculations to as large an extent as possible.

The aim of this paper is to show how one can derive embedded pairs of commutator-free Runge–Kutta methods with the smallest additional computational cost measured in terms of stages and flow calculations (exponentials) per step. Of course, it by no means clear that this way of measuring efficiency yields the smallest possible global error for a given computational cost. But we believe that our approach is a good starting point that illustrates in a clear way some of the challenges involved in constructing embedded pairs of commutator-free methods.
The rest of this paper is organised as follows: in section 2 we begin by briefly describing the class of schemes we consider, and in particular we review the order conditions which were derived in <cit.>. We then give particular examples of embedded pairs of type 3(2) and 4(3). These schemes are constructed such that the pairs share common stages, and employ the well-known “First same as last” (FSAL) property known from classical Runge-Kutta methods. We also design the pair to reuse flow calculations to a largest possible extent.
In the third section, we apply the CF43 method to three examples; first the Euler free rigid body, a standard test case for Lie group integrators. 
As a second example, we apply the method to the stiff Van der Pol equation, where we formulate the problem via the standard matrix-vector action by GL(2) on ℝ^2. Finally, we consider the heavy top, formulated via the coadjoint action of the Lie group SE(3) on the dual 𝔰𝔢(3)^* of its  Lie algebra.



§ COMMUTATOR FREE METHODS

We consider the numerical solution of autonomous initial value problems on manifolds. For a given manifold M and vector field F, that is a section of the tangent bundle TM, a smooth curve y(t) is an integral curve of F if it satisfies

    ẏ = F(y(t)),    t∈(a,b)

For any given starting location y(a)=p, the integral curve through p exists for sufficiently small time intervals. In the special case that M is a finite-dimensional vector space (which is identified with ℝ^n), integral curves may be approximated numerically using Runge-Kutta methods, i.e.

    k_i    =    f(y_n + h∑_j=1^s a^j_i k_j),   i=1,…,s 
    
    y_n+1   =    y_n + h ∑_i=1^s b^i k_i

It is not so apparent how to extend such methods to manifolds not possessing a linear structure. We begin by observing that the Euler scheme may be interpreted as a method which first `freezes' the vector field f at the point p=y_n so that f_p(y)=f_p(y_n) for all y, and then flows a distance h along the integral curve of the frozen vector field f_p, i.e.

    y_n+1 = exp(hf_p)y_n,
 
where exp(tf_p) is the one-parameter group of diffeomorphisms of ℝ^n induced by the vector field f_p. However, our notion of freezing the vector field at p relies on the identification of the tangent spaces T_p(M) and T_y(M) for any y. In the general setting, this requires a notion of parallel transport, which in turn requires a choice of linear connection. In general, if we freeze a vector field f by parallel transporting f_p(p)∈ T_p(M) onto T_y(M), the result is not independent of the path connecting p and y along which the transport is performed. For Lie group integrators, the connection is typically chosen to be flat so we can ignore this issue, but more generally assuming y is sufficiently close to p, the path could be taken to be the geodesic connecting y and p. In any case, the consequence is that exp(tf_p).p will coincide with the geodesic γ_t in the direction of f(p) from p.

Two strategies emerge from the above comments when generalizing Runge-Kutta methods to manifolds. The first <cit.> consists of considering the k_i as vector fields frozen at different points, parallel transporting these to a base tangent space T_p(M), taking a linear combination of the resulting tangent vectors and flowing along the geodesic in this direction. This strategy, when implemented using a Lie group action on a homogeneous space, results in the RKMK family of integrators, and allows the use of standard Runge-Kutta tableaux.

In some situations, it is advantageous to avoid using the linear structure on T_p(M) by moving entirely by composition of exponentials of frozen vector fields. This leads to the commutator-free Runge-Kutta methods:

    g_r    =   exp(h ∑_k α^k_r,J f_g_k)⋯exp(h∑_kα^k_r,1 f_g_k) p,   r=1,…,s 
    
    y_n+1   =   exp(h∑_k β^k_J  f_g_k)⋯exp(h∑_kβ^k_1  f_g_k) p




 §.§ Lie group integrators

The above methods are in principle very general, but require the computation of geodesics and possibly parallel transport, which is typically impractical. The observation underlying Lie group integrators is that the geodesics exp(hf_p) on a Lie group or reductive homogeneous space equipped with the canonical connection can be computed using the Lie group exponential. In practice, this typically means a matrix exponential, which can be approximated to machine accuracy with tolerable computational effort. Indeed, we suppose that the ODE we wish to solve may be written in the form

    y' = ( λ_* f(y))(y),   y(0)=p,

where f:M→𝔤, and λ_*:𝔤→𝒳(M) is the infinitesimal action arising from a group action Λ:G× M→ M as

    λ_*(u)(p) = .d/dt|_t=0Λ(exp(tu),p)

The vector field λ_* f is then frozen straightforwardly at a point p by taking

    λ_*f_p(y) = .d/dt|_t=0Λ(exp(tf(p)),y)

The exponentials of these frozen vector fields are seen to obey

    exp(tλ_* f_p)y = Λ(exp(tf(p)),y)

As a consequence, commutator-free RK Lie-group integrators take the form

    g_r    =   exp(h ∑_k α^k_r,J f_k)⋯exp(h∑_kα^k_r,1 f_k) 
    
    f_r    =    f(Λ(g_r,p)) 
    
    y_n+1   =   Λ(exp(h∑_k β^k_J  f_k)⋯exp(h∑_kβ^k_1  f_k),p)




 §.§ Order conditions


Let φ̂_h:M→ M be the mapping corresponding to taking a single timestep of a given commutator-free RK method, and φ_h the mapping which flows for time h along the solution curve through the initial point. The method is said to be of order p if, for all C^∞-functions ψ on M, we have

    ψ(φ_h(p))-ψ(φ̂_h(p)) = 𝒪(h^p+1)


A commutator-free RK scheme is in general specified by the collection of coefficients α^k_r,j,β^k_j. An order theory comprising systems of algebraic equations in the coefficients to be satisfied to attain a given order was derived in <cit.>, analogous to the theory of order conditions of standard RK methods. We will summarize the results without proof here. First, define the coefficients

    a^k_r = ∑_j α^k_r,j,   b^k = ∑_j β^k_j,   c_k = ∑_i a^i_k


A necessary condition for a commutator-free RK scheme with coefficients α^k_r,j,β^k_j to have order p is that the associated a^k_r, b^k are the coefficients of a standard Runge-Kutta method of order p or greater.


To attain a method of order 1 or 2, the above conditions are also sufficient. In particular, we can take J=1, such that there is only one exponential computed at each stage. On the other hand, we must satisfy an additional condition to attain order 3:


There are no order 3 methods employing exclusively J=1. Given that J=2, a sufficient condition to attain order 3 is that the coefficients a and b form a classical RK3 method, and in addition

    ∑_k β_1^k c_k + 1/2β_2^k = 1/3,
 


We see that we only require two exponentials for computation of y_n+1; one exponential suffices for the computation of g_r and f_r, i.e. we can take α^k_r,1=a^k_r. The situation changes for methods of order 4:


Given J=2, the order 4 conditions are the classical conditions together with the non-classical order 3 condition and the following:

    ∑_k β_1^k c_k + 1/3β_2^k    =   1/4
    ∑_k β_1^k c_k^2 + 1/3β_2^k    =   1/6
    ∑_j,kβ_1^k a_k^j c_j + 1/6∑_k β_2^k    =   1/12
    ∑_i,j b^i c_i α^i_j,1 c_j +  ∑_i,j,k b^i a_i^j c_i α^i_j,2 c_j    =   1/12


To satisfy these conditions, it satisfies to take two exponentials in one of the four intermediate stages, two in the final stage, and one exponential in other stages. We illustrate this by giving a sample tableau, displaying the coefficients of a method given in <cit.> which extends the classical RK4 method:

    [     0      ;   1/2   1/2      ;   1/2     0   1/2      ;   2*1   1/2     0     0      ;        -1/2     0     1      ;         1/4   1/6   1/6 -1/12;       -1/12   1/6   1/6   1/4 ]

This is a four-stage method, using two exponentials for the fourth stage, where the tableau displayed above is such that α_4,1=(1/2,0,0), and α_4,2=(-1/2,0,1).

To obtain a method of order 5, at least three exponentials are required for the final stage. There results a large nonlinear system of algebraic equations which has so far proved resistant to all attempts to attain a solution. Indeed, the construction of commutator-free methods of order 5 or higher remains an open problem.



 §.§ Embedded pairs reusing exponentials


Classical Runge-Kutta schemes are typically implemented as an embedded pair, i.e. two different sets of coefficients b^i, b̂^i are given such that the associated approximations y_n,ŷ_n are of different orders, typically p̂=p-1. In this section, we extend this idea to commutator free RK methods, devising schemes with coefficients β^i_j,β̂^i_j to create embedded pairs of orders differing by 1.

In practice, it is important to note that the dominant computational cost in the implementation of a commutator free method is likely to be the evaluation of Lie group (matrix) exponentials. Each unique horizontal row in the tableau generally requires the evaluation of one exponential. The observation of Owren <cit.> was that it is possible to design tableaux where certain rows coincide, such as the commutator free RK4 method above where the second and fourth rows are identical. This consideration looms large in the implementation of embedded pairs, where there is even greater opportunity for reuse of exponentials. A general CF3 scheme requires at least 3 exponentials, whilst a CF4 scheme requires 5. Our main achievement is the construction of CF32 and CF43 schemes using only one additional exponential and function evaluation compared to the constant stepsize scheme of the same order. 



  §.§.§ CF32

The order conditions to obtain a commutator free method of second order coincide with the classical RK2 conditions. As in the classical case the construction of a RK32 pair is always possible from an RK3 scheme, there is no difficulty in constructing a commutator free RK32 pair. In general, this will require one extra exponential and function evaluation. Suppose the tableau is constructed as follows

    [     0      ;   c_2  a_21      ;   c_3  a_31  a_32      ;   2*y β^1_1 β^2_1 β^3_1;       β^1_2 β^2_2 β^3_2;     ŷ  a_31  a_32    0, ]  or  [     0      ;   c_2  a_21      ;   c_3  a_31  a_32      ;   2*1 β^1_1 β^2_1 β^3_1      ;       β^1_2 β^2_2 β^3_2      ;   2*y β^1_1 β^2_1 β^3_1      ;       β^1_2 β^2_2 β^3_2      ;     ŷ  a_31  a_32     0    0, ]

implying a reuse of the third stage and hence saving an exponential. The two schemes above are identical, but the right hand tableau is of a more general type that allows ŷ to depend on the value of f(y). Schemes of this type have the FSAL (first same as last) property, as the first function evaluation at the next step coincides with the function evaluation f(y) assuming the step is accepted. We will henceforth write all tableaux in this form, and abbreviate the row for y as FSAL.

The order 2 conditions become a_31+a_32=1 and a_21a_32=1/2, hence any 3-stage commutator free RK3 scheme obeying

    a_32=1/2a_21,   a_31=2a_21-1/2a_21

would admit an embedded pair without the need for an extra exponential. Unfortunately, such a scheme does not exist, as is shown readily using symbolic computation software. We therefore focus on constructing commutator free RK3 schemes reusing an exponential, a topic as yet unexplored in the literature. All of the results stated in the remainder of this section were proven using symbolic computation software.

Let ω be a root of 36z^2 + (9a-30)z + 3a + 1. The general tableau of an RK32 scheme reusing the third stage in the second row of the fourth is given in Table <ref>.


Once a is chosen and a consistent choice of root for ω(a) is made, there then follows a system of two linear equations for the four coefficients of ŷ, in general resulting in a two-parameter family once the rest of the tableau is set. In general, the values will be irrational, but we highlight that some choices of a give rational coefficients, for instance a=-1/3,1/3,2/3,7/9,101/9,34/3,35/3 etc. In general these can be found by setting the discriminant 81a^2 - 972a + 756 to be a square. We give two samples of rational tableaux from the case a=1/3 in the second row of Table <ref>.






























We can also reuse the third stage in the first row of the fourth; indeed this may be preferable for some integrators as it permits the storage of the action of the exponential on y_0 rather than the exponential itself. Let ν be the root of 36z^2 + (9a-6)z -3a + 1. The general tableau for such a scheme is found in Table <ref>.

There are also one-parameter families of RK32 schemes reusing the second stage, either in the first or second row of y. Fix a and let γ be a root of 4a(3a-1)z^2 +4(3a-1)z + 3 and δ be a root of 4az^2 + (12a-2)z + 9a+6. The respective tableaux are given in Table <ref>.

The discriminants of the equations for γ and δ are 16(1-3a) and 4(1-36a) respectively, which allows for easy generation of rational versions of the above schemes should this be desired, for instance a=0 works in both cases.


  §.§.§ CF43

In the classical case, there are no RK43 pairs with only 4 stages, so it is essential to use a 5-stage FSAL scheme to attain RK43 at minimal cost. This is no longer true for commutator free methods, as it is possible to take a single CF4 method together with β,β̂ such that the underlying b=b̂, but that only β fulfills all of the non-classical order four conditions. Nonetheless, we will preserve the greater generality afforded by the FSAL schemes in our presentation.

An order three method requires two exponentials in the final computation, so in the generic case a CF43 pair requires two extra exponentials compared to the CF4 case. It is possible to reduce the extra cost to one by reusing an exponential, but the full tableau must be generated with this in mind, as a generic CF4 scheme does not permit a CF43 pair which reuses an exponential. In general, it is possible to reuse two exponentials (but not more) in a CF43 tableau, one above the lines (in the αs) and one below the lines (i.e., generating ŷ). Once the reuse pattern has been specified, the top of the tableau is fixed, and a one parameter family of β̂ coefficients is typically admitted. Note that patterns involving reuse of the second stage do not give CF43 pairs. In contrast to the CF32 case, there are no rational CF43 pairs reusing the optimal number of exponentials. We give the exact form of one of the tableaux in Table <ref>, but typically print only floating point forms; more accurate descriptions of the coefficients are available from the authors on request.





The same tableau is given below as a decimal approximation:




















    [              0               ;    4.785707347    4.785707347               ;    .8093268944    .7701000600   .03922683443               ;            2*1    .7701000600   .03922683443              0               ;                   .6195164818   .06934556872   -.4981889449               ;            2*1    .4211354919 -.005776103764   -.1381183969    .2227590088               ;                  -.1403784973  .006491728470    1.302163795   -.6682770264               ;              y           FSAL               ;              ŷ    .6195164818  0.06934556872   -.4981889449               ;                 -.158427746 a  0.008517658 a              a  -.618653791 a  -.231436127 a;                 - 0.075415454   -0.082788288                 + .5828295568  + .3847010797 ]




Another possible choice of reuse pattern leads to the next tableau, where the one parameter family for β_2 has been set such that β^3_2=0.

    [            0             ;                 .67104050             ;               2.547687640 -1.355037274             ;               2.547687640 -1.355037274            0             ;                -.21944181   -0.0735967     .1003880             ;            1   .324015249    .15832891   -.21057643  .2282322824             ;               -.108005081    .84426683    .44843513 -.6846968472             ;            y         FSAL             ;                -.21944181   -0.0735967     .1003880            0            0;                 .45603817    .93310478            0    -.2660264   0.06953371;              ]


Perhaps the optimal CF43 method in terms of reuse pattern is the following, in which the reuse always occurs in the first part of a split stage, and the method is a true 4-stage method not using FSAL. It is the only CF43 scheme with these properties:

    [              0               ;                   1.351207192               ;                           0.5    0.097900176               ;                           0.5    0.097900176               ;                   7.900943678    2.989500877   -10.48834473               ;              y     .301574869   -0.054881885     .238291289  0.01501572796;                  -.1005249562    .1005249562    .5450471839 -0.04504718389;              ŷ            0.5    0.097900176              0              0;                  -.2989500877  -0.0522571042     .783338473 -0.03003145592;                ]




§ PRACTICAL IMPLEMENTATION

We give an outline of how the above methods are implemented in practice. First we show how an embedded pair allows for automatic step size control, following <cit.>. Indeed, suppose we have chosen an initial step size h, and obtain approximate solutions y_1 and ŷ_1. If y takes values in a normed space, we can consider y-ŷ as an estimate of the error, and aim to ensure

    ||y_1-ŷ_1||<sc,   sc=Atol + max(||y_0||,||y_1||)· Rtol,

for some user-specified absolute and relative tolerances Atol,Rtol. If y takes values in a manifold equipped with only a metric, the notion of relative error does not make sense, but we can still aim to ensure that d(y_1,ŷ_1)<sc=Atol. In either case, let err=d(y_1,ŷ_1)/sc. In general, for an embedded pair of order p(p-1), the optimal step size is h_opt=h· err^-1/p. In practice, it is usual to dampen the fluctuations in h, a typically procedure is to let

    h_new = h·min(facmax,max(facmin,fac· err^-1/p)),

for some fac<1 and appropriate facmin,facmax. The above discussion is in some sense only rigorous when the lower order approximation ŷ is used to continue the integration, as only then is d(y,ŷ) a good measure of the local error. Nonetheless, practice has shown that it is usually better to continue with the higher order integrator, a procedure known as local extrapolation.

We have not addressed automatic selection of the initial step size. For problems in vector spaces, this is typically performed by constructing estimates of the derivatives of vector field F, using some combination of function evaluations and taking small step(s) with the Euler scheme, see <cit.>. Similar ideas may be employed for Lie group integrators, but require modification as derivatives cannot be approximated so simply using evaluations of vector fields at different points due to the difficulty of identifying nearby tangent spaces.



 §.§ Free rigid body


We now give a simple illustrative example of Lie group integration. Euler's equations for the body angular momentum of a free rigid body in a reference frame parallel to the principal axes of inertia take the form

    dξ/dt = -mI^-1ξ̂.ξ,

where m is the body mass, I is the diagonal inertia tensor, and we use the standard convention of the hat-map

    ξ=([ ξ_1; ξ_2; ξ_3 ])  ⇒  ξ̂ =([    0 -ξ_3  ξ_2;  ξ_3    0 -ξ_1; -ξ_2  ξ_1    0 ])

The quantity |ξ|^2 is readily shown to be conserved, i.e. equation <ref> describes the evolution of ξ∈ S^2. The equations allow for the immediate application of Lie group integrators, indeed the right hand side is of the form λ_* f(ξ)ξ, where the manifold M=S^2, f:S^2→𝔰𝔬(3) is the hat map, and the infinitesimal action is that of 𝔰𝔬(3) on S^2 by matrix multiplication. The associated group action Λ is SO(3) acting by matrix multiplication. In this context, the commutator-free RK Lie-group integrator given by the coefficients α,β is therefore

    g_r    =   exp(h ∑_k α^k_r,J f_k)⋯exp(h∑_kα^k_r,1 f_k) 
    
    f_r    =   g_r·ξ_n
    ξ_n+1   =   exp(h∑_k β^k_J  f_k)⋯exp(h∑_kβ^k_1  f_k)·ξ_n,

where g_r∈ SO(3), f_r∈𝔰𝔬(3), and ξ_n∈ S^2.




 §.§ The Van der Pol oscillator

We consider the following non-conservative oscillator with non-linear damping, formulated as a scalar second order differential equation

    ẍ - μ (1-x^2) ẋ + x = 0

where μ is a parameter that affects the stiffness of the system. One can rephrase this problem in the form (<ref>) using the simple matrix-times-vector action of the Lie group GL(2) on ℝ^2\{0}, 

    d/dt([ x; ẋ ])
    =([        0        1;       -1 μ(1-x^2) ]) ([ x; ẋ ])

The corresponding Lie group integrator can be interpreted as an exponential integrator in the sense defined for instance in
<cit.>. An explicit Lie group integrator cannot be expected to work well for stiff problems, but we believe it is still of interest to observe how the new embedded pair of commutator-free schemes behaves through the “needle” of the Van der Pol oscillator.

We have implemented the CF32 scheme given in the bottom left tableau of Table <ref> and applied it in a variable stepsize fashion to the Van der Pol oscillator (<ref>) setting μ=60 and y_0=(1,1)^T in all the experiments. In Figure <ref> we show the two components of the solution computed by CF32 in the top graph, and note in particular the sharp downward spike (“needle")
for the second component in the approximate interval t∈[1.4,1.56]. The relative and absolute tolerances were both set to 10^-3 in this experiment.
The stepsizes chosen by the CF32 schemes are shown in the bottom graph of Figure <ref>, and one can see how the stepsizes are reduced through the needle. For comparison, we also show the stepsizes used by the builtin Matlab solver ODE45 which is based on the Dormand-Prince embedded Runge–Kutta pair <cit.>. We observe that the new CF32 solver behaves similarly to the Dormand-Prince scheme except that the former takes larger steps, this might be expected due to the fact that it computes (exact) matrix exponentials. For matrices with eigenvalues whose real parts tend to -∞, such exact exponentials are bounded as opposed to their explicit Runge–Kutta counterparts which use polynomial approximations to the exponential map.

In Figure <ref> we visualise the difference between constant and variable stepsize for the problem (<ref>).
By running the variable stepsize code for a number of different choices for tolerances, we compute the global error at a fixed point t=1.6, just after passing through the needle. The exact solution has been approximated various different ways, one was to run various builtin Matlab solvers with strict tolerances.  As a measure for the efficiency of the integrator, we have computed  the number of exponential calcualtions used by the code to obtain a prescribed global error. The dashed line in Figure <ref> shows the result for the constant stepsize version of the third order method in our embedded pair, and the solid line shows the result for the variable stepsize solver. The cost ratio between the constant and variable stepsize methods depend on the chosen global error, but for instance to obtain a global error of 10^-5 the constant stepsize integrator needs approximately 6.5 times as many exponential calculations as the variable stepsize method.



 §.§ The heavy top

Mathematical models for the heavy top can be found in many text books, see for instance
<cit.>. Lie group integrators were applied to this problem in <cit.>. The heavy top is a rigid body, but because of the gravitational forces, it is not invariant under the action of SO(3) and the dynamics can therefore not be reduced to 𝔰𝔬(3)^* as the free rigid body.
There is a smaller symmetry group S^1 corresponding to rotation about the vertical axis and it turns out that the system can be formulated on 𝔰𝔢(3)^*, the dual of the Lie algebra of the special Euclidean group SE(3) consisting of translations and rotations in 3-space. According to <cit.> the equations can be written in the form

    μ̇   = μ×𝕀^-1μ + β× mgχ
    β̇   = β×𝕀^-1μ

Here μ is the body angular momentum and β is the vertical direction as seen from the rotating body, more precisely β=R^Te_3 where R is the attitude matrix of the top.
χ is the unit vector in the direction from the fixed point to the center of mass of the heavy top, m is the mass, g the constant of gravity, and 𝕀 is the inertia tensor.


In this example, the coadjoint orbits are preserved, and for this reason it is natural to invoke Lie group integrators via the (right) coadjoint action of SE(3) on 𝔰𝔢(3)^*. It is convenient to take elements of both 𝔰𝔢(3) and 𝔰𝔢(3)^* to be vectors in 
ℝ^3×ℝ^3. Similarly, elements of SE(3) are represented as pairs (g,u) where g could be an orthogonal 3× 3-matrix an u∈ℝ^3. 

The coadjoint action is the map Λ: SE(3)×𝔰𝔢^*(3)→𝔰𝔢(3)^*

    Λ((g,u),(μ,β))=^*_(g,u)(μ,β)=
     (g^T(μ-u×β),g^Tβ)

whose infinitesimal generator is the map λ_*:𝔰𝔢(3)→𝒳(𝔰𝔢(3)^*)

    λ_*(ξ,u)(μ,β)=^*_(ξ,u)(μ,β)
    =(-ξ×μ-u×β,-ξ×β)

For the heavy top equations (<ref>) we now have

    (μ̇,β̇)=λ_*(𝕀^-1μ,mgχ)(μ,β)

so that with reference to (<ref>) we have

    f(μ,β)  = (𝕀^-1μ,mgχ)



The key properties of the group action are summarised in Table <ref>.



plain


