Numerical integrators for motion under a strong constraining force
This paper deals with the numerical integration of Hamiltonian systems in which a stiff anharmonic potential causes highly oscillatory solution behavior with solution-dependent frequencies. The impulse method, which uses micro- and macro-steps for the integration of fast and slow parts, respectively, does not work satisfactorily on such problems. Here it is shown that variants of the impulse method with suitable projection preserve the actions as adiabatic invariants and yield accurate approximations, with macro-stepsizes that are not restricted by the stiffness parameter.
Keywords. Oscillatory Hamiltonian systems, varying high frequencies, impulse method, mollified impulse method, projected impulse method.
AMS subject classifications. 34E13, 65L11, 65P10, 70H11, 70H15, 70H45
We are interested in the efficient numerical integration of Hamiltonian systems in which a stiff anharmonic potential causes highly oscillatory solution behavior with state-dependent slowly varying high frequencies.
1.1 The highly oscillatory Hamiltonian system
We consider Hamiltonians as studied, in varying degrees of generality and with different analytical techniques, by Rubin & Ungar , Takens , Bornemann , Lorenz  and Hairer, Lubich & Wanner [11, Section XIV.3]:
depending on positions and momenta . The mass matrix is symmetric and positive definite and depends smoothly on . The slow potential is smooth, and the stiff potential with a smooth function attains its minimum value on a -dimensional manifold
We assume that the potential is strongly convex along directions non-tangential to . More precisely, there exists such that for the Hessian satisfies
for all vectors in the -orthogonal complement of the tangent space .
Furthermore, we assume that a constraint function , with , is known such that
and the derivative matrix is of full rank on .
The corresponding system of Hamiltonian differential equations reads
A simple, yet nontrivial model example is the stiff spring double pendulum. The Hamiltonian reads where is the kinetic energy, , and is the stiff potential depending on the small parameter . The parameters denote the lengths of the springs, are the large spring constants.
Example 1 helps to fix ideas on a simple toy model. Obviously it extends to chains of stiff springs, which describe the dynamics of chains of atoms in a molecule with almost rigid bonds, cf., e.g., .
1.2 The effective Hamiltonian system
It has been known since Rubin & Ungar  that the motion of the system in the limit differs from the Hamiltonian dynamics constrained to the manifold (the rigid double pendulum in the above example) for general initial values that have an energy bounded independently of ,
Note that the set of admissible initial values satisfying (6) depends on . The effective constrained Hamiltonian has a correction potential ,
The correction potential depends on the frequencies , i.e., the square roots of the nonzero generalized eigenvalues of the pencil , and on parameters , known as the actions, which are determined by the initial values of (5). The actions vanish for consistent initial data that satisfy and .
As is outlined in  and will be recapitulated in Section 3, the effective Hamiltonian can be found by transforming the system to separate slow and fast variables as in  and [11, Section XIV.3], and transforming the obtained slow system () back via the effective dynamics of the fast variables.
The effective constrained Hamiltonian system is then given by
with Lagrange multipliers . This differs from the usual constrained equations of motion through the correction force .
To initial values of system (5) with bounded energy (6), we associate consistent initial values for the effective system (8). These are chosen by projecting -orthogonally onto the manifold of consistent values:
With the projection
the second equation can be rewritten as , and along the solution of (8) we note .
The effective Hamiltonian system describes the limit dynamics on the constraint manifold as long as the solution-dependent frequencies remain separated and are non-resonant: for some ,
where the constants in the -notation depend on and deteriorate as ; see [19, 11, 3]. In the case of a single frequency (), where no separation and non-resonance conditions appear, the approximation of the full system by the effective system was already studied by Rubin & Ungar . Note that the above estimates also imply
which is equivalent to saying that the tangential component of the velocity error is . The normal component of the velocity is, however, disregarded in the constrained effective equation.
Conditions (11) and (12) may appear rather severe at first sight, but in fact conditions of this type are needed for the above approximation result for the effective dynamics. Using the techniques of [11, Chap. XIV] it can be shown that the order of this approximation is still if have zeros of multiplicity . However, the separation cannot be omitted. If the distance of two frequencies becomes smaller than , then the slow motion can depend very sensitively on the initial values, and it is no longer approximated by the dynamics of the effective Hamiltonian system; see Takens . The indeterminacy of the slow motion in the case of non-separated frequencies is termed Takens chaos in .
1.3 Outline of the paper and relation to the literature
The objective of this paper is to devise and analyze a two-scale integrator for the highly oscillatory Hamiltonian system (5), such that for a macro-stepsize that is not restricted by , the method yields an error in the positions and the projected momenta over time intervals .
This paper is part of the vast literature on the numerical solution of highly oscillatory differential equations; see, e.g., the reviews [5, 15]. Recent work on the numerical integration of highly oscillatory mechanical systems includes [1, 4, 17, 18, 20].
While much work has been done on systems with constant high frequencies, the numerical analysis of the present case of solution-dependent high frequencies or even just the case of explicitly time-dependent high frequencies is scarce; see [11, Chapter XIV]. An important aspect here is to preserve the adiabatic invariants (see, e.g.,  for this notion) in the numerical discretization.
In this paper we study two-scale time integrators for (1) which aim at solving the effective system (8) over the time scale without, however, explicitly evaluating the correction force . This additional force is, in general, directly accessible only via a series of computationally expensive, nonlinear implicit coordinate transformations. Moreover, even in cases where the correction force is computationally accessible, it is of interest to have a numerical method that is able to monitor the possible breakdown of the validity of the effective equation due to the loss of adiabatic invariance of the actions in cases where frequencies come close or become resonant.
Heterogeneous multiscale methods (HMM) [6, 8, 7] have been developed for the very purpose to handle situations where the underlying effective dynamics is not known. In Brumm & Weiss  an HMM-approach for highly oscillatory mechanical systems with solution-dependent frequencies is analyzed. This approach shows, however, major drawbacks because of difficulties in initializing the micro-simulation.
In this article, we follow the alternative idea of the impulse method where the Hamiltonian is split into the slow potential and the fast part including the kinetic and stiff potential energy. The slow part is integrated in macro-steps, the fast part uses micro-steps. As it turns out, this must be complemented with a suitable projection to lead to a method with satisfactory error behavior.
We proceed as follows: In Section 2 we formulate the impulse method, a mollified impulse method, and a novel projected impulse method for highly oscillatory mechanical systems with solution-dependent frequencies. We state the main convergence theorem and show results of numerical experiments that highlight different behavior of the various methods. In Section 3 we transform the system, following  and [11, Section XIV.3] to variables that are appropriate for the further analysis. Moreover, a further mollified impulse method with a projection mollifier in the transformed variables is introduced, which is computationally impractical but serves as a theoretical reference method for the error analysis. This method is studied in Section 4. Using the obtained results, the analysis of the mollified and projected impulse methods of Section 2 is done in Section 5.
2 Numerical methods and statement of the main result
2.1 Impulse method
The impulse method was introduced in the context of the numerical treatment of molecular dynamics (Grubmüller, Heller, Windemuth & Schulten , Tuckerman, Berne & Martyna ). A mathematical study of this method is given by García-Archilla, Sanz-Serna & Skeel . The idea is to split the Hamiltonian
and to approximate the exact flow by the following symmetric decomposition:
Since the flow of the slow part can be trivially solved, one step is equivalent to
oscillate: solve system (5) with and initial values over a time step to obtain ,
another kick: .
Step 2. is solved approximately using, e.g., the Störmer–Verlet method with micro-stepsizes, or alternatively using a large-timestep method in suitably transformed variables (an adiabatic integrator) as in .
Compared to a direct numerical integration of the full system (5) with small stepsizes, this method saves many evaluations of the slow force , which is often the computationally most expensive part.
For our numerical experiments we consider the stiff spring double pendulum with initial values
and the parameters , over the time interval . In this situation the frequencies remain well-separated.
We observe unsatisfactory behavior of the impulse method in Figure 1. Here and in all following figures the dash-dotted straight line has slope 2, corresponding to the desired error behavior. We used the Störmer–Verlet method with very small stepsize ( for the impulse method and for the following methods) for the micro-integration in order to avoid any significant influence on the overall error. Throughout all computations, as a reference, the effective system (8) is approximated in transformed variables (see (29)) by the Störmer–Verlet method with small stepsize, the results being translated back into cartesian coordinates.
2.2 Mollified impulse method
García-Archilla, Sanz-Serna & Skeel  and Izaguirre, Reich & Skeel  improve the impulse method by replacing the slow potential by a mollified potential , where is an averaged or suitably projected value of . The mollified force then reads
The mollification considered in the present paper is given by the -orthogonal projection onto the configuration manifold , i.e., with
As it turns out in the analysis, the unsatisfactory behavior of the impulse method is due to -orthogonal components of the slow forces . The mollification reduces those -orthogonal components. Indeed, we observe the following.
Proof. The condition is equivalent to . Noting that is invertible in view of the full rank of , the implicit function theorem then yields such that , and hence . Differentiating both equations in (15) yields
Inserting the first into the second equation permits us to compute
Reinserting this expression into the first equation yields the stated result on recalling the definition of .
2.3 Projected impulse method
The preceding lemma motivates us to simplify the method by projecting the slow force:
oscillate: solve system (5) with and initial values over a time step obtaining ,
Using this new simplified scheme, we observe convergence behavior as in the case of the mollified impulse method, see Figure 3.
We have sacrificed the symplecticity of the method which is not of main interest here, but have nevertheless maintained the time-reversal symmetry.
2.4 Main Result
The idea of a projection as a mollification is proposed in . There, the use of this idea is shown experimentally but no analysis is given. On the other hand, in [14, 11] the adiabatic nature of the systems of interest is revealed by applying a series of canonical transformations. Combining the different ideas and techniques, we are now able to formulate and prove the result about the global error of the mollified impulse method with the projection mollifier. Additionally, we prove the same result for the computationally simpler projected impulse method. The proof is based on a further, different mollification introduced in Section 3.
Let the initial values satisfy the energy bound (6) and assume that the frequencies remain separated and non-resonant (see conditions (11)-(12)) along the solution of (5) for . Then, the errors of the mollified impulse method and of the projected impulse method after steps with stepsize satisfy
Combined with (13) this also yields the error bounds with respect to the solution of the highly oscillatory problem
We note, however, that and . Moreover, the method does not converge to the solution of the highly oscillatory system for a fixed as . This causes no problems since the interest of the method lies in the use of large step sizes . Note that in Theorem 2.2 there is no restriction of the step size by the small parameter .
Theorem 2.2 explains the error behavior observed in Figures 2 and 3.
3 Transformed variables and another mollified impulse method
where and and the appearing functions have all their partial derivatives bounded independently of and are as follows:
is a symmetric positive definite matrix;
is a diagonal matrix with positive entries, the frequencies ;
is a symmetric matrix with ;
The assumption (6) of bounded energy now becomes
We define the actions
see , p. 562.
As will become clear from our analysis, this is a key property that should be transfered to the numerical method. However, if we express the impulse method in the transformed variables, then the kick step becomes
and we see that the actions are not approximately preserved. This is illustrated in Figure 4. The non-preservation of the actions is at the base of the disappointing numerical behavior observed in Figure 1.
On the other hand, if we use a mollified impulse method for which the modified potential is chosen, in the transformed variables, as
then , and hence the actions are exactly preserved in the kick step. While this method is not practical in that it would require performing the coordinate transformation from to , it gives much theoretical insight into the error propagation behavior. We will therefore study its error in the next section. Subsequently we will interpret the mollified and projected impulse methods of Section 2, which work in the original variables, as perturbations of this theoretically interesting method.
As a numerical illustration, in Figure 6 we use the transformed-variable mollified impulse method with the initial value of Section 2 for the stiff spring double pendulum. We observe similar results as in Figures 2 and 3.
4 Error analysis of the transformed-variable method
We will show almost-conservation of the actions along the numerical solution, with .
Assume the energy bound (6). Furthermore, assume that the frequencies of the transformed-variable mollified impulse method with modified potential (22) remain separated and non-resonant (see conditions (11)-(12)) for . Then, this method approximately preserves the actions:
The constant symbolized by is independent of , , and .
To prove this result, we first need to look in more detail into the differential equations in the transformed variables. As is shown in , p. 560, the Hamiltonian equations of motion take the form
with functions , and , . Moreover we have (omitting the arguments in , which are all )
where is an matrix and the functions and are bilinear.
and introduce the diagonal phase matrix by
Following , p. 561, we transform the oscillatory part of the solution to adiabatic variables
We further introduce the matrices and as
so that and . In adiabatic variables, the differential equation for the oscillatory part becomes
The functions are those appearing in the remainder terms and in (24).
Proof. (of Theorem 4.1) We rewrite the mollified impulse method in adiabatic variables and note that the kick steps do not change the adiabatic variables: in the th time step, and . We thus obtain
with initial value on the interval . All terms with superscript are defined with respect to the solution of the oscillation step of the mollified impulse method on .
In the remaining part of the proof we show . The techniques are more or less the same as for the exact solution presented in . Therefore, we just consider the first term of the righthand side in (27). For partial integration gives
where the latter term is of size in the case of separated frequencies. Taking into account the -jumps from to and noting that and , we prove the same bound for the first term. We have thus shown
and since , the result follows.
We are now in the situation to prove an error bound.
Assume the energy bound (6). Furthermore, assume that the frequencies remain separated and non-resonant (see conditions (11)-(12)) on a fixed time interval . Then, the error of the transformed-variable mollified impulse method of Section 3 after steps with stepsize satisfies
The constants symbolized by do not depend on , and with .
We note, however, that the normal components of the momenta are not approximated correctly: we only have .
Proof. We consider the method in the slow components as a perturbed variant of the Störmer–Verlet scheme applied to the slow system
More precisely, if we write the Störmer–Verlet scheme for (29) in one-step form as
then the slow components of the mollified impulse method for (23) fulfill
with a local error , because
and by Theorem 4.1. Application of the discrete Gronwall Lemma gives the desired result for the slow components in the variables : for ,
For the fast variables we have, using (28),
but in view of the phase difference this does not yield an approximation estimate. Transforming back to the coordinates and considering the rescaling and the Lipschitz-continuity of the transformations then gives the result.
5 Error analysis of the mollified and projected impulse methods
In order to analyze the mollified and projected impulse methods of Section 2, we have to derive an appropriate expression of the kick-step in the transformed variables . We show that both methods are -perturbations of the transformed-variable mollified impulse method of Section 3, which uses the modified potential (22), that is, in the original variables, the modified potential with
The following result is essential in relating the various methods.
For the mollifier we have, under the bounded-energy condition ,
where the projection is defined in (10).
Proof. As the construction in [11, Chapter XIV.3] shows, the transformation is composed as with an -independent transformation and a rescaling . Since for the bounded-energy condition is equivalent to , we obtain
The proof of the result for is then obtained from the identity
used for . This identity is obtained from the transformation laws as follows: Under a change of coordinates , the constraint function changes to and its derivative to . The inverse mass matrix changes to . Consequently, the projection transforms as
For , the transposed derivative transforms in the same way for with :
Now, in the variables we have and a block diagonal mass matrix , and on the other hand . This gives us
and hence the result follows.
Using Lemma 5.1 and the corresponding result Lemma 2.1 for the mollified impulse method of Section 2, we find for the kick step of the mollified - and projected impulse methods expressed in the variables of Section 3
Here, the additional factor is due to the rescaling of the fast positions and momenta in the transformation. For the actions, we then obtain the estimate in the kick step, and as in Section 4 we obtain the near-invariance of the actions along the numerical solution.