Numerical integrators for motion under a strong constraining force
NNumerical integrators formotion under a strong constraining force
Christian Lubich † and Daniel Weiss ‡ October 24, 2018
Abstract.
This paper deals with the numerical integration of Hamiltonian systems in whicha stiff anharmonic potential causes highly oscillatory solution behavior with solution-dependentfrequencies. The impulse method, which uses micro- and macro-steps for the integration of fastand slow parts, respectively, does not work satisfactorily on such problems. Here it is shownthat variants of the impulse method with suitable projection preserve the actions as adiabaticinvariants and yield accurate approximations, with macro-stepsizes that are not restricted by thestiffness parameter.
Keywords.
Oscillatory Hamiltonian systems, varying high frequencies, impulse method,mollified impulse method, projected impulse method.
AMS subject classifications.
We are interested in the efficient numerical integration of Hamiltonian systems in which a stiffanharmonic potential causes highly oscillatory solution behavior with state-dependent slowlyvarying high frequencies.
We consider Hamiltonians as studied, in varying degrees of generality and with different analyticaltechniques, by Rubin & Ungar [16], Takens [19], Bornemann [2], Lorenz [14] and Hairer, Lubich& Wanner [11, Section XIV.3]: H ( x, y ) = y T M ( x ) − y + U ( x ) + 1 ε V ( x ) , < ε (cid:28) , (1) † Mathematisches Institut, Univ. T¨ubingen, Auf der Morgenstelle, D-72076 T¨ubingen, Germany. E-mail: [email protected] ‡ Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, Kaiserstr. 89-93, D-76133 Karlsruhe, Germany. E-mail: [email protected] a r X i v : . [ m a t h . NA ] J u l epending on positions x ∈ R n and momenta y ∈ R n . The mass matrix M ( x ) is symmetricand positive definite and depends smoothly on x . The slow potential U is smooth, and the stiffpotential ε V with a smooth function V : D ⊂ R n → R attains its minimum value 0 on a d -dimensional manifold V = { x ∈ D : V ( x ) = min V = 0 } . (2)We assume that the potential V is strongly convex along directions non-tangential to V . Moreprecisely, there exists α > x ∈ V the Hessian ∇ V ( x ) satisfies v T ∇ V ( x ) v ≥ α · v T M ( x ) v (3)for all vectors v in the M ( x )-orthogonal complement of the tangent space T x V .Furthermore, we assume that a constraint function g : D → R m , with m = n − d , is knownsuch that V = { x ∈ D : g ( x ) = 0 } (4)and the derivative matrix G ( x ) = g (cid:48) ( x ) is of full rank on V .The corresponding system of Hamiltonian differential equations reads˙ x = M ( x ) − y ˙ y = −∇ x (cid:16) y T M ( x ) − y (cid:17) − ∇ U ( x ) − ε ∇ V ( x ) . (5) Example 1
A simple, yet nontrivial model example is the stiff spring doublependulum. The Hamiltonian reads H ( x, y ) = y T y + U ( x ) + 1 ε V ( x ) , where y T y is the kinetic energy, U ( x ) = x + x , and ε V ( x ) = (cid:104) α ε (cid:105) ( (cid:107) x (cid:107) − l ) + (cid:104) α ε (cid:105) ( (cid:107) x − x (cid:107) − l ) is the stiff potential depending on the small parameter ε . Theparameters l i denote the lengths of the springs, α i /ε are the largespring constants. m ( x , x ) m ( x , x ) Example 1 helps to fix ideas on a simple toy model. Obviously it extends to chains of stiffsprings, which describe the dynamics of chains of atoms in a molecule with almost rigid bonds,cf., e.g., [13].
It has been known since Rubin & Ungar [16] that the motion of the system in the limit ε → V (the rigid double pendulum in2he above example) for general initial values ( x , y ) that have an energy bounded independentlyof ε , H ( x , y ) ≤ Const . (6)Note that the set of admissible initial values ( x , y ) satisfying (6) depends on ε . The effectiveconstrained Hamiltonian has a correction potential W , H eff ( X, Y ) = Y T M ( X ) − Y + U ( X ) + W ( I, X ) , g ( X ) . (7)The correction potential W ( I, X ) = (cid:80) mk =1 I k ω k ( X ) depends on the m frequencies ω k ( X ), i.e.,the square roots of the nonzero generalized eigenvalues of the pencil λM ( X ) − ∇ V ( X ), andon m parameters I = ( I , . . . , I m ), known as the actions , which are determined by the initialvalues ( x , y ) of (5). The actions vanish for consistent initial data that satisfy g ( x ) = 0 and G ( x ) M ( x ) − y = 0.As is outlined in [3] and will be recapitulated in Section 3, the effective Hamiltonian can befound by transforming the system to separate slow and fast variables as in [14] and [11, SectionXIV.3], and transforming the obtained slow system ( ε = 0) back via the effective dynamics ofthe fast variables.The effective constrained Hamiltonian system is then given by˙ X = M ( X ) − Y ˙ Y = −∇ X (cid:16) Y T M ( X ) − Y (cid:17) − ∇ U ( X ) − ∇ X W ( I, X ) − G ( X ) T Λ0 = g ( X ) , (8)with Lagrange multipliers Λ( t ) ∈ R m . This differs from the usual constrained equations of motionthrough the correction force F ( I, X ) = −∇ X W ( I, X ).To initial values ( x , y ) of system (5) with bounded energy (6), we associate consistent initialvalues ( X , Y ) for the effective system (8). These are chosen by projecting M -orthogonally ontothe manifold of consistent values: X = x + M ( x ) − G ( x ) T λ, g ( X ) ,Y = y + G ( X ) T µ, G ( X ) M ( X ) − Y . (9)With the projection P ( x ) = I − Q ( x ) , Q ( x ) = [ G T ( GM − G T ) − GM − ]( x ) , (10)the second equation can be rewritten as Y = P ( X ) y , and along the solution of (8) we note Y ( t ) = P ( X ( t )) Y ( t ).The effective Hamiltonian system describes the limit dynamics on the constraint manifold aslong as the solution-dependent frequencies ω k ( X ( t )) remain separated and are non-resonant: forsome δ > | ω j ( X ( t )) − ω k ( X ( t )) | ≥ δ for j (cid:54) = k (11) | ω j ( X ( t )) ± ω k ( X ( t )) ± ω l ( X ( t )) | ≥ δ for all j, k, l . (12)3f these conditions are satisfied for t ≤ t , 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 O -notation depend on δ and deteriorate as δ →
0; see [19, 11, 3]. Inthe case of a single frequency ( m = 1), 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 O ( ε ). Thenormal 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 conditionsof this type are needed for the above approximation result for the effective dynamics. Usingthe techniques of [11, Chap. XIV] it can be shown that the order of this approximation is still O ( ε / ( m +1) ) if ω j ± ω k ± ω l have zeros of multiplicity m . However, the separation cannot beomitted. If the distance of two frequencies becomes smaller than √ ε , then the slow motion candepend very sensitively on the initial values, and it is no longer approximated by the dynamicsof the effective Hamiltonian system; see Takens [19]. The indeterminacy of the slow motion inthe case of non-separated frequencies is termed Takens chaos in [2].
The objective of this paper is to devise and analyze a two-scale integrator for the highly oscillatoryHamiltonian system (5), such that for a macro-stepsize h that is not restricted by ε , the methodyields an O ( h ) + O ( ε ) error in the positions x ( t ) and the projected momenta P ( x ( t )) y ( t ) overtime intervals t = O (1).This paper is part of the vast literature on the numerical solution of highly oscillatory differ-ential equations; see, e.g., the reviews [5, 15]. Recent work on the numerical integration of highlyoscillatory mechanical systems includes [1, 4, 17, 18, 20].While much work has been done on systems with constant high frequencies, the numericalanalysis of the present case of solution-dependent high frequencies or even just the case of explic-itly time-dependent high frequencies is scarce; see [11, Chapter XIV]. An important aspect here isto 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 effectivesystem (8) over the time scale t = O (1) without, however, explicitly evaluating the correctionforce F ( I, X ) = −∇ X W ( I, X ). This additional force is, in general, directly accessible only viaa 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 anumerical method that is able to monitor the possible breakdown of the validity of the effectiveequation due to the loss of adiabatic invariance of the actions in cases where frequencies comeclose or become resonant.Heterogeneous multiscale methods (HMM) [6, 8, 7] have been developed for the very purposeto handle situations where the underlying effective dynamics is not known. In Brumm & Weiss [3]4n HMM-approach for highly oscillatory mechanical systems with solution-dependent frequenciesis analyzed. This approach shows, however, major drawbacks because of difficulties in initializingthe micro-simulation.In this article, we follow the alternative idea of the impulse method where the Hamiltonian issplit into the slow potential U 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, thismust be complemented with a suitable projection to lead to a method with satisfactory errorbehavior.We proceed as follows: In Section 2 we formulate the impulse method, a mollified impulsemethod, and a novel projected impulse method for highly oscillatory mechanical systems withsolution-dependent frequencies. We state the main convergence theorem and show results ofnumerical experiments that highlight different behavior of the various methods. In Section 3 wetransform the system, following [14] and [11, Section XIV.3] to variables that are appropriate forthe further analysis. Moreover, a further mollified impulse method with a projection mollifierin the transformed variables is introduced, which is computationally impractical but serves as atheoretical reference method for the error analysis. This method is studied in Section 4. Usingthe obtained results, the analysis of the mollified and projected impulse methods of Section 2 isdone in Section 5. The impulse method was introduced in the context of the numerical treatment of moleculardynamics (Grubm¨uller, Heller, Windemuth & Schulten [10], Tuckerman, Berne & Martyna [21]).A mathematical study of this method is given by Garc´ıa-Archilla, Sanz-Serna & Skeel [9]. Theidea is to split the Hamiltonian H ( x, y ) = H fast ( x, y ) + U ( x )and to approximate the exact flow ϕ Hh by the following symmetric decomposition: ϕ Hh ≈ ϕ slow h/ ◦ ϕ fast h ◦ ϕ slow h/ . Since the flow of the slow part can be trivially solved, one step is equivalent to1. kick: y + n = y n − h/ · ∇ U ( x n ),2. oscillate: solve system (5) with U = 0 and initial values ( y + n , x n ) over a time step h toobtain ( y − n +1 , x n +1 ),3. another kick: y n +1 = y − n +1 − h/ · ∇ U ( x n +1 ).Step 2. is solved approximately using, e.g., the St¨ormer–Verlet method with micro-stepsizes,or alternatively using a large-timestep method in suitably transformed variables (an adiabaticintegrator) as in [14].Compared to a direct numerical integration of the full system (5) with small stepsizes, thismethod saves many evaluations of the slow force −∇ U ( x ), which is often the computationallymost expensive part. 5 −2 −1 −3 −2 −1 stepsize m a x i m a l e rr o r i n po s i t i on s ε = 10 − ε = 10 − −2 −1 −3 −2 −1 stepsize m a x i m a l e rr o r : p r o j . m o m en t a ε = 10 − ε = 10 − Figure 1: Double logarithmic plots: stepsize versus maximal error (maximum norm, maximumover all discrete times) of the impulse method applied to the stiff spring double pendulum withinitial value (14). Left: Error in positions. Right: Error of the projected momenta.For our numerical experiments we consider the stiff spring double pendulum with initial values x (0) = ( √ . , −√ . , √ , ε ) T , y (0) = (0 , , , T (14)and the parameters α = α = 1, over the time interval 0 ≤ t ≤
10. In this situation thefrequencies remain well-separated.We observe unsatisfactory behavior of the impulse method in Figure 1. Here and in all fol-lowing figures the dash-dotted straight line has slope 2, corresponding to the desired h errorbehavior. We used the St¨ormer–Verlet method with very small stepsize ( ε/ ε/
100 for the following methods) for the micro-integration in order to avoid any sig-nificant influence on the overall error. Throughout all computations, as a reference, the effectivesystem (8) is approximated in transformed variables (see (29)) by the St¨ormer–Verlet methodwith small stepsize, the results being translated back into cartesian coordinates.
Garc´ıa-Archilla, Sanz-Serna & Skeel [9] and Izaguirre, Reich & Skeel [13] improve the impulsemethod by replacing the slow potential U ( x ) by a mollified potential ¯ U ( x ) = U ( α ( x )), where α ( x ) is an averaged or suitably projected value of x . The mollified force then reads −∇ U ( x ) = − α (cid:48) ( x n ) T ∇ U ( α ( x n )) . The mollification considered in the present paper is given by the M ( x )-orthogonal projectiononto the configuration manifold { X : g ( X ) = 0 } , i.e., α ( x ) = X with X = x + M ( x ) − G ( x ) T λ, g ( X ) . (15)Using the same initial value (14) as in the case of the impulse method, we observe better con-vergence behavior of the positions and the projected momenta, see Figure 2.6 −2 −1 −5 −4 −3 −2 stepsize m a x i m a l e rr o r i n po s i t i on s ε = 10 − ε = 10 − ε = 10 − −2 −1 −5 −4 −3 −2 stepsize m a x i m a l e rr o r : p r o j . m o m en t a ε = 10 − ε = 10 − ε = 10 − Figure 2: Double logarithmic plots: stepsize versus maximal error of the mollified impulse methodapplied to the stiff spring double pendulum with initial value (14). Left: Error in positions. Right:Error of the projected momenta.As it turns out in the analysis, the unsatisfactory behavior of the impulse method is due to M ( x )-orthogonal components of the slow forces ∇ U ( x ). The mollification reduces those M ( x )-orthogonal components. Indeed, we observe the following. Lemma 2.1
Under the bounded-energy condition V ( x ) = O ( ε ) , the mollifier α ( x ) of (15) satisfies α ( x ) = x + O ( ε ) ,α (cid:48) ( x ) T = P ( x ) + O ( ε ) , where the projection P ( x ) is defined in (10) .Proof. The condition V ( x ) = O ( ε ) is equivalent to g ( x ) = O ( ε ). Noting that ( GM − G T )( x )is invertible in view of the full rank of G , the implicit function theorem then yields λ = O ( ε )such that g ( x + M ( x ) − G ( x ) T λ ) = 0, and hence α ( x ) = x + O ( ε ). Differentiating both equationsin (15) yields α (cid:48) ( x ) = I + M − ( x ) G ( x ) T λ (cid:48) ( x ) + O ( ε )0 = G ( α ( x )) α (cid:48) ( x ) . Inserting the first into the second equation permits us to compute λ (cid:48) ( x ) = − ( GM − G T ) − G ( x ) + O ( ε ) . Reinserting this expression into the first equation yields the stated result on recalling the defini-tion of P ( x ). The preceding lemma motivates us to simplify the method by projecting the slow force:7. kick: y + n = y n − h/ · P ( x n ) ∇ U ( x n ),2. oscillate: solve system (5) with U = 0 and initial values ( y + n , x n ) over a time step h obtaining( y − n +1 , x n +1 ),3. kick: y n +1 = y − n +1 − h/ · P ( x n +1 ) ∇ U ( x n +1 ).Using this new simplified scheme, we observe convergence behavior as in the case of themollified impulse method, see Figure 3.We have sacrificed the symplecticity of the method which is not of main interest here, buthave nevertheless maintained the time-reversal symmetry. −2 −1 −5 −4 −3 −2 stepsize m a x i m a l e rr o r i n po s i t i on s ε = 10 − ε = 10 − ε = 10 − −2 −1 −5 −4 −3 −2 −1 stepsize m a x i m a l e rr o r : p r o j . m o m en t a ε = 10 − ε = 10 − ε = 10 − Figure 3: Double logarithmic plots: stepsize versus maximal error of the projected impulsemethod applied to the stiff spring double pendulum with initial value (14). Left: Error inpositions. Right: Error of the projected momenta.
The idea of a projection as a mollification is proposed in [13]. There, the use of this idea is shownexperimentally but no analysis is given. On the other hand, in [14, 11] the adiabatic nature ofthe systems of interest is revealed by applying a series of canonical transformations. Combiningthe different ideas and techniques, we are now able to formulate and prove the result about theglobal error of the mollified impulse method with the projection mollifier. Additionally, we provethe same result for the computationally simpler projected impulse method. The proof is basedon a further, different mollification introduced in Section 3.
Theorem 2.2
Let the initial values satisfy the energy bound (6) and assume that the frequenciesremain separated and non-resonant (see conditions (11) - (12) ) along the solution of (5) for ≤ t ≤ t . Then, the errors of the mollified impulse method and of the projected impulse method after n steps with stepsize h satisfy x n − X ( t n ) = O ( h ) + O ( ε ) P ( x n ) y n − Y ( t n ) = O ( h ) + O ( ε ) (16)8 here ( X ( t ) , Y ( t )) is the solution of the effective Hamiltonian system (8) with initial valuesdefined by (9) . The constants symbolized by O do not depend on ε , h and n with nh ≤ t . Combined with (13) this also yields the error bounds with respect to the solution ( x ( t ) , y ( t ))of the highly oscillatory problem x n − x ( t n ) = O ( h ) + O ( ε ) P ( x n )( y n − y ( t n )) = O ( h ) + O ( ε ) . (17)We note, however, that y n − Y ( t n ) = O (1) and y n − y ( t n ) = O (1). Moreover, the method doesnot converge to the solution ( x ( t ) , y ( t )) of the highly oscillatory system for a fixed ε as h → h > ε .Note that in Theorem 2.2 there is no restriction of the step size h by the small parameter ε .Theorem 2.2 explains the error behavior observed in Figures 2 and 3. Under conditions (2)-(4), [11, Section XIV.3] and [14] show that there exists a canonical changeof coordinates ( x, y ) = ψ ( q, p ) of the separated form x = χ ( q ), y = χ (cid:48) ( q ) − T p , which transformsthe Hamiltonian (1) into the form H ( q, p ) = p T M ( q ) − p + 12 ε p T Ω( q ) p + 12 ε q T Ω( q ) q + (cid:18) p ε − / p (cid:19) T R ( q , ε / q ) (cid:18) p ε − / p (cid:19) + ˇ U ( q , ε / q ) , (18)where q = ( q , q ) ∈ R d × R m and p = ( p , p ) ∈ R d × R m and the appearing functions have alltheir partial derivatives bounded independently of ε and are as follows: • M ( q ) is a symmetric positive definite d × d matrix; • Ω( q ) is a diagonal m × m matrix with positive entries, the frequencies ω k ( q ); • R ( q , ε / q ) is a symmetric n × n matrix with R ( q ,
0) = 0; • ˇ U ( q , ε / q ) = U ( x ) for x = χ ( q ).The assumption (6) of bounded energy now becomes q = O ( ε / ) , p = O ( ε / ) . (19)We define the actions I k = 12 ε (cid:0) q ,k + p ,k (cid:1) , k = 1 , . . . , m. (20)Under the separation and non-resonance conditions (11)–(12), the actions are adiabatic invari-ants : they remain nearly constant along solutions of the Hamiltonian system with boundedenergy, I k ( t ) = I k (0) + O ( ε ) , (21)9ee [11], p. 562.As will become clear from our analysis, this is a key property that should be transfered tothe numerical method. However, if we express the impulse method in the transformed variables,then the kick step becomes (cid:18) p + n, p + n, (cid:19) = (cid:18) p n, p n, (cid:19) − h (cid:18) ∇ ˇ U ( q n, , ε / q n, ) ε / ∇ ˇ U ( q n, , ε / q n, ) (cid:19) , and we see that the actions are not approximately preserved. This is illustrated in Figure 4. Thenon-preservation of the actions is at the base of the disappointing numerical behavior observedin Figure 1. Figure 4: Non-preservation of the actions for the impulse method. −6 time I − I ( ) −4 time I − I ( ) Figure 5: Near-preservation of the actions for the projected impulse method.10n the other hand, if we use a mollified impulse method for which the modified potential ischosen, in the transformed variables, as U ( q ) = ˇ U ( q , , (22)then p + n, = p n, , and hence the actions are exactly preserved in the kick step. While this methodis not practical in that it would require performing the coordinate transformation from ( x, y ) to( q, p ), it gives much theoretical insight into the error propagation behavior. We will thereforestudy its error in the next section. Subsequently we will interpret the mollified and projectedimpulse methods of Section 2, which work in the original variables, as perturbations of thistheoretically interesting method.As a numerical illustration, in Figure 6 we use the transformed-variable mollified impulsemethod with the initial value of Section 2 for the stiff spring double pendulum. We observesimilar results as in Figures 2 and 3. −2 −1 −5 −4 −3 −2 step size m a x i m a l e rr o r i n po s i t i on s ε = 10 − ε = 10 − ε = 10 − −2 −1 −5 −4 −3 −2 step size m a x i m a l e rr o r : p r o j . m o m en t a ε = 10 − ε = 10 − ε = 10 − Figure 6: Double logarithmic plots: stepsize versus maximal error of the mollified impulse methodwith (22) applied to the stiff spring double pendulum with initial value (14). Left: Error inpositions. Right: Error of the projected momenta.
We will show almost-conservation of the actions along the numerical solution, I n = ( I n, , . . . , I n,m )with I n,k = ε ( q n, ,k + p n, ,k ). Theorem 4.1
Assume the energy bound (6) . Furthermore, assume that the frequencies of thetransformed-variable mollified impulse method with modified potential (22) remain separated andnon-resonant (see conditions (11) - (12) ) for ≤ t ≤ t . Then, this method approximately preservesthe actions: I n = I + O ( ε ) for nh ≤ t. The constant symbolized by O is independent of n , h , and ε .
11o prove this result, we first need to look in more detail into the differential equations in thetransformed variables. As is shown in [11], p. 560, the Hamiltonian equations of motion take theform ˙ p = − ∇ q (cid:16) p T M ( q ) − p + U ( q , (cid:17) − ∇ q (cid:16) ε p T Ω( q ) p + 12 ε q T Ω( q ) q (cid:17) + f ( p, q )˙ q = M ( q ) − p + g ( p, q ) (23) (cid:18) ˙ p ˙ q (cid:19) = 1 ε (cid:18) − Ω( q )Ω( q ) 0 (cid:19) (cid:18) p q (cid:19) + (cid:18) f ( p, q ) g ( p, q ) (cid:19) with functions f = O ( ε ), g = O ( ε ) and f = O ( ε / ), g = O ( ε / ). Moreover we have(omitting the arguments p , q in a, b, c, L , which are all O (1)) f = − ε / c − Lp + ε − / a ( p , p ) − ε / ∇ ˇ U ( q ,
0) + O ( ε / ) g = L T q + ε − / b ( p , q ) + O ( ε / ) (24)where L is an m × m matrix and the functions a and b are bilinear.We diagonalize Γ ∗ (cid:18) − Ω( q )Ω( q ) 0 (cid:19) Γ = i (cid:18) Ω( q ) 00 − Ω( q ) (cid:19) =: i Λ( q ) , with Γ = 1 √ (cid:18) I I − iI iI (cid:19) , and introduce the diagonal phase matrix Φ( t ) byΦ( t ) = (cid:90) t Λ( q ( s )) ds. Following [11], p. 561, we transform the oscillatory part of the solution to adiabatic variables η ( t ) = ε − / exp (cid:18) − iε Φ( t ) (cid:19) Γ ∗ (cid:18) p ( t ) q ( t ) (cid:19) . (25)We further introduce the m × m matrices P ( t ) and Q ( t ) as (cid:18) P Q (cid:19) = Γ exp (cid:18) iε Φ (cid:19) , so that p = ε / P η and q = ε / Q η . In adiabatic variables, the differential equation for theoscillatory part becomes ˙ η = exp (cid:16) − iε Φ (cid:17) W ( p , q ) exp (cid:16) iε Φ (cid:17) η + exp (cid:16) − iε Φ (cid:17) Γ ∗ (cid:18) a ( P η, P η ; p , q ) b ( P η, Q η ; p , q ) (cid:19) (26) − P ∗ (cid:16) c ( p , q ) + ∇ ˇ U ( q , (cid:17) + O ( ε )12ith W = − (cid:18) L − L T L + L T L + L T L − L T (cid:19) . The functions
L, a, b, c are those appearing in the remainder terms f and g in (24). Proof. (of Theorem 4.1) We rewrite the mollified impulse method in adiabatic variables andnote that the kick steps do not change the adiabatic variables: in the j th time step, η + j = η j and η j +1 = η − j +1 . We thus obtain η j +1 = η j + (cid:90) t j +1 t j ˙ η j ( s ) ds, where η j ( t ) solves˙ η = exp (cid:18) − iε Φ j (cid:19) W ( p j , q j ) exp (cid:18) iε Φ j (cid:19) η + exp (cid:18) − iε Φ j (cid:19) Γ ∗ (cid:18) a ( P j η, P j η ; p j , q j ) b ( P j η, Q j η ; p j , q j ) (cid:19) + ( P j ) ∗ c ( p j , q j ) + O ( ε ) (27)with initial value η j on the interval [ t j , t j +1 ]. All terms with superscript j are defined with respectto the solution of the oscillation step of the mollified impulse method on [ t j , t j +1 ].In the remaining part of the proof we show (cid:80) n − j =0 (cid:82) t j +1 t j ˙ η j ( s ) ds = O ( ε ). The techniques aremore or less the same as for the exact solution presented in [11]. Therefore, we just consider thefirst term of the righthand side in (27). For l (cid:54) = k partial integration gives n − (cid:88) j =0 (cid:90) t j +1 t j exp (cid:18) − iε (Φ jl ( s ) − Φ jk ( s )) (cid:19) w lk ( p j ( s ) , q j ( s )) η jk ( s ) ds = iε n − (cid:88) j =0 exp (cid:18) − iε (Φ jl ( s ) − Φ jk ( s )) (cid:19) w lk ( p j ( s ) , q j ( s )) η jk ( s ) ω jl ( q ( s )) − ω jk ( q ( s )) (cid:12)(cid:12)(cid:12) t j +1 t j − iε n − (cid:88) j =0 (cid:90) t j +1 t j exp (cid:18) − iε (Φ jl ( s ) − Φ jk ( s )) (cid:19) dds w lk ( p j ( s ) , q j ( s )) η jk ( s ) ω jl ( q ( s )) − ω jk ( q ( s )) ds, where the latter term is of size O ( ε ) in the case of separated frequencies. Taking into account the O ( h )-jumps from p j ( t j +1 ) to p j +10 ( t j +1 ) and noting that η j ( t j +1 ) = η j +1 ( t j +1 ) and Φ j ( t j +1 ) =Φ j +1 ( t j +1 ), we prove the same bound for the first term. We have thus shown η n = η + O ( ε ) , (28)and since I n,k = | η n,k | , the result follows.We are now in the situation to prove an error bound.13 heorem 4.2 Assume the energy bound (6) . Furthermore, assume that the frequencies remainseparated and non-resonant (see conditions (11) - (12) ) on a fixed time interval ≤ t ≤ t . Then,the error of the transformed-variable mollified impulse method of Section 3 after n steps withstepsize h satisfies x n − X ( t n ) = O ( h ) + O ( ε ) P ( x n ) y n − Y ( t n ) = O ( h ) + O ( ε ) . The constants symbolized by O do not depend on ε , h and n with nh ≤ t . We note, however, that the normal components of the momenta are not approximated cor-rectly: we only have y n − y ( t n ) = O (1). Proof.
We consider the method in the slow components p , q as a perturbed variant of theSt¨ormer–Verlet scheme applied to the slow system˙ p = −∇ q (cid:0) p T M ( q ) p + U ( q , (cid:1) − m (cid:88) k =1 I k (0) ∇ q ω k ( q ) , ˙ q = M ( q ) − p . (29)More precisely, if we write the St¨ormer–Verlet scheme for (29) in one-step form as (cid:18) p n +1 , q n +1 , (cid:19) = Ψ h ( p n, , q n, ) , then the slow components of the mollified impulse method for (23) fulfill (cid:18) p n +1 , q n +1 , (cid:19) = Ψ h ( p n, , q n, ) + d n with a local error d n = O ( h ) + O ( hε ), because ∇ q (cid:16) ε p T Ω( q ) p + 12 ε q T Ω( q ) q (cid:17) = m (cid:88) k =1 I k ∇ q ω k ( q )and I n,k = I ,k + O ( ε ) by Theorem 4.1. Application of the discrete Gronwall Lemma gives thedesired result for the slow components in the variables ( q, p ): for nh ≤ t , q n, = q ( t n ) + O ( h ) + O ( ε ) , p n, = p ( t n ) + O ( h ) + O ( ε ) . For the fast variables we have, using (28), (cid:18) p ( t n ) q ( t n ) (cid:19) = ε / Γ exp (cid:18) iε Φ( t n ) (cid:19) η ( t n )= ε / Γ exp (cid:18) iε Φ( t n ) (cid:19) [ η (0) + O ( ε )]= ε / Γ exp (cid:18) iε Φ( t n ) (cid:19) [ η n + O ( ε )]= Γ exp (cid:18) iε [Φ( t n ) − Φ n ( t n )] (cid:19) Γ ∗ (cid:18) p n, q n, (cid:19) + O ( ε / ) , q n, = O ( ε / ) , p n, = O ( ε / ) , but in view of the phase difference Φ( t n ) − Φ n ( t n ) = O ( h ) + O ( ε ) this does not yield an approx-imation estimate. Transforming back to the coordinates ( x, y ) and considering the rescaling andthe Lipschitz-continuity of the transformations then gives the result. In order to analyze the mollified and projected impulse methods of Section 2, we have to derivean appropriate expression of the kick-step in the transformed variables ( q, p ). We show thatboth methods are O ( ε )-perturbations of the transformed-variable mollified impulse method ofSection 3, which uses the modified potential (22), that is, in the original variables, the modifiedpotential U ( π ( x )) with π ( x ) = χ ( q ,
0) for x = χ ( q ) with q = ( q , q ) . The following result is essential in relating the various methods.
Lemma 5.1
For the mollifier π ( x ) we have, under the bounded-energy condition V ( x ) = O ( ε ) , π ( x ) = x + O ( ε ) ,π (cid:48) ( x ) T = P ( x ) + O ( ε ) , where the projection P ( x ) is defined in (10) .Proof. As the construction in [11, Chapter XIV.3] shows, the transformation x = χ ( q ) iscomposed as χ = ξ ◦ φ ε with an ε -independent transformation ξ and a rescaling φ ε ( q , q ) =( q , ε / q ). Since for x = χ ( q ) the bounded-energy condition V ( x ) = O ( ε ) is equivalent to q = O ( ε / ), we obtain x = ξ ( q , ε / q ) = ξ ( q ,
0) + O ( ε ) = π ( x ) + O ( ε ) . The proof of the result for π (cid:48) ( x ) T is then obtained from the identity π (cid:48) ( X ) T = P ( X ) for X with g ( X ) = 0used for X = π ( x ). This identity is obtained from the transformation laws as follows: Un-der a change of coordinates x = χ ( q ), the constraint function changes to ˇ g ( q ) = g ( χ ( q )) andits derivative ˇ G = ˇ g (cid:48) to ˇ G ( q ) = G ( x ) χ (cid:48) ( q ). The inverse mass matrix changes to ˇ M ( q ) − = χ (cid:48) ( q ) − M ( x ) − χ (cid:48) ( q ) − T . Consequently, the projection ˇ P ( q ) = I − [ ˇ G T ( ˇ G ˇ M − ˇ G T ) − ˇ G ˇ M − ]( q )transforms as ˇ P ( q ) = χ (cid:48) ( q ) T P ( x ) χ (cid:48) ( q ) − T . For ˇ π ( q ) = χ − ( π ( χ ( q ))), the transposed derivative transforms in the same way for x = χ ( q )with x = π ( x ): ˇ π (cid:48) ( q ) T = χ (cid:48) ( q ) T π (cid:48) ( x ) T χ (cid:48) ( q ) − T . q we have ˇ g ( q , q ) = q and a block diagonal mass matrix ˇ M ( q ), and onthe other hand ˇ π ( q , q ) = ( q , P ( q ) = (cid:18) (cid:19) = ˇ π (cid:48) ( q ) T , and hence the result follows.Using Lemma 5.1 and the corresponding result Lemma 2.1 for the mollified impulse methodof Section 2, we find for the kick step of the mollified - and projected impulse methods expressedin the variables ( q, p ) of Section 3 (cid:18) p + n, p + n, (cid:19) = (cid:18) p n, p n, (cid:19) − h (cid:18) ∇ q ˇ U ( q n, , (cid:19) + (cid:18) O ( hε ) O ( hε / ) (cid:19) . Here, the additional factor ε / is due to the rescaling of the fast positions and momenta in thetransformation. For the actions, we then obtain the estimate I + n = I n + O ( hε ) 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 Sec-tion 2. Assume the energy bound (6) . Furthermore, suppose that the frequencies remain separatedand non-resonant (see conditions (11) - (12) ) for ≤ t ≤ t . Then, the method approximately pre-serves the actions: I n = I + O ( ε ) for nh ≤ t. The constant symbolized by O is independent of n , h , and ε . Therefore, Theorem 2.2 holds true, which can be proven in exactly the same way as Theorem4.2.
We devised and analyzed numerical integrators that capture the effective dynamics of stiff me-chanical systems where a strong constraining force leads to highly oscillatory solution behaviorwith state-dependent frequencies. The integrators are projected variants of the impulse method,which splits the Hamiltonian into the slow potential and the fast part and integrates the latterwith micro-steps. A key aspect is the preservation of the actions as adiabatic invariants in thenumerical method. It is shown that this can be ensured by approximately projecting out the nor-mal components of the slow force. Numerical experiments and theoretical error bounds illustratethe favorable properties of the proposed methods.
References [1]
G. Ariel., J.M. Sanz-Serna, and R. Tsai , A multiscale technique for finding slowmanifolds of stiff mechanical systems , Multiscale Model. Simul., 10 (2012), pp. 1180–1203.162]
F. Bornemann , Homogenization in Time of Singularly Perturbed Mechanical Systems ,vol. 1687 of Lecture Notes in Mathematics, Springer, 1998.[3]
B. Brumm and D. Weiss , Heterogeneous multiscale methods for highly-oscillatory mechan-ical systems with solution-dependent frequencies , IMA J. Numer. Anal., 34 (2014), pp. 55–82.[4]
M.P. Calvo and J.M. Sanz-Serna , Heterogeneous multiscale methods for mechanicalsystems with vibrations , SIAM J. Sci. Comput., 32 (2010), pp. 2029–2046.[5]
D. Cohen, T. Jahnke, K. Lorenz, and C. Lubich , Numerical integrators for highly os-cillatory Hamiltonian systems: a review , in Analysis, Modeling and Simulation of MultiscaleProblems, A. Mielke, ed., Springer, Berlin, 2006, pp. 553–576.[6]
W. E , Analysis of the heterogeneous multiscale method for ordinary differential equations ,Comm. Math. Sci., 1 (2003), pp. 423–436.[7]
W. E, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden , Heterogeneous multiscalemethods: a review , Commun. Comput. Phys., 2 (2007), pp. 367–450.[8]
B. Engquist and R. Tsai , Heterogeneous multiscale methods for stiff ordinary differentialequations , Math. Comput., 74 (2005), pp. 1707–1742.[9]
B. Garc´ıa-Archilla, J.M. Sanz-Serna, and R.D. Skeel , Long-time-step methods foroscillatory differential equations , SIAM J. Sci. Comput., 20 (1998), pp. 930–963.[10]
H. Grubm¨uller, H. Heller, A. Windemuth, and K. Schulten , Generalized Verletalgorithm for efficient molecular dynamics simulations with long-range interactions , Mol.Sim., 6 (1991), pp. 121–142.[11]
E. Hairer, C. Lubich, and G. Wanner , Geometric Numerical Integration. Structure–Preserving Algorithms for Ordinary Differential Equations , Springer, Berlin, 2006.[12]
J. Henrard , The adiabatic invariant in classical mechanics , in Expositions in DynamicalSystems, vol. 2 of Dynamics reported. New series, Springer, Berlin, 1993, pp. 117–235.[13]
J.A. Izaguirre, S. Reich, and R.D. Skeel , Longer time steps for molecular dynamics ,J. Chem. Phys., 110 (1999), pp. 9853–9864.[14]
K. Lorenz , Adiabatische Integratoren f¨ur hochoszillatorische Hamilton-Systeme , PhD the-sis, Universit¨at T¨ubingen, 2006. http://nbn-resolving.de/urn:nbn:de:bsz:21-opus-24060.[15]
L.R. Petzold, L.O. Jay, and J. Yen , Numerical solution of highly oscillatory ordinarydifferential equations , Acta Numerica, 7 (1997), pp. 437–483.[16]
H. Rubin and P. Ungar , Motion under a strong constraining force , Comm. Pure Appl.Math., 10 (1957), pp. 65–87.[17]
J.M. Sanz-Serna , Mollified impulse methods for highly oscillatory differential equations ,SIAM J. Numer. Anal., 46 (2008), pp. 1040–1059.1718]
R. Sharp, R. Tsai, and B. Engquist , Multiple time scale numerical methods for theinverted pendulum problem , in Multiscale Methods in Science and Engineering, vol. 44 ofLecture Notes in Computational Science and Engineering, Springer, 2005, pp. 241–261.[19]
F. Takens , Motion under the influence of a strong constraining force , in Global theoryof dynamical systems, Proc. Int. Conf., Evanston/Ill. 1979, vol. 819 of Lecture Notes inMathematics, Springer, 1980, pp. 425–445.[20]
M. Tao, H. Owhadi, and J.E. Marsden , Nonintrusive and structure preserving multi-scale integration of stiff ODEs, SDEs, and Hamiltonian systems with hidden slow dynamicsvia flow averaging , Multiscale Model. Simul., 8 (2010), pp. 1269–1324.[21]