¿Cuál es la solución analítica de forma cerrada del problema de dos cuerpos para verificar sus resultados de integración numérica?

Esta es una respuesta complementaria por ahora porque si bien sabemos que una órbita de dos cuerpos se puede reducir a una órbita de un cuerpo alrededor de un potencial central, hacer eso aquí será un poco molesto y creo que el resultado para el cuerpo en el potencial central parece limpiador. Ver también las respuestas a ¿ Se pueden resolver las oscilaciones radiales de una órbita elíptica usando un potencial centrífugo ficticio?

Según este comentario , sé que tuve una discusión en algún lugar de este sitio (o en Astronomy SE ) donde se me explicó por primera vez que las órbitas de Kepler tienen soluciones analíticas que puede escribir para el tiempo en función de la posición , aunque nosotros todavía es necesario utilizar técnicas numéricas (p. ej., el método de Newton) para resolver la posición en función del tiempo. (ver también ¿Cómo lo hicieron Newton y Kepler (en realidad)? )

Si alguien lo encuentra antes que yo, no dude en agregar un enlace aquí, ¡gracias!

Ecuación 27 en la órbita de Kepler de Wikipedia; Propiedades de la ecuación de trayectoria es

dónde$a$es el semieje mayor,$\mu$es el parámetro gravitacional estándar también conocido como el producto$GM$,$e$es la excentricidad y$E$es la anomalía excéntrica .

La relación entre$E$y la verdadera anomalía $\theta = \arctan2(y, x)$es

y resolviendo para$E$:

conectando de nuevo a la primera ecuación (pero sin escribirlo todo):

Probemos una verificación numérica de este increíble resultado. Tenga en cuenta que con$a=1$y$\mu=1$el periodo es$2 \pi$.

El último gráfico en la parte inferior izquierda muestra que el análisis$t(\theta)$Residencia en$\theta$desde una órbita integrada numéricamente coincide con el tiempo utilizado en el cálculo numérico para un$e=0.8$órbita elíptica. Habrá fallas numéricas o singularidades en los puntos finales y para$e=1$pero parece comprobar hacia fuera muy bien!

Escritura de Python:

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

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

e = 0.8
a = 1.0
mu = 1.0
r_peri, r_apo = a*(1.-e), a*(1.+e)
v_peri, v_apo = [np.sqrt(2./r - 1./a) for r in (r_peri, r_apo)]
T = twopi * np.sqrt(a**3/mu)

X0 = np.array([r_peri, 0, 0, v_peri])
X0 = np.array([-r_apo, 0, 0, -v_apo])
times = np.linspace(-T/2., T/2., 1001)

answer, info = ODEint(deriv, X0, times, full_output=True)

x, y = answer[1:-1].T[:2]
theta = np.arctan2(y, x)
E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
t = a * np.sqrt(a/mu) * (E - e * np.sin(E))

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(x, y)
    plt.plot([0], [0], 'ok')
    plt.gca().set_aspect('equal')
    plt.title('y vs. x numerical')
    plt.subplot(2, 1, 2)
    plt.plot(times[1:-1], x)
    plt.plot(times[1:-1], y)
    plt.xlim(-pi, pi)
    plt.title('x(t) and y(t) numerical')
    plt.show()

    plt.subplot(2, 2, 1)
    plt.title('theta(t_numerical)')
    plt.plot(times[1:-1], theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 2)
    plt.title('E_analytic(theta_numerical)')
    plt.plot(E, theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 3)
    plt.title('theta(t_analytic)')
    plt.plot(t, theta)
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.subplot(2, 2, 4)
    plt.title('t_analytic(t_numerical)')
    plt.plot(t, times[1:-1])
    plt.xlim(-pi, pi)
    plt.ylim(-pi, pi)
    plt.gca().set_aspect('equal')
    plt.show()

La distancia desde el foco de atracción de una órbita se puede expresar como una función de la verdadera anomalía (ángulo) dada por$r(\theta)=a\frac{1-e^2}{1+ecos(\theta)}$, dónde$a$es el semieje mayor y$e$es la excentricidad.