This paper is concerned with the numerical integration of initial value problems for Newton’s second order equations \(x''=F(x)\), where \(F=F(x)\) is a conservative force, i.e. \(F(x)=-\text{grad }U(x)\) and the dimensionality \(N\) of the space variable \(x\) is very high. In typical applications which arise in astrophysics and molecular dynamics models, \(N={\mathcal O}(10^4)\) and for simple potentials, the dominating cost for evaluating \(F\) has \({\mathcal O}(N^2)\) complexity. As a consequence of the high dimensionality, numerical integrators for these types of problems must limit as far as possible the number of function evaluations per step and therefore the leap-frog formula given by \(v(t + \Delta t/2)=v(t-\Delta t/2)+\Delta t F(x(t))\), \(x(t+\Delta t/2) =x(t)+\Delta tv(t+\Delta t/2)\) is standard for such a kind of problem.
The aim of the paper is to improve the performance of the leap-frog method without increasing the number of function evaluations per step by modifying appropriately the advancing formula of the velocity. This means that instead of using a constant interpolation for \(F(x)\) to advance the velocity vector, the authors propose an interpolation based on three predicted values with a symmetry that preserves reversibility. In this way since the new formula is based on a three point approximation of \(F(x)\), it integrates exactly one-dimensional, harmonic oscillators and for multidimensional oscillators it integrates exactly some quadratic approximation of their potential.
Finally, the paper presents the results of some numerical experiments comparing the performance of the new method with the standard leap-frog and the fourth order Runge-Kutta method for some harmonic and anharmonic oscillators.