4
$\begingroup$

I'm very new to computational Physics and am finding conflicting statements on whether the velocity Verlet algorithm, defined as:

$\begin{align} x_{n+1} &= x_n + v_n \Delta t + \frac{1}{2} a_n \Delta t^2 \tag{1}\\ v_{n+1} &= v_n + \frac{1}{2}(a_n + a_{n+1}) \Delta t \tag{2} \end{align}$

is time-reversible.

From my understanding, a method is time-reversible if it is invariant under $\Delta t \to -\Delta t$. Applying this to the velocity Verlet algorithm, equation 1 becomes:

$\begin{align} x_{n-1} &= x_n - v_n \Delta t + \frac{1}{2} a_n \Delta t^2 \\ x_n &= x_{n-1} + v_n \Delta t - \frac{1}{2} a_n \Delta t^2 \\ x_{n+1} &= x_{n} + v_{n+1} \Delta t - \frac{1}{2} a_{n+1} \Delta t^2 \tag{1'} \end{align}$

and equation 2 becomes:

$\begin{align} v_{n-1} &= v_n - \frac{1}{2}(a_n + a_{n-1}) \Delta t \\ v_n &= v_{n-1} + \frac{1}{2}(a_n + a_{n-1}) \Delta t \\ v_{n+1} &= v_{n} + \frac{1}{2}(a_{n+1} + a_{n}) \Delta t \tag{2'} \end{align}$

So far, we see equation 2' is equivalent to equation 2 (so the velocity part is time reversible for step $n$), but equation 1' has the wrong times for acceleration and velocity, as well as wrong sign for the acceleration part. However, if we substitute in equation 2 (or 2') into equation 1', this recovers back equation 1 again. Therefore, I would have thought velocity Verlet is time reversible.

On the contrary, in Basic Concepts in Computational Physics by Stickler and Schachinger, they state on page 108:

The Stormer-Verlet algorithm$^{\dagger}$ is time-reversal symmetric (invariant under the transformation $\Delta t \to - \Delta t$), hence reversible. This is a direct consequence of its relation to the symplectic Euler method. [...] The leap-frog algorithm or the velocity Verlet algorithm [methods] are not time-reversal invariant.

They repeatedly state that leapfrog and velocity verlet are not time-reversible throughout the rest of the book. They define leapfrog as:

$\begin{align} x_{n+1} &= x_{n} + v_{n+\frac{1}{2}} \Delta t \tag{3} \\ v_{n+\frac{1}{2}} &= v_{n-\frac{1}{2}} + a_n \Delta t \tag{4} \\ v_{\frac{1}{2}} &= v_0 + \frac{1}{2}a_0 \Delta t \tag{5} \end{align}$

Using the same procedure as for velocity Verlet, I found that equations 3 and 4 are invariant under time reversal, but the initialising condition does change to:

$$v_{\frac{1}{2}} = v_0 + \frac{1}{2}a_{\frac{1}{2}} \Delta t \tag{5'}$$

where 5' differs from 5 since the acceleration at $t = \frac{1}{2} \Delta t$ is used instead of at $t = 0$.

Looking at the half time step implementation of velocity Verlet, defined as:

$\begin{align} v_{n + \frac{1}{2}} &= v_n + \frac{1}{2} a_n \Delta t \tag{6} \\ x_{n+1} &= x_n + v_{n + \frac{1}{2}} \Delta t \tag{7} \\ v_{n+1} &= v_{n + \frac{1}{2}} + \frac{1}{2} a_{n+1} \Delta t \tag{8} \end{align}$

I have found under time-reversal, these transform to:

$\begin{align} v_{n + \frac{1}{2}} &= v_n + \frac{1}{2} a_{n+\frac{1}{2}} \Delta t \tag{6'} \\ x_{n+1} &= x_n + v_{n + \frac{1}{2}} \Delta t \tag{7'} \\ v_{n+1} &= v_{n + \frac{1}{2}} + \frac{1}{2} a_{n+\frac{1}{2}} \Delta t \tag{8'} \end{align}$

So equation 7 is invariant, but the times for the accelerations at equations 6 and 8 are changed. Does this suggest the substitution of 2' into 1' was invalid and velocity Verlet isn't time reversible? Actually, equation 6' is equivalent to equation 8 and equation 8' is equivalent to equation 6 (as seen by adding or subtracting $\frac{1}{2}\Delta t$ to the times the equations are evaluated at). So does this mean velocity Verlet is still time reversible?

Lastly, in Computational Physics by Thijssen, it's stated that:

There exist two alternative formulations of the Verlet algorithm$^{\dagger}$, which are exactly equivalent to it in exact arithmetic but which are less susceptible to errors resulting from finite numerical precision in the computer than the original version. The first of these [is] the leap-frog form [and the second is] the so-called velocity-Verlet algorithm which is also more stable than the original Verlet form.

To summarise, my questions are:

  1. Are Stickler and Schachinger incorrect in stating that the velocity Verlet and leapfrog algorithms are not time reversible?
  2. Is the substitution of 2' into 1' valid to show that velocity Verlet is time reversible?
  3. Does the fact that the initialising condition of leapfrog is not time reversible mean that the leapfrog method itself is not time reversible? Other resources, including answers on this stack exchange, state leapfrog is time reversible.
  4. Since Thijssen states the velocity verlet and leapfrog algorithm are "exactly equivalent" in terms of the arithmetic, does this mean they inherit the time-reversibility (and possibly even the symplectic nature) of the Stormer-Verlet algorithm?

$^{\dagger}$For reference, the Stormer-Verlet algorithm is outlined here. This Wikipedia page also entails how to obtain velocities within this framework.

https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6496 and https://math.stackexchange.com/questions/1448005/what-does-the-time-reversibility-of-verlet-or-other-integration-mean also both state velocity Verlet is time reversible

https://www5.in.tum.de/lehre/vorlesungen/sci_compII/ss13/uebungen/blatt9solution.pdf uses the substitution outlined earlier to show velocity Verlet is time reversible

$\endgroup$
1
  • $\begingroup$ If you are interested in MD, you might find the new-ish Matter Modeling SE site useful as well. $\endgroup$
    – Tyberius
    Nov 11, 2021 at 16:46

2 Answers 2

2
$\begingroup$

I can share my perspective on this topic and can try to prove that the Leapfrog (velocity Verlet) is time-reversible, according to an appropriate definition of time reversibility (very nice property that increases this integration method's accuracy, alongside its $\Delta t^2$ convergence and symplecticity).

The idea behind Leapfrog:

Assume you have a system of differential equations, that can be written, in vector form, as follows: $$\frac{dx}{dt} \,=\, f(x) \, + \, g(x)$$ You can split it into two systems \begin{align} &\frac{dx}{dt} \,=\, f(x)\\ &\frac{dx}{dt} \,=\, g(x) \end{align} Denote by $\Phi^t(x)$ the phase flow of the original system and by $\varphi^t(x)$ and $\psi^t(x)$ the phase flows of the two separated systems. Now, if the vector fields, associated to the two separated systems commute, i.e. $[f(x), g(x)] = 0$, then one can prove the following theorem: $$\Phi^t(x) = \varphi^t\big(\, \psi^t(x) \,\big) = \psi^t\big(\, \varphi^t(x) \,\big) = \varphi^{\frac{t}{2}}\Big(\, \psi^t\big(\, \varphi^{\frac{t}{2}}(x) \,\big) \,\Big) $$ Furthermore, loosely speaking, if the vector fields above do not commute, but their commutator is small enough, i.e. $[f(x), g(x)] \sim 0$, then we could expect that $$\Phi^t(x) \sim \varphi^{\frac{t}{2}}\Big(\, \psi^t\big(\, \varphi^{\frac{t}{2}}(x) \,\big) \,\Big) $$ and in particular, there is a theorem that asserts that $$\Phi^t(x) = \varphi^{\frac{t}{2}}\Big(\, \psi^t\big(\, \varphi^{\frac{t}{2}}(x) \,\big) \,\Big) \, + \, O(t^2) $$ or something like that.

Now, in classical mechanics, the Newtonian equations of motion very often look like this: \begin{align} &\frac{dx}{dt} = v\\ &\frac{dv}{dt} = F(x) \end{align} where $x, \, v \, \in\, \mathbb{R}^n$. Such a system is time reversible in the following sense. Assume you start from any point $(x, v)$ and you following the solution of the differential equation, moving over time $t$ to the point $(x_t, v_t)$. Now, invert the velocity at the point of arrival, i.e. take $(x_t, -\,v_t)$. Starting from this new point $(x_t, -\,v_t)$, follow the solutions to the differential equations over the same time interval $t$. You will arrive at the state $(x, -\,v)$.

Important: The fact that $F(x)$ depends exclusively on $x$ and does not depend on $v$ is crucial for the time-reversibility! In cases where the vector function $F(x, v)$ depends on the velocity, time-reversibility my no longer hold!

One can express this symmetry by saying that the system is invariant under the $2n \times 2n$ linear transformation $$\begin{bmatrix} x \\ v \end{bmatrix} \,\mapsto \,\begin{bmatrix} I_n & 0 \\ 0 & -\,I_n \end{bmatrix} \begin{bmatrix} x \\ v \end{bmatrix}$$ or symbolically, $$(x, v) \, \mapsto \, T(x, v) = (x, -v)$$ combined with the inversion of time $t \mapsto -t$ (hence the name time-reversibility).

On the level of the vector field corresponding to the system of differential equations above, one has to check two things. First, one has to change the variables $T : (x, v) \mapsto (x, -v)$ and plugs them in the vector field, obtaining: $$\begin{bmatrix} v \\ F(x) \end{bmatrix} \,\mapsto \,\begin{bmatrix} -\,v \\ F(x) \end{bmatrix}$$ Second, on needs to check how the derivative $DT$ of $T$ changes the vector field: $$\begin{bmatrix} v \\ F(x) \end{bmatrix} \,\mapsto \,\begin{bmatrix} I_n & 0 \\ 0 & -\,I_n \end{bmatrix} \begin{bmatrix} v \\ F(x) \end{bmatrix} \, = \, \begin{bmatrix} v \\ -\,F(x) \end{bmatrix}$$ After that, inverting time $t \mapsto -t$ results in the multiplication of the latter vector field by $-1$ and thus $$\begin{bmatrix} v \\ -\, F(x) \end{bmatrix} \,\mapsto \, (-1)\begin{bmatrix} v \\ -\, F(x) \end{bmatrix} \, = \, \begin{bmatrix} -\,v \\ F(x) \end{bmatrix}$$ equals the variable transformed vector field form the first step.

Consequently, if $\Phi^t(x,v)$ is the phase flow of the position-velocity Newtonian system above, then time-reversibility can be expressed as the time-inversion commutativity relation $$\Phi^{(-t)} \circ T(x, v) = T \circ \Phi^t(x,v)$$

The Leapfrog method comes from the splitting of the system as follows: \begin{align} &\frac{dx}{dt} = 0 \, + \, v\\ &\frac{dv}{dt} = F(x) \, + \, 0 \end{align} or in a vector form $$\frac{d}{dt}\begin{bmatrix} x \\ v \end{bmatrix} \, = \, \begin{bmatrix} 0 \\ F(x) \end{bmatrix} \, + \, \begin{bmatrix} v \\ 0 \end{bmatrix}$$ so here $f = \begin{bmatrix} 0 \\ F(x) \end{bmatrix}$ and $g = \begin{bmatrix} v \\ 0 \end{bmatrix} $ which means we can form the two separated systems \begin{align} &\frac{dx}{dt} = 0 \\ &\frac{dv}{dt} = F(x) \end{align} and \begin{align} &\frac{dx}{dt} = v\\ &\frac{dv}{dt} = 0 \end{align} The good news is that these two are exactly solvable (no numerical approximation is needed!) $$\varphi^t(x,v)\, = \big(\, x\,, \,\, v \, + \,t\,F(x)\,\big)$$ $$ \psi^t(x,v) \, = \big(\, x \, + \, t\, v\,, \,\, v\,\big)$$ You can directly check that both of these systems are time reversible, i.e. $$\varphi^{(-t)} \circ T(x,v) = T \circ \varphi^t(x,v)$$ and $$\psi^{(-t)} \circ T(x,v) = T \circ \psi^t(x,v)$$ By meticulously following of the composition of the three transformation (actually it is pretty straightforward), one can show that Leapfrog (velocity Verlet) is obtained as the composition of phase flows $$\varphi^{\frac{\Delta t}{2}}\Big(\, \psi^{\Delta t}\big(\, \varphi^{\frac{\Delta t}{2}}(x, v) \,\big) \,\Big) $$ As each of them is time-invertible $T-$invariant (commutes with $T$ combined with time-inversion) the total composition is also $T-$invariant, i.e. $$\varphi^{\frac{-\,\Delta t}{2}}\Big(\, \psi^{-\,\Delta t}\big(\, \varphi^{\frac{-\,\Delta t}{2}} \circ T (x, v) \,\big) \,\Big) = T \circ \varphi^{\frac{\Delta t}{2}}\Big(\, \psi^{\Delta t}\big(\, \varphi^{\frac{\Delta t}{2}}(x, v) \,\big) \,\Big) $$

Ok, for the record, let's show how the derivation of the velocity Verlet works. Start with point $(x_n, v_n)$: \begin{align} &x_{n+\frac{1}{3}} \, = \, x_n\\ &v_{n+\frac{1}{3}} \ = \, v_n \, + \, \frac{\Delta t}{2}\, F(x_n) \end{align} $$$$ \begin{align} &x_{n+\frac{2}{3}} \, = \, x_{n+\frac{1}{3}} \, + \, \Delta t \, v_{n+\frac{1}{3}} \, =\, x_n \, + \, \Delta t \, v_n \, + \, \frac{\Delta t^2}{2}\, F(x_n)\\ &v_{n+\frac{2}{3}} \ = \, v_{n+\frac{1}{3}} \, = \, v_n \, + \, \frac{\Delta t}{2}\, F(x_n) \end{align} $$$$ \begin{align} &x_{n+1} \, = \, x_{n+\frac{2}{3}} \, =\, x_n \, + \, \Delta t \, v_n \, + \, \frac{\Delta t^2}{2}\, F(x_n)\\ &v_{n+1} \ = \, v_{n+\frac{2}{3}} \, + \, \frac{\Delta t}{2}\, F\Big(x_{n+\frac{2}{3}}\Big) \, = \, v_n \, + \, \frac{\Delta t}{2}\, F(x_n) \,+\, \frac{\Delta t}{2}\, F(x_{n+1}) \end{align} And there you finally have it: \begin{align} &x_{n+1} \, =\, x_n \, + \, \Delta t \, v_n \, + \, \frac{\Delta t^2}{2}\, F(x_n)\\ &v_{n+1} \ = \, v_n \, + \, \frac{\Delta t}{2}\Big(\, F(x_n) \,+\, F(x_{n+1}) \, \Big) \end{align}

$\endgroup$
3
  • $\begingroup$ You have not changed time, but reversed $v$. This is not enough, you have to reverse the force field $F$ also, for the system to retrace its path. The other way would be to just reverse $t$, which has the same effect as reversing both $v,F$. $\endgroup$
    – cfdlab
    Dec 21, 2020 at 4:29
  • $\begingroup$ @cfdlab Yes, you are right that I did not invert time and I should have. Sloppy on my part (that's what happens when one skips a proof). I have fixed the mistake, included the relevant arguments and now the argumentation should be correct. $\endgroup$ Dec 21, 2020 at 13:57
  • $\begingroup$ If you have a general interest in MD and molecular modeling, you may also be interested in the new-ish Matter Modeling SE $\endgroup$
    – Tyberius
    Nov 11, 2021 at 16:48
2
$\begingroup$

I think for time reversibility we need to see if we can predict velocity and position at time $t$ if the values at $t+\Delta t$ is known.

Verlet Algorithm

The integrator equation is given as:

$r_{t+\Delta t} = 2r_t -r_{t-\Delta t} + \Delta t^2 a_t$

Given we know the positions at $t-\Delta t$ and $t$ we can predict the position at $t+\Delta t$. So if this is time reversible we should be able to predict $t-\Delta t$ if we know the positions at $t+\Delta t$ and $t$. By simply moving variables you can see that,

$r_{t-\Delta t} = 2r_t -r_{t+\Delta t} + \Delta t^2 a_t$

So Verlet is time time reversible

Leap-Frog

The integrator equations:

(1) $v_{t+\Delta t/2} = v_{t-\Delta t/2} + \Delta ta_t $

(2) $r_{t+\Delta t} = r_{t} + \Delta tv_{t+\Delta t/2}$

For time reversibility you first find,

$r_t = r_{t+\Delta t} -\Delta t v_{t+\Delta t/2} $

Note we know the value of at $v_{t+\Delta t/2}$

Then, we can compute $a_t$ for $r_t$. Followed by,

$v_{t-\Delta t/2} = v_{t+\Delta t/2} - \Delta t a_t$

Velocity Verlet

You can also prove the same for Velocity Verlet. I will summarize the key equations.

$\alpha = v_{t} + \frac{1}{2}\Delta t a_t = v_{t+\Delta t} -\frac{1}{2}\Delta t a_{t+\Delta t}$

Using $\alpha$ you can find $r_t$ as,

$r_t = r_{t+\Delta t} -\Delta t \left( v_t + \frac{1}{2}\Delta t a_t \right)$ $= r_{t+\Delta t} - \Delta t \alpha$

Once you know $r_t$, $a_t$ is known and,

$v_t = \alpha - \frac{1}{2}\Delta t a_t $

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.