Hacker News new | ask | show | jobs
by Xcelerate 4925 days ago
This same thing is used in molecular dynamics simulations. For instance, there is an algorithm called RESPA that is used to break integrations of different types of particle interactions into appropriate timestep intervals. Bond vibrations must be calculated much more frequently than non-bonded interactions.

The algorithm (reversible RESPA) is formally derived from the Liouville operator (which governs the time evolution of any property):

    A(t) = exp(iLt) * A(0)
For instance, A(t) can be position or momentum. The Liouville operator must be symmetric in order to generate a reversible numerical integration algorithm.

The result of all this is basically that:

    p(t + ∆t/2) = p(t) + ∆t/2 F(r(t))
    r(t + ∆t) = r(t) + ∆t p(t + ∆t/2)
    p(t + ∆t) = p(t + ∆t/2) + ∆t/2 F(r, t + ∆t)
where p is momentum, r is position, and F is force.
1 comments
zbyszek 4925 days ago
This was also colloquially known as the "leapfrog" algorithm and is the simplest of a class of integrators that are symmetric (in simulation time) and symplectic which are crucial properties for some simulations. I take it that RESPA is the partitioning of the time evolution operator into components with different force gradients, and the application of different integration schemes to those components. One can also generalise leapfrog to integrate the momentum (or some part) with n steps of dt/n.

The metric used to describe their accuracy is the degree to which they violate the conservation of energy, which can be shown to be an odd integral power of the timestep dt per step [0]. The error in leapfrog goes like the 3rd power.

Higher order integration schemes can be derived e.g. [1]. They may not be useful in practice, depending on the cost of computing the individual terms and the accumulation of finite precision errors. But the known scaling behaviour provides a nice way of verifying the calculation of the evolution operators. Another nice thing to do is to compute the "round trip", i.e. integrate forwards in time, and then backwards. With a symmetric integrator you should end up where you started in terms of position, momentum and energy, regardless of step size, so computing a suitable difference and seeing how it scales with trajectory length can be informative. (e.g. one can compute a Lyapunov exponent from such round-trips to see if the underlying dynamics are chaotic).

[0] McLachlan and Atela, Nonlinearity 5 (1992) http://www.massey.ac.nz/~rmclachl/si.ps (PS) [1] M. Creutz, A. Gocksch; Phys. Rev. Lett. 63 (1989) http://thy.phy.bnl.gov/~creutz/mypubs/pub106.pdf (PDF)

link