Всемирное тяготение и круговая орбита

Планета не падает на звезду только потому, что всё время «промахивается» мимо неё — это и есть орбита.

Закон всемирного тяготения: два тела притягиваются с силой, пропорциональной массам и обратно пропорциональной квадрату расстояния: $F = G\dfrac{m_1 m_2}{r^2}$.

Сила, обратная квадрату расстояния

Гравитация связывает каждую пару тел во Вселенной. Её величина:

$$F = G\frac{m_1 m_2}{r^2},$$

где $G \approx 6.674\cdot10^{-11}$ — гравитационная постоянная, а $r$ — расстояние между центрами. Сила направлена вдоль линии, соединяющей тела, и всегда притягивающая. Зависимость $\frac{1}{r^2}$ означает, что при удвоении расстояния сила слабеет вчетверо. Для движка важно, что это центральная сила: она всегда смотрит на притягивающий центр, и её направление задаётся единичным вектором от тела к центру.

Условие круговой орбиты

Чтобы планета двигалась по кругу радиуса $r$, гравитация должна обеспечивать ровно ту центростремительную силу, что удерживает на окружности: $\frac{m v^2}{r} = G\frac{M m}{r^2}$. Отсюда круговая скорость:

$$v = \sqrt{\frac{G M}{r}}.$$

Слишком медленно — упадёт по спирали внутрь; слишком быстро — улетит по вытянутой траектории. Чтобы числа были удобными, перейдём к «астрономическим» единицам: расстояние в а.е. (астрономических единицах), время в годах, массу — в массах Солнца. В этих единицах $G M_{\odot} = 4\pi^2$, и Земля на $r = 1$ обращается за $1$ год со скоростью $2\pi$ а.е./год. Очень аккуратно.

Симулируем орбиту

Запустим планету по круговой орбите вокруг неподвижной звезды (полу-неявный Эйлер) и проследим её положение в течение года.

import math
GM = 4*math.pi**2            # G·M_солнца в единицах а.е.³/год²
x, y = 1.0, 0.0              # старт на 1 а.е.
vx, vy = 0.0, math.sqrt(GM)  # круговая скорость 2π а.е./год
dt = 0.005
for step in range(201):
    if step % 25 == 0:
        t = step*dt
        r = math.hypot(x, y)
        print(f"t={t:.3f} год  x={x:+.3f}  y={y:+.3f}  r={r:.4f}")
    r = math.hypot(x, y)
    ax = -GM*x/r**3          # компоненты центральной силы /m
    ay = -GM*y/r**3
    vx += ax*dt; vy += ay*dt
    x += vx*dt; y += vy*dt

Вывод:

t=0.000 год  x=+1.000  y=+0.000  r=1.0000
t=0.125 год  x=+0.693  y=+0.706  r=0.9890
t=0.250 год  x=-0.031  y=+0.984  r=0.9848
t=0.375 год  x=-0.736  y=+0.662  r=0.9900
t=0.500 год  x=-1.000  y=-0.061  r=1.0015
t=0.625 год  x=-0.679  y=-0.750  r=1.0122
t=0.750 год  x=+0.027  y=-1.016  r=1.0162
t=0.875 год  x=+0.718  y=-0.712  r=1.0112
t=1.000 год  x=+1.000  y=-0.004  r=1.0000

Планета прошла полный круг: стартовав в $(1, 0)$, она через год вернулась почти точно туда же. Расстояние $r$ держится у $1.0$ с колебанием в пределах процента — это численный артефакт грубого шага, а не реальный дрейф (симплектический метод не даёт орбите раскручиваться). Мы воспроизвели орбиту Земли в миниатюре.

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

Ключ к коду — разложение центральной силы на компоненты. Единичный вектор от планеты к звезде (в начале координат) равен $\left(-\frac{x}{r}, -\frac{y}{r}\right)$, а величина ускорения — $\frac{GM}{r^2}$. Перемножив, получаем $a_x = -\frac{GM\,x}{r^3}$, $a_y = -\frac{GM\,y}{r^3}$ — куб в знаменателе появляется из-за совмещения нормировки ($\frac{1}{r}$) и закона $\frac{1}{r^2}$. Выбор единиц с $GM = 4\pi^2$ — стандартный трюк астрономов, избавляющий от гигантских и крошечных чисел, при которых страдает точность float. Симплектический интегратор здесь критичен: на явном Эйлере орбита за несколько витков раскрутилась бы и планета улетела.

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

• Делить на $r^2$ вместо $r^3$ при покомпонентной силе. Лишняя степень $r$ берётся из нормировки направления; забыв её, получите неверную силу.
• Брать явный Эйлер для орбит. Он раскручивает энергию — орбита «распухает» и тело улетает. Нужен симплектический метод.
• Работать в СИ с огромными числами. Массы $10^{30}$ и расстояния $10^{11}$ режут точность; удобнее астрономические единицы.

Итог

• Тяготение: $F=G\frac{m_1 m_2}{r^2}$ — центральная сила, обратная квадрату расстояния.
• Круговая скорость: $v=\sqrt{GM/r}$; быстрее — улетит, медленнее — упадёт.
• Покомпонентно $a_x=-\frac{GM\,x}{r^3}$ (куб из-за нормировки направления).
• Удобные единицы ($GM=4\pi^2$) и симплектический метод дают стабильную орбиту.