Symplectic Euler

SpUtting our Hamiltonian using H = Tip) = p M p/2 and H2 = U q) gives the symplectic Euler method. The BCH expansion gives the perturbed Hamiltonian 

We can think of the velocity Verlet method as being defined by a splitting into three parts  

This is a sequence of two types of operations it consists of a kick exhibited by a jump in the momentum, linear drift with the resulting momentum, followed by a final kick. The symmetry of the method is one of its important features. Switching the order of the operations, i.e. drift-kick-drift, gives the position Verlet method. 

In general for the Verlet method applied to Hamiltonians of the form H(q,p) = Tip) + U(q), the BCH lemma gives 

The extra power of h in the leading perturbation reflects the fact that this method is 2nd order accurate, rather than first order, as Symplectic Euler. Note that only even order terms (h, . ..) appear in the perturbative expansion this is a consequence 

---

The method (2.18)-(2.19) is called the Symplectic Euler method. Its adjoint method has a similar structure ... 

Example 2.5 (Verlet is Symplectic) For the Symplectic Euler method /, and its adjoint method consider the composition... 

Fig. 3.2 We compare the symplectic Euler, Verlet and Yoshida schemes in application to a Lennard-Jones oscillator. The plot shows the absolute deviation in the computed Hamiltonian (left) as a function of time, for each scheme at a fixed timestep h = 0.005. Moreover we simulate the system using different stepsizes and compute for each stepsize the mtiximum deviation in the Hamiltonian (right), comparing the results with guide lines associated to various powers of the step size. The scheme used as the base method for the Yoshida composition methods is denoted in the parenthesis (either position or velocity Verlet)...

Let us perform a similar analysis for the Symplectic Euler method. Because Sym-plectic Euler cannot be applied to scalar equations (due to its partitioned structure) we must work directly with the second order system, but this is straightforward for the harmonic oscillator. The timestep map is defined by... 

This is the so-called linear stability condition of the Symplectic Euler method if hQ < 2 the integrator is stable. When hQ > 2, the eigenvalues of the discretization method are both real, with one strictly inside and one strictly outside the unit circle. This implies that the method will exhibit exponentially growing solutions. We say that the stability threshold of the Symplectic Euler method is 2/f2. 

Fig. 4.1 Stability regions for Euler s method (l ) and Symplectic Euler/Verlet (right). When a harmonic oscillator is treated using these methods, the origin is unstable for Euler s method, regardless of stepsize—this means that there is no choice of scaling h which will allow us to ensure that 11 + ftA, < 1. On the other hand, the Verlet method has an interval of stability on the imaginary axis, and it is always possible to find a value of h which guarantees that hQ < 2...

Actually, many methods in common use for molecular dynamics cannot be seen as maps of T M since the hidden constraint is allowed to be violated. A simple example of such a method is the constrained Symplectic Euler-Uke method... 

The projected symplectic constrained method (4.20)-(4.24) is only first order accurate. We forego providing a detailed proof of this fact, but note that it could be demonstrated using standard methods [164]. Note that (4.20)-(4.24) reduces to the symplectic Euler method in the absence of constraints, and the projection of the momenta would not alter this fact. There are several constraint-preserving, second-order alternatives which generalize the Stormer-Verlet scheme. One of these is the SHAKE method [322]. The original derivation of the SHAKE method began from the position-only, two-step form of the Stormer rule for q = F(q)... 

In practice, the implementation of this method is essentially the same as that of the constrained Symplectic Euler method (including the cotangency projection). Moreover, the symplecticness of (4.30)-<4.34) may be demonstrated using a very similar approach to that used to demonstrate the symplectic condition for (4.20)-(4.24). 

Recall that the Stormer-Verlet method could be constructed by composing steps using Symplectic Euler and its adjoint method. Using more complicated methods it is possible to build higher order schemes. It seems natural that a similar procedure should be possible in the constrained setting. But what, precisely, is the adjoint method in the case of (4.20)-(4.24) ... 

For Symplectic Euler with constraints, we view Qh as a mapping of the co-tangent bundle and write Q,P) = Qh(s,p) where... 

Now consider the application of the Symplectic Euler method, q + = q + hpn+i, Pn+1 = Pn - Mn- As wc know, the solution through any given initial point remains confined to an ellipse for all n defined by Hh(q,p) = Hh(qo,Po), where Hh is the perturbed Hamiltonian defined by the backward error analysis. This implies that a density which is a function of Hh will remain invariant under the associated Liouvilie equation. There are many solutions of Clp = 0, none of which is attractive. 

This corresponds to using the adjoint symplectic Euler scheme to solve the Newtonian part of the Langevin dynamics SDE, followed by an exact OU solve. An alternative is to use velocity Verlet for the Hamiltonian part, resulting in a scheme denoted BABO] with associated map... 

A simple first order example of a Langevin dynamics integrator is the method obtained by composing one of the Symplectic Euler variants with an Ornstein-Uhlenbeck step. To get a feel for how the expansion goes, let us work out its terms for such a splitting scheme in the case of a one degree of freedom model with unit mass H = + f/(g). To be explicit, let us say our numerical method first solves... 

---