Greenspan's method

Update for positions:

x_i^m+1 = x_i^m + h/2 (v_i^m+1 + v_i^m)

Update for velocities:

v_i^m+1 = v_i^m + h A_i^m

where we use a different acceleration than before. That is, in other methods the acceleration a_i^m was the sum over j≠i of:

G (x_j^m - x_i^m) / (r_ij^m)³

where r_ij=||x_j-x_i||. On the other hand, Greenspan's method uses the sum over j≠i of:

G ((x_j^m+1 - x_i^m+1) + (x_j^m - x_i^m)) / ( r_ij^m+1 r_ij^m (r_ij^m+1 + r_ij^m))

for the acceleration A_i^m.

This method is energy conserving. The design comes from approximating the force as a discrete derivative of the potential energy (rather than, as before, keeping the force exact while effectively approximating the potential energy). The energy conservation comes at the cost of the method being implicit.

The method also conserves both linear and angular momentum. It is time reversible.

References

Chapter 2 from "N-body problems and models" by Donald Greenspan, 2004
page 173 in Section V.5 of the second edition of "Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations" by Hairer, Lubich, and Wanner, 2006
Section 2.3 from "Numerical Solutions of the N-body Problem" by Andrzej Marciniak, 1985
in Hairer, Norsett, and Wanner's "Solving Ordinary Differential Equations I: Nonstiff Problems": ???
not in the Burden and Faires "Numerical Analysis" textbook: Section 5.2
not in the Silly-pedia?