¿Cómo resolver el problema de los dos cuerpos en el marco ECI a través de la integración numérica?

Hay varias formas de hacer esto. La más fácil y directa es dividirlo en dos conjuntos incluyendo la velocidad como variable y resolver juntos.

En lugar de una única ecuación diferencial de segundo orden

Podemos resolver el siguiente par de ecuaciones diferenciales de primer orden en paralelo

usando varios métodos simples que incluyen bibliotecas estándar o implementaciones propias de Runge-Kutta, incluido mi RK4/5 simple favorito con tamaño de paso variable .

Es maravilloso y realmente educativo codificar eso una vez por sí mismo y apreciar la tarea primero antes de usar bibliotecas estándar a partir de ese momento.

Para obtener más información sobre los errores que usan los solucionadores de ecuaciones diferenciales y una forma de probarlos usted mismo, consulte mi pregunta en Math SE: necesito comprender mejor la precisión de la solución ODE frente a la precisión numérica (que debería haber preguntado en Computational Science SE )

Para fines de programación se puede escribir$\mathbf{r} / r^3$comor * (r**2).sum()**-1.5

A partir de esta respuesta, puede ver una implementación 2D utilizando no solo el monopolo$1/r^2$término de gravedad pero el cuadrupolo adicional$J_2$término para la forma achatada y el campo de la Tierra. Para obtener más información, consulte esta respuesta a Problemas para derivar componentes rectangulares de aceleración del satélite en órbita alrededor de la Tierra con consideración J2.

Consulte también una solución similar en esta respuesta a ¿Cómo determinar la órbita de un satélite para la detección de colisiones? .

Además, puede considerar ir sin unidades usando$\mu = 1$y periodo$T = 2 \pi$y semi-eje mayor$a=1$. Algunos integradores manejarán eso un poco mejor.

La línea azul es correcta, la línea roja tiene J2 diez veces su valor real para mostrar una precesión absidal exagerada por diversión.

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