1 Introduction

Implicit-Explicit Variational Integration of Highly Oscillatory Problems


In this paper, we derive a variational integrator for certain highly oscillatory problems in mechanics. To do this, we take a new approach to the splitting of fast and slow potential forces: rather than splitting these forces at the level of the differential equations or the Hamiltonian, we split the two potentials with respect to the Lagrangian action integral. By using a different quadrature rule to approximate the contribution of each potential to the action, we arrive at a geometric integrator that is implicit in the fast force and explicit in the slow force. This can allow for significantly longer time steps to be taken (compared to standard explicit methods, such as Störmer/Verlet) at the cost of only a linear solve rather than a full nonlinear solve. We also analyze the stability of this method, in particular proving that it eliminates the linear resonance instabilities that can arise with explicit multiple-time-stepping methods. Next, we perform some numerical experiments, studying the behavior of this integrator for two test problems: a system of coupled linear oscillators, for which we compare against the resonance behavior of the r-RESPA method; and slow energy exchange in the Fermi–Pasta–Ulam problem, which couples fast linear oscillators with slow nonlinear oscillators. Finally, we prove that this integrator accurately preserves the slow energy exchange between the fast oscillatory components, which explains the numerical behavior observed for the Fermi–Pasta–Ulam problem.

2000 Mathematics Subject Classification:
65P10, 70K70
First author’s research partially supported by a Gordon and Betty Moore Foundation fellowship at Caltech, and by NSF grant CCF-0528101.
Second author’s research partially supported by the NSF (MSPA Award No. IIS-05-28402, CSR Award No. CNS-06-14770, CAREER Award No. CCF-06-43268).


1. Introduction

1.1. Problem Background

Many systems in Lagrangian mechanics have components acting on different time scales, posing a challenge for traditional numerical integrators. Examples include:

  1. Elasticity: Several spatial elements of varying stiffness, resulting from irregular meshes and/or inhomogeneous materials (Lew et al., 2003).

  2. Planetary Dynamics: -body problem with nonlinear gravitational forces, arising from pairwise inverse-square potentials. Multiple time scales result from the different distances between the bodies (Farr and Bertschinger, 2007).

  3. Highly Oscillatory Problems: Potential energy can be split into a “fast” linear oscillatory component and a “slow” nonlinear component. These problems are widely encountered in modeling molecular dynamics (Leimkuhler et al., 1996), but have also been used to model other diverse applications, for example, in computer animation (Eberhardt et al., 2000; Boxerman and Ascher, 2004).

Because these systems each satisfy a Lagrangian variational principle, they lend themselves readily to variational integrators: a class of geometric numerical integrators designed for simulating Lagrangian mechanical systems. By construction, variational integrators preserve a discrete version of this Lagrangian variational structure; consequently, they are automatically symplectic and momentum-conserving, with good long-time energy behavior (Marsden and West, 2001).

Explicit Methods, Multiple Time Stepping, and Resonance Instability.

The Störmer/Verlet (or leapfrog) method is one of the canonical examples of a geometric (and variational) numerical integrator (see Hairer et al., 2003). Yet, it and other simple, explicit time stepping methods do not perform well for problems with multiple time scales. The maximum stable time step for these methods is dictated by the stiffest mode of the underlying system; therefore, the fastest force dictates the number of evaluations that must be taken for all forces, despite the fact that the slow-scale forces may be (and often are) much more expensive to evaluate.

To reduce the number of costly function evaluations associated to the slow force, several explicit variational integrators use multiple time stepping, whereby different time step sizes are used to advance the fast and slow degrees of freedom. These include substepping methods, such as Verlet-I/r-RESPA and mollified impulse, where for each slow time step, an integer number of fast substeps are taken (Izaguirre et al., 2002). More recently, asynchronous variational integrators (AVIs) have been developed, removing the restriction for fast and slow time steps to be integer (or even rational) multiples of one another (Lew et al., 2003). Multiple-time-stepping methods can be more efficient than single-time-stepping explicit methods, like Störmer/Verlet, since one can fully resolve the fast oscillations while taking many fewer evaluations of the slow forces. This is especially advantageous for highly oscillatory problems, where the slow forces are nonlinear and hence more computationally expensive to evaluate.

One drawback of multiple-time-stepping methods, however, is that they can exhibit linear resonance instability. This phenomenon occurs when the slow impulses are nearly synchronized, in phase, with the the fast oscillations. These impulses artificially drive the system at a resonant frequency, causing the energy (and hence the numerical error) to increase without bound. The problem of numerical resonance is well known for substepping methods (Biesiadecki and Skeel, 1993), and has also recently been shown for AVIs as well—in fact, the subset of fast and slow time step size pairs leading to resonance instability is dense in the space of all possible parameters (Fong et al., 2007). Resonance instability can therefore be difficult to avoid, particularly in highly oscillatory systems with many degrees of freedom, as in molecular dynamics applications.

Implicit Methods for Single Time Stepping with Longer Step Sizes

Because multiple-time-stepping methods have these resonance problems, a number of single-time-stepping methods have been developed specifically for highly oscillatory problems. As noted earlier, single-time-stepping methods cannot fully resolve the fast oscillations without serious losses in efficiency. Therefore, the goal of these methods is to take long time steps, without actually resolving the fast oscillations, while still accurately capturing the macroscopic behavior that emerges from the coupling between fast and slow scales. The challenge is to design methods that allow for these longer time steps, without destroying either numerical stability or geometric structure.

One obvious candidate integrator is the implicit midpoint method, which is (linearly) unconditionally stable, as well as variational (hence symplectic) and symmetric. Unfortunately, the stability of the method comes at a cost: because the integrator is implicit in the slow force, which is generally nonlinear, a nonlinear system of equations must be solved at every time step. Therefore, just like the fully-resolved Störmer/Verlet method, this means that the implicit midpoint method requires an excessive number of function evaluations.

Implicit-Explicit Integration

For highly oscillatory problems, implicit-explicit (IMEX) integrators have been proposed as a potentially attractive alternative to either explicit, multiple-time-stepping methods or implicit, single-time-stepping methods. Rather than using separate fast and slow time step sizes, IMEX methods combine implicit integration (e.g., backward Euler) for the fast force with explicit integration (e.g., forward Euler) for the slow force. Because the fast force is linear, this semi-implicit approach requires only a linear solve for the implicit portion, as opposed to the expensive nonlinear solve that would be required for a fully implicit integrator, like the implicit midpoint method.

IMEX methods were developed by Crouzeix (1980), and have continued to progress, including the introduction of IMEX Runge–Kutta schemes for PDEs by Ascher et al. (1997). However, in all of these methods, the splitting is done at the level of the Euler–Lagrange differential equations, rather than at the variational level of the Lagrangian. Consequently, a wide variety of IMEX schemes have been created, both geometric and non-geometric, but in general they cannot guarantee properties such as symplecticity, momentum conservation, or good long-time energy behavior, which automatically hold for variational integrators. As an example of an IMEX integrator that is not “geometric” in the usual sense, consider the LI and LIN methods of Zhang and Schlick (1993), which combine the backward Euler method with explicit Langevin dynamics for molecular systems. In particular, to ameliorate the artificial numerical dissipation introduced by using backward Euler, these methods rely on stochastic forcing to inject the missing energy back into the system.

In this paper, we develop IMEX numerical integration from a Lagrangian, variational point of view. We do this by splitting the fast and slow potentials at the level of the Lagrangian action integral, rather than with respect to the differential equations or the Hamiltonian. From this viewpoint, implicit-explicit integration is an automatic consequence of discretizing the action integral using two distinct quadrature rules for the slow and fast potentials. The resulting discrete Euler–Lagrange equations coincide with a semi-implicit algorithm that was originally introduced by Zhang and Skeel (1997) as a “cheaper” alternative to the implicit midpoint method; Ascher and Reich (1999b) also studied a variant of this method for certain problems in molecular dynamics, replacing the implicit midpoint step by the energy-conserving (but non-symplectic) Simo–Gonzales method.

We also show that this variational IMEX method is free of resonance instabilities; the proof of this fact is naturally developed at the level of the Lagrangian, and does not require an examination of the associated Euler–Lagrange equations. We then compare the resonance-free behavior of variational IMEX to the multiple-time-stepping method r-RESPA in a numerical simulation of coupled slow and fast oscillators. Next, we evaluate the stability of the variational IMEX method, for large time steps, in a computation of slow energy exchange in the Fermi–Pasta–Ulam problem. Finally, we prove that the variational IMEX method accurately preserves this slow energy exchange behavior (as observed in the numerical experiments) by showing that it corresponds to a modified impulse method.

1.2. A Brief Review of Variational Integrators

The idea of variational integrators was studied by Suris (1990) and Moser and Veselov (1991), among others, and a general theory was developed over the subsequent decade (see Marsden and West, 2001, for a comprehensive survey).

Suppose we have a mechanical system on a configuration manifold , specified by a Lagrangian . Given a set of discrete time points with uniform step size , we wish to compute a numerical approximation , to the continuous trajectory . To construct a variational integrator for this problem, we define a discrete Lagrangian , replacing tangent vectors by pairs of consecutive configuration points, so that with respect to some interpolation method and numerical quadrature rule we have

Then the action integral over the whole time interval is approximated by the discrete action sum

If we apply Hamilton’s principle to this action sum, so that when variations are taken over paths with fixed endpoints, then this yields the discrete Euler–Lagrange equations

where and denote partial differentiation in the first and second arguments, respectively. This defines a two-step numerical method on , mapping . The equivalent one-step method on the cotangent bundle , mapping , is defined by the discrete Legendre transform

where the first equation updates , and the second updates .


Consider a Lagrangian of the form , where , is a constant mass matrix, and is a potential. If we use linear interpolation of with trapezoidal quadrature to approximate the contribution of to the action integral, we get

which we call the trapezoidal discrete Lagrangian. It is straightforward to see that the discrete Euler–Lagrange equations for correspond to the explicit Störmer/Verlet method. Alternatively, if we use midpoint quadrature to approximate the integral of the potential, this yields the midpoint discrete Lagrangian,

for which the resulting integrator is the implicit midpoint method.

2. A Variational IMEX Method

In this section, we show how to develop a variational integrator that combines aspects of the Störmer/Verlet and implicit midpoint methods mentioned above. The main idea is that, given a splitting of the potential energy into fast and slow components, we define the discrete Lagrangian by applying the midpoint quadrature rule to the fast potential and the trapezoidal quadrature rule to the slow potential. The resulting variational integrator is implicit in the fast force and explicit in the slow force. After this, we focus on the specific case of highly oscillatory problems, where the fast potential is quadratic (corresponding to a linear fast force). In this case, we show that the IMEX integrator can be understood as Störmer/Verlet with a modified mass matrix.

It should be noted that the results in this section can also be verified directly in terms of the numerical algorithm, and do not strictly require making use of the Lagrangian variational structure. However, we find the variational perspective to be useful and illustrative, both in arriving at this particular IMEX algorithm and in interpreting its numerical features.

2.1. The IMEX Discrete Lagrangian and Equations of Motion

Suppose that we have a Lagrangian of the form , where is a slow potential and is a fast potential, for the configuration space . Then define the IMEX discrete Lagrangian

using (explicit) trapezoidal approximation for the slow potential and (implicit) midpoint approximation for the fast potential. The discrete Euler–Lagrange equations give the two-step variational integrator on

and the corresponding discrete Legendre transform is given by

To see how this translates into an algorithm for a one-step integrator on , it is helpful to introduce the intermediate stages

Substituting these into the previous expression and rearranging yields the algorithm

Step 1:
Step 2:
Step 3:

where Step 2 corresponds to a step of the implicit midpoint method.

This can be summarized, in the style of impulse methods, as:

  1. kick: explicit kick from advances ,

  2. oscillate: implicit midpoint method with advances ,

  3. kick: explicit kick from advances .

In particular, notice that this reduces to the Störmer/Verlet method when and to the implicit midpoint method when . Also, if the momentum does not actually need to be recorded at the full time step (i.e., collocated with the position ), then Step 3 can be combined with Step 1 of the next iteration to create a staggered “leapfrog” method.

Interpretation as a Hamiltonian Splitting Method

This algorithm on can also be interpreted as a fast-slow splitting method (McLachlan and Quispel, 2002; Hairer et al., 2006, II.5 and VIII.4.1) for the separable Hamiltonian , where is the kinetic energy, as follows. Let denote the numerical flow of the implicit midpoint method with time step size , applied to the fast portion of the Hamiltonian , and let be the exact Hamiltonian flow for the slow potential (i.e., constant acceleration without displacement). Then the variational IMEX method has the flow map , which can be written as the following composition of exact and numerical flows:

This formulation highlights the fact that variational IMEX is symmetric (since it is a symmetric composition of symmetric methods) as well as symplectic (since it can be written as a composition of symplectic maps).

2.2. Application to Highly Oscillatory Problems

For highly oscillatory problems on , we start by taking a quadratic fast potential

A prototypical is given by the block-diagonal matrix , where some of the degrees of freedom are subjected to an oscillatory force with constant fast frequency . We also denote the slow force and assume, without loss of generality, that the constant mass matrix is given by . Therefore, the nonlinear system we wish to approximate numerically is

This is the conventional setup for highly oscillatory problems, used by Hairer et al. (2006, XIII) and others.

Applying the IMEX method to this example, we get the discrete Lagrangian

and so the two-step IMEX scheme is given by the discrete Euler–Lagrange equations

Combining terms, we can rewrite this as

which is equivalent to Störmer/Verlet with a modified mass matrix . This equivalence can similarly be shown to hold for the one-step formulation of the IMEX scheme on —that is, the two methods also produce the same , as well as the same .

In fact, this correspondence between IMEX and a modified Störmer/Verlet method is true not just in the discrete Euler–Lagrange equations, but in the discrete Lagrangian itself. This follows immediately from the following proposition.

Proposition 2.1.

Suppose we have a Lagrangian and its corresponding midpoint discrete Lagrangian . Next, define the modified Lagrangian , having the same quadratic potential but a different mass matrix , and take its trapezoidal discrete Lagrangian . Then when .


The midpoint discrete Lagrangian is given by

Now, notice that we can rearrange the terms

Therefore the discrete Lagrangian can be written in the trapezoidal form

which is precisely when . ∎

Corollary 2.2.

Consider a highly oscillatory system with an arbitrary slow potential , quadratic fast potential , and constant mass matrix , so that the Lagrangian and IMEX discrete Lagrangian are defined as above. Next, take the modified Lagrangian with the same potentials but different mass matrix . Then when .

2.3. Analysis of Linear Resonance Stability

To study the linear resonance stability of this IMEX integrator, we consider a model problem where and both correspond to linear oscillators. Let and , where for some , and again let the mass matrix . Although this is something of a “toy problem”—obviously, one could simply combine and into a single quadratic potential —it is illustrative for studying the numerical resonance of multiple-time-stepping methods, since the system has no external forcing terms and hence no real physical resonance.

To prove that the IMEX method does not exhibit linear resonance instability, we show that the stability condition only requires that the time step be stable for the explicit slow force, and is independent of the fast frequency . The idea of the proof is to use the results from 2.2, showing that the IMEX method is equivalent to Störmer/Verlet with a modified mass matrix, and then to apply the well-known stability criteria for Störmer/Verlet.

In particular, for a harmonic oscillator with unit mass and frequency , the Störmer/Verlet method is linearly stable if and only if , as can be shown by a straightforward calculation of the eigenvalues of the propagation matrix (Hairer et al., 2006, p. 23). For a system with constant mass and spring constant , this condition generalizes to .

Theorem 2.3.

The IMEX method is linearly stable, for the system described above, if and only if (i.e., if and only if is a stable time step size for the slow oscillator alone).


As proved in the previous section, the IMEX method for this system is equivalent to Störmer/Verlet with the modified mass matrix . Now, this modified oscillatory system has constant mass and spring constant . Therefore, the necessary and sufficient condition for linear stability is

and since the terms cancel on both sides, this is equivalent to , or . ∎

This shows that, in contrast to multiple-time-stepping methods, the IMEX method does not exhibit linear resonance instability. In particular, one can interpret the modified mass matrix as giving the system an effective frequency of , which attenuates the destabilizing high frequencies in the original system. It should be noted that nonlinear instability is known to be possible for the implicit midpoint method, although even that can be avoided with a time step size restriction that is considerably weaker than that required for explicit methods (see Ascher and Reich, 1999a).

3. Numerical Experiments

3.1. Coupled Linear Oscillators

To illustrate the numerical resonance behavior of the variational IMEX scheme, as compared with a multiple-time-stepping method, we consider the model problem from 2.3 for dimension (i.e., ). Figure 1 shows a log plot of the maximum absolute error in total energy (i.e., the Hamiltonian) for both r-RESPA and the variational IMEX method, for a range of frequencies . MATLAB simulations were performed over the time interval , with fixed time step size , and with the normalized frequency ranging over . Additionally, to fully resolve the fast oscillations, r-RESPA took 100 fast substeps of size for each full time step of size .

Figure 1. Maximum energy error of r-RESPA and variational IMEX, integrated over the time interval for a range of parameters . The r-RESPA method exhibits resonance instability near integer values of , while the variational IMEX method remains stable.

The r-RESPA method exhibits “spikes” in the total energy error near integer values of , corresponding to the parameters where resonance instability develops and the numerical solution becomes unbounded. (The finite size of these spikes is due to the fact that the numerical simulation was run only for a finite interval of time. Interestingly, one also sees “negative spikes,” where the fast and slow oscillations are exactly out-of-phase and cancel one another.) It should be noted that the small substep size of r-RESPA is sufficient for stable integration of the fast force alone; it is only the introduction of the slow force that makes things unstable. By contrast, the maximum energy error for the variational IMEX method is nearly constant for all values of , showing no sign of resonance. This is fully consistent with the theoretical result obtained in Theorem 2.3.

3.2. The Fermi–Pasta–Ulam Problem

As an example of a nontrivial highly oscillatory problem with nonlinear slow potential, we chose the modified Fermi–Pasta–Ulam (FPU) problem considered by Hairer et al. (2006, I.5 and XIII), whose treatment we will now briefly review. The FPU problem consists of unit point masses, which are chained together, in series, by alternating weak nonlinear springs and stiff linear springs. (This particular setup is due to Galgani et al., 1992, and is a variant of the problem originally introduced by Fermi et al., 1955.) Clearly, this system becomes rather trivial if we make the stiff springs “infinitely stiff,” replacing them by rigid constraints (as done by some numerical methods, such as SHAKE/RATTLE). However, for finite stiffness, the FPU system exhibits interesting dynamics due to the coupling between fast and slow springs.

Let us denote the displacements of the point masses by (where the endpoints are taken to be fixed), and their conjugate momenta by for . In these variables, the FPU system has the Hamiltonian

which contains a quadratic potential for the stiff linear springs, each with frequency , and a quartic potential for the soft nonlinear (cubic) springs. However, it is helpful to perform the coordinate transformation (following Hairer et al., 2006, p. 22)

so that (modulo rescaling) corresponds to the location of the th stiff spring’s center, corresponds to its length, and are the respective conjugate momenta. Writing the Hamiltonian in these new variables, we have

which considerably simplifies the form of the fast quadratic potential.

Following the example treated numerically by Hairer et al. (2006); McLachlan and O’Neale (2007), we consider an instance of the FPU problem, integrated over the time interval , with parameters , , whose initial conditions are

with zero for all other initial values. This displays an interesting and complex property of the FPU problem, called slow energy exchange, which results from the slow nonlinear coupling between the stiff springs. If we consider only the energies in the stiff springs, written as

then the initial conditions start with all of the energy in and none in . Over the course of the time interval, this energy is transferred in a characteristic way from to , gradually transitioning through the middle spring . Furthermore, the total stiff energy remains nearly constant, i.e., is an adiabatic invariant of the system.

(a) Reference solution: Störmer/Verlet with time step size .
(b) Störmer/Verlet with .
(c) Störmer/Verlet with .
(d) IMEX with .
(e) IMEX with .
(f) IMEX with .
(g) IMEX with .
(h) IMEX with .
(i) IMEX with .
Figure 2. The IMEX method robustly captures slow energy exchange in the Fermi–Pasta–Ulam problem with , even for large time steps. Because the fast force is integrated implicitly, IMEX remains stable and degrades gradually as the time step size increases—unlike the fully explicit Störmer/Verlet method, which rapidly becomes unstable.

Figure 2 shows several numerical simulations of this FPU energy exchange, computed both with Störmer/Verlet and with the variational IMEX method, for different choices of time step size. The first plot is a reference solution, computed using Störmer/Verlet with , fully resolving the fast oscillations. However, we see that the Störmer/Verlet solution’s quality and stability degrade rapidly as we increase the step size (for , we have , which is near the upper end of the stability region ). By contrast, the variational IMEX method performs extremely well for , degrading gradually as the time step size increases. Even as the numerical solution begins to undergo serious degradation for , the qualitative structure of the energy exchange behavior between is still maintained. Compare Hairer et al., 2006, p. 24, Figure 5.3; see also McLachlan and O’Neale, 2007, who examine a wide variety of geometric integrators, particularly trigonometric integrators, for the FPU problem, with respect to both resonance stability and slow energy exchange. In particular, these authors found that the existing trigonometric integrators exhibit a trade-off between correct energy exchange behavior and resonance stability, and that these features tend to be mutually exclusive.

(a) IMEX with
(b) IMEX with
Figure 3. Numerical simulation of the FPU problem for , which shows the behavior of the IMEX method on the scale. For , we already have , yet the oscillatory behavior and adiabatic invariant are qualitatively correct. By contrast, for , the method has begun to blow up; oscillatory coupling is a drawback of implicit midpoint methods for large time steps.

In Figure 3 we show the numerical behavior of the variational IMEX method, applied to this same FPU problem, on a longer time scale () and for large time steps (). At , the IMEX simulation still displays the correct qualitative energy behavior, with respect to both the slow energy exchange and the adiabatic invariant , and the numerical solution remains bounded. However, by , numerical stability has broken down, as oscillatory coupling in the fast modes leads to unbounded amplitude growth. This illustrates one of the drawbacks of implicit midpoint-type methods: despite the lack of linear resonances, numerical instability can still result for very large time steps due to nonlinear coupling (Ascher and Reich, 1999a, b).

This example was chosen to demonstrate that the variational IMEX method does not attain its stability merely by “smoothing out” the fast frequencies, in a way that might destroy the structure of any fast-slow nonlinear coupling. Rather, despite the fact that it does not resolve the fast frequencies, the method is still capable of capturing the complex multiscale interactions seen in the FPU problem.

4. Analysis of Slow Energy Exchange in the IMEX Method

In the previous section, the numerical experiments for the Fermi–Pasta–Ulam problem seemed to suggest that the variational IMEX method preserves the slow energy exchange between the fast oscillatory modes. This is somewhat surprising, since the method does not actually resolve these fast oscillations. However, in this section, we will prove that, in fact, this method does accurately reproduce the slow energy exchange behavior, as long as the numerical solutions remain bounded. This is demonstrated by showing that the variational IMEX method can be understood as a modified impulse method; that is, the midpoint step exactly resolves the oscillations of some modified differential equation. We can then apply some of the existing theory about numerical energy exchange for impulse methods. (It should be noted that impulse methods, which originated with the work of Deuflhard, 1979, can be understood as a special case of trigonometric integrators when applied to highly oscillatory problems.)

First, let us rewrite the fast oscillatory system as the first-order system

so it follows that the exact solution satisfies

We will now show that the implicit midpoint method effectively replaces this rotation matrix for by the rotation matrix corresponding to a modified . In the transformed coordinates just introduced, the implicit midpoint method has the expression

Therefore, if we take the skew matrix

it follows that

Notice that the expression is the Cayley transform, which maps skew matrices to special orthogonal matrices (and can be seen as an approximation to the matrix exponential map, which gives the exact solution). Hence the stability matrix is special orthogonal, so we can write

for some modified frequency . Therefore, the stability matrix for the implicit midpoint method corresponds to the exact flow matrix for a modified oscillatory system.

As an example, suppose we have for some constant frequency . Applying the Cayley transform, it can be seen that the modified frequency satisfies

Squaring both sides, this becomes

which finally gives the solution for the modified frequency,


This perspective provides another explanation as to why the variational IMEX method does not exhibit resonance: we always have . In fact, the Cayley transform does not map to a rotation by , except in the limit as . Therefore, for any finite and , we will never encounter the resonance points corresponding to integer multiples of .

As an aside, this also leads to another possible interpretation for the onset of nonlinear instability, if the time step size becomes too large (as we saw in Figure 3). Since , the modified frequency must shrink as grows. Informally, then, if is very small, this can be seen as leading to amplitude growth in the fast modes, since it requires less energy to induce this amplification.

Since the implicit midpoint method has now been seen as the exact solution of a modified system, we can write the variational IMEX method as the following modified impulse scheme:

Step 1:
Step 2:
Step 3:

Suppose again that for some constant frequency , and likewise . (This includes the case of the FPU problem.) The slow energy exchange behavior of this system was analyzed in detail by Hairer et al. (2006, XIII, see especially p. 495) using the so-called modulated Fourier expansion; we now give a brief, high-level summary of this work. In the notation of Hairer et al., the exact solution is asymptotically expanded as , where

Here, is a smooth real-valued function and is a smooth complex-valued function, and these can be partitioned according to the blocks of as and . They show that, assuming the exact solution has energy bounded independent of , this implies , so describes the slow-scale evolution of the system. Plugging in this ansatz for a highly oscillatory system, and eliminating the variables and , the slow evolution turns out to be described by

where here the slow force has also been block-decomposed as .

Hairer et al. compare the above with the numerical solution for a trigonometric integrator, which is similarly expanded as

with and . For the unmodified Deuflhard/impulse method, in particular, the slow-scale numerical evolution is given by

which implies that the equation for the impulse solution is consistent with that for .

We now finally have what we need to prove our main result on the slow energy exchange behavior of the variational IMEX method.

Theorem 4.1.

Let the variational IMEX method be applied to the highly oscillatory problem above, and suppose the numerical solution remains bounded. Then the ordinary differential equation for , describing the slow energy exchange in the numerical solution, is consistent with that for in the exact solution; this holds up to order .


As we have previously shown, the IMEX scheme corresponds to a modified impulse method with frequency . Therefore, to get the equation for , we must modify the equation above for by replacing with on the left hand side. However, notice that Step 2 of the modified method advances the original state vector , rather than the modified . Changing from to also introduces a scaling factor of on the right hand side. Therefore, the variational IMEX solution satisfies the slow-scale equation

Finally, cancelling the factors and multiplying by , we once again get

which is the same as that for the coefficient in the original, unmodified impulse method. This completes the proof. ∎


The authors would like to thank Will Fong and Adrian Lew for helpful conversations about the stability of AVI methods, Houman Owhadi for suggesting that we examine the FPU problem, Dion O’Neale for advice regarding the FPU simulation plots, Teng Zhang for pointing out a bug in an earlier version of our FPU simulation code, and Jerry Marsden and Reinout Quispel for their valuable suggestions and feedback about this work. We especially wish to acknowledge Marlis Hochbruck, Arieh Iserles, and Christian Lubich for suggesting that we try to understand the variational IMEX scheme as a modified impulse method; this perspective led us directly to the results in 4. Finally, this paper benefited greatly from the thoughtful critiques and suggestions of the anonymous referees, whom we also wish to thank.


  1. Ascher, U. M., and S. Reich (1999a), The midpoint scheme and variants for Hamiltonian systems: advantages and pitfalls. SIAM J. Sci. Comput., 21 (3), 1045–1065 (electronic).
  2. Ascher, U. M., and S. Reich (1999b), On some difficulties in integrating highly oscillatory Hamiltonian systems. In Computational molecular dynamics: challenges, methods, ideas (Berlin, 1997), volume 4 of Lecture Notes in Computational Science and Engineering, pages 281–296. Springer, Berlin.
  3. Ascher, U. M., S. J. Ruuth, and R. J. Spiteri (1997), Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Appl. Numer. Math., 25 (2-3), 151–167. Special issue on time integration.
  4. Biesiadecki, J. J., and R. D. Skeel (1993), Dangers of multiple time step methods. J. Comput. Phys., 109 (2), 318–328.
  5. Boxerman, E., and U. Ascher (2004), Decomposing cloth. In SCA ’04: Proceedings of the 2004 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pages 153–161. Eurographics Association, Aire-la-Ville, Switzerland. doi:10.1145/1028523.1028543.
  6. Crouzeix, M. (1980), Une méthode multipas implicite-explicite pour l’approximation des équations d’évolution paraboliques. Numer. Math., 35 (3), 257–276.
  7. Deuflhard, P. (1979), A study of extrapolation methods based on multistep schemes without parasitic solutions. Z. Angew. Math. Phys., 30 (2), 177–189.
  8. Eberhardt, B., O. Etzmuß, and M. Hauth (2000), Implicit-explicit schemes for fast animation with particle systems. In Proceedings of the 11th Eurographics Workshop on Computer Animation and Simulation (EGCAS), pages 137–151. Springer, Interlaken, Switzerland.
  9. Farr, W. M., and E. Bertschinger (2007), Variational integrators for the gravitational -body problem. Astrophys. J., 663 (2), 1420–1433. arXiv:astro-ph/0611416.
  10. Fermi, E., J. Pasta, and S. Ulam (1955), Studies of nonlinear problems. Report LA-1940. Los Alamos National Laboratory, Los Alamos, NM.
  11. Fong, W., E. Darve, and A. Lew (2007), Stability of asynchronous variational integrators. In PADS ’07: Proceedings of the 21st International Workshop on Principles of Advanced and Distributed Simulation, pages 38–44. IEEE Computer Society Press, Washington, DC. doi:10.1109/PADS.2007.29.
  12. Galgani, L., A. Giorgilli, A. Martinoli, and S. Vanzini (1992), On the problem of energy equipartition for large systems of the Fermi-Pasta-Ulam type: analytical and numerical estimates. Phys. D, 59 (4), 334–348.
  13. Hairer, E., C. Lubich, and G. Wanner (2003), Geometric numerical integration illustrated by the Störmer-Verlet method. Acta Numer., 12, 399–450.
  14. Hairer, E., C. Lubich, and G. Wanner (2006), Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition. Structure-preserving algorithms for ordinary differential equations.
  15. Izaguirre, J. A., Q. Ma, T. Matthey, J. Willcock, T. Slabach, B. Moore, and G. Viamontes (2002), Overcoming instabilities in Verlet-I/r-RESPA with the mollified impulse method. In T. Schlick and H. H. Gan, editors, Computational Methods for Macromolecules: Challenges and Applications—Proceedings of the 3rd International Workshop on Algorithms for Macromolecular Modeling, volume 24 of Lecture Notes in Computational Science and Engineering (LNCSE), pages 146–174. Springer, Berlin.
  16. Leimkuhler, B. J., S. Reich, and R. D. Skeel (1996), Integration methods for molecular dynamics. In Mathematical approaches to biomolecular structure and dynamics (Minneapolis, MN, 1994), volume 82 of IMA Vol. Math. Appl., pages 161–185. Springer, New York.
  17. Lew, A., J. E. Marsden, M. Ortiz, and M. West (2003), Asynchronous variational integrators. Arch. Ration. Mech. Anal., 167 (2), 85–146.
  18. Marsden, J. E., and M. West (2001), Discrete mechanics and variational integrators. Acta Numer., 10, 357–514.
  19. McLachlan, R. I., and D. R. J. O’Neale (2007), Comparison of integrators for the Fermi-Pasta-Ulam problem. Preprint NI07052-HOP. Isaac Newton Institute for Mathematical Sciences, Cambridge, UK. Available from: http://www.newton.ac.uk/preprints/NI07052.pdf.
  20. McLachlan, R. I., and G. R. W. Quispel (2002), Splitting methods. Acta Numer., 11, 341–434.
  21. Moser, J., and A. P. Veselov (1991), Discrete versions of some classical integrable systems and factorization of matrix polynomials. Comm. Math. Phys., 139 (2), 217–243.
  22. Suris, Y. B. (1990), Hamiltonian methods of Runge-Kutta type and their variational interpretation. Mat. Model., 2 (4), 78–87.
  23. Zhang, G., and T. Schlick (1993), LIN: a new algorithm to simulate the dynamics of biomolecules by combining implicit-integration and normal mode techniques. J. Comput. Chem., 14 (10), 1212–1233. doi:10.1002/jcc.540141011.
  24. Zhang, M., and R. D. Skeel (1997), Cheap implicit symplectic integrators. Appl. Numer. Math., 25 (2-3), 297–302. Special issue on time integration.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description