Как решить задачу двух тел в системе ECI с помощью численного интегрирования?

Есть несколько способов сделать это. Самый простой и простой способ — разбить его на два набора, включив скорость в качестве переменной и решить вместе.

Вместо одного дифференциального уравнения второго порядка

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

Мы можем решить следующую пару дифференциальных уравнений первого порядка параллельно

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

используя различные простые методы, включая стандартные библиотеки или самодельные реализации Рунге-Кутта, включая мой любимый простой RK4/5 с переменным размером шага .

Замечательно и по-настоящему познавательно написать код один раз для себя и сначала оценить задачу, прежде чем использовать стандартные библиотеки.

Чтобы узнать больше об ошибках с использованием решателей дифференциальных уравнений и о том, как проверить их самостоятельно, см. мой вопрос в Math SE: нужно понять лучшую точность решения ODE по сравнению с численной точностью (что я должен был задать в Computational Science SE )

Для целей программирования вы можете написать$\mathbf{r} / r^3$какr * (r**2).sum()**-1.5

Из этого ответа вы можете увидеть 2D-реализацию, использующую не только монополь$1/r^2$гравитационный член, но дополнительный квадруполь$J_2$термин для обозначения сплюснутой формы и поля Земли. Подробнее об этом см. в этом ответе на проблему с получением прямоугольных составляющих ускорения спутника на орбите вокруг Земли с учетом J2.

См. также аналогичное решение в этом ответе на Как определить орбиту спутника для обнаружения столкновения? .

Кроме того, вы можете подумать о том, чтобы обойтись без единиц измерения, используя$\mu = 1$и период$T = 2 \pi$и большая полуось$a=1$. Некоторые интеграторы справятся с этим немного лучше.

Синяя линия верна, красная линия имеет J2 в десять раз больше реального значения, чтобы показать преувеличенную прецессию апсид для развлечения.

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()