Задача двух тел численно

Аналитически орбита — эллипс, но тот же результат можно получить и численно, шагая по времени маленькими шажками.

Численное интегрирование орбиты — пошаговое вычисление положения и скорости тела по закону тяготения, без аналитических формул эллипса.

Идея метода

Аналитика даёт красивый ответ — эллипс — но только для двух тел. Для трёх и более тел общего решения нет, и единственный путь — численное моделирование. Принцип тот же, что в курсах «Дифференциальные уравнения» и «Численные методы»: разбиваем время на маленькие шаги $\Delta t$, на каждом вычисляем ускорение от гравитации, обновляем скорость, затем положение. Повторяем тысячи раз.

На каждом шаге: ускорение $\vec{a} = -\dfrac{GM}{r^3}\vec{r}$ (направлено к центру), затем $\vec{v} \mathrel{+}= \vec{a}\,\Delta t$, затем $\vec{r} \mathrel{+}= \vec{v}\,\Delta t$.

import math

# Единицы: а.е., годы, массы Солнца. Тогда GM = 4*pi^2.
GM = 4 * math.pi ** 2

def simulate_orbit(a, e, dt, years):
    """Симуляция орбиты, старт в перигелии. Возвращает финальное (x, y, r)."""
    x = a * (1 - e)          # перигелийное расстояние
    y = 0.0
    r = x
    v = math.sqrt(GM * (2 / r - 1 / a))  # vis-viva в перигелии
    vx, vy = 0.0, v
    steps = int(years / dt)
    for _ in range(steps):
        r = math.sqrt(x * x + y * y)
        ax = -GM * x / r ** 3
        ay = -GM * y / r ** 3
        vx += ax * dt
        vy += ay * dt
        x += vx * dt
        y += vy * dt
    return x, y, math.sqrt(x * x + y * y)

# Орбита a=1 а.е., e=0.2 — за 1 год должна вернуться в перигелий
x, y, r = simulate_orbit(a=1.0, e=0.2, dt=0.0001, years=1.0)
print("После 1 года: x =", round(x, 3), " y =", round(y, 3))
print("Расстояние r =", round(r, 3), "а.е. (перигелий = 0.8)")

Вывод:

После 1 года: x = 0.8  y = -0.0
Расстояние r = 0.8 а.е. (перигелий = 0.8)

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

Выбор единиц делает код чистым: в а.е., годах и массах Солнца гравитационный параметр $GM = 4\pi^2$ (потому что для Земли $T = 1$ год, $a = 1$ а.е., и третий закон $T^2 = \dfrac{4\pi^2}{GM}a^3$ даёт ровно это). Орбита с $a = 1$ имеет период ровно $1$ год по третьему закону, поэтому за год тело возвращается в стартовую точку — перигелий на расстоянии $a(1 - e) = 0.8$ а.е. Симуляция это и подтверждает: $r$ вернулось к $0.8$. Это отличная проверка корректности интегратора.

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

• Слишком большой шаг $\Delta t$ — орбита «раскручивается» или схлопывается из-за накопления ошибки.
• Перепутать знак ускорения — гравитация тянет к центру, поэтому минус.
• Забыть про деление на $r^3$ (а не $r^2$): мы умножаем на компоненту $\vec{r}$, и одна степень $r$ уходит на нормировку направления.

Итог

• Орбиту можно получить численно: ускорение → скорость → положение, тысячи маленьких шагов.
• В единицах а.е./год/масса Солнца $GM = 4\pi^2$ — формулы становятся компактными.
• Возврат в перигелий за период — простой тест правильности симулятора.