Time Reversible Algorithms

From QED

Time-Reversible Algorithms for Differential Equations

A while back I learned about a cute class of numerical algorithms that are exactly time-reversible at the bit level. These are also called "bit-reversible" algorithms.

Background and Motivation

There is a class of algorithms for integrating ODE's for Hamiltonian systems that insure conservation of a surrogate Hamiltonian (a discrete approximation to the true Hamiltonian), to within round off error. These are called "symplectic algorithms" (the simplest example is a leap-frog type of algorithm). What this means is that the solution is guaranteed to lie (within round off error) on a surface corresponding to conservation of this surrogate Hamiltonian (which therefore lies close to a surface of constant true energy). Thus there may be phase-like errors, but there is no numerical dissipation in the solution. However, there can still be round-off errors that can build up over many time steps. While an exact Hamiltonian system is time-reversible, round-off errors will cause a standard numerical solution (even a symplectic one) to not be exactly time reversible. Particularly for chaotic systems where a large Lyapunov exponent implies exponential separation of solutions (and thus exponential growth of errors), this might significantly limit the length of time for which one could try to simulate a Hamiltonian system before it was no longer accurately time-reversible because of the build up of errors. Of course this is to be expected, and, as usual for simulations of chaotic systems, one can still get meaningful and important information from the statistics of the long-time solutions.

A solution to achieve exact time-reversibility

However, there is a very simple way to modify a symplectic algorithm to make it time-reversible. To see this, consider the simple system

dx / dt = v

dv / dt = F(x)

Write this as

d²x / dt² = F(x)

and do the equivalent of a leap-frog discretization:

x(n+1) - 2 x(n) + x(n-1) = dt**2 * F(x)

This is symplectic. To make it fully bit reversible, all you have to do is add 5 characters, so that it becomes

x(n+1) - 2 x(n) + x(n-1) = Fix(dt**2 * F(x))           Eq.(1)

where Fix(...) is a function that rounds its argument to a fixed-point precision (and x(n) is in the same fixed-point precision). This equation can then be advanced in time to determine x(n+1) given x(n) and x(n-1). Or given x(n+1) and x(n), it can be solved backwards in time to determine x(n-1), and it will exactly reproduce the same x(n-1) that was needed to give that x(n+1). Thus it is exactly time-reversible at the bit level. Of course there are still round-off errors introduced when the term dt**2*F(x) is truncated to a fixed resolution, and there are still truncation errors from the finite order of the algorithm, but these are both time-reversible errors and numerically preserve the time-reversible property of a Hamiltonian system.

The essential idea here is that there are no round off errors in integer arithmetic. One doesn't need to literally restrict the calculations to integers to do this, as there are also no round off errors in fixed-point arithmetic. By rounding off certain key terms to a given fixed-point precision before adding, one can mimic fixed-point arithmetic with floating-point hardware.

Consider the case where x is an integer, scaled to the appropriate units so that a 16-digit integer could accurately represent the position in those units, and the time step and force are such that a 16-bit integer could also accurately represent dt^2*F(x). Then one could convert x to a real number to do floating point calculation of the force and dt^2*F(x), and then Fix(dt^2*F(x)) would convert back to an integer.

But it isn't necessary to be restricted to an integer, one can work with real numbers as long as one uses fixed-point arithmetic for the final addition/subtraction operations. The position of the fixed-point can be chosen to be whatever is convenient to accurately represent the scales in which x is measured. In any case, this fixed-point representation then avoids the irreversible round-off errors of floating-point representations. Consider the example where x is the position in meters of a planet that is going around a star, so that the values of x never exceeded the range [-xmin,xmax], where xmax = 10¹³. If the precision of floating point numbers on your computer is N_prec bits (N_prec=53 for double precision in the IEEE standard, about 16 decimal digits of accuracy), then one could define the Fix() function by

N_max = int(log(xmax)/log(2.0))+1 ! so 2^N_max >= xmax.

scale = 2**(N_max - N_prec)

Fix(x) = scale * nint( x/scale )

(where nint() returns the nearest integer). (I might be off by 1 in the definition of N_max or N_Prec for an optimal result that guarantees avoidance of overflows for all |x|<=xmax? Should test this before using it.) (In Fortran90 one might use the anint() function, or one could use bit manipulations to avoid the floating point multiplications, since scale is a power of 2.) Thus one can convert a symplectic code into a time-reversible code with just a few lines of code changed. Using this Fix() function is a simple way of mimicking fixed-point arithmetic with floating-point hardware.

Instead of using this fixed-point arithmetic trick, the other option is to convert your code to using integer arithmetic in most places, which allows you to use the full 64 bits of accuracy of 8 byte integers (if one has optimally scaled the variables in the problem), instead of the 53 bits of accuracy of 8 byte double-precision floating point numbers in the above method that used the Fix() function.

One can extend this idea to turn other symplectic algorithms into bit-level time-reversible algorithms.

While an algorithm that is exactly time/bit-reversible is a neat trick, in many cases it might not be important. However, there might be some cases where it really is helpful to ensure numerically-exact time reversibility, see some of the references below. One can use such an algorithm to verify that the key properties of a simulation of chaotic system didn't depend on whether the numerical algorithm was exactly time-reversible or not. This is what we expect, as we understand from statistical mechanics and chaos theory how a deterministic time-reversible Hamiltonian system can exhibit behavior that appears to be (and effectively is) random and irreversible.

Symplectic integrator properties

As mentioned above, symplectic integrators don't exactly conserve the original Hamiltonian, but they conserve a "surrogate Hamiltonian" (a.k.a. a "shadow Hamiltonian" or "modified Hamiltonian"), which is a numerical or discrete variation of the original Hamiltonian. If the integrator is Nth order, the difference between the modified Hamiltonian and the exact Hamiltonian will also be Nth order.

Consider a simple leapfrog algorithm written in a "drift-kick-drift" form as

x_{n+1/2} = x_n + (dt/2) * v_n

v_{n+1} = v_n + dt*F(x_{n+1/2})/m

x_{n+1} = x_{n+1/2} + (dt/2) * v_{n+1}

where the force F(x) = - dV/dx with a quadratic potential V(x) = k*x^2/2. One can show that the modified Hamiltonian that is conserved by this algorithm is

H = (m/2) v_n^2 (1-dt^2 Omega^2 / 4) + k x_n^2 / 2

where Omega^2 = k/m. (For more complicated systems, one usually can't find an exact analytic expression for the shadow Hamiltonian, though one can sometimes work out next order asymptotic expressions for it. It is still nice to know that this conserved quantity exists.)

(The simple symplectic algorithm used here may work only for Hamiltonians of a particular form, such as a separable Hamiltonian. Symplectic algorithms for some more complicated Hamiltonians may require implicit solution techniques.)

References

I learned about this idea from my colleague Charles Karney, who did ground-breaking work in chaos theory (including applications to plasma physics and fusion energy research) in the 70's and 80's. The fundamental idea of using integer arithmetic to preserve time-reversibility in mappings goes back to work by Miller and Predergast (1968) and by Rannou (1972). Karney doesn't claim to to have invented this idea and his 1986 paper cites these two papers. But these earlier works were apparently on a coarse lattice and it seems to me that Karney might deserve the credit for showing how to extend this idea to obtain resolution comparable to floating point numbers in a simple way. These ideas have been applied and rediscovered in various fields. I might not have the history completely right, but the following references (and references therein) are a good place to start:

Charles F. F. Karney (1983), "Long-time correlations in the stochastic regime", Physica D 8, 360

Charles F. F. Karney (1986) "Numerical Techniques for the Study of Long-Time Correlations", Particle Accelerators 19, 213

R. H. Miller (1996), "Some comments on numerical methods for chaos problems", Celestial Mechanics and Dynamical Astronomy 64, 33.

David J. D. Earn and Scott Tremain (1992), "Exact numerical studies of Hamiltonian maps: Iterating without roundoff error", Physica D 56, 1

F. Rannou (1974), "Numerical Study of Discrete Plane Area-preserving Mappings", Astronomy and Astrophysics 31, 289 (Ph.D. Dissertation, 1972).

The next two papers and book are very interesting, and I got the above leapfrog example from them (but it seems they weren't aware of the earlier work above):

D. Levesque and L. Verlet (1993), "Molecular Dynamics and Time Reversibility", J. of Statistical Physics, 72, 519

Oyeon Kum and William G. Hoover (1994), "Time-Reversible Continuum Mechanics", J. Stat. Phys. 76, 1075 (1994)

William Graham Hoover (1999), Time Reversibility, Computer Simulation, and Chaos.

The form of the algorithm in Eq.1 makes clear the reversibility property, but to minimize round-off errors it is better to solve it in another form, like the drift-kick-drift version of leapfrog shown above. See

Ross A. Lippert, Kevin J. Bowers, Ron. O. Dror, et al. (2007) "A Common, Avoidable Source of Error in Molecular Dynamics Integrators", Journal of Chemical Physics 126, 046101.

Other related articles:

Nobuyoshi Komatsu, Takahior Kiwata, and Shigeo Kimura (2008), "Numerical irreversibility in self-gravitating small N-body systems" Physica A 37, 2267

Miguel Preto and Scott Tremaine (1999), "A Class of Symplectic Integrators with Adaptive Time Step for Separable Hamiltonian Systems", The Astronomical Journal 118, 2532

J. Candy and W. Rozmus (1991), "A symplectic integration algorithm for separable Hamiltonian functions", JCP 92, 230.

Wayne B. Hayes (2008), "Surfing on the edge: chaos versus near-integrability in the system of Jovian planets", MNRAS 386, 295.

M. Tuckerman, B. J. Berne, and G. J. Martyna (1992), "Reversible multiple time scale molecular dynamics", J. Chem. Phys. 97, 1990.

N. F. Loureiro, G. W. Hammett (2008), "An Iterative Semi-Implicit Scheme with Robust Damping", J. Comp. Phys.227, 4518. This paper can be thought of as an iterative implicit method with a physics-based preconditioner (for a particular class of hyperbolic wave problems). One interesting feature is that each iteration of the algorithm remains symplectic, one does not have to iterate to round-off error accuracy.

Other applications:

Time-reversible ODE algorithms can be thought of as invertible maps, which can have other widely useful applications, such as a lossless (and fast) algorithm for image rotation, based on composing 3 shear operations. This is actually used in the pnm package (see Karney 86, and his PPPL-2218 report which came out in 1985, and the pnmrotate man page, which cites a paper by Alan Paeth, 1986).

--Hammett 03:18, 22 June 2009 (UTC)