How to solve the two-body problem in the ECI frame through numerical integration?

Question

I need to know how to solve two-body problem by solving a system of first order equation derived from the equation below.

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

How do I go about this, and how would I then use this to output a trajectory in MATLAB?

Answer 1

There are several ways to do this. The easiest and most straightforward is to break it into two sets by including velocity as a variable, and solve together.

Instead of a single second order differential equation

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

We can solve the following pair of first order differential equations in parallel

$$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$$ $$\dot{\mathbf{r}} = \mathbf{v}$$

using various simple methods including standard libraries or homegrown implementations of Runge-Kutta including my favorite simple RK4/5 with variable step size.

It's wonderful and really educational to code that once for yourself and appreciate the task first before using standard libraries from then on.

For more on errors using differential equation solvers and for a way to test them yourself, see my question in Math SE: need to understand better ODE solution accuracy vs numerical precision (which I should have asked in Computational Science SE)

For programming purposes you can write $\mathbf{r} / r^3 $ as r * (r**2).sum()**-1.5

From this answer you can see a 2D implementation using not only the monopole $1/r^2$ gravity term but the additional quadrupole $J_2$ term for Earth's oblate shape and field. For more on that see this answer to Trouble deriving rectangular components of acceleration of satellite in orbit around Earth with J2 consideration.

See also a similar solution in this answer to How to determine an orbit of a satellite for a collision detection?.

Also, you can consider going unit-less by using $\mu = 1$ and period $T = 2 \pi$ and semi-major axis $a=1$. Some integrators will handle that a little better.

Blue line is correct, red line has J2 ten times its real value to show exaggerated apsidal precession for fun.

def deriv(X, t):

    x, v = X.reshape(2, -1)

    acc0  = -GMe * x * ((x**2).sum())**-1.5
    acc2  = -1.5 * GMe * J2 * Re**2 * x * ((x**2).sum())**-2.5

    return np.hstack([v, acc0 + acc2])


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# David Hammen's nice table https://physics.stackexchange.com/a/141981/83380
# See http://www.iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/51-groten.pdf
# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere

GMe = 3.98600418E+14  # m^3 s^-2
J2e = 1.08262545E-03  # unitless
Re  = 6378136.3 # meters

X0 = np.hstack([6778000.0, 0.0, 0.0, 10000.])  # x, y, vx, vy

time = np.arange(0, 300001, 100)

J2 = J2e  # correct J2
answerJ2, info = ODEint(deriv, X0, time, full_output=True)

J2 = 10*J2e # 10x larger J2
answer10xJ2, info = ODEint(deriv, X0, time, full_output=True)

if True:
    plt.figure()
    x, y = answerJ2.T[:2]
    plt.plot(x, y, '-b')
    x, y = answer10xJ2.T[:2]
    plt.plot(x, y, '-r')
    plt.plot([0], [0], 'or')
    plt.gca().set_aspect('equal')
    plt.show()