Context: I'm implementing an engine that simulates $n$ physics bodies, and I am trying to understand and refine my $2$-body-problem engine first.
If $\bf r$ is my position vector, my ODE is
$$ \ddot{{\bf r}} = f({\bf r}) $$
where $f$ is the acceleration due to gravity at given location.
As I understand it, though euler's method is generally applied to first order ODEs, a simple extension of the method to second order ODEs might look like:
$$ \begin{aligned} \textbf{r}_0 &= \text{<initial value>} \\ \dot{\textbf{r}}_0 &= \text{<initial value>} \\ \textbf{r}_{n+1} &= \textbf{r}_n + h\dot{\textbf{r}}_n \\ \dot{\textbf{r}}_{n+1} &= \dot{\textbf{r}}_n + hf(\textbf{r}_n) \\ \end{aligned} $$
Where $h$ is the step size.
I've found that this method does a really poor job of approximating an elliptical orbit, even for extremely small $h$. I've also implemented the standard Runge-Kutta 4th order method, which does a much much better job, but has consistent error relative to what Kepler's equations tell us an elliptical orbit should look like. In both cases, bodies whose trajectories are modelled with these methods gain energy, gradually reaching ever more eccentric orbits.
I have however stumbled on a different method, of which I don't know the name:
$$ \begin{aligned} \textbf{r}_0 &= \text{<initial value>} \\ \dot{\textbf{r}}_0 &= \text{<initial value>} \\ \dot{\textbf{r}}_{n+1} &= \dot{\textbf{r}}_n + hf(\textbf{r}_n) \\ \textbf{r}_{n+1} &= \textbf{r}_n + h\dot{\textbf{r}}_{n+1} \\ \end{aligned} $$
Crucially, the next position vector is calculated using the next velocity vector. Does this method have a name? For this particular use case of approximating an elliptical orbit, it has the interesting property of doing a really good job at keeping the total energy of the system constant (also the angular momentum, the eccentricity, the period, etc.), much better than Runge-Kutta-4, despite having much larger error (several orders of magnitude for my choice of $h$) on the position and velocity, relative to proper keplerian motion. Oddly enough, the one thing it doesn't seem to do a good job at (which RK4 does) is keeping the "argument of periapsis" (the orientation of the ellipse) constant. (Caveat: it's possible I have a bug in my RK4 implementation...)
I have a hunch that when I throw in a third body, this method will start performing worse than RK4, but who knows?
Things I'm trying to understand:
what is this method called. I haven't been able to find mention of it online
why is it so good at keeping the energy of the system constant. Is it some property of the two-body problem?
Is there some way I can adapt Runge-Kutta to have this energy-invariant property?
