Which numerical methods preserve time reversal symmetry?

Question

If I have a physical system which contains a time reversal symmetry (for example a Hamiltonian $H(x,p)=p^2/2m + V(x)$ with $V(x)$ real) and I want to solve the differential equations which describe this system, which solver for ODEs should I use in order to keep the time reversal symmetry (for example in mathematica)? Which solvers break this symmetry?

EDIT: I want to extend this question. Let us consider a system of coupled first order differential equations $$\dot{a}_1 (t) = f_1(a_1,a_2,a_3,\ldots,a_n;t) \\ \dot{a}_2(t) = f_2(a_1,a_2,a_3,\ldots,a_n;t) \\ \dot{a}_3(t) = f_3(a_1,a_2,a_3,\ldots,a_n;t) \\ \vdots$$ What integration method is best used if the underlying system contains a time reversal symmetry?

Answer 1

What one usually wants in this situation is to preserve a discrete analog of time symmetry: namely, if the time discretization is applied to solve first forward and then backward in time, the initial condition is recovered. This is true if the method is invariant under the following substitutions:

$$ \Delta t \to -\Delta t$$ $$ a^{n+j} \to a^{n-j}$$

(here $a^n$ is the numerical approximation of the solution $a(t^n)$, so the second substitution is implied by the first).

I will give two examples to illustrate. The explicit Euler method $$a^{n+1} = a^n + \Delta t f(a^n)$$ does not preserve time symmetry; applied backward in time it becomes the implicit Euler method: $$a^n = a^{n-1} + \Delta t f(a^n).$$ On the other hand the midpoint (or leapfrog) method $$a^{n+1} = a^{n-1} + 2\Delta t f(a^n)$$ does preserve time-reversal symmetry. Other well-known methods that preserve time-reversal symmetry include the trapezoidal method and (as mentioned in the comments) the Verlet method.