I'm writing code simulating the 2-dimensional motion of two massive bodies with gravitational fields.

The bodies' masses are known and I have a gravitational force equation. I know from that force I can get a differential equation for coordinates. I know that I once I solve this equation I will get the coordinates. I will need to make up some initial position and some initial velocity.

I'd like to end up with a numeric solver for the ordinal differential equation for coordinates to get the formulas that I can write in code. Could someone break down how from laws and initial conditions we get to the formulas that calculate x and y at time t?

  • \$\begingroup\$ This is exactly the gravitational 2-body problem. The derivations on Wikipedia's page on general two-body problems seem to be what you're looking for. \$\endgroup\$ – Anko Apr 26 '14 at 16:34
  • \$\begingroup\$ How accurate do you want to get? I wrote some code recently to follow Keplerian Orbits given a set of orbital elements (which can be calculated from initial position & velocity), but it's not entirely simple. Calculating elliptical/parabolic/hyberbolic anomaly requires some iterative solving with Newton's method. This will correctly show the acceleration of the bodies at periapsis and deceleration at apoapsis, and is more stable than the integration-based approaches in the answers below - it's just a matter of whether this makes it worth the heftier code for your particular needs. \$\endgroup\$ – DMGregory May 27 '14 at 22:12

To a first approximation, this behaves just like any other mass-under-a-force problem: since the gravitational force is F=G*m1*m2/r^2, the acceleration on body 1 is a1=G*m2/r^2 and that on body 2 is a2=G*m1/r^2 — note that these accellerations are in the direction of the other body.

To keep this pseudo-code relatively minimal, I'm going to presume that there's already a well-defined Vector class with all of the relevant operations (addition, subtraction, scalar multiplication and normalization) available, and all of my vector variables will be prefaced with a 'v'; e.g., vVel1 and vPos2.

Given this setup, then a first approximation is pretty straightforward — but there's one huge caveat that I'll cover at the end.

First of all, there's a straightforward initialization step:

vPos1 = vInitPos1; vVel1 = vInitVel1;
vPos2 = vInitPos2; vVel2 = vInitVel2;

Then on every update, we first find the vector from the first body to the second — of course, the vector from the second body to the first is just the negative of this. While we're at it, we pre-scale it by the gravitational constant and the inverse squared distance between the two vectors; these values will be the same for both accelleration computations, so now's a good time to get them out of the way. Note that since we'll normalize vDeltaPos by dividing it by its length, and then we divide it again by its squared length, we could actually just divide by the length cubed - but vector normalization routines often have special code for numerical stability, so it's generally best to do this 'the hard way'.

vDeltaPos = vPos2-vPos1;
radSq = vDeltaPos.LengthSquared();
vScaledDeltaPos = vDeltaPos*(GRAVITATIONAL_CONSTANT/radSq);

Now with this in hand, the acceleration for the two bodies is easy; body 1's accel just points towards body two, scaled by body 2's mass, while body 2's accel points towards body 1, scaled by that body's mass:

vAccel1 = mass2*vScaledDeltaPos;
vAccel2 = (-mass1)*vScaledDeltaPos;

With the accelerations in hand, then the update step is just two iterations of the Euler method, updating velocity from acceleration and updating position from velocity:

vVel1 += deltaTime*vAccel1;
vVel2 += deltaTime*vAccel2;
vPos1 += deltaTime*vVel1;
vPos2 += deltaTime*vVel2;

And that's all there is to it! But... there's a huge caveat here. Actually, two huge caveats. The first is that Euler's method is as inaccurate as a differential equation solver can get while still being useful; there are methods that give you much better accuracy for a given delta-Time, but they're a lot more complicated. (If you're curious about these, the magic names to start with are Runge-Kutta).

The second, though, is a little more subtle, and the code above luckily accounts for it. You may notice that I update the velocity before I update the position; this means that the position update uses the updated velocity rather than the 'old' velocity. This minor tweak helps to correct for a problem with the naive Euler methods: they don't necessarily conserve the energy of the system, so that even a configuration that looks like it should be a stable mutual orbit can wind up with the two bodies flying off towards infinity. Using the updated velocity is a very simple version of a so-called semi-implicit Euler method, which helps with energy conservation (at least to the approximation provided by the timestep), but if you want guaranteed long-term conservation - i.e., stable orbits over long timescales - then you'll want to look even deeper down the rabbit hole, into methods that are guaranteed to be energy-preserving.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.