How to best modify Leap Frog algorithm to account for binary motion of particles.

I am trying to write a program that can output data that can be plotted to show the orbital motion of two bodies. I am starting with the LeapFrog algorithm for a single particle orbiting about a fixed mass at the origin, so I am only having to calculate parameters (position, velocity, force) for one particle at present. I am using the following order of operations:

Increment the position vector by $0.5 \delta t \times \textbf{v}$

Increment t by $o.5\delta t$

Find the force $\textbf{f}$

Increment $\textbf{v}$ by $0.5\delta t \times \frac{1}{m}\textbf{f}$

Increment $\textbf{x}$ by $0.5\delta t \times \textbf{v}$

Increment t by $\delta t$

Now I would like to add motion of the central body, presently fixed at the origin, as well. This will mean computing the force on this body for certain time steps and incrementing its position and velocity, as above.

However O am struggling to see how best to interleave the calculations for the second body (i.e. computing the force on this body, increasing its position and/or velocity) with computing these quantities for the first body, as above.

Theoretically, I suppose I could apply one iteration of the leapfrog to the planet, and then once to the star, but I feel that interleaving these steps will provide a result with less error.

Any ideas would be very welcome!

• I suggest the Runge-Kutta method . It is based on the Simpson-rule which already gives quite precise results and is not too difficult to be implemented. – Peter Jan 4 '18 at 19:09
• @Peter Thank you for your reply, however unfortunately my question was slightly different (I will edit the title of my post to remove ambiguity!) I was meaning about the best way to add the effects of the motion of the second body within the framework of the algorithm I have! – Meep Jan 4 '18 at 19:27
• @Peter Actually something along these lines is very important in something like celestial mechanics where numerically preserving conservation laws is essential. – Ian Jan 4 '18 at 19:36
• I think your scheme is not consistent: you only increment $v$ by half a force the whole time. Is that a typo? Anyway, in general you want to do everything in lockstep: update the vector position then the vector velocity, doing each of those steps somewhere between 1 and 3 times depending on the integrator you're using. See en.wikipedia.org/wiki/Symplectic_integrator – Ian Jan 4 '18 at 19:39
• I explained a non-Leapfrogged Verlet variant in stackoverflow.com/a/23671054/3088138. In the question and in the answer with jsfiddle you can see how the simultaneous treatment of the coupled system has to work. – LutzL Jan 5 '18 at 13:53

If you have some kind of time scale separation that you want to use as a numerical optimization (e.g. a star that moves but much more slowly than its planets, motivating you to only move the star every $N$ time steps), you will need to adjust the scheme as a whole in order to preserve your numerical conservation law(s).
Incidentally, your scheme is wrong because it doesn't do a full time step update on $v$, but that's easy enough to fix.