The physics and equations of the two- and three-body problem

This section is optional and covers the physics and mathematical details of the two-body problem and three-body problem.

The two-body problem

Consider two point masses, m₁ and m₂, moving within a plane (that is, the 2-dimensional case). Let the position of m₁ be denoted by (x₁, y₁) and m₂ by (x₂, y₂). Because the masses are moving with respect to time, their positional coordinates (x, y) are actually functions of time t, that is (x(t), y(t)). Figure 1 shows such a situation for a particular value of t (a moment frozen in time):

Figure 1 showing values of two point masses

Figure 1

In figure 1, we see that m₁ and m₂ are a distance r apart and are mutually attracted to one another through the force of gravity F₁ and F₂ (whose magnitudes are the same).

This force of gravity (F) is given by Newton's law of universal gravitation:

Equation 1 applying Newton's law of universal gravitation

The force vector F₁ acting on the first point mass m₁ has the following x- and y-components (Fₓ and Fy, respectively), as can be seen here:

Figure 2 shows force vector acting on point masses

Figure 2

Using trigonometry, we see that:

Equations 2 and 3

And, using the larger triangle, shows that:

Equations 4 and 5

If we focus on Fₓ and substitute for F and cos θ, we have:

Equations 6 through 8

Now, Newton’s second law states:

Equation 9

Where F(t) and a(t) are both vector functions of time. If, for now, we consider just the x-direction of m₁, we see that:

Equation 10

Where, in the last expression, we drop the implied function of time variable t.

We now have two expressions for Fₓ, the first from Newton’s second law and the second from Newton’s law of universal gravitation:

Equations 11 and 12

By using the Pythagorean theorem (see figure 2), the distance r between m₁ and m₂ is given by:

Equation 13

Therefore (for m₁ in the x-direction), we have:

Equation 14

Now given the initial positions (x₁, y₁) and (x₂, y₂) at t = 0, we can calculate the acceleration of m₁ due to m₂ in the x-direction using equation 14.

We can use the same reasoning as above to show that the y-component of m₁’s acceleration (due to m₂) is given by:

Equation 15

From m₂’s perspective, we see that the acceleration equations are identical except for a change in sign (indicating opposite force directions) as provided by the swap in positional terms in the numerator:

Equations 16 and 17

Thus, we have the following set of equations, where the numeric subscripts indicate the mass m₁ or m₂:

Equations 18 through 21

Where: Expression for alpha used in equations 18 through 21 (when α = 0, an infinite acceleration is imparted when the two masses exactly collide, which is physically impossible).

Now we have acceleration as a function of each mass’s current position but how do we approximate its velocity and position at some small amount of time h (that is, Δt) later? This is where numerical integration comes in. In particular, we’ll use leapfrog integration (the position Verlet algorithm) in that it’s relatively simple but reasonably accurate and stable.

Consider N celestial masses. Let i indicate one of the masses (i = 1, …, N) and h a small time interval. For the position Verlet algorithm, the ith mass’s next position and velocity values are calculated as follows:

Equations 22 through 24

Where, for a given directional component, xⁱ is the position, vⁱ the velocity, and aⁱ the acceleration of the ith mass.

In order to improve accuracy, on first use, equation 22 is replaced with the following:

Equation 25

After first use, equation 25 is not used and equation 22 is used instead (see Leapfrog Integration).

This information leads to the following algorithm:

  1. Choose a small value for ht).
  2. Choose initial (component-wise) position and velocity values for all masses.
  3. Using equations 18 through 21, calculate initial acceleration values.
  4. Using the a subscript 0 superscript i values from the prior step, use equation 25 to calculate the initial x subscript one-half superscript i values.
  5. Using the previously calculated x subscript one-half superscript i values, use equations 18 through 21 to calculate the a subscript n + one-half superscript i values.
  6. Using the a subscript n + one-half superscript i values from the prior step, calculate the v subscript n+1 superscript i values using equation 23
  7. Using the v subscript n+1 superscript i values from the prior step, calculate the x subscript n+1 superscript i values using equation 24.
  8. Using the x subscript n+1 superscript i values from the prior step, use equation 22 to calculate the x subscript n+one-half superscript i values.
  9. Go to step 5 until the user chooses to stop.

The three-body problem

When there are three masses, any one mass has two forces acting on it from the other two masses. For example, m₁ has the following forces (F₂ and F₃) acting on it:

Figure 3 showing how any one mass has two forces acting on it from two other masses

Figure 3

To start, we note that the net force F₁ acting on m₁ is the sum of F₂ and F₃. That is, F₁ = ma₁ = F₂ + F

Now applying the same trigonometry used in figure 2 to figure 3, the magnitude of the net force F₁ on m₁ can be broken into its x- and y-components as follows:

Equations 26 and 27

Using the red and green triangles in figure 3, we see that:

Equations 28 through 31

From Newton’s law of universal gravitation, F₂ and F₃ can be expressed as:

Equations 32 and 33

Substituting equation 28 through 33 into equations 26 and 27 results in:

Equations 34 and 35

Simplifying 34 and 35 yields:

Equations 36 and 37

Replacing r₂ and r₃ with their Pythagorean formulae yields the following x- and y-component acceleration equation for m₁:

Equations 38 and 39

Where α and β are as follows (for r, see equation 13):

Equations 40 and 41

Applying the same procedure to m₂ and m₃ produces the following system of acceleration equations:

Equations 42 through 47


Equations 48 through 50

As was done in the two-body problem, leapfrog integration (equations 22 through 25) can be used to numerically solve this system of equations (42 through 47) given a set of initial conditions (mass, position, and velocity values for each body) with reasonable accuracy and stability. To allow for speedy high precision, a web worker can be utilized to perform the numeric integration in a thread separate from the main page UI thread.

Related topics

The leapfrog method and other "symplectic" algorithms for integrating Newton’s laws of motion
Leapfrog Integration
Basic 3D graphics using Three.js
The one-body problem
The two-body problem
The three-body problem