Задача двух тел и орбита планеты

Закон всемирного тяготения превращаем в орбиту и проверяем сохранением энергии.

Задача двух тел — движение двух масс под взаимным притяжением; единственный случай гравитации, имеющий точное решение (конические сечения: эллипс, парабола, гипербола).

Сила тяготения

Закон Ньютона: две массы притягиваются с силой F = G·M·m/r², направленной вдоль линии между ними. Для планеты вокруг тяжёлого Солнца удобно считать Солнце неподвижным в начале координат, а ускорение планеты записать векторно: a = -G·M·r⃗/|r|³ (множитель 1/|r|³ даёт и величину 1/r², и направление к центру).

Используем астрономические единицы: расстояние в а.е. (астрономических единицах), время в годах. Тогда G·M_солнца = 4π² — это удобство сильно упрощает числа. Круговая скорость на радиусе 1 а.е.: v = √(GM/r) = 2π а.е./год.

Орбита скоростным Верле

Запустим планету по круговой орбите радиуса 1 а.е. на один год и проверим, замкнётся ли орбита и сохранится ли энергия:

import math
GM = 4*math.pi**2
x, y = 1.0, 0.0
vx, vy = 0.0, 2*math.pi
dt = 0.001
def accel(x, y):
    r = math.hypot(x, y)
    a = -GM/r**3
    return a*x, a*y
ax, ay = accel(x, y)
def energy(x,y,vx,vy):
    r = math.hypot(x,y)
    return 0.5*(vx*vx+vy*vy) - GM/r
E0 = energy(x,y,vx,vy)
for _ in range(int(1.0/dt)):
    x += vx*dt + 0.5*ax*dt*dt
    y += vy*dt + 0.5*ay*dt*dt
    ax2, ay2 = accel(x, y)
    vx += 0.5*(ax+ax2)*dt
    vy += 0.5*(ay+ay2)*dt
    ax, ay = ax2, ay2
print(f"После 1 года: позиция ({x:+.3f}, {y:+.3f}) а.е.")
print(f"(старт был (1.000, 0.000) — орбита замкнулась)")
E1 = energy(x,y,vx,vy)
print(f"Энергия: старт {E0:.4f}, конец {E1:.4f} (дрейф {100*(E1-E0)/abs(E0):+.3f}%)")

Вывод:

После 1 года: позиция (+1.000, -0.000) а.е.
(старт был (1.000, 0.000) — орбита замкнулась)
Энергия: старт -19.7392, конец -19.7392 (дрейф -0.000%)

Земная орбита замкнулась ровно за один год (планета вернулась в (1.000, 0.000)), а энергия не дрейфовала вовсе. Это идеальная проверка: правильная физика плюс симплектический интегратор дают стабильную, замкнутую орбиту. Отрицательная энергия означает связанное состояние — планета не улетит.

Типы орбит и энергия

Знак полной энергии определяет судьбу тела. E < 0 — связанная орбита (эллипс или окружность), тело вечно кружит. E = 0 — параболическая траектория, тело улетает «на грани». E > 0 — гиперболическая, тело пролетает мимо и уходит навсегда. Скорость, разделяющая связанные и свободные орбиты, — вторая космическая (escape velocity), при которой E=0.

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

Почему орбита — эллипс, а не спираль или розетка? Это следствие точного закона 1/r²: только обратный квадрат (и ещё линейная сила) даёт замкнутые орбиты — теорема Бертрана. Малейшее отклонение силы от 1/r² заставляет эллипс медленно поворачиваться (прецессия). Именно крошечная прецессия орбиты Меркурия, не объяснимая ньютоновской гравитацией, стала одним из первых триумфов общей теории относительности. В симуляции же эллипс держится строго замкнутым ровно потому, что мы заложили точный 1/r².

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

Использовать 1/r² для вектора силы напрямую. Нужен множитель r⃗/|r|³, чтобы получить и величину, и направление.
Брать Эйлер для орбиты. Он раскрутит планету по спирали; нужен Верле или leapfrog.
Деление на ноль при r=0. Если тело подходит вплотную к центру, сила взрывается; в реальных кодах добавляют смягчение (softening).

Итоги

Сила тяготения 1/r² даёт ускорение -GM·r⃗/|r|³.
В единицах а.е./год удобно положить GM=4π².
Симплектический Верле даёт замкнутую орбиту без дрейфа энергии.
Знак энергии решает: эллипс (E<0), парабола (E=0) или гипербола (E>0).