Two-Body Numerical Solution in an Inertial Frame

Two-Body Numerical Solution in an Inertial Frame#

The equations of motion for the two-body problem in an inertial frame are given by (21), repeated here:

\[\begin{split}\begin{aligned} \ddot{X}_1 &= G m_2 \frac{X_2 - X_1}{r^3} & \ddot{Y}_1 &= G m_2 \frac{Y_2 - Y_1}{r^3} & \ddot{Z}_1 &= G m_2 \frac{Z_2 - Z_1}{r^3}\\ \ddot{X}_2 &= G m_1 \frac{X_1 - X_2}{r^3} & \ddot{Y}_2 &= G m_1 \frac{Y_1 - Y_2}{r^3} & \ddot{Z}_2 &= G m_1 \frac{Z_1 - Z_2}{r^3} \end{aligned}\end{split}\]

The initial positions of the masses are:

\[\begin{split}\begin{aligned} \vector{R}_{1,0} &= \left(0 \uvec{I} + 0 \uvec{J} + 0 \uvec{K}\right) \mathrm{km} \\ \vector{R}_{2,0} &= \left(3000 \uvec{I} + 0 \uvec{J} + 0 \uvec{K}\right) \mathrm{km} \end{aligned}\end{split}\]

and the initial velocities are:

\[\begin{split}\begin{aligned} \dot{\vector{R}}_{1,0} &= \left(10 \uvec{I} + 20 \uvec{J} + 30 \uvec{K}\right) \mathrm{km/s} \\ \dot{\vector{R}}_{2,0} &= \left(0 \uvec{I} + 40 \uvec{J} + 0 \uvec{K}\right) \mathrm{km/s} \end{aligned}\end{split}\]

where the subscript \(0\) indicates the initial position and velocity. These are the 12 components of the initial state vector, which we will give the symbol \(y_0\).

The two masses are equal, \(m_1 = m_2 = 1.0\times 10^{26} \text{ kg}\) and Newton’s gravitational constant is \(G = 6.67430\times 10^{-2} \text{ km}^3 \text{ kg}^{-1} \text{ s}^{-2}\).

The State Vector#

In the numerical solution, the state vector is stored in an array. Arrays in programming are data structures that are optimized to store numbers. The order that the components appear in the state vector isn’t important, only that we keep track of the order and don’t forget it. For the sake of choosing something, we will set the positions followed by the velocities:

(22)#\[y = \left[X_1\ Y_1\ Z_1\ X_2\ Y_2\ Z_2\ \dot{X}_1\ \dot{Y}_1\ \dot{Z_1}\ \dot{X}_2\ \dot{Y}_2\ \dot{Z}_2\right]\]

The array is indicated by using square brackets, \([\dots]\), and lists all the components of the array inside the brackets. Depending on the programming language you’re using, there are different ways to create arrays.

In the following code samples we use arrays to store the initial positions and velocities of both masses and then construct the state vector.

import numpy as np

G = 6.67430e-20  # km**3/(kg * s**2)
m_1 = m_2 = 1.0e26  # kg

R_1_0 = np.array((0, 0, 0))  # km
R_2_0 = np.array((3000, 0, 0))  # km
dotR_1_0 = np.array((10, 20, 30))  # km/s
dotR_2_0 = np.array((0, 40, 0))  # km/s

y_0 = np.hstack((R_1_0, R_2_0, dotR_1_0, dotR_2_0))


G = 6.67430E-20;  % km^3/(kg * s^2)
m1 = 1.0E26;  % kg
m2 = 1.0E26;  % kg

R10 = [0 0 0];  % km
R20 = [3000 0 0];  % km
dotR10 = [10 20 30];  % km/s
dotR20 = [0 40 0];  % km/s

y0 = [R10 R20 dotR10 dotR20];

These code samples first set the constants in the problem, \(G\) and \(m_1 = m_2\). Then, they create the initial position and velocity arrays. Finally, the arrays are stuck together into the initial state vector, y_0.

Transforming the System of Ordinary Differential Equations#

The system of ordinary differential equations (ODEs) in (20) or (21) is a second order system. However, it can be transformed to a system of first-order ODEs by recognizing the acceleration as the first derivative of the velocity. Solving a first-order system is much simpler than solving a second-order system.

We know that the system of equations requires the

(23)#\[\dot{\vector{y}} = f\left(t, \vector{y}\right)\]

where \(\vector{y}\) is the state vector, \(t\) is time, and \(f\) is a function that depends on the current time and the current state vector.

The left hand side of (23) is the time-derivative of the state vector, which we find by taking the time derivative of each element of the state vector:

(24)#\[\dot{\vector{y}} = \left[\dot{X}_1\ \dot{Y}_1\ \dot{Z}_1\ \dot{X}_2\ \dot{Y}_2\ \dot{Z}_2\ \ddot{X}_1\ \ddot{Y}_1\ \ddot{Z}_1\ \ddot{X}_2\ \ddot{Y}_2\ \ddot{Z}_2\right]\]

In words, the left side of (23) is another array, where the first six elements are the velocity components (the first derivative of position) and the second six are the acceleration (the first derivative of velocity).

Solution Algorithm#

We now have enough information to start to solve the problem. The first step is to calculate the initial acceleration using (21). This can be done for one direction at a time, as shown in the following code:

X_1 = y_0[0]
Y_1 = y_0[1]
Z_1 = y_0[2]
X_2 = y_0[3]
Y_2 = y_0[4]
Z_2 = y_0[5]

r = np.sqrt((X_2 - X_1) ** 2 + (Y_2 - Y_1) ** 2 + (Z_2 - Z_1) ** 2)

ddotX_1 = G * m_2 * (X_2 - X_1) / r**3
ddotY_1 = G * m_2 * (Y_2 - Y_1) / r**3
ddotZ_1 = G * m_2 * (Z_2 - Z_1) / r**3
ddotX_2 = -G * m_1 * (X_2 - X_1) / r**3
ddotY_2 = -G * m_1 * (Y_2 - Y_1) / r**3
ddotZ_2 = -G * m_1 * (Z_2 - Z_1) / r**3


X_1 = y0(1);
Y_1 = y0(2);
Z_1 = y0(3);
X_2 = y0(4);
Y_2 = y0(5);
Z_2 = y0(6);

r = sqrt((X_2 - X_1).^2 + (Y_2 - Y_1).^2 + (Z_2 - Z_1).^2);

ddotX_1 = G .* m2 .* (X_2 - X_1) ./ r.^3;
ddotY_1 = G .* m2 .* (Y_2 - Y_1) ./ r.^3;
ddotZ_1 = G .* m2 .* (Z_2 - Z_1) ./ r.^3;
ddotX_2 = -G .* m1 .* (X_2 - X_1) ./ r.^3;
ddotY_2 = -G .* m1 .* (Y_2 - Y_1) ./ r.^3;
ddotZ_2 = -G .* m1 .* (Z_2 - Z_1) ./ r.^3;

This code first retrieves the position components from the state vector. Then it calculates each component of acceleration for the two masses.

This code is pretty long, and we’ve created a bunch of variables to keep track of. Fortunately, there are simpler ways to approach this calculation. The next code samples show how to take advantage of the nature of array computations:

R_1 = y_0[:3]
R_2 = y_0[3:6]

r = np.sqrt(np.sum(np.square(R_2 - R_1)))
ddot = G * (R_2 - R_1) / r**3
ddotR_1_0 = m_2 * ddot
ddotR_2_0 = -m_1 * ddot


R1 = y0(1:3);
R2 = y0(4:6);

r = norm(R2 - R1);
ddot = G .* (R2 - R1) ./ r.^3;
ddotR10 = m2 .* ddot;
ddotR20 = -m1 .* ddot;

In this code, you retrieve the position of each mass as an array instead of into a single variable. Then, using array functions, you compute the distance and the accelerations.

Now that you’ve calculated the acceleration, the velocity and position at the next time can be found by integration:

(25)#\[\begin{split}\begin{aligned} \dot{R}_1 &= \int_0^t \ddot{R}_1 dt + \dot{R}_{1,0}\\ \dot{R}_2 &= \int_0^t \ddot{R}_2 dt + \dot{R}_{2,0}\\ R_1 &= \int_0^t \dot{R}_1 dt + R_{1,0}\\ R_2 &= \int_0^t \dot{R}_2 dt + R_{2,0}\\ \end{aligned}\end{split}\]

Assuming that the acceleration and velocity are constant over some time interval \(\Delta t\), the integrals simplify to the product of the term and \(\Delta t\). The choice of \(\Delta t\) is up to you. Smaller values of \(\Delta t\) will be more accurate but take longer to compute.

Let’s choose \(\Delta t = 1\text{ s}\). Then, to compute the state vector at the next time step:

Delta_t = 1  # s
dotR_1_1 = ddotR_1_0 * Delta_t + dotR_1_0
dotR_2_1 = ddotR_2_0 * Delta_t + dotR_2_0

R_1_1 = dotR_1_0 * Delta_t + R_1_0
R_2_1 = dotR_2_0 * Delta_t + R_2_0


Delta_t = 1;  % s
dotR11 = ddotR10 .* Delta_t + dotR10;
dotR21 = ddotR20 .* Delta_t + dotR20;

R11 = dotR10 .* Delta_t + R10;
R21 = dotR20 .* Delta_t + R20;

However, it would be inefficient to do this by hand and there are more accurate methods available. I don’t see a reason to re-implement standard functions, so we are going to use the functions built-in to SciPy or Matlab, depending on which software you’re using.

Numerical Solution Using Pre-Built Libraries#

In SciPy, the function is called solve_ivp. In Matlab, the function is called ode45.

First, you need to start by importing the appropriate Python libraries. In Matlab, all the functions you need are built-in.

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

The two imports from matplotlib will plot the solution of the problem.

Then, we need to define the function that describes the motion of our system. This function needs to compute the left hand side of the equation above (\(\dot{\vector{y}}\)) and return it to the solver, so that the solver can calculate the value of the state vector at time \(t + \Delta t\).

The function takes the time \(t\) and current value of the state vector \(y\) as inputs. In this way, we can use the current state vector to compute the current acceleration.

Inside the function, we use the values in the state vector to fill the ydot vector. Then, we return ydot back to the solver.

def absolute_motion(t, y):
    """Calculate the motion of a two-body system in an inertial reference frame.

    The state vector ``y`` should be in the order:

    1. Coordinates of $m_1$
    2. Coordinates of $m_2$
    3. Velocity components of $m_1$
    4. Velocity components of $m_2$
    """
    # Get the six coordinates for m_1 and m_2 from the state vector
    R_1 = y[:3]
    R_2 = y[3:6]

    # Fill the derivative vector with zeros
    ydot = np.zeros_like(y)

    # Set the first 6 elements of the derivative equal to the last
    # 6 elements of the state vector, which are the velocities
    ydot[:6] = y[6:]

    # Calculate the acceleration terms and fill them in to the rest
    # of the derivative array
    r = np.sqrt(np.sum(np.square(R_2 - R_1)))
    ddot = G * (R_2 - R_1) / r**3
    ddotR_1 = m_2 * ddot
    ddotR_2 = -m_1 * ddot

    ydot[6:9] = ddotR_1
    ydot[9:] = ddotR_2
    return ydot



function dydt = absolute_motion(~,y)
    % Calculate the motion of a two-body system in an inertial reference frame.
    %
    % The state vector ``y`` should be in the order:
    %
    % 1. Coordinates of $m_1$
    % 2. Coordinates of $m_2$
    % 3. Velocity components of $m_1$
    % 4. Velocity components of $m_2$

    R1 = y(1:3);
    R2 = y(4:6);

    dotR1 = y(7:9);
    dotR2 = y(10:12);

    r = norm(R2 - R1);

    ddotR1 = G .* m2 .* (R2 - R1) ./ r.^3;
    ddotR2 = G .* m1 .* (R1 - R2) ./ r.^3;

    dydt = [dotR1; dotR2; ddotR1; ddotR2];

end % absolute_motion

With the function defined, we can call solve_ivp() or ode45(). We need to tell it the function it should solve, the beginning and end times, the initial state vector, and then some information to help control the output.

Once the solver finishes, the solution is stored in sol.y in Python or just y in Matlab. Each column of sol.y corresponds to a single timestep and each row corresponds to one of the state variables. It is more convenient to work with the transpose of this array, so we do that and define y. In Matlab, the solution already has each timestep in a row.

Then we extract the position and velocity of each mass as a function of time, and compute the barycenter (the center of gravity of the system).

t_0 = 0  # seconds
t_f = 480  # seconds
t_points = np.linspace(t_0, t_f, 1000)

sol = solve_ivp(absolute_motion, [t_0, t_f], y_0, t_eval=t_points)

y = sol.y.T
R_1 = y[:, :3]  # km
R_2 = y[:, 3:6]  # km
V_1 = y[:, 6:9]  # km/s
V_2 = y[:, 9:]  # km/s
barycenter = (m_1 * R_1 + m_2 * R_2) / (m_1 + m_2)  # km


t0 = 0;
tf = 480;
[~,y] = ode45(@absolute_motion, [t0 tf], y0);

R1 = y(:, 1:3); % km
R2 = y(:, 4:6); % km
barycenter = (m1 .* R1 + m2 .* R2) ./ (m1 + m2);

Finally, we construct some plots of the situation. In Fig. 9, we are plotting the absolute motion of each of the two masses as well as the barycenter. Notice that the barycenter moves in a straight line. We will discuss this further in the next section.

The two masses spiral around the barycenter. One way to imagine this system is as the Earth and the Moon viewed as though you were sitting on the Sun. The Earth and Moon would move through space, and they would appear to be orbiting around each other. If you observed them for a short enough time, their motion would appear to be in a straight line.

Fig. 9 The motion of two bodies subject to mutual gravitational attraction, viewed from an external inertial frame.#

Another way to view this system is by setting the barycenter to be the origin of the coordinate system, as shown in Fig. 10. Remember that since the barycenter is moving with constant velocity, it is allowed to be used as an inertial reference frame. This is kind of like sitting above the barycenter of the Earth-Moon system. You would see them orbit around the barycenter, and the orbits would be ellipses.

Fig. 10 The motion of two bodies subject to mutual gravitational attraction, viewed from an inertial frame attached to the system barycenter. In this reference frame, the orbits of \(m_1\) and \(m_2\) appear to be ellipses with the barycenter at one of the foci.#

Fig. 11 fixes the coordinate system on the first mass and plots the motion of the barycenter and the second mass relative to the position of the first mass. This is kind of like sitting on the Earth and watching the Moon go around. Notice that the barycenter of the system also orbits around the first mass in this reference frame.

Fig. 11 The motion of two bodies subject to mutual gravitational attraction, viewed from a non-inertial frame attached to \(m_1\). In this reference frame, the orbits of the barycenter and \(m_2\) appear to be ellipses with \(m_1\) at one of the foci.#

Interestingly, the equations for this solution are symmetric. We can reverse the roles of \(m_1\) and \(m_2\) and have exactly the same plot as Fig. 11. This means that sitting on the Moon watching the Earth orbit is the same as sitting on the Earth watching the Moon orbit. Just like the Moon has phases when viewed from Earth, the Earth has phases when viewed from the Moon!

The code to generate the plots is shown below.

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(R_1[:, 0], R_1[:, 1], R_1[:, 2], label="m_1")
ax.plot(R_2[:, 0], R_2[:, 1], R_2[:, 2], label="m_2")
ax.plot(barycenter[:, 0], barycenter[:, 1], barycenter[:, 2], label="COG")
ax.legend()

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
R1_rel_COG = R_1 - barycenter
R2_rel_COG = R_2 - barycenter
ax.plot(R1_rel_COG[:, 0], R1_rel_COG[:, 1], R1_rel_COG[:, 2], label="m_1")
ax.plot(R2_rel_COG[:, 0], R2_rel_COG[:, 1], R2_rel_COG[:, 2], label="m_2")
ax.plot(0, 0, 0, "ro", label="COG")
ax.legend()

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
R2_rel_R1 = R_2 - R_1
COG_rel_R1 = barycenter - R_1
ax.plot(R2_rel_R1[:, 0], R2_rel_R1[:, 1], R2_rel_R1[:, 2], label="m_2")
ax.plot(COG_rel_R1[:, 0], COG_rel_R1[:, 1], COG_rel_R1[:, 2], label="COG")
ax.plot(0, 0, 0, "ro", label="m_1")
ax.legend()

output

return

function output
    % Plot the output in three graphs
    figure(1)
    title('Motion of two masses in an inertial frame')
    hold on
    plot3(R1(:, 1), R1(:, 2), R1(:, 3), '-r')
    plot3(R2(:, 1), R2(:, 2), R2(:, 3), '-b')
    plot3(barycenter(:, 1), barycenter(:, 2), barycenter(:, 3), '-g')
    axis_settings

    figure(2)
    title('Motion of two masses relative to their barycenter')
    hold on
    R1_rel_COG = R1 - barycenter;
    R2_rel_COG = R2 - barycenter;
    plot3(R1_rel_COG(:, 1), R1_rel_COG(:, 2), R1_rel_COG(:, 3), '-r')
    plot3(R2_rel_COG(:, 1), R2_rel_COG(:, 2), R2_rel_COG(:, 3), '-b')
    plot3(0, 0, 0, 'og')
    axis_settings

    figure(3)
    title('Motion of mass 2 relative to mass 1')
    hold on
    R2_rel_R1 = R2 - R1;
    COG_rel_R1 = barycenter - R1;
    plot3(R2_rel_R1(:, 1), R2_rel_R1(:, 2), R2_rel_R1(:, 3), '-b')
    plot3(COG_rel_R1(:, 1), COG_rel_R1(:, 2), COG_rel_R1(:, 3), '-g')
    plot3(0, 0, 0, 'or')
    axis_settings

    function axis_settings
        view([2,4,1.2])
        grid on
        axis equal
        xlabel('X (km)')
        ylabel('Y (km)')
        zlabel('Z (km)')
    end % axis_settings

end % output
Contents