# Motion of two gravitational bodies?

I'm doing a basic realtime simulation of two bodies, but the orbits are unstable for some reason.

This code is run at every timestep:

rsq = (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y)

a.vx = a.vx - (a.x - b.x) * g * b.m * (1/rsq) * dt
a.vy = a.vy - (a.y - b.y) * g * b.m * (1/rsq) * dt

b.vx = b.vx - (b.x - a.x) * g * a.m * (1/rsq) * dt
b.vy = b.vy - (b.y - a.y) * g * a.m * (1/rsq) * dt

a.x = a.x + a.vx * dt
a.y = a.y + a.vy * dt

b.x = b.x + b.vx * dt
b.y = b.y + b.vy * dt


Where m is mass, x&y are position, vx&vy are velocity, g is the gravitational constant, and dt is the elapsed time.

My bodies orbit eachother but the orbit is unstable and acts as though R is to an exponent other than 2.

Have I missed something in my understanding of this problem or somewhere in my code?

• your acceleration is wrong. Do a unit analysis and you'll get $m^2/s^2$. Essentially, you want the component of acceleration along the x or y direction. Try replacing 1/rsq with (1/rsq)^(3/2)
– Jim
Apr 1, 2014 at 21:18
• Thank you, that seems to have been the problem with my maths. However, my orbit still changes over time, is my timestep too long? Apr 1, 2014 at 21:33
• The other issue is the algorithm is only first order accurate. Do a half-move before and after the accelerate to get to second order.
– webb
Apr 1, 2014 at 21:47
• Would Computational Science be a better home for this question? Apr 1, 2014 at 21:55

## 2 Answers

I'm not sure if it would necessarilly lead to such an instability, as Joce says for a first order method you'll need a really small time step to maintain any accuracy, but at the moment your acceleration is very wrong.

What you effectively have here is ${\bf a} = - \frac{GM}{r^2} {\bf r}$ when you want ${\bf a} = - \frac{GM}{r^2} {\bf \hat{r}}$

(or ${\bf a} = - \frac{GM}{r^3} {\bf r}$ might be easier for you here)

The time discretisation you have chosen is an explicit Euler scheme. In order for it to be stable, you need the time step to be low enough, see e.g. wikipedia.

You could use an implicit method, or increase the order of the method, but in any case there will always be a numerical drift in the orbits, proportional to the numerical accuracy. This can be considered an issue in systems without damping like yours, as it does not enforce the energy conservation principle from which the equations you discretise are derived.

In order to get rid of that drift, you need to use a method that uses the energy of the system as a primary variable, thus the numerical error will be on other quantities (e.g. the period) but the energy will be conserved up to machine accuracy.

Note that some algorithms employ a (dirty) workaround that is commonly accepted: they iterate the state of the system using any (stable) scheme, and then "project" the approximate solution to the closest one with conserved energy, that is they simply correct e.g. the distance between the bodies so that the energy at time $n+1$ is the same as initially. Because this correction is of the order of the numerical accuracy, it is numerically acceptable.

• Thank you. This is very helpful and opens up a lot of room for further exploration! Apr 1, 2014 at 22:09
• You don't necessarily need to compute something that conserves the energy exactly -- you can still use a second order explicit symplectic integrator and get perfectly stable orbits that don't exactly conserve the energy. Implicit methods are going to be a lot trickier to code up.
– webb
Apr 1, 2014 at 22:10
• @webb: You don't need it for stability, as explained. But it is the way you go when you want to have no drift at all in one specific observable. Implicit methods are trickier if fully implicit, here you may linearize the attraction term and have a code about twice the length of the above.
– Joce
Apr 2, 2014 at 7:13
• @Joce that's very true, and the issue if whether you want exact energy conservation or symplecticity (or both) is a discussion I've spent too many hours on. Also, I think Zephyr's answer is relevant as well.
– webb
Apr 2, 2014 at 16:38