Solving the Two Body Problem Numerically in Python

numerical methodsnumerical-calculusordinary differential equationspython

I am following Example 2.2 in Orbital Mechanics for Engineering Students.

The equations of motion of the two bodies in inertial space are given as follows:
$$\ddot{X_1} = Gm_2\frac{X_2-X_1}{r^3},\space \ddot{Y_1} = Gm_2\frac{Y_2-Y_1}{r^3},\space \ddot{Z_1} = Gm_2\frac{Z_2-Z_1}{r^3}$$
$$\ddot{X_2} = Gm_1\frac{X_1-X_2}{r^3},\space \ddot{Y_2} = Gm_1\frac{Y_1-Y_2}{r^3},\space \ddot{Z_2} = Gm_1\frac{Z_1-Z_2}{r^3}$$

The objective is to use specified initial conditions and a numerical integrator to solve the equations numerically.

The initial conditions are:
$$y_0=[X_1^0, Y_1^0, Z_1^0, X_2^0, Y_2^0, Z_2^0, \dot{X}_1^0, \dot{Y}_1^0, \dot{Z}_1^0, \dot{X}_2^0, \dot{Y}_2^0, \dot{Z}_2^0]$$

I am using odeint from scipy in Python to integrate.
The solver requires a system of first-order ordinary differential equations. To do this, I am introducing several auxiliary variables.

$$a_1=X_1,\space a_2=Y_1,\space a_3=Z_1,\space a_4=X_2,\space a_5=Y_2,\space a_6=Z_3$$
$$a_7=\dot{X}_1, a_8=\dot{Y}_1, a_9=\dot{Z}_1, a_{10}=\dot{X}_2, a_{11}=\dot{Y}_2, a_{12}=\dot{Z}_2$$

There are now 12 first-order ODEs, two of which are:
$$\dot{a}_1=a_7,\space \dot{a}_2=a_8,\space \dot{a}_3=a_9,\space \dot{a}_4=a_{10},\space \dot{a}_5=a_{11},\space \dot{a}_6=a_{12},\space $$

$$\dot{a}_7=Gm_2\frac{a_4-a_1}{r^3},\space \dot{a}_8=Gm_2\frac{a_5-a_2}{r^3},\space \dot{a}_9=Gm_2\frac{a_6-a_3}{r^3}$$

$$\dot{a}_{10}=Gm_1\frac{a_1-a_4}{r^3},\space \dot{a}_{11}=Gm_1\frac{a_2-a_5}{r^3},\space \dot{a}_{12}=Gm_1\frac{a_3-a_6}{r^3}$$

Implementing this in Python:

# time array
time = np.arange(0, 480, 0.1)

# body m1 initial conditions
m1 = 10e26  # mass (kg)
r10 = np.array([0, 0, 0])  # initial position (m)
v10 = np.array([10e3, 20e3, 30e3])  # initial velocity (m/s)

# body m2 initial conditions
m2 = 10e26  # mass (kg)
r20 = np.array([3000e3, 0, 0])  # initial position (m)
v20 = np.array([0, 40e3, 0])  # initial velocity (m/s)

# magnitude of position vector from r1 to r2
r_mag = np.linalg.norm(r20 - r10)

# initial conditions
# [X1 (0), Y1 (1), Z1 (2), X2 (3), Y2 (4), Z2 (5), VX1 (6), VY1 (7), VZ1 (8), VX2 (9), VY2 (10), VZ2 (11)]
y0 = np.concatenate((r10, r20, v10, v20))

def two_body_eqm(_y, _t, _G, _m1, _m2, _r_mag):
    c0 = _y[6:12]
    c1 = _G * _m2 * ((_y[3:6] - _y[:3]) / np.power(_r_mag, 3))
    c2 = _G * _m1 * ((_y[:3] - _y[3:6]) / np.power(_r_mag, 3))
    return np.concatenate((c0, c1, c2))

y = odeint(two_body_eqm, y0, time, args=(G, m1, m2, r_mag))

When I plot the results the motion does not appear to be correct as indicated by the figures in the textbook.

Plot result

Textbook result

Can someone give me an indication of what is wrong with my solution?

Additionally, to solve for the relative motion, I am using the following:

for yk in y:
    # approximate y at time t

    # extract inertial positions of body 1 and body 2
    r1 = yk[:3]
    r2 = yk[3:6]

    # determine position of centre of mass
    rg = ((m1 * r1) + (m2 * r2)) / (m1 + m2)

    # position vector from m1 to m2
    r12 = r2 - r1

    # position vector from m1 to g
    r1g = rg - r1

    # position vector from g to m1
    rg1 = r1 - rg

    # position vector from g to m2
    rg2 = r2 - rg

Thanks for your time.

Updated plot

Updated code:

def two_body_eqm(_y, _t, _G, _m1, _m2):
    r_mag = np.linalg.norm(_y[3:6] - _y[:3])
    c0 = _y[6:12]
    c1 = _G * _m2 * ((_y[3:6] - _y[:3]) / np.power(r_mag, 3))
    c2 = _G * _m1 * ((_y[:3] - _y[3:6]) / np.power(r_mag, 3))
    return np.concatenate((c0, c1, c2))

y = odeint(two_body_eqm, y0, time, args=(G, m1, m2))

Best Answer

The most glaring problem is that you are not solving the gravitational equations for the 2-body problem. For that to be the case, you would have to compute r_mag at the start of two_body_eqm as the distance norm(_y[3:6] - _y[:3]) between the two bodies.

Related Question