Numerical integrators for motion under a strong constraining force

# Numerical integrators for motion under a strong constraining force

Abstract. 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

## 1 Introduction

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 [16], Takens [19], Bornemann [2], Lorenz [14] and Hairer, Lubich & Wanner [11, Section XIV.3]:

 H(x,y) =12yTM(x)−1y+U(x)+1ε2V(x),0<ε≪1, (1)

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

 V ={x∈D:V(x)=minV=0}. (2)

We assume that the potential is strongly convex along directions non-tangential to . More precisely, there exists such that for the Hessian satisfies

 vT∇2V(x)v ≥α⋅vTM(x)v (3)

for all vectors in the -orthogonal complement of the tangent space .

Furthermore, we assume that a constraint function , with , is known such that

 V ={x∈D:g(x)=0} (4)

and the derivative matrix is of full rank on .

The corresponding system of Hamiltonian differential equations reads

 (5)
###### Example 1

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., [13].

### 1.2 The effective Hamiltonian system

It has been known since Rubin & Ungar [16] 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 ,

 H(x0,y0) ≤Const. (6)

Note that the set of admissible initial values satisfying (6) depends on . The effective constrained Hamiltonian has a correction potential ,

 Heff(X,Y)=12YTM(X)−1Y+U(X)+W(I,X),0=g(X). (7)

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 [3] 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 [14] 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

 (8)

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:

 X0=x0+M(x0)−1G(x0)Tλ,0=g(X0),Y0=y0+G(X0)Tμ,0=G(X0)M(X0)−1Y0. (9)

With the projection

 P(x)=I−Q(x),Q(x)=[GT(GM−1GT)−1GM−1](x), (10)

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 ,

 ∣∣ωj(X(t))−ωk(X(t))∣∣ ≥δ for j≠k (11) ∣∣ωj(X(t))±ωk(X(t))±ωl(X(t))∣∣ ≥δ for all j,k,l. (12)

If these conditions are satisfied for , then we have for the corresponding solutions of (5) and (8) over this time interval

 X(t)−x(t)=O(ε)Y(t)−P(x(t))y(t)=O(ε), (13)

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 [16]. Note that the above estimates also imply

 P(X(t))(y(t)−Y(t))=O(ε),

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 [19]. The indeterminacy of the slow motion in the case of non-separated frequencies is termed Takens chaos in [2].

### 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., [12] 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 [3] 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 [14] 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 [10], Tuckerman, Berne & Martyna [21]). A mathematical study of this method is given by García-Archilla, Sanz-Serna & Skeel [9]. The idea is to split the Hamiltonian

 H(x,y) =Hfast(x,y)+U(x)

and to approximate the exact flow by the following symmetric decomposition:

 φHh≈φslowh/2∘φfasth∘φslowh/2.

Since the flow of the slow part can be trivially solved, one step is equivalent to

1. kick: ,

2. oscillate: solve system (5) with and initial values over a time step to obtain ,

3. 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 [14].

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

 x(0) y(0) =(0,0,0,0)T (14)

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 [9] and Izaguirre, Reich & Skeel [13] 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

 −∇¯¯¯¯U(x) =−α′(xn)T∇U(α(xn)).

The mollification considered in the present paper is given by the -orthogonal projection onto the configuration manifold , i.e., with

 X=x+M(x)−1G(x)Tλ,0=g(X). (15)

Using the same initial value (14) as in the case of the impulse method, we observe better convergence behavior of the positions and the projected momenta, see Figure 2.

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.

###### Lemma 2.1

Under the bounded-energy condition , the mollifier of (15) satisfies

 α(x) =x+O(ε), α′(x)T =P(x)+O(ε),

where the projection is defined in (10).

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

 α′(x) =I+M−1(x)G(x)Tλ′(x)+O(ε) 0 =G(α(x))α′(x).

Inserting the first into the second equation permits us to compute

 λ′(x)=−(GM−1GT)−1G(x)+O(ε).

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:

1. kick: ,

2. oscillate: solve system (5) with and initial values over a time step obtaining ,

3. kick: .

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 [13]. 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.

###### Theorem 2.2

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

 xn−X(tn)=O(h2)+O(ε)P(xn)yn−Y(tn)=O(h2)+O(ε) (16)

where is the solution of the effective Hamiltonian system (8) with initial values defined by (9). The constants symbolized by do not depend on , and with .

Combined with (13) this also yields the error bounds with respect to the solution of the highly oscillatory problem

 xn−x(tn)=O(h2)+O(ε)P(xn)(yn−y(tn))=O(h2)+O(ε). (17)

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

Under conditions (2)-(4), [11, Section XIV.3] and [14] show that there exists a canonical change of coordinates of the separated form , , which transforms the Hamiltonian (1) into the form

 H(q,p) = 12pT0M0(q0)−1p0+12εpT1Ω(q0)p1+12εqT1Ω(q0)q1 (18) + 12(p0ε−1/2p1)TR(q0,ε1/2q1)(p0ε−1/2p1)+ˇU(q0,ε1/2q1),

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 ;

• for .

The assumption (6) of bounded energy now becomes

 q=O(ε1/2), p=O(ε1/2). (19)

We define the actions

 Ik=12ε(q21,k+p21,k),k=1,…,m. (20)

Under the separation and non-resonance conditions (11)–(12), the actions are adiabatic invariants: they remain nearly constant along solutions of the Hamiltonian system with bounded energy,

 Ik(t)=Ik(0)+O(ε), (21)

see [11], 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

 ¯¯¯¯U(q)=ˇU(q0,0), (22)

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 .

###### Theorem 4.1

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:

 In =I0+O(ε)  for nh≤¯t.

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 [11], p. 560, the Hamiltonian equations of motion take the form

 ˙p0 = −∇q0(12pT0M0(q0)−1p0+U(q0,0)) ˙q0 = M0(q0)−1p0+g0(p,q) (23) (˙p1˙q1) = 1ε(0−Ω(q0)Ω(q0)0)(p1q1)+(f1(p,q)g1(p,q))

with functions , and , . Moreover we have (omitting the arguments in , which are all )

 f1 = −ε1/2c−Lp1+ε−1/2a(p1,p1)−ε1/2∇1ˇU(q0,0)+O(ε3/2) g1 = LTq1+ε−1/2b(p1,q1)+O(ε3/2) (24)

where is an matrix and the functions and are bilinear.

We diagonalize

 Γ∗(0−Ω(q0)Ω(q0)0)Γ=i(Ω(q0)00−Ω(q0))=:iΛ(q0),

with

 Γ=1√2(II−iIiI),

and introduce the diagonal phase matrix by

 Φ(t) =∫t0Λ(q0(s))ds.

Following [11], p. 561, we transform the oscillatory part of the solution to adiabatic variables

 η(t) =ε−1/2exp(−iεΦ(t))Γ∗(p1(t)q1(t)). (25)

We further introduce the matrices and as

 (P1Q1) =Γexp(iεΦ),

so that and . In adiabatic variables, the differential equation for the oscillatory part becomes

 ˙η = exp(−iεΦ)W(p0,q0)exp(iεΦ)η + exp(−iεΦ)Γ∗(a(P1η,P1η;p0,q0)b(P1η,Q1η;p0,q0)) − P∗1(c(p0,q0)+∇1ˇU(q0,0))+O(ε)

with

 W =−12(L−LTL+LTL+LTL−LT).

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

 ηj+1 =ηj+∫tj+1tj˙ηj(s)ds,

where solves

 ˙η=exp(−iεΦj)W(pj0,qj0)exp(iεΦj)η+exp(−iεΦj)Γ∗(a(Pj1η,Pj1η;pj0,qj0)b(Pj1η,Qj1η;pj0,qj0))+(Pj1)∗c(pj0,qj0)+O(ε) (27)

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 [11]. Therefore, we just consider the first term of the righthand side in (27). For partial integration gives

 n−1∑j=0∫tj+1tjexp(−iε(Φjl(s)−Φjk(s)))wlk(pj0(s),qj0(s))ηjk(s)ds = iεn−1∑j=0exp(−iε(Φjl(s)−Φjk(s)))wlk(pj0(s),qj0(s))ηjk(s)ωjl(q0(s))−ωjk(q0(s))∣∣tj+1tj −iεn−1∑j=0∫tj+1tjexp(−iε(Φjl(s)−Φjk(s)))ddswlk(pj0(s),qj0(s))ηjk(s)ωjl(q0(s))−ωjk(q0(s))ds,

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

 ηn=η0+O(ε), (28)

and since , the result follows.

We are now in the situation to prove an error bound.

###### Theorem 4.2

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

 xn−X(tn) =O(h2)+O(ε) P(xn)yn−Y(tn) =O(h2)+O(ε).

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

 ˙p0=−∇q0(12pT0M0(q0)p0+U(q0,0))−m∑k=1Ik(0)∇q0ωk(q0),˙q0=M0(q0)−1p0. (29)

More precisely, if we write the Störmer–Verlet scheme for (29) in one-step form as

 (pn+1,0qn+1,0) =Ψh(pn,0,qn,0),

then the slow components of the mollified impulse method for (23) fulfill

 (pn+1,0qn+1,0) =Ψh(pn,0,qn,0)+dn

with a local error , because

 ∇q0(12εpT1Ω(q0)p1+12εqT1Ω(q0)q1)=m∑k=1Ik∇q0ωk(q0)

and by Theorem 4.1. Application of the discrete Gronwall Lemma gives the desired result for the slow components in the variables : for ,

 qn,0=q0(tn)+O(h2)+O(ε), pn,0=p0(tn)+O(h2)+O(ε).

For the fast variables we have, using (28),

 (p1(tn)q1(tn)) =ε1/2Γexp(iεΦ(tn))η(tn) =ε1/2Γexp(iεΦ(tn))[η(0)+O(ε)] =ε1/2Γexp(iεΦ(tn))[ηn+O(ε)] =Γexp(iε[Φ(tn)−Φn(tn)])Γ∗(pn,1qn,1)+O(ε3/2),

which shows

 qn,1=O(ε1/2), pn,1=O(ε1/2),

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

 π(x)=χ(q0,0)for x=χ(q) with q=(q0,q1).

The following result is essential in relating the various methods.

###### Lemma 5.1

For the mollifier we have, under the bounded-energy condition ,

 π(x) =x+O(ε), π′(x)T =P(x)+O(ε),

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

 x=ξ(q0,ε1/2q1)=ξ(q0,0)+O(ε)=π(x)+O(ε).

The proof of the result for is then obtained from the identity

 π′(X)T=P(X) for X with g(X)=0

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

 ˇP(q)=χ′(q)TP(x)χ′(q)−T.

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

 ˇP(q)=(1000)=ˇπ′(q)T,

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

 (p+n,0p+n,1) =(pn,0pn,1)−h2(∇q0ˇU(qn,0,0)0)+(O(hε)O(hε3/2)).

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.

###### Theorem 5.2

Consider the mollified impulse method or the projected impulse method of Section 2. Assume the energy bound (6). Furthermore, suppose that the frequencies remain separated and non-resonant (see conditions (11)-(12)) for . Then, the method approximately preserves the actions:

 In