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

Это дополнительный ответ на данный момент, потому что, хотя мы знаем , что орбита двух тел может быть уменьшена до орбиты одного тела вокруг центрального потенциала, выполнение этого здесь будет немного отвлекать, и я думаю, что результат для одного тела в центральном потенциале выглядит очиститель. См. также ответы на вопрос Можно ли решить радиальные колебания эллиптической орбиты с помощью фиктивного центробежного потенциала?

Согласно этому комментарию, я знаю, что где-то на этом сайте (или в Astronomy SE ) у меня была дискуссия, где мне впервые объяснили, что орбиты Кеплера имеют аналитические решения, которые вы можете записать для времени как функции положения , хотя мы по-прежнему необходимо использовать численные методы (например, метод Ньютона) для определения положения как функции времени. (см. также Как Ньютон и Кеплер (на самом деле) сделали это? )

Если кто-то найдет его раньше меня, пожалуйста, не стесняйтесь добавить сюда ссылку, спасибо!

Уравнение 27 орбиты Кеплера из Википедии; Свойства уравнения траектории

где$a$большая полуось,$\mu$стандартный гравитационный параметр, также известный как произведение$GM$,$e$это эксцентриситет и$E$эксцентрическая аномалия .

Отношения между$E$и настоящая аномалия $\theta = \arctan2(y, x)$является

и решение для$E$:

подключаясь к первому уравнению (но не записывая все):

Давайте попробуем численно проверить этот удивительный результат. Обратите внимание, что с$a=1$и$\mu=1$период$2 \pi$.

Последний график внизу слева показывает, что аналитический$t(\theta)$на основе$\theta$с численно интегрированной орбиты соответствует времени, используемому в численном расчете для$e=0.8$эллиптическая орбита. Будут числовые сбои или сингулярности на конечных точках и для$e=1$но, кажется, проверить красиво!

Скрипт 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()

Расстояние от центра притяжения орбиты может быть выражено как функция истинной аномалии (угла), определяемой выражением$r(\theta)=a\frac{1-e^2}{1+ecos(\theta)}$, где$a$является большой полуосью и$e$является эксцентриситет.