Article
Reference
Long-time energy conservation of numerical methods for oscillatory
differential equations
HAIRER, Ernst, LUBICH, Christian
Abstract
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 determines its coefficients. 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 %of order 2 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.
HAIRER, Ernst, LUBICH, Christian. Long-time energy conservation of numerical methods for
oscillatory differential equations. SIAM journal on numerical analysis, 2000, vol. 38, no. 2, p.
414-441
Available at:
http://archive-ouverte.unige.ch/unige:12329
Disclaimer: layout of this document may differ from the published version.
1 / 1
LONG-TIME ENERGY CONSERVATION OF NUMERICAL
METHODS FOR OSCILLATORY DIFFERENTIAL EQUATIONS
∗
ERNST HAIRER
†
AND CHRISTIAN LUBICH
‡
SIAM J. N
UMER. ANAL.
c
2000 Society for Industrial and Applied Mathematics
Vol. 38, No. 2, pp. 414–441
Abstract. 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.
Key words. oscillatory differential equations, long-time energy conservation, second-order sym-
metric methods, frequency expansion, backward error analysis, Fermi–Pasta–Ulam problem
AMS subject classifications. 65L05, 65P10
PII. S0036142999353594
1. Introduction. Long-time near-conservation of the total energy and of adia-
batic invariants in numerical solutions to Hamiltonian differential equations is impor-
tant 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 adia-
batic 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 symplec-
tic 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
hω → 0 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 +Ω
2
x = g(x),(1.1)
∗
Received by the editors March 17, 1999; accepted for publication (in revised form) January 13,
2000; published electronically July 19, 2000.
http://www.siam.org/journals/sinum/38-2/35359.html
†
Dept. de math´ematiques, Universit´e de Gen`eve, CH–1211 Gen`eve 24, Switzerland (Ernst.Hairer
@math.unige.ch).
‡
Mathematisches Institut, Universit¨at T¨ubingen, Auf der Morgenstelle 10, D–72076 T¨ubingen,
Germany (lubich@na.uni-tuebingen.de).
414
LONG-TIME ENERGY CONSERVATION 415
where
Ω=
00
0 ωI
,ω 1(1.2)
(with blocks of arbitrary dimensions), and where the nonlinearity
g(x)=−∇U(x)(1.3)
has a Lipschitz constant bounded independently of ω. This situation arises in the cel-
ebrated 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 fre-
quencies are confined to the linear part of the problem, and that the linear part has
a single high frequency. The diagonal form of Ω 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 ω = 0. 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),
x(t)=y(t)+
k=0
e
ikωt
z
k
(t),(1.4)
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
H(x, ˙x)=
1
2
(|˙x
1
|
2
+ |˙x
2
|
2
)+
1
2
ω
2
|x
2
|
2
+ U(x),(1.5)
where x =(x
1
,x
2
) according to the partitioning of Ω. The other invariant is close to
I(x, ˙x)=
1
2
|˙x
2
|
2
+
1
2
ω
2
|x
2
|
2
,(1.6)
which represents the oscillatory energy of the system.
For the numerical solution we derive a similar frequency expansion which is valid
on grid points t = nh. Under a nonresonance assumption on hω,
|sin(
1
2
khω)|≥c
√
h for k =1,...,N (N ≥ 2),(1.7)
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 C
N
h
−N
. The only restriction on N comes from the
above nonresonance condition. The analysis uses only the symmetry of the methods
and does not require symplecticness.
416 ERNST HAIRER AND CHRISTIAN LUBICH
In section 2 we describe the numerical methods and we present numerical ex-
periments with the Fermi–Pasta–Ulam problem. These experiments illustrate the
long-time conservation of the total and the oscillatory energy in nonresonance sit-
uations, 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.
Section 3 gives a complete analysis of the two-dimensional linear case of (1.1)
over the whole range of nonresonant, near-resonant, and exactly resonant cases. This
already gives 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 section
4, that of the numerical solution in section 5. The numerical invariants are derived
in section 6. The main result on the numerical long-time conservation of energy for
(1.1) is formulated and proved in section 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), where Ω
2
is a symmetric and positive semidefinite (not necessarily diagonal) real matrix, and we
assume that initial values x
0
and ˙x
0
are given at t
0
= 0. By the variation-of-constants
formula, the exact solution of (1.1) satisfies
x(t)
˙x(t)
=
cos tΩΩ
−1
sin tΩ
−Ω sin tΩ cos tΩ
x
0
˙x
0
+
t
0
Ω
−1
sin(t − s)Ω
cos(t − s)Ω
g
x(s)
ds.(2.1)
(Observe that Ω
−1
sin tΩ is well-defined also for singular Ω.) It is therefore natural to
consider, for a fixed step size h, the explicit discretization
x
n+1
= cos hΩ x
n
+Ω
−1
sin hΩ˙x
n
+
1
2
h
2
Ψ g
n
,(2.2)
˙x
n+1
= −Ω sin hΩ x
n
+ cos hΩ˙x
n
+
1
2
h
Ψ
0
g
n
+Ψ
1
g
n+1
,(2.3)
where g
n
= g(Φx
n
) and Φ = φ(hΩ), Ψ = ψ(hΩ), Ψ
0
= ψ
0
(hΩ), Ψ
1
= ψ
1
(hΩ) with
real functions φ(ξ), ψ(ξ), ψ
0
(ξ), ψ
1
(ξ) depending smoothly on ξ
2
.Forg(x) ≡ 0 this
method integrates the problem (1.1) without error.
If φ(0) = ψ(0) = ψ
0
(0) = ψ
1
(0) = 1, the method is consistent of order 2. For
fixed Ω and for h → 0, second-order convergence follows from classical results. In this
article we are mainly interested in the situation where hΩ can take large values.
For long-time integrations, symmetric and/or symplectic methods are expected
to have favorable properties. By exchanging n ↔ n + 1 and h ↔−h in (2.2)-(2.3), it
is seen that the method is symmetric for all g(x) if and only if
ψ(ξ) = sinc ξ ·ψ
1
(ξ),ψ
0
(ξ) = cos ξ ·ψ
1
(ξ)(2.4)
(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), φ(ξ)=ψ
1
(ξ) also holds. This
condition will not be required for the analysis of our paper.
LONG-TIME ENERGY CONSERVATION 417
stiff
harmonic
soft
nonlinear
x
1
x
2
x
2n−1
x
2n
···
Fig. 1. Alternating soft and stiff springs.
Since the nonlinearity g(x) in (1.1) does not depend on ˙x, we can eliminate ˙x
n
in
(2.2) with the help of (2.3). In the case of a symmetric discretization we thus get the
two-step recurrence
x
n+1
− 2 cos hΩ x
n
+ x
n−1
= h
2
Ψ g
n
.(2.5)
The starting value x
1
is obtained from (2.2) with n = 0. For the case Ω = 0 we
recognize the well-known St¨ormer method.
Methods of the type (2.5) or (2.2)–(2.3) have been proposed and studied by several
authors. Gautschi [Gau61] suggests taking ψ(ξ) = sinc
2
(ξ/2). With this choice, (1.1)
is integrated exactly even for g(x) = const. Deuflhard [Deu79] discretizes the integral
in (2.1) by the trapezoidal rule and thus arrives at (2.5) with ψ(ξ) = sinc ξ and
φ(ξ) = 1. More recently, Garc´ıa-Archilla, Sanz-Serna, and Skeel [GSS99] introduce
a function φ(ξ) in the argument of g, and they consider the case where the method
is symplectic, so that ψ(ξ) = sinc ξ · φ(ξ). Hochbruck and Lubich [HoL99] consider
ψ(ξ) = sinc
2
(ξ/2) and general φ(ξ). 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 φ, ψ with φ(0)=1,ψ(0) = 1, which
have no zeros except possibly at integral multiples of π. Since we are interested in
the energy conservation of the numerical solution, we also need 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
x
n+1
− x
n−1
=2h sinc hΩ˙x
n
.(2.6)
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 ˙x
n−1
by exchanging n ↔ n + 1 and h ↔−h in (2.3). Subtracting the
resulting formula from (2.3), we obtain the two-step recurrence
˙x
n+1
− ˙x
n−1
= −2Ω sin hΩ x
n
+
1
2
h
Ψ
1
g
n+1
+2Ψ
0
g
n
+Ψ
1
g
n−1
.(2.7)
Formulas (2.5) and (2.7) give a symmetric two-step method even if (2.4) is not satis-
fied. If ψ(ξ) = sinc
2
(ξ/2) and ψ
0
(ξ)+ψ
1
(ξ) = 2 sinc ξ, then this method is exact for
g(x) = const. The choice ψ
1
(ξ) = 0 has 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 Figure 1). The variables x
1
,...,x
2n
(and x
0
=0,x
2n+1
=0)
stand for the displacements of end-points of the springs. The movement is described
by a Hamiltonian system with
H(x, ˙x)=
1
2
n
i=1
˙x
2
2i−1
+˙x
2
2i
+
K
2
n
i=1
(x
2i
− x
2i−1
)
2
+
n
i=0
(x
2i+1
− x
2i
)
4
.
