Уравнение Кеплера и итерация Ньютона

Чтобы узнать, где планета на орбите в данный момент, надо решить трансцендентное уравнение Кеплера — итерацией.

Уравнение Кеплера: $M = E - e\sin E$, где $M$ — средняя аномалия (растёт равномерно со временем), $E$ — эксцентрическая аномалия, $e$ — эксцентриситет.

Проблема: время не равно углу

По второму закону Кеплера планета движется неравномерно: быстро у Солнца, медленно вдали. Поэтому угол, на который она повернулась, не растёт равномерно со временем. Зато равномерно растёт средняя аномалия $M$ — воображаемый «равномерный» угол, $M = \dfrac{2\pi}{T}(t - t_0)$. Чтобы найти настоящее положение, нужно сначала из $M$ получить эксцентрическую аномалию $E$ через уравнение Кеплера. Но это уравнение трансцендентное — $E$ входит и снаружи, и под синусом, аналитически не решается.

Решение: метод Ньютона

Уравнение $M = E - e\sin E$ решают итеративно. Метод Ньютона ищет корень функции $f(E) = E - e\sin E - M$, последовательно уточняя оценку: $E_{n+1} = E_n - \dfrac{f(E_n)}{f'(E_n)}$, где $f'(E) = 1 - e\cos E$. Сходится за несколько шагов.

import math

def solve_kepler(M, e, tol=1e-12, max_iter=100):
    """Решить M = E - e*sin(E) методом Ньютона. M, E в радианах."""
    E = M if e < 0.8 else math.pi  # начальное приближение
    for _ in range(max_iter):
        f = E - e * math.sin(E) - M
        f_prime = 1 - e * math.cos(E)
        dE = f / f_prime
        E -= dE
        if abs(dE) < tol:
            break
    return E

# Планета: средняя аномалия 60 градусов, эксцентриситет 0.2
M = math.radians(60.0)
e = 0.2
E = solve_kepler(M, e)
print("Эксцентрическая аномалия E =", round(math.degrees(E), 4), "град")

# Истинная аномалия (реальный угол от перигелия)
nu = 2 * math.atan2(math.sqrt(1 + e) * math.sin(E / 2),
                    math.sqrt(1 - e) * math.cos(E / 2))
print("Истинная аномалия nu =", round(math.degrees(nu), 4), "град")

Вывод:

Эксцентрическая аномалия E = 70.8233 град
Истинная аномалия nu = 82.0958 град

Как работает под капотом

Метод Ньютона имеет квадратичную сходимость: число верных цифр примерно удваивается на каждом шаге. Для $e = 0.2$ хватает $4$–$5$ итераций, чтобы дойти до точности $10^{-12}$. Заметьте логику начального приближения: при умеренном эксцентриситете берём $E_0 = M$ (хорошая стартовая точка), а при большом $e$ ($\gt 0.8$) — $E_0 = \pi$, потому что у сильно вытянутых орбит наивный старт может увести итерацию не туда. Это тот же метод Ньютона, что в курсе «Численные методы», только применённый к конкретному астрономическому уравнению.

Получив $E$, мы находим истинную аномалию $\nu$ — настоящий угол планеты от перигелия — через формулу с половинными углами. А расстояние до Солнца: $r = a(1 - e\cos E)$.

Частые ошибки

• Путать три аномалии: среднюю $M$ (равномерную), эксцентрическую $E$ и истинную $\nu$ (реальную).
• Забыть, что метод Ньютона работает в радианах, а не в градусах.
• Брать плохое начальное приближение для очень вытянутых орбит — итерация может разойтись.

Итог

• Уравнение Кеплера $M = E - e\sin E$ связывает время (через $M$) с положением (через $E$).
• Оно трансцендентно — решается итерацией Ньютона за несколько шагов.
• Из $E$ получают истинную аномалию $\nu$ и расстояние $r = a(1 - e\cos E)$.