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

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?

• +1 We don't have a Matlab tag but we do have a Python tag and they are similar looking enough languages (superficially) that hopefully the tag will attract programmers of both. I've added MathJax formatting standard for most Stack Exchange sites. – uhoh Nov 13 '20 at 13:08
• I've posted an answer using free, open source, friendly, popular, fun and well-documented Python. Perhaps it might help sway you into trying to move towards it and away from expensive and black-box MATLAB. – uhoh Nov 13 '20 at 13:43

## 1 Answer

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

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.

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 it's 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(, , 'or')
plt.gca().set_aspect('equal')
plt.show()