Numerical Integrators that Preserve Symmetries and Reversing Symmetries.
We consider properties of flows, the relationships between them, and whether numerical integrators can be made to preserve these properties. This is done in the context of automorphisms and antiautomorphisms of a certain group generated by maps associated to vector fields. This new framework unifies several known constructions. We also use the concept of "covariance" of a numerical method with respect to a group of coordinate transformations. The main application is to explore the relationship between spatial symmetries, reversing symmetries, and time symmetry of flows and numerical integrators.
Introduction
Recently there has been a lot of interest in constructing numerical integration schemes
for ordinary differential equations (ODEs) in such a way that some qualitative geometrical
property of the solution of the ODE is exactly preserved. This has resulted in much
work on integration schemes that can preserve the symplectic structure for Hamiltonian
ODEs [12, 13, 14, 23, 33]. Other authors have constructed volume-preserving integrators
for divergence-free vector fields [3, 19, 30]. Other authors again have concentrated on
preserving energy and other first integrals [14, 20, 21] or other mathematical properties [7].
In this paper we are interested in constructing integrators that preserve the symmetries
and reversing symmetries of a given ODE. One reason this is important is that nongeneric
bifurcations can become generic in the presence of symmetries, and vice versa. One can
also consider the time step of the integration scheme as a bifurcation parameter, showing
how vitally important it is to stay within the smallest possible class of systems. Reversing
symmetries are particularly important, as they give rise to the existence of invariant tori
and invariant cylinders [15, 22, 25, 28].
One of the first authors to study integrators that preserve a reversing symmetry was
Scovel [24]. His method, as it stands, can only preserve a single reversing symmetry and no
ordinary symmetries. In this paper we give a construction of integration schemes preserving
an arbitrary number of symmetries and reversing symmetries.
In exploring this question, it became clear that constructions were working because
of the effect of certain operators on compositions of maps. This led us to express these
effects in terms of automorphisms and antiautomorphisms of groups. In particular, it turns
out that many different compositions designed to give methods particular features are all
examples of the same general principle. We prepare the basic algebraic material in section
2.
A curious structure emerges several times. Often a subset of a group is not closed under
products AB, but is closed under symmetric products ABA. This is true for symmetric or
antisymmetric matrices, reversible maps, and time-symmetric maps. We have called such
a subset a 'pseudogroup', and it is very useful to us here in integration, although it remains
to be seen whether it is an interesting algebraic object in its own right.
Flows of vector fields have many nice properties, and many of these are desirable in an
integrator as well. Some of these properties are defined for all vector fields, such as covariance
under transformation groups, closure under differentiation [2], and time symmetry
(see (27) below). These are discussed in section 3. Others are defined only for the flows
of a certain class of vector fields, such as symplecticity for Hamiltonian systems, closure
under restriction to closed subsystems [2], and energy conservation [14]. We believe that
in the context of this new viewpoint of numerical methods, it is natural to try to 'lift'
these properties so that they are naturally defined for all vector fields. For example, Ge's
definition of invariance with respect to symplectic transformations [5] is only defined for
Hamiltonian systems, but it can also be seen as a special case of general covariance. In
section 4 we study this lift for symmetries and reversing symmetries, and also look at
methods that make explicit use of the symmetries of the vector field.
The automorphism point of view elucidates the relationship between different prop-
erties, such as time-symmetry and reversibility, which are in fact independent, although
related in some special cases considered previously (Runge-Kutta methods with a linear
reversing symmetry by Stoffer [27], and splitting into reversible pieces by McLachlan [10]).
For example, we show in section 4 how a map covariant under a group larger than strictly
necessary can be more flexible under composition.
Algebraic preliminaries
Let G be a group with elements ' and consider functions G. We only consider
functions which are either automorphisms of the group, that is, A+ is a bijection and
or antiautomorphisms, that is, A \Gamma is a bijection and
G: (2)
Examples of automorphisms are the inner automorphisms
G. Examples of antiautomorphisms are inverse
and, if G is a linear group (a group of matrices), the transpose T '. (This can be generalized
to the quantum case where ' is a linear operator on a Hilbert space.)
Note that the set of automorphisms and antiautomorphisms itself also forms a group,
where the group operation is composition. The automorphisms form a normal subgroup.
More specifically, every such group G \Sigma is homomorphic to Z 2
there is an onto map
We call the number oe(A) the grading of A. In the case at hand, the automorphisms have
grading +1 and the antiautomorphisms have grading \Gamma1. Each such group G \Sigma is generated
by, for example, its antiautomorphisms alone (if it has any), or by one antiautomorphism
together with all the automorphisms.
Define the fixed sets
and
Notice that
but not. However, we do have that the set of group elements fixed by a given
antiautomorphism is closed under the symmetric triple product:
This concept seems to be useful in, for example, constructing integrators, and we shall call
it a pseudogroup 2 .
Example. Let G be the linear group GL(n) and consider the antiautomorphism T ,
transpose. The fixed set of this antiautomorphism is the pseudogroup of nonsingular
symmetric matrices: if X and Y are symmetric matrices, then so is XY X.
1 Unless there are no antiautomorphisms in the group.
A discussion of (what we call) the pseudogroup of maps with time-reversal symmetry, is given in [9].
Example. Let G be any group and consider the inner automorphism N / . The fixed
set of this automorphism is the group of so-called /-equivariant elements: if ' 1 , ' 2 are
/-equivariant (i.e., ' then so is ' 1 ' 2 .
If the antiautomorphism A \Gamma is an involution (i.e., A 2
there is a 'projection' to
because
first introduced in a special case by
Benzel, Ge, and Scovel [1].) This can be generalized to a group G \Sigma of automorphisms and
(not necessarily involutory) antiautomorphisms. Unfortunately, it is difficult to project an
arbitrary group element to fix(G \Sigma ); but QA \Gamma
does map the fixed set of the subgroup G+ of
automorphisms to the fixed set of the entire group.
Proposition: (Generalized Scovel projection) Let G \Sigma be a group of automorphisms
and antiautomorphisms, and let G+ be its subgroup of automorphisms. Let A \Gamma 2 G \Sigma be an
antiautomorphism, so that G
(The proof is given below.) This is not a true projection because it does not satisfy
but it will have some useful properties in our applications.
The 'projection' is not always surjective. For example, let take a group
consisting of one automorphism (the identity) and one antiautomorphism (transpose), that
is, G . Projection to the symmetric matrices produces
symmetric matrices with nonnegative determinant.
The projection Q and the pseudogroup property are in fact both examples of the
following more general relationship, which is of central importance in this paper:
Proposition: (Composition property) Let G \Sigma be a group of automorphisms
fA j+ g and antiautomorphisms of the group G. Let ' 1 ,
for all A
Proof: We only need consider the effect of an antiautomorphism; membership of fix(G
will follow because the antiautomorphisms generate all of G \Sigma .
A
for all k.
From this proposition the pseudogroup property follows by taking
the generalized Scovel projection Q follows by taking ' 2 to be the identity.
As an example of the composition property (12), take G g. All matrices B are
in matrices A are in fix(G \Sigma ), and indeed, BAB T is symmetric. The
reader is invited to check the consequences of the proposition for the sets of matrices given
in table 1.
We note that in most of this work, we can weaken the antiautomorphism requirement
to
G. This allows, for example, the operator A \Gamma
the antisymmetric matrices. The composition property (12) is retained, but the projection
(11) is lost if the identity is not in
Projecting onto the fixed set of automorphisms is more difficult. However, one may
have (as in the case of maps on R n ) that G is a near-ring: (G; +; :) is a near-ring if (G; +)
is a group, (G; :) is a semigroup, and ('
that maps under composition do not form a ring because they are not left-distributive.) If
the automorphisms are linear with respect to addition in the near-ring, and G+ is a finite
group with jG then we can use the following, which is a true projection:
The half-sized group construction
As mentioned above, groups G \Sigma are always homomorphic to Z 2 . In our application
to reversing symmetries, we shall use groups G+ of automorphisms which themselves are
homomorphic to Z 2 . From such a group we can construct a group H \Sigma of the same size as
G+ in the following way. Let A+;oe denote an element of G+ of grading
be a fixed antiautomorphism whose square A 2
\Gamma is an element of G+ with grading +1 which
commutes with G+ .
The automorphisms in H \Sigma are the elements A+;1 . The antiautomorphisms in H \Sigma are
the elements . It is straightforward to check that this gives a group H \Sigma of the
same size as G+ , i.e., half the size of G \Sigma .
Example. Let G be GL(2) and G+ be the inner automorphisms corresponding to a
subgroup of G homomorphic to Z 2 , say rotations by 3-=2. The antiau-
tomorphism I (inverse) satisfies the above requirements: its square (the identity) is in
G+ , does have grading +1, and does commute with G+ . The inner automorphisms N g;+1
together with the antiautomorphisms IN g;\Gamma1 form the group
This group comprises two automorphisms and two antiautomorphisms-it is Z 2 \Theta Z 2 .
Then in order for Q IN -=2
to be a projection to fix(H \Sigma ) we need ' to be fixed under all the
automorphisms in H \Sigma . Here we require
operator sign fixed set type 'projection'
I anti
IT auto AA orthogonal

Table

1. Sets of matrices as fixed sets of (anti)automorphisms. Here I is the identity
matrix and
\GammaI
3 General properties of integration methods
For simplicity we restrict our attention to systems of first-order autonomous ODEs on R n ,
i.e., ODEs of the form
for some n. Let F be a set of vector fields on R n (each vector field need not be on a space
of the same dimension), and let \Phi be a set of functions
An element ' 2 \Phi is called an integrator of f if
When we want to emphasize the functional dependence of ' on - and f we shall call it
a "method;" at other times we may want to think of fixing - and f and looking at the
resulting map.
Examples. If all f 2 F are differentiable, then the exact flow ' -
method Runge-Kutta method, and any one-step method which is
functionally dependent only on f are possible members of \Phi. If we impose additional
conditions on the set of vector fields F , then there is larger family of possible 's; e.g., if
all f 2 F are C r , then the Taylor series method of order r may be in \Phi. Moreover,
functions which are not integrators at all can also be in \Phi. Examples are ' - (f) : y 7! y,
and, (if all f are C r ), elementary differentials such as - 2 f 0 f .
To apply the results of the preceding section, a group is needed. To construct it we let
the set of maps
generate a group G by composition. For example,
We will not specify the sets F and \Phi precisely, as we wish to make the constructions
apply as generally as possible. Minimum requirements will be clear in each application.
Not all the elements of G associate maps with vector fields. For example, let ' - (f) be
an integrator of f . Then even though is an integrator of all vector fields
not an operator on vector fields. In fact the group G may seem at first
to be much bigger than necessary: one could also let, e.g., \Phi generate a group, with group
operation
But we wish to allow flexibility in the choice of time steps (for composition methods) and
vector fields (for splitting methods).
The functions A are now operators on the space of methods. They may or may not
preserve the property of being an integrator, and if they do, they may not preserve the
order of the integrator.
The exact flows of differential equations are important elements of G. It is their properties
that we are trying to mimic in constructing actual numerical methods. Some properties
(e.g., the semigroup property ' - (f)' oe hold and are defined for all f ; others
only hold and are defined for some f (e.g., volume-preservation det(d' -
divergence-free f ). Where possible, we try and extend the definition of these restricted
properties to all f ; but even with the guiding principle that they should hold for flows,
this extension cannot be unique. In this section we consider properties defined for all f , in
particular, those shared by fixed sets of automorphisms and antiautomorphisms of G. 3
Our first such property is "covariance". It is often useful to consider whether a method
produces the same integrator in different coordinate systems. This concept was studied for
symplectic transformations by Ge [5] with the aim of preserving integrals of Hamiltonian
systems by symplectic integrators. It was also used (with general linear transformations
and Runge-Kutta methods) by one of the present authors in [11]. We extend this concept
to general integrators and arbitrary groups of transformations acting on phase space.
Let H be a group with a left action (h; y) 7! hy on phase space R n 3 y. A method
- (f) is covariant with respect to H if the following diagram commutes 8h 2 H:
f=h   f
y x=hy \Gamma! ' - ( ~
The top arrow pushes forward the vector field f by h; the bottom arrow conjugates the
Writing out the transformations through a clockwise loop shows that
where
3 Note that if the antiautomorphism A \Gamma is such that flows are in is an integrator of
f of a certain order, then so is A \Gamma ' - (f). This is useful when composing as in (12).
(Compare the coordinate transformation automorphism K h to the inner automorphism
whose fixed set is the equivariant maps.) Now that covariance is
defined for \Phi by (24-26) (i.e., it is defined for the generators of G), the definition extends
to all of G in the natural way.
If ' - (f) is the time- flow of f , then ' is covariant with respect to the group of
diffeomorphisms of R n . In general, however, discrete approximations will only be covariant
with respect to the group of affine transformations or one of its subgroups. Different
methods can then be classified by their covariance group. Many traditional methods, such
as Euler's method, the midpoint and trapezoidal rules, and any Runge-Kutta or Taylor
series method, are covariant with respect to the affine group (the semi-direct product
of GL(n) and the translations R n ). (Some cases have been treated in [5, 27, 31].) Most
methods are covariant with respect to translations, but some are not covariant with respect
to the entire linear group but only with respect to a subgroup. For example, the known
volume-preserving methods [3, 19] are covariant with respect to the subgroup of diagonal
transformations. Splitting methods which partition the coordinates are covariant with
respect to arbitrary transformations which preserve the partitioning.
Our second important property is time-symmetry. A method ' is time-symmetric if
This important property, also known as 'symmetry' or `self-adjointness' or 'reversibility',
is considered in [4, 6, 17, 24, 32, 33]. (We reserve the name 'symmetry' for the spatial
symmetries considered in the next section.) Because the parameter - is included in G, we
can introduce the time-reversal automorphism
so that time symmetry means being in the fixed set of the antiautomorphism IR.
(One can also consider inverse-negative 's with ' -
time-scaling
's with ' c- (f=c). However, since there is a simple projection onto the fixed
set of the time-scaling operator, namely, ' - (f) 7! ' 1 (-f ), and the inverse-negative and
time-symmetry properties are then equivalent, we will not emphasize these properties.)
If our favourite integration method is not time-symmetric, we can use the projection
\Gamma-
to make it so. Now so the group G \Sigma generated by IR is just f1; IRg. It has
no nontrivial automorphisms so all ' are in fix(G + ) and (29) gives a time-symmetric map
for any '. (This projection was first given in [1].)
We can also use the composition property (12) to get time-symmetric methods of the
any method and ' 2 is time-symmetric, which has been used by many authors
starting with Suzuki [29] and Yoshida [33] for constructing higher-order methods.
Now consider the problem of constructing methods which are both covariant and time-
symmetric. Take a group of coordinate transformation automorphisms
and the antiautomorphism IR which together generate G (1)
\Sigma . Firstly, if ' is covariant with
respect to H (so that ' 2 fix(G
for all h 2 H. Secondly, if H, and hence G+ is homomorphic to Z 2 , then the half-sized
group construction gives a subgroup of G (1)
G (2)
with projection
where g is an element with grading \Gamma1 and ' 2 fix(K h ) for all elements h with grading +1.
All exact flows are in fix(G (1)
are desirable properties for integrators
as well. The advantage of also considering G (2)
\Sigma is that then the projection requires less of
'.
4 Integrators that preserve symmetries and reversing
symmetries
In this section we are interested in properties of integrators that do not hold for all vector
fields, but only for those with some given property, such as possessing a symmetry or
reversing symmetry. We first give some definitions:
The ODE (17) is S-symmetric if the vector field f satisfies
The ODE (17) is R-reversible if the vector field f satisfies
Here S and R are arbitrary diffeomorphisms of phase space, i.e., R is not necessarily an
involution.
The map ' is S-symmetric if ' 2 fix(N S ), i.e.
The map ' is R-reversible if ' 2 fix(IN R ), i.e.
If f is S-symmetric, then S-symmetry of ' - (f) is equivalent to S-covariance (' 2 fix(K S )).
If f is R-reversible, then R-reversibility of ' - (f) is equivalent to ' 2 fix(IRKR ).
Of course, a dynamical system may possess more than one (reversing) symmetry. The
set of all symmetries and reversing symmetries of a given dynamical system form a group
under composition, the reversing symmetry group \Gamma [8, 22]. It is homomorphic to Z 2 : the
symmetries have grading +1 and the reversing symmetries have grading \Gamma1. The flow
exp(-f) of f has the same reversing symmetry group as f .
The question we now address is: Given a vector field f , possessing a reversing symmetry
group \Gamma, how does one construct integrators ' - (f) that possess the same reversing symmetry
group? For the case of a single reversing symmetry R, Scovel [24] introduced the following
projection: If R
- R is R-reversible. This is another instance of the
projection (11), in this case, Q IN R
. (Just as for the antiautomorphism IR, the only
generated automorphism is the identity.) For the case of a larger reversing symmetry
group \Gamma, the half-sized group construction means that we can use 's which are not covariant
with respect to all of \Gamma, but only with respect to its symmetries. Then for any reversing
symmetry R in the group, the projection Q IN R
R possesses the full reversing
symmetry group \Gamma.
There are two possible extensions of this projection to all vector fields f . The first,
', has the reversing symmetry group \Gamma, regardless of whether the vector field does. In
particular, it is not an integrator of f even if ' is. Thus we prefer the following extension:
The conceptual advantage is all flows are in fix(IRKR ). This is exactly the group and
projection used in (33).
A drawback of the method Q IRKR ' is that it is not time-symmetric, and that making
it so by a second projection would destroy the desired reversibility: Q IRQ IRK R
There are three ways out.
The first is to use the composition property (12) to increase the order, ignoring time-
symmetry completely. For example, if ' - is reversible and of order 1, then one can check
that
is reversible and of order 2.
The second is to argue that time-symmetry does not affect the dynamics of the method
as much as reversibility. So we could first construct a non-reversible high-order method
(which may well be time-symmetric), and then form Q IRKR ', preserving the order but
not time-symmetry.
The third way is the most interesting. Abstractly, the problem is that we have two
antiautomorphisms A \Gamma and B \Gamma and would like to construct an element invariant under
them both. For the projection to work, the initial element ' must be invariant under all
automorphisms generated by A \Gamma and B \Gamma . Here we have ' 2
not ' 2
does not map to the fixed set of the full group of anti- and
automorphisms. However, if also ' 2 is easy to see that ' is fixed by
all generated automorphisms, and that projection with any generated antiautomorphism
gives the same Q'.
In the case at hand, A . Thus we are interested
in maps covariant under the reversing symmetry R as well, instead of just the ordinary
symmetry R 2 as previously:
Note that if f is R-reversible, its flow is R-reversible, R-covariant, and time-symmetric.
If a map satisfies any two of time-symmetry, R-reversibility, and R-covariance, then it also
satisfies the third. This gives a very flexible approach for constructing integrators, for
R-covariance-invariance under an automorphism-is preserved under arbitrary composi-
tions. At any stage one may project using Q IR to gain all three properties.
The case of a reversing symmetry group \Gamma is similar. We would now require full \Gamma-
covariance: reversing symmetries R in the
group. (In practice, one only checks covariance under a set of generators of \Gamma, such as the
reversing symmetries.) For example, affinely-covariant methods such as Runge-Kutta are
covariant with respect to any affine reversing symmetry group. So by asking a little more
of ' - (f ), we can directly apply the symmetric compositions of Yoshida [33] to increase the
order while still preserving reversibility.
Thus there are many routes to structure-preserving integrators. The preferred path will
depend on what additional structure (e.g. volume-preservation, integrals) the ODE has,
and on whether the symmetries and/or the reversing symmetries are linear. First find a
first-order method ' - (f ). If it is \Gamma-covariant, we are done. If \Gamma is linear, (see (15)) is
\Gamma-covariant, although this may break some structure in '. (For example, there is no known
linearly-covariant volume-preserving integrator.) The case when ' is covariant only with
respect to the symmetries is handled by the generalized Scovel projection, (38). Finally,
splitting methods in which the constituent vector fields are not reversible introduce some
new possibilities, discussed in the example below.
Example 1: Affine method.
One possible way of integrating an ODE whose reversing symmetry group is affine is
well known. It is not difficult to show that the implicit midpoint rule
is covariant with respect to the group of affine transformations (as indeed all Runge-Kutta
methods are). It follows that if the vector field f has an affine reversing symmetry group
\Gamma, then so does ' - (f ). Since the midpoint rule also has time symmetry (27), ' - (f) has
reversing symmetry group \Gamma. (The midpoint rule is also symplectic when the ODE (17)
is Hamiltonian ([23]), but not all Runge-Kutta methods with the time-symmetry property
are symplectic, nor vice versa.)
Example 2:
Consider the divergence-free vector fields in R 4 with parameter ff
@x
(\Gammay 3
@y
@
@z
@
which possess the following reversing symmetry group
the cyclic group generated by R 1 ). We want a method preserving the divergence-free
structure, the symmetries, and the reversing symmetries. The midpoint rule would
preserve the reversing symmetry group but is not, in general, volume preserving. It is shown
in [17] that the midpoint rule is volume preserving if and only if P
In R 4 , this requires the coefficient of - in P (-) to be zero. It can be
checked that this coefficient is not zero here.
Instead, we shall preserve volume by using a splitting method. Suppose
each f i is divergence free and possesses all (nonreversing) symmetries (but not necessarily
all reversing symmetries), and the solution of each -
can be approximated by a
preserving both volume and the symmetries and any reversing symmetries which f i
may have. Then we only need to combine the approximations so as to recover the reversing
symmetries. In general, there are three cases.
Firstly, if each f i is reversible, then the ' - (f i ) are reversible by assumption, so the
composition
is reversible-and, if ' - (f) is time-symmetric, (44) is also time-symmetric.
Secondly, if there is no structure when the f i are reversed, the best we can do is let
and apply the generalized Scovel projection, giving / - R/ \Gamma1
is irrelevant here. (Volume is preserved because the reversing symmetry group preserves
volume in this case.)
Thirdly, and perhaps most efficiently, suppose we can split into three pieces, f 1 , f 2 , and
with the following structure:
(The reversible piece f 2 may be split further if desired.) Let ' - (f 2 ) be reversible, and
integrator. Then
\Gamma-
is reversible. This situation occurs frequently when there is a reversing symmetry which is
linear and satisfies R \Sigma1. For then let f 1 be any vector field and define f 3 := \GammaRf
(because f is reversible by
assumption.) We only need to find an f 1 such that this splitting is useful.
In the present case we split into the planar vector fields
@
@x
@
@x
@y
@y
@
@z
so that each f i is divergence-free and hence Hamiltonian in two dimensions, and has symmetries
(because their coefficients are odd functions). The f 's also satisfy (46)
where the R i are given in (43). We let ' - be the midpoint rule,
which is volume-preserving for these Hamiltonian systems, and get the reversible method
This highlights the point that reversible integrators need not be time-symmetric, nor need
the pieces in a splitting method possess the reversing symmetries.
we have the nice situation of a splitting into two pieces, with
the direct composition of their flows (or time-symmetric, \Gamma-covariant approximations of
being reversible. Exactly the same constructions apply to a system with integrals,
once basic maps preserving the integrals and symmetries have been found.

Acknowledgement

. G.S.T. is grateful to the Australian Research Council for partial support
during the time this paper was written.



--R

Elementary construction of higher order Lie-Poisson integrators

Dynamical systems and geometric construction of algorithms

Equivariant symplectic difference schemes and generating functions
Solving Ordinary Differential Equations I
Chaotic Numerics
Reversing symmetries in dynamical systems
Symmetries and reversing symmetries in kicked systems
On the numerical integration of ordinary differential equations by symmetric composition methods
"Poisson schemes for Hamiltonian systems on Poisson mani- folds"
The accuracy of symplectic integrators
A survey of open problems in symplectic integration
Lectures in Mechanics
Convergent series expansions for quasi-periodic motions
Coexistence of conservative and dissipative behavior in reversible dynamical systems



Solving ODE's numerically while preserving a first integral
Solving ODE's numerically while preserving all first integrals
Chaos and time-reversal symmetry: order and chaos in reversible dynamical systems

Symplectic numerical integration of Hamiltonian systems
Reversible Systems
Symmetric two-step algorithms for ordinary differential equations
Variable steps for reversible integration methods
Numerical analysis of dynamical systems
Fractal decomposition of exponential operators with applications to many-body theories and Monte-Carlo simulations
Representation of volume-preserving maps induced by solenoidal vector fields


Construction of higher order symplectic integrators
--TR

--CTR
G. R. W. Quispel , D. I. McLaren, Explicit volume-preserving and symplectic integrators for trigonometric polynomial flows, Journal of Computational Physics, v.186 n.1, p.308-316, 20 March
R. I. McLachlan , M. Perlmutter , G. R. W. Quispel, Lie group foliations: dynamical systems and integrators, Future Generation Computer Systems, v.19 n.7, p.1207-1219, October
symmetries;automorphisms;numerical integrators