Long-Time Energy Conservation of Numerical Methods for Oscillatory Differential Equations.
We consider second-order differential systems where high-frequency oscillations are generated by a linear part. We present a frequency expansion of the solution, and we discuss two invariants of the system that determine the coefficients of the frequency expansion. These invariants are related to the total energy and the oscillatory harmonic energy of the original system.For the numerical solution we study a class of symmetric methods that discretize the linear part without error. We are interested in the case where the product of the step size with the highest frequency can be large. In the sense of backward error analysis we represent the numerical solution by a frequency expansion where the coefficients are the solution of a modified system. This allows us to prove the near-conservation of the total and the oscillatory energy over very long time intervals.
Introduction
. Long-time near-conservation of the total energy and of adiabatic
invariants in numerical solutions to Hamiltonian differential equations is important
in a wide range of physical applications from molecular dynamics to nonlinear
wave propagation. Backward error analysis [BG94, HaL97, Rei98] has shown that
symplectic numerical integrators approximately conserve the total energy and adiabatic
invariants over times that are exponentially long in the step size; more precisely,
over times of length exp(c=h!) where ! is the highest frequency in the system. Such
a result is meaningful only for h! ! 0, which is often not a practical assumption.
For example, in spatially discretized wave equations h! is the CFL number, which is
not chosen small in actual computations. Recently, in [GSS99, HoL99] new symplectic
or symmetric time-stepping methods have been studied which admit second-order
error bounds on finite time intervals independently of the frequencies of the dominant
linear part of the system. In particular for such "long-time-step methods", the case
is of no computational interest. The situation is reminiscent of stiff versus
nonstiff differential equations, where stiff integrators are not appropriately analyzed
by considering only the limit behavior h ! 0. In the stiff case, much insight has been
gained by studying the behavior of numerical methods on well-chosen, rather simple
linear and nonlinear stiff model problems.
As a first step towards an understanding of the numerical energy behavior in
Hamiltonian systems when the product of the step size and the highest frequency is
not a small quantity, we consider in this article the nonlinear, highly oscillatory model
problem
x
where
Dept. de math'ematiques, Universit'e de Gen'eve, CH-1211 Gen'eve 24, Switzerland
y Mathematisches Institut, Universit?t T-ubingen, Auf der Morgenstelle 10, D-72076 T-ubingen,
Germany (E-mail: lubich@na.uni-tuebingen.de)
(with blocks of arbitrary dimensions), and where the nonlinearity
has a Lipschitz constant bounded independently of !. This situation arises in the celebrated
Fermi-Pasta-Ulam model, for which we will present numerical experiments.
Clearly, in the model problem (1.1) we take strong restrictions in that the high frequencies
are confined to the linear part of the problem, and that the linear part has
a single high frequency. The diagonal form
of\Omega is not essential, since the numerical
methods are invariant under a diagonalization of the matrix.
We study the long-time energy behavior of a class of symmetric numerical methods
which are used with step sizes h such that the product h! is bounded away from zero
and can be arbitrarily large. The methods integrate the linear part of (1.1) exactly
and reduce to the St-ormer/Verlet method for The class includes the methods
of [GSS99, HoL99]. Classical symmetric methods such as the St-ormer/Verlet method,
the trapezoidal rule, or Numerov's method are not considered in this article. However,
using the results of the present paper, their energy behavior on (1.1) for h! in the
range of linear stability is analyzed in [HaL99].
Our approach to the near-conservation of the energy is based on a frequency
expansion of the solution x(t) of (1.1),
which is an asymptotic series where the coefficient functions y(t) and z k (t) together
with their derivatives are bounded independently of !. It turns out that the system
determining the coefficient functions has two (formally exact) invariants. One of these
is close to the total energy
according to the partitioning of \Omega\Gamma The other invariant is close to
which represents the oscillatory energy of the system.
For the numerical solution we derive a similar frequency expansion which is valid
on grid points Under a non-resonance assumption on h!,
the equations determining the coefficient functions have a similar structure to those
of the continuous problem. This allows us to obtain two almost-invariants close to H
and I , and rigorous estimates for the near-conservation of the total and the oscillatory
energy over time intervals of size CN h \GammaN . The only restriction on N comes from the
above non-resonance condition. The analysis uses only the symmetry of the methods
and does not require symplecticness.
In Sect. 2 we describe the numerical methods and we present numerical experiments
with the Fermi-Pasta-Ulam problem. These experiments illustrate the long-time
conservation of the total and the oscillatory energy in non-resonance situations,
which will later be completely explained by the theory. We also show the energy behavior
of the methods near resonances. This behavior depends strongly on properties
of the filter functions that determine the numerical method. We identify conditions
that yield satisfactory energy conservation near resonances and the correct energy
exchange between highly oscillatory components.
Sect. 3 gives a complete analysis of the two-dimensional linear case of (1.1) over
the whole range of non-resonant, near-resonant and exactly resonant cases. This
gives already much insight into conditions determining the energy conservation in the
general situation.
The frequency expansion of the analytical solution of (1.1) is introduced in Sect. 4,
that of the numerical solution in Sect. 5. The numerical invariants are derived in
Sect. 6. The main result on the numerical long-time conservation of energy for (1.1)
is formulated and proved in Sect. 7.
2. Numerical methods and numerical experiments. In this section we
present the numerical methods and we illustrate the main results of this paper with
the Fermi-Pasta-Ulam problem.
2.1. The discretization. We consider the differential equation (1.1),
is a symmetric and positive semi-definite (not necessarily diagonal) real matrix, and
we assume that initial values x 0 and -
x 0 are given at t By the variation-of-
constants formula, the exact solution of (1.1) satisfies
t\Omega \Gamma\Omega sin
t\Omega cos
t\Omega
s)\Omega
ds
(observe
t\Omega is well-defined also for
singular\Omega\Gamman It is therefore natural to
consider, for a fixed step size h, the explicit discretization
h\Omega xn
\Gamma\Omega sin
real functions OE(-), depending smoothly on - 2 . For
method integrates the problem (1.1) without error.
If the method is consistent of order 2. For
fixed\Omega and for h ! 0, second-order convergence follows from classical results. In this
article we are mainly interested in the situation where
h\Omega can take large values.
For long-time integrations, symmetric and/or symplectic methods are expected
to have favorable properties. By exchanging n in (2.2)-(2.3), it
is seen that the method is symmetric for all g(x) if and only if
(where sinc sin -). It can be shown by direct verification that (2.2)-(2.3) is
a symplectic discretization if, in addition to (2.4), also holds. This
condition will not be required for the analysis of our paper.
Since the nonlinearity g(x) in (1.1) does not depend on -
x, we can eliminate -
xn in
(2.2) with the help of (2.3). In the case of a symmetric discretization we thus get the
two-step recurrence
The starting value x 1 is obtained from (2.2) with For the
recognize the well-known St-ormer method.
stiff
harmonic
soft
nonlinear
Fig. 1. Alternating soft and stiff springs
Methods of the type (2.5) or (2.2)-(2.3) have been proposed and studied by several
authors. Gautschi [Gau61] suggests to take (-=2). With this choice, (1.1)
is integrated exactly even for const . Deuflhard [Deu79] discretizes the integral
in (2.1) by the trapezoidal rule and thus arrives at (2.5) with
recently, Garc'ia-Archilla, Sanz-Serna and Skeel [GSS99] introduce
a function OE(-) in the argument of g, and they consider the case where the method
is symplectic, so that Hochbruck and Lubich [HoL99] consider
OE(-). The papers [GSS99] and [HoL99] derive error
bounds on finite time intervals which are independent of ! and of the smoothness of
the solution.
In this article, we consider general functions OE; / with
have no zeros except possibly at integral multiples of -. Since we are interested in
the energy conservation of the numerical solution, we need also an approximation to
the derivative if we use the two-term recurrence relation (2.5). This can be obtained
by the relation (2.3) or, in the case of a symmetric method, also by the formula
This is possible if h! is not a nonzero integral multiple of -. We obtain (2.6) by
subtracting (2.5) from twice the formula of (2.2). For a symmetric method we obtain
a formula for -
by exchanging n in (2.3). Subtracting the
resulting formula from (2.3), we obtain the two-step recurrence
\Gamma2\Omega sin
Formulas (2.5) and (2.7) give a symmetric two-step method even if (2.4) is not satis-
fied. If -, then this method is exact for
const . The choice / 1 been considered in [HoL99].
2.2. Experiments with the Fermi-Pasta-Ulam problem. We consider a
chain of springs, where soft nonlinear springs alternate with stiff harmonic springs
(see [GGMV92] and Fig. 1). The variables x
for the displacements of end-points of the springs. The movement is described by a
Hamiltonian system with
Using the symplectic change of variables u
2,
we get a new Hamiltonian system with
(a) H
I
I
I
(a) H
I
I
I
(a) H
I
I
(c) H
I
100 200 3001
(a) H
I
100 200 300
I
100 200 300
(c) H
I
Fig. 2. Energy exchange of stiff springs
where This is exactly of the form (1.1).
For our numerical experiments we consider the case shown in Fig. 1)
As initial values we take
and zero for the remaining initial values. We apply the method (2.2)-(2.3) with the
following data:
(a)
(c)
with / 0 (-) and / 1 (-) given by (2.4).
We study the total energy (2.8) and the oscillatory energy
along the numerical solution on the interval 0 - t - 400. With the chosen initial
values we have In Fig. 2 we have
plotted, for three different step sizes and for all three methods, the numerical values
for I 1 ; I 2 ; I 3 ; I and H \Gamma 0:8. We see that an exchange of energy takes place, going
from the first stiff spring with energy I 1 to the second stiff spring and later to the
third one. For the smallest step size we have also plotted in gray the numerical values
for perturbed initial values obtained by adding 10 \Gamma8 to u 1 (0), -
illustrates that the solution is very sensitive to perturbations.
6 E. HAIRER AND CH. LUBICH
In all cases we see that H and I are well preserved over the whole interval, even
for step sizes where the numerical solution is completely wrong. Further experiments
have shown that such a preservation holds for much longer intervals (we tested up to
An explanation of this phenomenon is the main objective of this paper.
2.3. Numerical experiments in near-resonant situations. When the product
of the step size and the frequency h! is close to a multiple of -, then the different
methods show widely different behavior. Energy is conserved only for some choices of
/ and OE. Satisfactory numerical behavior also in near-resonance situations is obtained
if the numerical method satisfies the following additional conditions:
These conditions yield long-time energy conservation for all values of h! with the
exception of h! in intervals of width O(h) near integral multiples of 2-. The total
energy appears to be conserved uniformly for arbitrary values of h! if
The necessity of these conditions is seen from an analysis of the linear case, which is
given in Sect. 3.
When h! is close to 2m- with a positive integer m, the condition (2.9) requires
a double zero of / at 2m-. Similarly, for h! close to an odd multiple of -, condition
there requires a simple zero of /. For the choice
HoL99] the condition (2.9) is obviously satisfied for all values of h!, but (2.10) is
violated near odd multiples of -. For condition (2.10) is trivially
satisfied for all h!, but condition (2.9) fails near even multiples of -. The choice
satisfies the three conditions (2.9)-
for all h!. Condition (2.12) is not satisfied by any of the methods previously
proposed in the literature.

Table
Methods used for the numerical experiments of Sect. 2.3
A sinc (-) 1 2k- p
Let us illustrate the effect of the conditions (2.9), (2.10) and (2.11) on the numerical
solution when h! is close to a multiple of -. We consider the Fermi-Pasta-Ulam
problem of Sect. 2.2 with the same initial values, and we apply six different methods.
Their characteristics are given in Table 1. The sign p
indicates that the corresponding
condition on / and OE is satisfied. If a condition is not satisfied for all values of
-22 (C)
out of scale
-22 (D)
Fig. 3. Energy conservation of different methods for
Fig. 4. Energy conservation of different methods for
Fig. 5. Energy conservation of different methods for
h!, we give the values close to which it is violated. For each of the methods (B), (C),
(D) only one of the conditions (2.9)-(2.11) is not fulfilled.
In Fig. 3 we show the errors of the Hamiltonian over the interval [0; 1000]. We
have used the step size such that
a maximal error of size 396497, because / 1 (-) given by (2.4) has a singularity at -.
.2
.2
.2
I
.2
I
Fig. 6. Error in the total and oscillatory energies as a function of h! for the FPU problem
Methods (A) and (D) show a clear drift from the constant value of the Hamiltonian.
Only the methods (B), (E), and (F), for which all three conditions are satisfied close to
-, conserve the Hamiltonian very well. For the pictures corresponding to these three
methods, we have changed the scale so that the small oscillations become visible.
Fig. 4 shows the same experiment, where this time ! is chosen such that
2:0000001 \Delta -. For both situations, we get the same qualitative behavior when we plot
the oscillatory energy instead of the Hamiltonian. The results of these experiments
confirm that the conditions (2.9), (2.10), (2.11) cannot be omitted if we are interested
in long-time energy estimates that hold uniformly in h!.
In Fig. 5 the numerical results of the nonresonant case are included.
All methods give satisfactory results. The most accurate results are obtained by the
method (C).
In the upper pictures of Fig. 6 we plot the maximal errors in the Hamiltonian
as a function of h!, and we take step sizes 0:025. The
picture to the right corresponds to the method (E) of Table 1. The picture to the left
is obtained with method (F) which satisfies (2.12). Uniform convergence of the error
can be observed only in this case. The lower pictures of Fig. 6 show the analogue
for the deviations in the oscillatory energy. Close to integral multiples of 2- this
deviation is large for both methods. The same phenomenon can be observed already
for linear problems (see Fig. 8), for which a complete analysis is given in Sect. 3.
2.4. Energy exchange. The energy exchange between stiff components takes
place on time intervals of length O(!). In Fig. 2, this is reproduced qualitatively
correctly for large h! only in the case where 1. The numerical
frequency expansion of Sect. 5 shows that the condition
.2
HH
.2
HH
.2
II
.2
II
Fig. 7. Same experiment as in Fig. 6 with method (2.5),
with method (2.16) (right pictures)
is needed for the approximation of the energy exchange between stiff components
when h! is bounded away from zero (compare the equations for the z 2 component in
(4.8) and (5.8)). This is a severe condition which excludes all methods considered so
far with the exception of the above-mentioned method
the other hand, we have seen that this method has rather poor energy conservation
properties. We therefore extend the class of methods (2.5) to
(h\Omega\Gamma2 For consistency, the functions / k , OE k must satisfy
Conditions (2.9) and (2.10) are now needed for
condition (2.13) is
replaced with
For example, the method
with (2.6), or equivalently in one-step form
h\Omega xn
\Gamma\Omega sin
shows the correct energy exchange
for large h!, as for method (b) in Fig. 2. In Fig. 7, we compare the energy
conservation for the method (2.5) with and the method
for the same step sizes as in Fig. 6. For (2.16), the total and oscillatory
energies are well conserved over long times except for h! in intervals of length O(h)
around integral multiples of -.
For ease of presentation, the following analysis will be done for methods of class
(2.5), but the arguments extend in an obvious way to the class (2.14).
3. Long-time energy conservation for linear problems. We start our analysis
with the case where
with a two-dimensional symmetric matrix A
satisfying a 11 ? 0, so that
is linear. This gives already a lot of insight and illustrates the importance of the
conditions (2.9)-(2.12). In this situation, the differential equation (1.1) becomes
x
The total energy H given by (1.5) is an invariant of the system. In the following we
assume that H
\Delta is bounded uniformly in !. This requires x 2
3.1. Analytical solution. The exact solution of (3.1) is given by
O(! \Gamma2 )'
i are the eigenvalues
so that
a 11 +O(! \Gamma2
For given initial values x(0), -
Consequently, we have
be
This implies
that the quantity
remains O(! \Gamma1 )-close to the constant value I
\Delta for all times t.
3.2. Numerical solution. We search for functions b
v, such that xn :=
x(nh) satisfies the numerical scheme (2.5) with This implies
so that has to be an eigenvalue of
and v a corresponding eigenvector. In this section we use the short notation
and OE(h!). Since the off-diagonal elements of (3.3) are small, the eigenvalues -
are close to ff
respectively. The corresponding eigenvectors are
'fl/
\Gammafl OE'
the general solution can thus be written as
where the complex coefficients are computed from the initial values b
obtained from (2.2).
Inserting
To study the long-time near-conservation of H
\Delta and I
\Delta we distinguish
two cases.
Case I: well-separated eigenvalues. We assume that one of the conditions
a 11 or
is satisfied. This covers nearly all choices of h!. Only values in intervals of length
O(h) are excluded.
Theorem 3.1. Consider the numerical method (2.2)-(2.3) applied to (3.1) with
a Under one of the restrictions (3.5) or (3.6) on the
step size, and under the conditions (2.9), (2.10), (2.11) on the numerical method, the
energies H and I of (1.5) and (1.6) along the numerical solution satisfy
for all n - 0. The constants symbolized by O(\Delta) are independent of !, h and n.
Proof. The characteristic polynomial of the matrix (3.3) is
Its zeros can be computed explicitly. Each of the conditions (3.5)
or (3.6) together with (2.11) implies that
for both eigenvalues - i of (3.3). Hence, for sufficiently small h, - 1 and - 2 are real and
in the interval [\Gamma1; 1]. The angles - i , defined by - are therefore also real
and satisfy
a
We next estimate the coefficients a; a; b; b in (3.4). Since the - i and v i are real, the
coefficients a; b are the complex conjugates of a; b. In both situations, (3.5) and (3.6),
we have under the assumption (2.9), we further have
Consequently, the condition b
of using b more convenient to work with b x(h) \Gamma cos
cos
From the estimates
sinc (h!=2)), which all follow from (2.9) and (2.10), we then get
be
The second relation is a consequence of the fact that
(which follows from -
and (2.11)). The statement (3.8) is now an immediate consequence of (3.9) and of the
because the modulus of c = be i-2 t is independent
of t.
The near-conservation of the Hamiltonian can be seen similarly. If one of the
conditions (3.5) and (3.6) is satisfied, we get flOE sin(- 2 so that
This implies jbx 0
Const
all t, this
together with (3.8) proves the statement (3.7) for the total energy.
Case II: nearly collapsing eigenvalues. We now consider the complementary case4 h 2 a
so that the two eigenvalues of (3.3) are very close.
Theorem 3.2. Consider the numerical method (2.2)-(2.3) applied to (3.1) with
a Under the condition (2.12) on the numerical method,
viz.
for all n - 0 and uniformly in h! - c ? 0. The constants symbolized by O(\Delta) are
independent of !, h and n.
Proof. Under the condition (2.12) the numerical method satisfies (2.9), (2.10),
so that Theorem 3.1 is applicable. It therefore remains to consider the situation
where h! is restricted by (3.10).
The condition (3.10) implies 1
2a 11 for sufficiently small
h. We now assume that OE/ - 0, which is satisfied by the choice (2.12). This guarantees
that the eigenvalues of (3.3) are real, that -
OE/), and that
a 11 +O(h 3
For the special choice (2.12) we further have jflj - 1= p
O(!), so that the coefficients of (3.4) satisfy a = O(1) and
(3.
.2
.4
.2
.4
.2
.4
I
.2
.4
I
Fig. 8. Error in the total and oscillatory energies as a function of h! for a linear problem
holds (recall that
This is a consequence of the estimate
OE/) and of the fact that jflj - 1=
OE/.
The relations (3.12) and (3.13) thus yield
a 11
Together with (3.4) this implies
To prove that the t-dependent term is small, we need the relation (2.12) between OE
and /. Using also (3.12) and (3.13) we obtain
a
a
sin
This completes the proof of Theorem 3.2, because
The proof above shows that in general
I
In the situation of Theorem 3.1 we have so that the first term in
the right-hand expression of (3.14) becomes negligible, and long-time conservation
of I can be concluded. If a 12 6= 0 and cos(h!) -
14 E. HAIRER AND CH. LUBICH
OE/ and fl - 1=
OE/, so that fl/ -
/=OE and flOE -
OE=/. The
initial condition b
so that none of the terms in (3.14) can be neglected. If - is different from an
integral multiple of 2-, the expression (3.14) cannot remain close to a constant value.
This result is intuitively clear, because in the situation (3.10) the two frequencies are
indistinguishably close for the numerical method.
The upper pictures of Fig. 8 show the maximal error in the Hamiltonian H(xn ; -
in dependence of h! for the problem (3.1) with a
initial values
. The three curves correspond
to the step sizes 0:05. The picture to the left it obtained
with a method satisfying (2.12). Uniform convergence of the error can nicely be
observed. The picture to the right corresponds to the method (E) of Table 1. The
lower pictures of Fig. 8 plot the maximal deviation of the oscillatory energy I(xn ; -
as a function of h!. It confirms the analysis above, which shows that for h! satisfying
a 11 the oscillatory energy cannot be well conserved.
4. Frequency expansion of the analytical solution. The main tool of our
analysis for nonlinear problems is a decomposition of the solution x(t) of (1.1) into
a smooth part and into highly oscillatory terms with smoothly varying amplitudes.
This decomposition is valid over finite time intervals. We show the existence of two
almost-invariants for the coefficients of this decomposition, which are related to the
total energy and the oscillatory energy of the system. A repeated use of these almost-
invariants then allows us to prove the long-time near-conservation of the oscillatory
energy.
4.1. The frequency expansion. We assume the nonlinearity g in (1.1) analytic
on an open set D, and we consider solutions of (1.1) which satisfy
where K is a compact subset of D. We assume further that the initial values have
limited harmonic energy:2 k -
where E is independent of !.
Theorem 4.1. Under the assumptions (4.2) and (4.1) for 0 - t - T , the solution
x(t) of Eq. (1.1) has for arbitrary N - 2 an expansion of the form
e ik!t z k (t) +RN (t);
where the remainder term and its derivative are bounded by
The real functions and the complex functions z
are bounded,
together with all their derivatives, by
and we have z . They are unique up to terms of size O(! \GammaN \Gamma2 ). The constants
symbolized by the O-notation are independent of ! and t with
on E, N , T , and on the order of the derivative).
Proof. To determine the smooth functions
we put
e ik!t z k (t);
insert this function into (1.1), expand the nonlinearity around y(t) and compare the
coefficients of e ik!t . With the notation g (m) (y)z
the following system of differential equations:
z k
z k
Here the sums range over all m - 1 and all multi-indices
integers ff j satisfying which have a given sum
For large !, the dominating terms in these differential equations are given by
the left-most expressions. However, since the central terms involve higher derivatives,
we are confronted with singular perturbation problems. We are interested in smooth
functions that satisfy the system up to a defect of size O(! \GammaN ). In the spirit
of Euler's derivation of the Euler-Maclaurin summation formula (see e.g. [HaW96])
we remove the disturbing higher derivatives by using iteratively the differentiated
equations (4.7)-(4.9). This leads to a system
z k
are formal series in powers of ! \Gamma1 . Since we get formal algebraic relations
for , we can further eliminate these variables in the functions F
.
We finally obtain for y the algebraic relations
z k
z k
and a system of real second-order differential equations for y 1 and complex first-order
differential equations for z
At this point we can forget the above derivation and we can take it as a motivation for
the ansatz (4.10)-(4.11), which we truncate after the O(! \GammaN ) terms. Inserting this
ansatz and its first and second derivatives into (4.7)-(4.9) and comparing like powers
recurrence relations for the functions F k
jl . This shows that these
functions together with their derivatives are all bounded on compact sets.
We determine initial values for (4.11) such that the function e x(t) of (4.6) satisfies
e
x(0). Because of the special structure of the ansatz (4.10)-
(4.11), this gives a system
which, by the implicit function theorem, yields (locally) unique initial values y 1 (0),
(0). The assumption (4.2) implies that z 2 It further follows
from the boundedness of F 2l that z 2 looking closer at
the structure of the function G k
jl it can be seen that it contains at least k times the
factor z 2 . This implies the stated bounds for all other functions.
We still have to estimate the remainder RN x(t). For this we consider
the solution of (4.10)-(4.11) with initial values (4.12). By construction, these functions
satisfy the system (4.7)-(4.9) up to a defect of O(! \GammaN ). This gives a defect of size
O(! \GammaN ), when the function e x(t) of (4.6) is inserted into (1.1). Hence on a finite time
To obtain the
slightly sharper bounds (4.4), we apply the above proof with N replaced by N 2.
4.2. The Hamiltonian of the frequency expansion. Consider now the situation
so that (1.1) is a Hamiltonian system
x
with Hamiltonian
where U(x) is assumed to be analytic. Let v k
that by (4.7)-(4.9) these functions satisfy
Here, the sum is again over all m - 1 and all multi-indices
integers ff which have a given sum
and we write
Further we denote
From the above it follows that the vector (y; V ) satisfies the system
y
which, neglecting the O(! \GammaN ) terms, is Hamiltonian with
Theorem 4.2. Under the assumptions (4.2) and (4.1) for
The constants symbolized by O(\Delta) are independent of ! and t with
depend on E, N and T .
Proof. Multiplying (4.17) and (4.18) with -
y T and ( -
respectively, gives
dt
Integrating from 0 to t and using v
By the bounds of Theorem 4.1, we have for
On the other hand, we have from (4.14) and (4.3) that
Using
it follows from
Inserted into (4.22) and (4.23) this yields the statement (4.21).
4.3. Another almost-invariant. Besides the Hamiltonian H(y; -
V ), the
coefficients of the frequency expansion have another almost-invariant. It only depends
on the oscillating part and it is given by
This almost-invariant turns out to be close to the energy of the harmonic oscillator,
Theorem 4.3. Under the assumptions (4.2) and (4.1) for
The constants symbolized by O(\Delta) are independent of ! and t with
depend on E, N and T .
Proof. With the vector holds that U(y;
Differentiating the identity
with respect to t yields
z \Gammak
because
The proof of Theorem 4.3 is now very similar to that of Theorem 4.2. We multiply
the relation (4.18) with \Gammai!k(v \Gammak ) T instead of ( -
Summing up yields, with the
use of (4.28),
\Gammai!
The derivative of I(V; -
by (4.24), is
d
dt
(v \Gammak
In the sums
the terms with k and \Gammak cancel.
Hence, the statement (4.26) follows from (4.29) and (4.30).
Using -
from the bounds of
Theorem 4.1 that
On the other hand, using the arguments of the proof of Theorem 4.2, we have
This proves the second statement of the theorem.
Corollary 4.4. If x(t) 2 K for
The constants symbolized by O(\Delta) are independent of ! and t with
depend on E and N .
Proof. With a fixed T ? 0, let V j denote the vector of frequency expansion terms
that correspond to starting values (x(jT ); -
x(jT )). For
we have by (4.27)
Vn ('T
Vn ('T
Vn
We note that I(V j+1 (0); -
by the uniqueness
statement of Theorem 4.1, we have V j+1
we have the bound (4.26) of Theorem 4.3. The same argument
applies to I(Vn ('T ); -
Vn ('T
Vn (0)). This yields the result.
Remark. It is already known from the article [BGG87] that the oscillatory energy
x(t)) is nearly preserved over long times. The proofs in [BGG87] are completely
different. They use coordinate transforms from Hamiltonian perturbation theory and
show that I is nearly preserved over time intervals which grow exponentially with !.
By carefully tracing the N-dependence of the constants in the O(! \GammaN )-terms, it is
possible to obtain near-conservation of I over exponentially long time intervals also
within the present framework of frequency expansions.
5. Frequency expansion of the numerical solution. In this section we show
that the numerical solution (2.2), (2.3) for nonlinear problems (1.1) has a frequency
expansion similar to that of the analytical solution. Following the idea of backward
analysis and motivated by the results of Sect. 4 we look for a function
e ik!t z k (t)
(with smooth y(t) and z k (t) depending on h) 1 such that, up to a small defect,
We assume throughout this section that
and that the numerical solution \Phix n remains in a compact subset of the region where
g(x) is analytic, i.e.,
5.1. Functional calculus. For the computation of the functions y(t) and z k (t)
the following functional calculus is convenient. Let f be an entire complex function
bounded by jf(i)j - C e fljij . Then,
converges for every function x which is analytic in a disk of radius r ? flh around t.
We note that (hD) k
are two such entire functions, then
whenever both sides exist. In particular, we have
To avoid an overloaded notation with hats, we use the same letters y and z k as for the analytical
solution. We hope that this does not cause confusion.
We therefore introduce the operator
h\Omega
sin
which, for h ! 0, is an approximation to h 2 (D 2
We next study the application of such an operator to functions of the form e i!t z(t).
By Leibniz' rule of calculus we have (hD) k e i!t z(t). After a
short calculation this also yields
f(hD)e i!t
5.2. Modified equations for the coefficient functions of the frequency
expansion. With the operator L(hD) of (5.5) the condition (5.2) becomes
Inserting the ansatz (5.1), expanding the right-hand side of (5.7) into a Taylor series
around \Phiy(t), and comparing the coefficients of e ik!t yields for the functions y(t) and
z k (t)
Here, multi-index as in the proof of Theorem 4.1,
ff is an abbreviation for the m-tupel (\Phiz ff To get
smooth functions y(t) and z k (t) which solve (5.8) up to a small defect, we look at the
dominating terms in the Taylor expansions of L(hD) and L(hD ik!h). With the
abbreviations
2 kh!) we have
(ihD)
The situation is now more complicated than in (4.7)-(4.9) for the frequency expansion
of the analytical solution, because several of the coefficients in (5.9) may vanish due
to numerical resonance. We here confine the discussion to the non-resonant case. We
assume that h and ! \Gamma1 lie in a subregion of the small parameters
for which there exists a positive constant c such that
2:
LONG-TIME ENERGY CONSERVATION 21
The condition excludes that h! is o(
close to integral multiples of -. For given h
and !, the condition imposes a restriction on N . In the following, N is a fixed integer
such that (5.10) holds.
Theorem 5.1. Under the limited-energy condition (4.2), under the non-resonance
condition (5.10), under the conditions (5.3), (5.4), and under the conditions (2.9) and
(2.10) on the numerical method (2.2)-(2.4), the numerical solution is of the form
uniformly for where the functions
satisfy (5.8) up to a defect of O(h N+2 ) in their first components, and O(/(h!)h N+2 )
in their second components. Together with all their derivatives these functions are
bounded by
, and the constants symbolized by the O-notation
are independent of ! and h, but depend on E, N , and T .
Proof. Under assumption (5.10), the first non-vanishing coefficients in (5.9) are
the dominant ones, and the derivation of the defining relations for y and z k is the
same as for the analytical solution in Theorem 4.1. We insert (5.9) into (5.8) and we
eliminate recursively the higher derivatives. This motivates the following ansatz for
the computation of the functions y and z k :
z k
z k
where the functions depend smoothly on the variables y 1 , -
and on the bounded
parameters
/(h!). Inserting this ansatz and its derivatives into
(5.8) and comparing like powers of
h yields recurrence relations for the functions
jl . The functions g k
jl (for k - 1) contain at least k times the factor OE(h!)z 2 , and
f 2l contains at least once this factor. Since the series in (5.12) need not converge, we
truncate them after the (
We next determine the initial values y 1 (0), -
x(h) of (5.1) coincide with the starting values x 0 and x 1 of the numerical scheme (x 1 is
computed from x 0 and -
x 0 via the formula (2.2) with Using the non-resonance
22 E. HAIRER AND CH. LUBICH
assumption (5.10), the condition b
The formula for the first component of (2.2), x
with b
implies that
For the second component we have x
x
(2.2), and b x 2
\Delta , which after division by h sinc h! yields
O
The four equations (5.13), (5.14), (5.15) constitute a nonlinear system for the four
quantities y 1 (0), -
z
. By the implicit function
theorem and using the limited-energy assumption (4.2), we get a locally unique
solution for sufficiently small h, if the conditions (2.9) and (2.10) are satisfied.
The initial value for z 2 satisfies z 2 it follows from (2.10) that
by (5.12). This implies z 2
for . The other estimates (5.11) are directly obtained from (5.12). Conse-
quently, the values b x(nh) inserted into the numerical scheme (2.5) yield a defect of
size O(h N+2
Standard convergence estimates then show that on bounded time intervals
is of size O(t 2 h N ) in the first component and of size O(/(h!)t 2 h N ) in the second
component. This completes the proof of Theorem 5.1.
5.3. Frequency expansion of the derivative approximation. Under the
condition (5.10) we have h! 6= k- for integer k, so that the derivative approximation
xn is given by (2.6). We now define b x 0 (t) by the continuous analogue
Using condition (2.10), Theorem 5.1 implies that
on bounded time intervals. We next write the function b x 0 (t) as
e ik!t z 0k (t):
Inserting the relation (5.1) into \Gammai
which is equivalent
to (5.17), and comparing the coefficients of e ik!t we obtain
sinc (ihD) -
In particular, we get for z 1
2 that
z 01
cos !h
sin !h
Theorem 5.2. Under the assumptions of Theorem 5.1, the numerical solution
xn , given by (2.6), satisfies
uniformly for where the functions y
together with all their derivatives are bounded by
z 01
The constants symbolized by the O-notation are independent of ! and h, but depend
on E, N and T .
Proof. The estimates follow from (5.19) and from Theorem 5.1. For y 0
1 and
z 01
2 we use the formulas (5.12) to get the sharper result.
5.4. Energy along the numerical solution. In the Hamiltonian case
\GammarU (y), the total energy H and the oscillatory energy I are related to the frequency
expansion coefficients as follows.
Lemma 5.3. If the coefficients of the frequency expansions for xn and -
xn satisfy
(5.11) and (5.21), respectively, then
Proof. By definition (4.25) we have I(bx; b
e i!t z 1
e i!t z 1
the statement (5.23) follows from the fact that jv . The formula
(5.22) can be proved in the same way.
6. Almost-invariants of the numerical frequency expansion. In this section
we show that, in the Hamiltonian case the coefficients of the
frequency expansion of the numerical solution have invariants that can be obtained as
in Sect. 4. We denote
z k are the coefficients of the frequency expansion (5.1). Similar to (4.16) we consider
the function
where the sum is taken over all m - 1 and all multi-indices
non-vanishing integral components for which
It
then follows from Theorem 5.1 that the coefficients y and v k satisfy
r y
The factor OE(h!) in the defect of (6.3) is due to the presence of the factor OE(h!)z 2 in
the relations (5.12) defining the z-components.
The similarity of these relations to (4.17), (4.18) allows us to obtain invariants
that are the analogues of H and I of Sect. 4.
6.1. First invariant. As in Sect. 4.2 we multiply (6.2) and (6.3) by -
y T and
respectively, and we thus obtain
dt
Since we know bounds on z k and on its derivatives (Theorem 5.1), we switch to the
quantities z k and we get the equivalent relation
dt
We shall show that the left-hand side is the total derivative of an expression that only
depends on y, z k and its derivatives. Indeed, the term -
y T y (2l) can be written as
dt
Similarly, we get for z \Gammak that
Re -
z T
z
dt
z T
z
Re z T z
dt
z T z
z T
z
dt
z T
z
Im z T z
dt
z T z
z T
Hence, there exists a function b
which depends on the values at t of the
functions y and of their first N derivatives, such
that (6.4) reads
d
dt
This yields immediately the first statement of the following result.
Theorem 6.1. Under the assumptions of Theorem 5.1, the coefficient functions
y and
for
Proof. The formula for b
H 0 is obtained from the formulas (5.9) for L(hD
together with the estimates of Theorem 5.1.
Remark. Symplectic discretizations have
6.2. Second invariant. As in the proof of Theorem 4.3 we have for the function
of (6.1) that
Consequently, it follows from (6.3) that
\Gammai!
Written in the z variables, this becomes
\Gammai!
As in Sect. 6.1, the left-hand expression can be written as the total derivative of a
function b I 0 [Z](t) which depends on the values at t of the function Z and its first N
derivatives:
d
dt
I
Theorem 6.2. Under the assumptions of Theorem 5.1, the coefficient functions
I
I
for
with -(h!) as in Theorem 6.1.
Proof. From (6.5) and the estimates of Theorem 5.1, we obtain
Because of condition (2.10), this yields the stated formula for b I 0 .
7. Long-time energy conservation of the numerical discretization. We
are now able to prove the main result of this paper. This shows that the total energy
H and the oscillatory energy I are nearly conserved over time intervals of length
CN h \GammaN , for any N for which the non-resonance condition (5.10) is satisfied.
For the convenience of the reader we restate our assumptions:
ffl the limited-energy condition
26 E. HAIRER AND CH. LUBICH
ffl the boundedness condition (5.4) for the numerical solution sequence: \Phix n stays in
a compact subset of the domain of analyticity of g;
ffl the condition (5.3): h! - d ? 0;
ffl the conditions (2.9) and (2.10) on the numerical method:
ffl the non-resonance condition (5.10): for some N - 2,
Theorem 7.1. Under the above conditions, the numerical solution of (1.1) obtained
by the method (2.2)-(2.4) satisfies
nh. The constants symbolized by O(\Delta) are independent of n and of h and !
satisfying the above conditions, but depend on N .
Proof. (a) If we consider the linear combinations b
I 0 and
I 0 =-, it follows from Theorem 6.1 and Theorem 6.2 that
Moreover, by Theorem 6.1 and Theorem 6.2 together with Lemma 5.3 we have
I
where again b x(t) is defined by the frequency expansion (5.1) with coefficients y(t) and
defined by (5.17). The relations (7.1) and (7.2) hold only on finite
time on which the frequency expansion is defined.
(b) We now apply the above relations repeatedly on intervals of length h, for
frequency expansions corresponding to different starting values. As long as
satisfies the limited-energy condition (4.2) (possibly with a larger constant E), Theorem
5.1 gives us frequency expansion coefficients yn (t); Zn (t) corresponding to starting
values
Because of the uniqueness (up to O(h N+1 )) of the coefficients of the
frequency expansion, the following diagram commutes up to terms of size O(h N+1
numerical
method
(up to O(h N+1
The construction of the coefficient functions via (5.12) shows that also higher derivatives
of (y n ; Zn ) at h and (y n+1 ; Zn+1 ) at 0 differ by only O(h N+1 ). We thus have
from (7.3) and (7.1)
Using this relation repeatedly, we obtain
Moreover, from (7.2) we have the following for the coefficient functions corresponding
to the starting values
xn ) and
by construction, and b x 0
Theorem 5.2, we
obtain
which gives the desired bound for the deviation of the total energy along the numerical
solution. The same argument applies to I(xn ; -

Acknowledgment

. We are grateful to Sebastian Reich for drawing our attention
to the Fermi-Pasta-Ulam problem.



--R

Realization of holonomic constraints and freezing of high frequency degrees of freedom in the light of classical perturbation theory.
On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms
A study of extrapolation methods based on multistep schemes without parasitic solutions
On the problem of energy equipartition for large systems of the Fermi-Pasta-Ulam type: analytical and numerical estimates

Numerical integration of ordinary differential equations based on trigonometric polynomials
The life-span of backward error analysis for numerical integrators
Energy conservation by St-ormer-type numerical integrators
Analysis by Its History
A Gautschi-type method for oscillatory second-order differential equations
Dynamical systems
--TR

--CTR
M. Van Daele , G. Vanden Berghe, Geometric numerical integration by means of exponentially-fitted methods, Applied Numerical Mathematics, v.57 n.4, p.415-435, April, 2007
Ernst Hairer, Important Aspects of Geometric Numerical Integration, Journal of Scientific Computing, v.25 n.1, p.67-81, October   2005
J. M. Franco, New methods for oscillatory systems based on ARKN methods, Applied Numerical Mathematics, v.56 n.8, p.1040-1053, August 2006
Paul J. Atzberger , Peter R. Kramer , Charles S. Peskin, A stochastic immersed boundary method for fluid-structure dynamics at microscopic length scales, Journal of Computational Physics, v.224 n.2, p.1255-1292, June, 2007
Fermi-Pasta-Ulam problem;frequency expansion;backward error analysis;second-order symmetric methods;oscillatory differential equations;long-time energy conservation