Законы Кеплера численно

Три закона Кеплера, выведенные не из формул, а из симуляции тяготения.

Третий закон Кеплера: квадрат периода обращения пропорционален кубу большой полуоси орбиты ($T^2 \propto a^3$); в астрономических единицах коэффициент равен единице.

Три закона из одной силы

Кеплер вывел свои законы из наблюдений Тихо Браге, а Ньютон показал, что все три следуют из закона тяготения $1/r^2$. В вычислительной физике мы можем пройти этот путь сами: заложить только силу притяжения и увидеть, как законы Кеплера возникают из симуляции, без всякой подсказки.

• Первый закон: орбита — эллипс, в одном из фокусов которого Солнце.
• Второй закон: радиус-вектор заметает равные площади за равные времена (планета быстрее у Солнца).
• Третий закон: $T^2 \propto a^3$ — связь периода и размера орбиты.

Проверяем третий закон

Запустим планеты на круговых орбитах разного радиуса, измерим период каждой и посчитаем отношение $T^2/a^3$. Если третий закон верен, оно должно быть одинаковым (и равным 1 в наших единицах):

import math
GM = 4*math.pi**2
def orbital_period(r):
    x, y = r, 0.0
    v = math.sqrt(GM/r)
    vx, vy = 0.0, v
    dt = 0.0005
    def acc(x,y):
        rr = math.hypot(x,y); a=-GM/rr**3; return a*x, a*y
    ax, ay = acc(x,y)
    t = 0.0; prev_y = y
    while True:
        x += vx*dt + 0.5*ax*dt*dt
        y += vy*dt + 0.5*ay*dt*dt
        ax2, ay2 = acc(x,y)
        vx += 0.5*(ax+ax2)*dt; vy += 0.5*(ay+ay2)*dt
        ax, ay = ax2, ay2
        t += dt
        if prev_y < 0 <= y and t > 0.1:
            return t
        prev_y = y
print("радиус a (а.е.)   период T (год)   T²/a³")
for r in (0.4, 1.0, 1.5, 5.2):
    T = orbital_period(r)
    print(f"   {r:4.1f}            {T:6.3f}        {T*T/r**3:.3f}")
print("T²/a³ = 1 для всех орбит — третий закон Кеплера.")

Вывод:

радиус a (а.е.)   период T (год)   T²/a³
    0.4             0.253        1.000
    1.0             1.000        1.001
    1.5             1.837        1.000
    5.2            11.858        1.000
T²/a³ = 1 для всех орбит — третий закон Кеплера.

Отношение $T^2/a^3$ вышло равным единице для всех четырёх орбит — третий закон Кеплера родился из симуляции сам собой. Более того, числа совпадают с реальностью: Меркурий (0.4 а.е.) обращается за 0.24 года, Юпитер (5.2 а.е.) — почти за 12 лет. Мы воспроизвели Солнечную систему из одной строчки про силу.

Второй закон и площади

Второй закон — следствие сохранения момента импульса $L = m(x\,v_y - y\,v_x)$. У эллиптической орбиты вблизи Солнца (перигелий) планета движется быстро, вдали (афелий) — медленно, но заметаемая площадь в единицу времени постоянна. В симуляции эллипса величина $L$ сохраняется так же строго, как энергия, — это ещё один обязательный тест.

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

Третий закон легко получить для круга: центростремительное ускорение $v^2/r$ равно гравитационному $GM/r^2$, откуда $v^2 = GM/r$. Период $T = 2\pi r/v$, и подстановка даёт $T^2 = (4\pi^2/GM)\,r^3$ — ровно $T^2 \propto r^3$. Для эллипса вместо $r$ стоит большая полуось $a$, а вывод требует интегрирования, но симуляция получает результат напрямую, просто измеряя период. В этом красота вычислительного подхода: закон, на доказательство которого Ньютону потребовалась новая математика, мы подтверждаем измерением в коде.

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

• Определять период по числу шагов, а не по возврату в исходную точку. Период надо измерять по пересечению старта.
• Путать радиус круга и большую полуось эллипса. В третьем законе фигурирует именно полуось $a$.
• Не проверять момент импульса. Его дрейф — признак ошибки в симуляции эллиптической орбиты.

Итоги

• Три закона Кеплера следуют из одной силы $1/r^2$ и видны в симуляции.
• Третий закон: $T^2/a^3$ одинаково для всех орбит (=1 в единицах а.е./год).
• Второй закон — следствие сохранения момента импульса.
• Симуляция воспроизводит реальные периоды планет Солнечной системы.