# How to simulate gravitational interaction of bodies in space?

I want to try simulating some bodies in space using the basic formulas I learned in my college mechanics class.

So I know Newton's Law and some basic mechanics equations:

$\vec{F} = G\frac{m_1m_2}{r^2}\hat{r}$

$\vec{v} = \vec{v}_0 + \vec{a}t$

$\vec{s} = \vec{s}_0 + \vec{v}t+\frac{1}{2}\vec{a}t^2$

Now let's say I have a large object with mass $M$ and a small object (say a satellite) with mass $m$. I can assume $m << M$ so that I only need to simulate the effects on the large object.

I also define $n$ to be the current "frame" and $\Delta t$ to be a timestep. I then end up with the following equations:

$\vec{a}_n = \frac{GM}{r_n^2}\hat{r}$

$\vec{v}_{n+1} = \vec{v}_n+\vec{a}_n\Delta t$

$\vec{s}_{n+1} = \vec{s}_n+\vec{v}_n\Delta t+\frac{1}{2}\vec{a}_n(\Delta t)^2$

To me, this seems like everything I would need to simulate gravitational interactions (according to Newtonian mechanics). Adding more bodies just means applying superposition and adding up the vectors accordingly.

My problem is that these equations are really meant to be continuous while I'm applying them discretely. From one frame to the next, the error won't be noticeable, but over say 1000 frames it will make the mechanics slightly off. Of course the smaller my $\Delta t$ the less noticeable everything is, but I'm not sure how small I need it to be. Is there a way to perform "error correction" every dozen frames or so so that everything is in its right place?

• Keep in mind that the two motion equations you have shown only count for constant acceleration $\vec a$ during a step. You timestep must be very small before these are correct to use. – Steeven Feb 13 '17 at 19:46
• @Steeven In that case, would I have to use integration/differential equations? – rcplusplus Feb 13 '17 at 19:51
• Would Computational Science be a better home for this question? – Qmechanic Feb 13 '17 at 19:54
• Note that rudimentary gravitational $n$-body codes are $\mathcal O(n^2)$ operations (I mentioned this recently in a similar problem). You'll need smarter algorithms to improve performance for larger $n$. – Kyle Kanos Feb 13 '17 at 20:22
• The differential versions of the formulae, $\vec v=d\vec s/dt$ and $\vec a=d\vec v/dt$ work for all situations, constant $a$ or not. – Steeven Feb 13 '17 at 20:33

Another method is to plot the mean-square distance $\sqrt{\sum_i |\textbf{x}_i(t)-\textbf{x}_i^\prime(t)|^2}$ over time for two simulations: one with a very small time-step and the other with a larger one. This error will grow exponentially with time $\sim\exp(\lambda t)$ (it's called Lyapunov instability). So you can run a short simulation with an extremely fine time-step and one with a much larger time-step, and measure the associated Lyapunov exponent $\lambda$, and then use that to work out what the overall error will be after some simulation time $t$ and with a time-step $\Delta t$.
• RK4 is not a symplectic integrator, so please do not advocate its use in $n$-body simulations. For some in-house mentions of velocity Verlet, see this search. – Kyle Kanos Feb 13 '17 at 20:18