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):
.png)
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:
.png)
The force vector F₁ acting on the first point massm₁ has the following x- and y-components (Fₓ and Fy, respectively), as can be seen here:
.png)
Figure 2
Using trigonometry, we see that:
.png)
And, using the larger triangle, shows that:
.png)
If we focus on Fₓ and substitute for F and cos θ, we have:
.png)
Now, Newton’s second law states:
.png)
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:
.png)
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:
.png)
By using the Pythagorean theorem (see figure 2), the distance r between m₁ and m₂ is given by:
.png)
Therefore (for m₁ in the x-direction), we have:
.png)
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:
.png)
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:
.png)
Thus, we have the following set of equations, where the numeric subscripts indicate the mass m₁ or m₂:
.png)
Where: .png)
(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:
.png)
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:
.png)
After first use, equation 25 is not used and equation 22 is used instead (see Leapfrog Integration).
This information leads to the following algorithm:
- Choose a small value for h (Δt).
- Choose initial (component-wise) position and velocity values for all masses.
- Using equations 18 through 21, calculate initial acceleration values.
- Using the
.png)
values from the prior step, use equation 25 to calculate the initial .png)
values. - Using the previously calculated values, use equations 18 through 21 to calculate the
.png)
values. - Using the values from the prior step, calculate the
.png)
values using equation 23 - Using the values from the prior step, calculate the
.png)
values using equation 24. - Using the values from the prior step, use equation 22 to calculate the
.png)
values. - 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:
.png)
Figure 3
To start, we note that the net forceF₁ acting on m₁ is the sum of F₂ and F₃. That is, F₁ = m₁a₁ = 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:
.png)
Using the red and green triangles in figure 3, we see that:
.png)
From Newton’s law of universal gravitation, F₂ and F₃ can be expressed as:
.png)
Substituting equation 28 through 33 into equations 26 and 27 results in:
.png)
Simplifying 34 and 35 yields:
.png)
Replacing r₂ and r₃ with their Pythagorean formulae yields the following x- and y-component acceleration equation for m₁:
.png)
Where α and β are as follows (for r, see equation 13):
.png)
Applying the same procedure to m₂ and m₃ produces the following system of acceleration equations:
.png)
Where:
.png)
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