1

Consider a general set of Hamilton's equations $$ \begin{align} \dot{q}(q, p) &= \partial_p H(q, p) \\ \dot{p}(q, p) &= -\partial_q H(q, p) \end{align} $$ A first-order integrator one could use is the sympletic-Euler. There are two second-order methods I know of that can be used to numerically integrate these equations: Leapfrog and position-Verlet.

Is there any difference at all between using Leapfrog or position-Verlet?

Leapfrog (second order, sympletic, reversible)

This is like a sympletic Euler but sandwiched between a half-momentum step at the beginning and at the end. $$ \begin{align} p_{t+\delta/2} &= p_t + \frac{\delta}{2}\dot{p}(q_t, p_t) \\ q_{t+\delta} &= q_t + \delta \dot{q}(q_t, p_{t+\delta/2}) \\ p_{t+\delta} &= p_{t+\delta/2} + \frac{\delta}{2}\dot{p}(q_{t+\delta}, p_{t+\delta/2}) \end{align} $$

Position-Verlet (second order, sympletic)

Similar to the Leapfrog method but starts with a half-position update rather than half-momentum update. \begin{align} q_{t+\delta/2} &= q_t + \frac{\delta}{2}\dot{q}(q_t, p_t) \\ p_{t + \delta} &= p_t + \delta \dot{p}(q_{t+\delta}, p_t) \\ q_{t+\delta} &= q_{+\delta/2} + \frac{\delta}{2}\dot{q}(q_{t+\delta/2}, p_{t+\delta}) \end{align}

Is there any difference between using Leapfrog or position-Verlet? Either numerically or theoretically.

Euler_Salter
  • 247

1 Answers1

1

The two methods are mathematically and numerically equivalent. The easiest way to see this is to consider the update as, \begin{align} p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ q &= q + \delta t \cdot \dot{q}(p,\,q) \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \end{align} This kind of update is done $n$ times in a sequential manner: \begin{align} p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ q &= q + \delta t \cdot \dot{q}(p,\,q) \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ q &= q + \delta t \cdot \dot{q}(p,\,q) \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ &\cdots \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \\ q &= q + \delta t \cdot \dot{q}(p,\,q) \\ p &= p + \delta t \cdot \dot{p}(p,\,q) / 2 \end{align} You can then combine the pairs of half-steps with each other while also breaking the individual $q$ updates into two halves. To rectify the fact that Leapfrog updates on a half-cycle offset from the position Verlet, we can add or subtract terms at the beginning & end as needed, which then would yield, \begin{align} q &= q + \delta t \cdot \dot{q}(p,\,q) / 2 \\ p &= p + \delta t \cdot \dot{p}(p,\,q) \\ q &= q + \delta t \cdot \dot{q}(p,\,q) / 2 \end{align} Hence, the two methods are mathematically equivalent. For a programmatic viewpoint, see this blog post.

The only difference is that in one case, the force is computed once while in the second case it is computed twice. It likely would be preferable, for computation efficiency reasons, to compute the more expensive function just once per step. If it happens that $\dot{p}$ is more expensive than $\dot{q}$, use position Verlet; if not, use Leapfrog/velocity Verlet.

In the comments, you indicated that "a few papers [...] show that the velocity-Verlet is numerically more accurate". Note that the thesis in question is discussing the original Verlet integrator, $$ q(t + \Delta t) = 2q(t) − q(t − \Delta t) + \Delta t^2\frac{F\left(q(t)\right)}{M}+\mathcal{O}\left(\Delta t^4\right) $$ with Leapfrog/velocity Verlet method, \begin{align} q\left(t+\Delta t\right) &= q(t)+\Delta tv(t)+\frac{\Delta t^2}{2}\frac{F\left(q(t)\right)}{M} \\ v\left(t+\Delta t\right) &= v(t) + \frac{\Delta t}{2}\frac{F\left(q(t+\Delta t)\right)+F\left(q(t)\right)}{M} \end{align} The reference Swope et al 1982, for instance, indicates that it is due to the rounding errors in finite precision on a computer that the velocity Verlet method is superior. Their appendix doesn't seem to cover examples. Tuckerman, Berne and Martyna 1992 doesn't actually seem to discuss numerical accuracy between the two methods.

Kyle Kanos
  • 29,127