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!
In general in this class of integrators you want to do vector updates, meaning that in each position update you adjust the entire position vector and in each velocity update you adjust the entire velocity vector, simultaneous for all particles. Of course in the program this doesn't necessarily need to be done in "vectorized" form, but in any case you should update all components of the full "position" vector prior to making any "velocity" updates.
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.