Backward Error Analysis for Numerical Integrators.
Backward error analysis has become an important tool for understanding the long time behavior of numerical integration methods. This is true in particular for the integration of Hamiltonian systems where backward error analysis can be used to show that a symplectic method will conserve energy over exponentially long periods of time. Such results are typically based on two aspects of backward error analysis: (i) It can be shown that the modified vector fields have some qualitative properties which they share with the given problem and (ii) an estimate is given for the difference between the best interpolating vector field and the numerical method. These aspects have been investigated recently, for example, by Benettin and Giorgilli in [ J. Statist. Phys., 74 (1994), pp. 1117--1143], by Hairer in [Ann. Numer. Math., 1 (1994), pp. 107--132], and by Hairer and Lubich in [ Numer. Math., 76 (1997), pp. 441--462]. In this paper we aim at providing a unifying framework and a simplification of the existing results and corresponding proofs. Our approach to backward error analysis is based on a simple recursive definition of the modified vector fields that does not require explicit Taylor series expansion of the numerical method and the corresponding flow maps as in the above-cited works. As an application we discuss the long time integration of chaotic Hamiltonian systems and the approximation of time averages along numerically computed trajectories.
Introduction
. In this paper, we consider the relationship between solutions
to a given system of ordinary differential equations (vector fields)
d
dt
numerical approximations
to them, and solutions to associated modified equations
d
dt
The vector fields ~
are formulated in terms of an asymptotic expansion in the
step-size ffit, i.e., are chosen such that the numerical solution can formally be inter-
preted, with increasing index i, as the more and more accurate solution of the modified
equation. Previous papers on backward error analysis for differential equations include
those by Warming & Hyett [38], Griffiths & Sanz-Serna [12], Beyn [6],
Feng [10], Fiedler & Scheurle [11], and Sanz-Serna [31]. Another early reference
to related ideas is Moser [24] who discusses the approximation of a symplectic
map near an equilibrium by the flow map of a Hamiltonian vector field.
More recently, general formulas for the computation of the modified vector fields
~
derived by Hairer [16], Calvo, Murua & Sanz-Serna [7], Be-
nettin & Giorgilli [5], and Reich [26]. In papers by Neishtadt [25], Benettin
& Giorgilli [5], and Hairer & Lubich [18], the question of closeness of the numerical
approximations and the solutions of the modified equations has been addressed.
Konrad-Zuse-Zentrum, Takustr. 7, D-14195 Berlin, current address: Department of Mathematics
and Statistics, University of Surrey, Guildford, Surrey GU2 5XH, UK (S.Reich@surrey.ac.uk).
In particular, it has been shown in these papers that the difference can be made
exponentially small in the step-size ffit, i.e.
provided the vector field Z and the numerical one step method \Psi ffit are real analytic
are appropriate constants, \Phi ffit; ~
denotes the time-ffit-flow map
of the vector field ~
, and the index i   (ffit) has been chosen such that the difference
is minimized.
Backward error analysis is of utmost importance for an understanding of the
qualitative behavior of symplectic methods [32] for Hamiltonian problems. It has
been shown by Hairer [16], Calvo, Murua & Sanz-Serna [7], Reich [26], and
Benettin & Giorgilli [5] that for symplectic discretizations, the modified vector
fields ~
are Hamiltonian. For special cases see also the papers by Auerbach &
Friedman [4] and Yoshida [40]. The Hamiltonian structure of the modified equations
implies that a symplectic integrator almost preserves the total energy over an
exponentially long period of time [25, 23, 5, 18]. Similarly, the adiabatic invariant of
a Hamiltonian system with rapidly rotating phase is also preserved over an exponentially
long periods of time provided a symplectic method is used [29].
The fact that symplectic methods lead to modified equations that are Hamiltonian
is a special instance of the so called geometric properties of backward error analysis.
By this we mean the following: If the vector field Z belongs to a certain class of
vector fields, like integral preserving or divergence-free vector fields, and the numerical
approximation \Phi ffit also preserves the corresponding quantities, then the modified
vector fields ~
will be in the same class as Z. Besides symplectic methods, special
instances of these geometric aspects have been discussed before. See, for example, the
papers by Reich [26], Hairer & Stoffer [20] and Gonzalez, Higham & Stuart
[13].
In this paper, we revisit backward error analysis by using a simple recursive
scheme for the definition of the modified vector fields ~
first proposed by the
author in the unpublished report [26]. The main advantage of this formulations is
that it does not require Taylor series expansions of the numerical one step method
\Psi ffit and the flow maps \Phi ffit; ~
in terms of the step size ffit as it used in the papers
by Benettin & Giorgilli [5], Hairer [16], and Hairer & Lubich [18]. This in
turn allows for a simple characterization 1 of the geometric properties of the modified
vector fields ~
and a rather simple proof for the exponentially small truncation
error (1.1). Our approach is close to the one discussed by Benettin & Giorgilli [5]
in the sense that we consider general one step methods 2 and that we use a "direct"
approach 3 . However, different techniques are used and we will discuss this in more
detail in x4.
In x5, we consider the numerical integration of a "chaotic" Hamiltonian system
by a symplectic method and discuss the approximation of time-averages along numerically
computed trajectories. We assume that a Poincar'e section [14] can be defined
and that the corresponding Poincar'e section is uniformly hyperbolic. Backward error
1 The general idea can already been found in the report [26].
Hairer & Lubich [18] consider methods that can be represented by P-series [19]. Note that
Runge-Kutta and partitioned Runge-Kutta methods fall under this category.
3 This is in contrast to the "indirect" approach used by Neishtadt [25] where the one step method
is first interpolated by the flow of a time-dependent vector field and averaging in time is then used
to obtain an optimal approximating time-independent vector field.
analysis and the shadowing lemma [33] will be used to show that a numerically computed
trajectory stays close to an exact solution over exponentially long periods of
time. This and a large deviation theorem [36] allow us to discuss the convergence of
long time averages along numerically computed trajectories. The anisotropic Kepler
problem [15] will serve us as a numerical illustration. This problem requires the application
of a symplectic variable step-size method as first discussed by the author in
the technical report [27] and independently by Hairer in [17].
2. The Modified Vector Field Recursion. Let us consider a smooth vector
field
d
dt
its discretization by a one step method [19]
We assume that \Psi ffit : U ae R n ! R n is a smooth map and a method of order p  1,
i.e.
for all x 2 U where \Phi ffit;Z is the time-ffit-flow map of the differential equation (2.1).
As in described in the Introduction, we look for a family of vector fields ~
such that
\Phi ffit; ~
or, equivalently,
for all ffit sufficiently small. Here \Phi 1;X denotes the time-one-flow map of the vector
field X(ffit). The family of modified vector fields X(ffit), ffit  0, is formally defined in
terms of an asymptotic expansion in the step-size ffit, i.e.
The formally infinite sequence of vector fields f\DeltaX i g i=1;::: ;1 can be obtained by
Taylor series expansion of the one step method \Psi ffit , i.e.
x the identity map, and comparison of this series with the expansion of the
time-one-flow map \Phi 1;X(ffit) in terms of ffit. The vector fields \DeltaX i are chosen such
that these two series coincide term by term. This is the general approach followed
by Benetttin & Giorgilli [5] and Hairer [16]. The two papers differ in the
way the Taylor series expansions are written down. But they lead to exactly the
same sequence of vector fields f\DeltaX i g i=1;::: ;1 . We obviously have \DeltaX
for a method of order p.
We now give a recursive definition of the modified vector field X(ffit) that does
not require an explicit Taylor series expansion. This recursion was introduced by the
author in the unpublished report [26]. First we formally introduce the "truncated"
expansions by means of
We obviously have
Let us assume that X i (ffit) has been chosen such that the difference between the time-
one-flow map of X i (ffit) and the numerical one step method \Psi ffit is O(ffit i+1 ). This
suggests to consider the following recursion:
\DeltaX
Indeed, this definition of \DeltaX i+1 implies that X i+1 (ffit), defined by (2.3), generates a
time-one-flow map that is O(ffit i+2 ) away from the numerical method \Psi ffit . This can
be seen from
Thus eqs. (2.3) and (2.4) recursively define the modified vector fields X i
1. The recursion is started with X 1 ffitZ. The generated sequence
is, of course, equivalent to the sequences obtained by using Taylor
series expansions as described in [16, 5].
Throughout this paper, we will exclusively work with the recusion (2.3)-(2.4).
In x3, it will be shown that this leads to a simple characterization of the geometric
properties of the modified vector fields and, in x4, explicit estimates for the difference
between the time-one-flow map of the modified vector field X i (ffit) and the numerical
method will be given. We like point out that these results can also be (and have
been [16, 5, 18, 20, 13]) derived using explicit Taylor series expansion of the flow map
and the numerical method. However, we feel that the application of the recursion
leads to a simplification in the presentation of these results.
3. Geometric Properties of Backward Error Analysis. In this section, we
consider differential equations (2.1) whose corresponding vector field Z belongs to a
certain linear subspace g of the infinite dimensional Lie algebra 4 of smooth vector
fields on R n [21],[1].
Assumption. Given a linear subspace g of the infinite dimensional Lie algebra
of smooth vector fields on R n , let us assume that there is a corresponding subset G of
the infinite dimensional Frechet manifold [21] of diffeomorphisms on R n such that
4 The algebraic operation is the Lie bracket [X;Y ] of two vector fields X and Y [3].
Here T id G is defined as the set of all vector fields X := @  [\Psi  ] =0 for which the
one-parametric family of diffeomorphisms \Psi  2 G is smooth in  and \Psi
For the linear space (Lie algebra) of Hamiltonian vector fields on R n this is, for
example, the subset of canonical transformations [1]. An important aspect of those
differential equations is that the corresponding flow map \Phi t;Z forms a one-parametric
subgroup in G [21],[1]. Especially in the context of long term integration, it is desirable
to discretize differential equations of this type in such a way that the corresponding
iteration map \Psi ffit belongs to the same subset G as \Phi t;Z . We will call those integrators
geometric integrators.
The following result concerning the backward error analysis of geometric integrators
has been first stated in the unpublished report [26]:
Theorem 3.1. Let us assume that the vector field Z in
d
dt
belongs to a linear subspace g of the Lie algebra of all smooth vector fields on R n . Let
us assume furthermore that
is a geometric integrator for this subspace g, i.e., \Psi ffit 2 G for all ffit  0 sufficiently
small. Then the perturbed vector fields X i defined through the
recursion (2.3)-(2.4) belong to g, i.e.
Proof. The statement is certainly true for us assume that it
also holds for for all ffit  0 sufficiently small. Since
and
for all ffit  0 sufficiently small as well as
we have
and \DeltaX g. This implies X i+1 (ffit) 2 g as required.
Remark. Often the linear subspace g is, in fact, a subalgebra under the Lie
bracket [3]
@
g. But this property is not needed in Theorem 3.1.
Let us discuss five examples:
Example 1. Consider the subspace g of all vector fields that preserve a particular
first integral F : R n ! R. In fact, this space is a subalgebra under the Lie bracket
(3.1). In other words
and
imply that
To show this we differentiate (3.2) w.r.t. x which gives
The same procedure is applied to (3.3). Using these identities and the definition
(3.1) in (3.4) yield the desired result. The corresponding subset G is given by the
F -preserving diffeomorphisms \Psi, i.e.
In fact, let \Psi  be a smooth family of F -preserving diffeomorphisms with \Psi
then X := @  [\Psi  since
Thus, T id and we can apply Theorem 3.1. In particular, if a numerical method
then the modified vector fields X i possess F as a first integral.
The same result was recently derived by Gonzalez, Higham & Stuart [13] using
a contradiction argument. \Pi
Example 2. Consider the Lie subalgebra of all divergence-free vector fields Z, i.e.
The corresponding subset G are the volume preserving diffeomorphisms,
i.e.
det
@
Again we have T id g. Namely:
@
@x \Psi  (x)
trace [@ x@  \Psi  (x)] =0
Thus, if the numerical method \Psi ffit is volume conserving, then
the modified vector fields X i are divergence-free. Again, the
same result has been formulated by Gonzalez, Higham & Stuart [13] using a
contradiction argument. \Pi
Example 3. Let an involution 5 S 2 R n\Thetan be given and consider the subspace g
of vector fields Z on R n that satisfy the time-reversal symmetry
This subspace is not a subalgebra under the Lie bracket (3.1). The corresponding
subset G is given by the time-reversible diffeomorphisms \Psi, i.e. \Psi
Let \Psi  2 G be smooth in  with \Psi
\Theta S\Psi
which implies that X := @  [\Psi  g. It follows that T id and we can apply
Theorem 3.1. Thus, if a numerical method \Psi ffit satisfies the time-reversal symmetry,
then the modified vector fields X i are time-reversible. This result
has been first proven by Hairer & Stoffer in [20] using ideas from [26]. \Pi
Example 4. Let f:; :g denote the Poisson bracket of a (linear) Poisson manifold
. Then the Lie algebra of Hamiltonian vector fields on P is given by
d
dt
R is a smooth function. The corresponding subset G is given by the
set of smooth diffeomorphisms on P that preserve the Poisson bracket f:; :g [1]. Let
\Psi  be a family of maps in G with \Psi
for all smooth functions F; . This is the condition for
a vector field X to be locally Hamiltonian. Since P is simply connected, the vector
field is also globally Hamiltonian [3].
If the discrete evolution (2.2) satisfies \Psi ffit 2 G for all ffit ? 0, then \Psi ffit is called
a symplectic method and it follows from Theorem 3.1 that the modified vector fields
are Hamiltonian vector fields on P . This result can also be
found in [5, 16, 26].
If a symplectic method can be expanded as a P-series, then the vector fields X i (ffit)
are globally Hamiltonian even if the phase space P ae R n is not simply connected
[16]. This result applies to all symplectic Runge-Kutta and partitioned Runge-Kutta
methods. Furthermore, symplectic methods defined by a generating function of the
third kind [32] are also always globally Hamiltonian [5]. The same statement is true
for symplectic methods based on the composition of exact flow maps [40]. \Pi
Example 5. Let us now consider differential equations on a matrix Lie group
G ae R n\Thetan [35]. In general, time independent differential equations on G can be
written in the form
d
dt
5 An involution is a non-singular matrix that satisfies S
ae R n\Thetan the Lie algebra of G. Many recent papers (see [8]
and references therein) have been devoted to methods that preserve the Lie group
structure, i.e.
and Y n 2 G implies Y n+1 2 G. Thus, \Psi ffit is a diffeomorphism defined on the
submanifold G ae R n\Thetan . In fact, this submanifold can be characterized, at least
locally, by a set of nonlinear equations which we denote by F (Y
is an F integral preserving map, i.e. F ffi \Psi G. Following Example 1, we
know then that the modified vector fields (as well as the given vector field) satisfy
Hence the modified vector fields X i (ffit) are vector fields on G and
give rise to modified differential equations of type
d
dt
with
and ~
g. See [28] for further results on backward error analysis for
numerical methods on manifolds. \Pi
4. Truncation Error of Backward Error Analysis. We like to derive an
explicit estimate for the norm of the vector fields \DeltaX i+1 and the difference between
the time-one-flow map \Phi 1;X i (ffit) and the numerical approximation \Psi ffit ,
To do so we assume from now on that the vector field Z in (2.1) is real analytic. We
also introduce the following notations: Let B r denote the complex ball of
Let us consider a compact subset K ae R n of phase space and a constant r ? 0 such
that a given real analytic vector field Y is bounded on B r
we define
with
We also define B
To find an estimate for \DeltaX i+1 , as defined in (2.4), we need estimates for the mappings
appearing on the right hand side of (2.4). We start with an estimate for the map \Psi ffit .
Lemma 4.1. Let us assume that the vector field Z in (2.1) is real analytic and
that there is a compact subset K of phase space and constants K; R ? 0 such that
We also assume that the numerical method \Psi ffit is real analytic. Then there exists a
constant M  K such that
Proof. Under the given assumptions, the flow map
Z Z(\Phi t;Z (x)) dt
is defined for complex valued  2 C where the integral on the right hand side is
independent of the path from zero to  . The complexified flow map satisfies
Z jjZ(\Phi t;Z (x))jj jdtj
1). Consistency of the numerical method implies that there exists a constant
\DeltaK ? 0 such that
for the (complexified) map \Psi  . Take M := K
Remark. Let us consider a s-stage Runge-Kutta method with coefficients
i;j=1;::: ;s and fb i g i=1;::: ;s [19] satisfying
s
s
d  1, and assume that the Runge-Kutta method uniquely 6 defines a real analytic
map \Psi  for all step-sizes  2 C with j j  R=K. Then we have
4.1. This follows from the fact that, under the stated assumptions, all stage variables
will be in BRK, where the vector field Z is bounded by the constant K. A similar
statement holds for partitioned Runge-Kutta methods. \Pi
Theorem 4.2. Let the assumptions of Lemma 4.1 be satisfied. Then there exists
a family of real analytic vector fields ~
, such that
with the order of the method. The family
of modified vector fields ~
X(ffit) satisfies the estimate
R
6 For an implicit method, the solution can be obtained by fixpoint iteration if j j is sufficiently
small.
with d p  1 a constant depending on the order p of the method. For example,
Proof. We know that X 1 and that \Psi ffit is a method of order p  1.
Thus jjX 1 ()jj R  j jM and \DeltaX Next we find an estimate for
the difference between the time-one-flow map \Phi 1;X 1 () and the map \Psi  . Using (4.1)
and (4.2), we obtain
Since the mappings are real analytic and their difference
is O(ffit p+1 ), we obtain the estimate [5]
1). Using this in (2.4) with
Next we show that
for 1). The estimate is true for
(compare (4.3)). We proceed by induction. First note that, for
We replace the parameter ff 2 [0; 1) in this formula by ff+
where
(j
0:067
This yields
for all ff 2 [0; 1) and all  2 C with
(j
Here we have used that
0:891 (4.6)
for 2. In particular, substitute  j for j j and ff
in (4.5). Then use the and the inequality
(4.6) to derive
1). Next we introduce the vector-valued, real analytic function
and observe that
1). Here we have used that x 2 B ffR K implies \Phi t;X j () (x) 2
Now we can find an estimate for
the difference between the time-one-flow map \Phi 1;X j () and the map \Psi  . Using (4.1)
and (4.7), we obtain
. Since the mappings are real analytic and their difference is O(ffit j+1 ), we
obtain the estimate [5]
1). Using this in (2.4) with finally obtain
which verifies (4.4) for
Next we need an estimate for the difference between the time-one-flow map
and the map \Psi ffit on the compact set K. Using (4.8) with
immediately have
R
The family of vector fields ~
X(ffit) is now defined by taking an optimal number i   (ffit)
of iterations. We take i   (ffit) as the integer part of
Thus
(ffit). This completes the first part of the
proof.
According to (4.5), the difference between the modified vector fields ~
X(ffit) and
Z is given by
R
R
Next we use
ffit  R
c (i
to obtain
R
e
R
R
Here d p  1 is chosen such that
for all j  2.
Remark. The proof of Theorem 4.2 is similar in spirit to the one given by
Benettin & Giorgilli [5] on the exponentially small difference between an optimal
interpolating vector field and a near-to-the-identity map. However, there are a couple
of important differences: (i) We explicitly take the order of a method into account.
(ii) We directly derive estimates on the difference between the flow maps \Phi 1;X i () and
\Psi ffit instead of using Taylor series expansions of \Phi 1;X i () and \Psi ffit and corresponding
estimates for the elements in the series. We believe that this simplifies the proof of
Theorem 4.2. (iii) By introducing the parameter ff 2 [0; 1), we do not have to shrink
the domain of definition of the vector fields X i () as the iteration index i increases.
Again we feel that this simplifies the proof. (iv) As in [18], we work directly with an
estimate for the given vector field Z instead of making assumptions on the map \Psi ffit .
The rather pessimistic constants entering the estimates seem the
main disadvantage of our approach. \Pi
A more elaborate version of the proof of Theorem 4.2 can be found in [30].
5. An Application: Ergodic Hamiltonian Systems. Let us consider a (real
analytic) Hamiltonian system
d
dt
d
dt
together with a smooth function A : R 2n ! R. We are interested in
evaluating the time-average of A along a trajectory (q(t); p(t)) of the Hamiltonian
system (5.1)-(5.2), i.e
We assume that
exists and is equal to the micro-canonical ensemble average corresponding to the
Hamiltonian
i.e., we assume that the system (5.1)-(5.2) is ergodic 7 (or even mixing) [37]. Thus
R A(q; p)
Z
and
Z
the inner product of A and ffi(E \Gamma H).
Let us write the equations (5.1)-(5.2) in more compact form as
d
dt
. The Hamiltonian H is preserved under the flow map \Phi t;H .
Let us assume that the hypersurface M 0 of constant energy
7 To be more precise: Ergodicity of a system implies that the time average is equivalent to the
ensemble average except for, at most, a set of initial conditions of measure zero.
is a compact subset of R 2n . We also assume that there is a constant
that jjrxH(x)jj ? fl 1 for all x 2 M 0 . This implies that M 0 is a smooth
dimensional compact submanifold. Furthermore, the family of hypersurfaces
sufficiently small, are smooth and compact as well (in fact diffeomorphic to
We define the open subset U of phase space by
U :=
E2(\Gamma\DeltaE;+\DeltaE)
So far we have made fairly generic assumptions. In the sequel, we become more
specific to ensure that the Hamiltonian system (5.1)-(5.2) is ergodic/mixing.
In a first step we construct a Poincar'e return map [14]. Let R be a
smooth function and positive constant such that jf/; Hg(x)j ? fl 2 on the
level sets
\Deltas ? 0 sufficiently small. Let us assume that S s defines a Poincar'e section for each
in the following way: For all x 2 S s , there is a positive number
such that the solution x(t), t  0, with initial condition
s and there is no
such that x(t 0
. The positive number
is called the Poincar'e return time of the point x Knowing the Poincar'e
return time for each x 2 S s , we define the "global" Poincar'e map
and
We assume that the Poincar'e return times t p (x), x 2 V , are bounded by some constant
We are interested in the solutions on a particular level set of constant energy. For
simplicity, we take the level set M 0 . Then it is sufficient to consider the "restricted"
which is defined as the restriction of \Pi to
Thus we have reduced the study of the dynamical properties of the Hamiltonian system
(5.1)-(5.2) on the energy shell M 0 to the study of the properties of the Poincar'e
is an ergodic (mixing) map, then the Hamiltonian system is ergodic
(mixing) on M 0 . Note that \Pi 0 is volume preserving, i.e. det @x \Pi 0
From now on we assume that \Pi 0 is a uniformly hyperbolic map, i.e., for each
D, the linearization @x \Pi 0 (x) at x possesses strictly expanding and contracting
directions only [14],[36]. The "stochastic" behavior of such a (deterministic) map has
been investigated, for example, in [36]. Here we only point out the four main results:
ffl There is a unique invariant density  0 on D that is invariant under \Pi 0 .
Furthermore,  0 is given by the Lebesgue measure on D.
ffl The autocorrelation function hA of a Holder continuous function
exponentially fast, i.e.
an appropriate constant.
ffl The time-averages
of A along trajectories fx i g i=1;::: ;N of \Pi 0 satisfy a central limit theorem.
ffl The time-average hAiN of A along trajectories of \Pi 0 with initial value x
satisfy a large deviation theorem. To be more specific [39]: Given any c ? 0
there is a h(c) ? 0 such that
for all large N  1.
These results can be proven (see, for example, [36]) by carefully studying the properties
of the corresponding Frobenius-Perron operator defined
by
Definition 5.1. We call a Hamiltonian system (5.1)-(5.2) with the above introduced
properties Poincar'e hyperbolic. In particular, we assume (i) that the level
sets of constant energy are compact submanifolds, (ii) that
there is a constant that a global
Poincar'e map \Pi can be defined on
which is uniformly hyperbolic as a map restricted to that the
Poincar'e return times t p (x), x 2 V, are bounded by some constant K ? 0, and (v)
that there is a constant
Let e
H be a perturbation of H such that
for all x 2 U and some ffl ? 0. Then we call e
H an ffl-perturbation of H.
Lemma 5.2. The property of being Poincar'e hyperbolic is stable under ffl-
perturbations of the Hamiltonian H provided ffl is sufficiently small.
Proof. The assumption jjrxH(x)jj ? fl 1 on the level sets ME implies that these
sets are persistent under small perturbations. Furthermore, there exists a constant
~
perturbed Hamiltonian ~
H and x 2 V . Thus a
Poincar'e map is also defined for the perturbed Hamiltonian ~
H . Uniform hyperbolicity
is also stable under small perturbations of the Poincar'e map [2].
Let us discretize (5.1)-(5.2) by a symplectic (real analytic) integrator \Psi ffit of order
1.
Assumption. We assume that backward error analysis can be applied on a compact
subset K with U ae K. The corresponding perturbed Hamiltonian is denoted by
e
i.e., for all x 2 K,
constants. Let the step-size ffit be sufficiently small such that
the perturbed Hamiltonian system is also Poincar'e hyperbolic. For simplicity, we shift
the modified Hamiltonian e
H(ffit) such that H(x 0
Let us introduce a couple of notations for the perturbed system. As for the unperturbed
system, we define the compact level sets f
ME and the open set e
U (replacing
H by e
H in the definition). Without loss of generality, we can assume that e
U ae K.
Furthermore,
e
+\Deltas). The corresponding sets e
V and e
D are now defined in the obvious way.
Finally, the global Poincar'e map e
\Pi and the reduced Poincar'e map e
\Pi 0 are introduced
as for the unperturbed system. Again, without loss of the generality, we can assume
that D ae e
.
We extend the discrete time map \Psi ffit to a map \Psi t , t 2 [0; ffit], by using the exact
flow map \Phi t; e
H of the modified problem as an interpolation for t 2 [0; ffit). The map is
then extended to t  ffit in the obvious way 8 as the composition of k steps with \Psi ffit
and one step with \Phi dt; e
Thus, in correspondence
with the definition of the global Poincar'e map
e
we define
for all x 2 e
. Here the Poincar'e return times e t p (x), x 2 e
V , are the same as in the
definition of e
\Pi. Lemma 5.2 implies that there is a constant e
sup
It follows from backward and forward error analysis [18] that there is a constant
for all x 2 e
V and for all ffit sufficiently small. More importantly, let fx i g i=1;::: ;N be a
"numerically" computed sequence of points with x
the corresponding sequence under the map e
\Pi with x
D. Each sequence fx i g,
respectively, generates two sequences of real numbers fE i g and fs i g, f e
fes i g respectively, which are defined by E
8 The resulting map is discontinuous at multiples of the step-size ffit. A smooth interpolation
could be defined. But this is not needed in the context of this paper.
and e s obviously have e
while The "drift"
in the values of E i and s i per step away from zero is exponentially small and sums up
linearly with the number of steps. This energy conserving property of a symplectic
method has been discussed by Benettin & Giorgilli [5] and Hairer & Lubich
[18]. The same exponentially slow drift follows for the sequence fs i g from
the Lipschitz constant of / on e
.
In other words, if we start initially on e
D, then the points computed "numerically"
with the Poincar'e map b
\Pi will stay in an exponentially small neighborhood of e
D over
exponentially many iterates of b
\Pi. Now, since our numerical method is of order p  1,
the compact manifolds f
and ME are O(ffit p ) away from each other, i.e., the
modified Hamiltonian e
H is an ffl-perturbation of the Hamiltonian H with ffl  ffit p .
Thus the sequence fx i g will also stay in a O(ffit p ) neighborhood of D as long as the
number of iterates N satisfies
an appropriate constant.
Now the Shadowing Lemma [33] is applied to the sequence fx i g i=1;::: ;N .
Proposition 5.3. There exists an exact trajectory f"x i g i=1;::: ;N of the Poincar'e
map \Pi 0 on D such that the "numerically" computed sequence fx i g i=1;::: ;N stays in
a O(ffit p ) neighborhood of the (shadowing) exact trajectory provided the number of
iterates N satisfies (5.4).
Proof. We first project the sequence fx i g i=1;::: ;N down onto D using a "search
direction" orthogonal to the manifold D. Denote the result by fx i g. The projected
sequence fx i g and the sequence fx i g are O(ffit p ) close to each other provided N
satisfies (5.4). The "local" error per step between the "exact" Poincar'e map \Pi and
the "numerical" Poincar'e map b
\Pi is also of order p in the step-size ffit. This follows
from standard forward error analysis. Thus the Shadowing Lemma [33] for uniformly
hyperbolic maps can be applied to the Poincar'e map and the projected
sequence fx i g on D. The shadowing distance is O(ffit p ). This shadowing result also
applies to the sequence fx i g.
Let us now assume that we want to compute the ensemble average of a smooth
R up to a certain accuracy c ? 0. The large deviation theorem
for hyperbolic maps tells us that the probability to obtain the ensemble average
in the desired accuracy as the time average along a single trajectory goes to one
exponentially fast as the length N of the trajectory is increased. If we numerically
compute an approximative trajectory for the system (5.1)-(5.2), then we know from
Proposition 5.3 that this trajectory is O(ffit p ) close to some exact trajectory over
exponentially many integration steps N . Let us denote the time average of A along
this exact trajectory by hAi e
N and the numerically computed time average by hAiN ,
then
for all N satisfying a bound of type (4.4). Thus we obtain the following:
Proposition 5.4. Let (5.1)-(5.2) be a Poincar'e hyperbolic (real-analytic) system
which we discretize by a symplectic method of order p  1 in the step-size ffit. Then
the time-average hAiN of an observable A along a "numerically" computed trajectory
fxn gn=1;::: ;N ,
satisfies (5.5) where hAi e
N is the time-average along some exact trajectory and the
number of steps N satisfies a bound of type (5.4). Furthermore, assume we want to
compute the ensemble average of A within a given accuracy c ? 0. We assume, for
simplicity, that the constant c is larger than the difference between the time averages
(5.5) which is always true for sufficiently small step-sizes ffit. Then the probability to
obtain the average in the desired accuracy as the time average along a numerically
computed trajectory goes to one exponentially fast as the number of integration steps
N is increased. Taking the maximum number (5.4) of steps, the probability can be
made double exponentially close to one in (5.3) as ffit ! 0.
Example. As a numerical example, we look at the following planar anisotropic
Kepler problem [15]:
d
dt
d
dt
and
The initial conditions are chosen such that
and Note that angular momentum L is not conserved.
We define the Poincar'e section S 0 by record the sequence of points
Conservation of energy implies that the thus defined sequence is restricted
to the subset
This subset has an awkward shape. But it can be transformed
into a rectangle by means of the area preserving transformation
10=2. The corresponding Poincar'e map is hyperbolic,
i.e, stable and unstable manifolds intersect transversally, and the dynamics can be
encoded in a binary Bernoulli shift [15].
The main computational difficulty consists in the existence of a weak singularity
at To remove this singularity, we have to scale the equations of motion by
introducing the time transformation
which implies that, in the new time  , the norm of the vector field remains bounded
at time transformation has to be introduced such that the transformed
equations of motion are still Hamiltonian. A constant step-size symplectic method
can then be used to integrate the transformed system. Let us describe the general
Assume that a Hamiltonian function H(q;
scaling function ae(q) are given. Following Zare & Szebehely [41], we introduce the
modified Hamiltonian function
with corresponding Hamiltonian equations of motion
d
d
d
d
d
d
d
d
in extended phase space R 2n \Theta R 2 . In particular, let us consider the case
p(0)). Then (5.9)-(5.11) can be simplified to
d
d
d
d
d
d
on the hypersurface of constant energy This is a scaled vector field as desired
but which is not Hamiltonian anymore. Therefore, as suggested by the author in
[27] and independently by Hairer [17], the Hamiltonian equations (5.9)-(5.11) are
discretized by a symplectic method and For example, the equations
can be discretized by the symplectic Euler method, i.e.
The method is explicit in the variable q. Unfortunately this implies that the method
is only first order in ffi . However, the method is symplectic and, therefore, the
time t
time
average
time
time
average
time
average
time t
time
average
(a)
time
actual
step-size
-22time t
angular
momentum
time
total
energy
(b)
Fig. 5.1. (a) The time evolution of the average hrin (mean distance) is shown for four different
initial conditions with equal initial energy evolution of the error in energy,
angular momentum, and the actual step-size. The bottom-right figure shows the intersections of the
trajectory with the Poincar'e section in the (X1 ; X2 ) coordinates. One thousand intersections are
plotted.
Hamiltonian is conserved to O(ffi ) over exponentially long
periods of time. A second-order symplectic discretization can be obtained by using
the second-order Lobatto IIIa-b partitioned Runge-Kutta formula [34], i.e.
The resulting scheme is implicit in ae(q).
This approach is applied to the anisotropic Kepler problem with a scaling function
We chose initial values such that The
equations of motion are integrated using the second order symplectic method (5.12)-
0:05. The time average of an observable A(q) along a trajectory
computed according to the recursive formula
hAin =t n
The time average of
was computed for four
different initial conditions and the evolution of the corresponding time averages hri n
can be found in Fig. 5.1(a). The different lengths of the time intervals are due to
the fact that the same number of steps with step-size ffi were taken which leads to
different actual step-sizes ffit n . Within a tolerance of c = 0:04, these averages converge
to the same value  1:33. In Fig. 5.1(b), the total energy H , the angular momentum
L, and the variation in the actual step-size can be found for a particular
trajectory. We also plotted the intersections of the trajectory with the (q x
in the (X Theoretically, these points should fill the
rectangle in a uniform way (the invariant meassure is the Lebesgue meassure). With
one thousand points plotted, the uniform distribution is satisfied quite well. \Pi

Acknowledgements

. I like to thank Ernst Hairer, Chus Sanz-Serna, Andrew
Stuart, and Claudia Wulff for comments on an earlier version of this paper and the
referees for making valuable suggestions.



--R

The Lie group structure of diffeomorphism groups and invertible Fourier integrals operators with Applications
Geometrical Methods in the Theory of Ordinary Differential Equations
Mathematical Methods of Classical Mechanics

On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms
Numerical methods for dynamical systems
Modified Equations for ODEs

Aspects of backward error analysis of numerical ODEs
Formal power series and numerical algorithms for dynamical systems
Discretization of homoclinic orbits and invisible chaos
On the scope of the modified equations
Qualitative properties of modified equations
Dynamical Systems and Bifurcations of Vector Fields
Chaos in Classical and Quantum Mechanics
Backward analysis of numerical integrators and symplectic methods
Variable time step integration with symplectic methods
The life-span of backward error analysis for numerical integrat- ors
Solving Ordinary Differential Equations
Reversible long-term integration with variable step sizes
The inverse function theorem of Nash and Moser


Lectures on Hamiltonian systems
The separation of motions in systems with rapidly rotating phase
Numerical integration of generalized Euler equations
Backward error analysis for numerical integrators
On higher order semi-explicit symplectic partitioned Runge-Kutta methods for constrained Hamiltonian systems
Preservation of adiabatic invariants under symplectic discretization
Dynamical Systems
Symplectic integrators for Hamiltonian problems: an overview

Rigorous verification of trajectories for the computer simulation of dynamical systems
Symplectic partitioned Runge-Kutta methods

Stochastic Dynamics of Deterministic Systems
Introduction to Ergodicity Theory
The modified equation approach to the stability and accuracy of finite-difference methods
Large deviations in dynamical systems
Construction of higher order symplectic integrators
Time transformations for the extended phase space
--TR

--CTR
P. F. Tupper, Computing statistics for Hamiltonian systems: A case study, Journal of Computational and Applied Mathematics, v.205 n.2, p.826-834, August, 2007
Ernst Hairer, Important Aspects of Geometric Numerical Integration, Journal of Scientific Computing, v.25 n.1, p.67-81, October   2005
Y.-K. Zou , W.-J. Beyn, On manifolds of connecting orbits in discretizations of dynamical systems, Nonlinear Analysis: Theory, Methods & Applications, v.52 n.5, p.1499-1520, February
B. Cano , A. Durn, Analysis of variable-stepsize linear multistep methods with special emphasis on symmetric ones, Mathematics of Computation, v.72 n.244, p.1769-1801, October
M. P. Calvo , A. Portillo, Are high order variable step equistage initializers better than standard starting algorithms?, Journal of Computational and Applied Mathematics, v.169 n.2, p.333-344, 15 August 2004
Brian E. Moore , Sebastian Reich, Multi-symplectic integration methods for Hamiltonian PDEs, Future Generation Computer Systems, v.19 n.3, p.395-402, April
Fong , Eric Darve , Adrian Lew, Stability of Asynchronous Variational Integrators, Proceedings of the 21st International Workshop on Principles of Advanced and Distributed Simulation, p.38-44, June 12-15, 2007
Jess A. Izaguirre , Scott S. Hampton, Shadow hybrid Monte Carlo: an efficient propagator in phase space of macromolecules, Journal of Computational Physics, v.200 n.2, p.581-604, November 2004
B. Cano , A. Durn, A technique to construct symmetric variable-stepsize linear multistep methods for second-order systems, Mathematics of Computation, v.72 n.244, p.1803-1816, October
hamiltonian systems;error analysis;long time dynamics;differential equations;numerical integrators