Введение в проблему численного моделирования гравитационных систем

Задача трех тел — классическая и одна из самых красивых проблем в физике и математике. Суть ее заключается в определении движения трех точечных масс, взаимодействующих согласно закону всемирного тяготения Ньютона. В отличие от задачи двух тел, которая имеет аналитическое решение и описывается простыми кеплеровыми орбитами, задача трех тел не имеет общего аналитического решения в квадратурах. Для ее исследования ученые, инженеры и разработчики игр вынуждены прибегать к методам численного интегрирования дифференциальных уравнений.

Среди множества численных алгоритмов адаптивные методы семейства Рунге-Кутты, и в частности метод Дорманда-Принца (более известный как RK45 или scipy.integrate.solve_ivp с методом 'DOP853' или 'RK45'), считаются золотым стандартом. Они автоматически подбирают шаг интегрирования, обеспечивая высокую точность и экономию вычислительных ресурсов. Однако разработчики, реализующие симуляции гравитационных систем с нуля или использующие стандартные библиотеки, часто сталкиваются с одной и той же неприятной аномалией: график движения со временем начинает «смещаться», орбиты «плывут», а полная энергия системы начинает дрейфовать, что приводит к физически некорректным результатам (например, планеты улетают в бесконечность или падают на центральное светило).

В этой статье мы подробно разберем, почему стандартный солвер RK45 при наивной реализации выдает «смещенный» или деградирующий график для задачи трех тел, какие скрытые математические ловушки подстерегают программиста и как превратить нестабильную симуляцию в надежный физический движок.

Анатомия задачи: уравнения движения и вектор состояния

Прежде чем искать баги в солвере, давайте вспомним математическую формулу, которую мы пытаемся решить. Система из трех тел с массами $m_1, m_2, m_3$ и радиус-векторами $\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3$ описывается вторым законом Ньютона в сочетании с законом гравитации:

m_i * d^2(r_i)/dt^2 = sum_{j \neq i} (G * m_i * m_j * (r_j - r_i)) / ||r_j - r_i||^3

Для применения численных методов вроде RK45 систему дифференциальных уравнений второго порядка необходимо привести к системе уравнений первого порядка. Для этого мы вводим вектор состояния, состоящий из координат и скоростей всех трех тел. Размерность нашего вектора состояния равна $3 \text{ тела} \times 3 \text{ измерения} \times 2 (\text{координата и скорость}) = 18$ переменных.

Функция правых частей (в терминологии Python часто называемая fun(t, y)) возвращает производные: для координат производной является скорость, а для скорости — ускорение, вычисленное по формуле гравитационного притяжения. На первый взгляд код выглядит элементарно:

import numpy as np

def three_body_ode(t, y, masses):
    G = 1.0
    # y содержит [r1, r2, r3, v1, v2, v3], где каждый ri и vi — это векторы длины 3
    r = y[:9].reshape((3, 3))
    v = y[9:].reshape((3, 3))
    
    dr = v
    dv = np.zeros_scipy_like_or_zero_matrix((3, 3)) if False else np.zeros((3, 3))
    
    for i in range(3):
        for j in range(3):
            if i != j:
                diff = r[j] - r[i]
                dist_cubed = np.linalg.norm(diff)**3
                dv[i] += G * masses[j] * diff / dist_cubed
                
    return np.concatenate([dr.flatten(), dv.flatten()])

Имея эту функцию, программист подключает scipy.integrate.solve_ivp, задает начальные условия (например, для знаменитой фигуры «восьмерки» Мурукаме или ретроградной стабильной орбиты Сухого) и запускает расчет. И тут график начинает расходиться.

Ловушка №1: Ошибки округления и жесткость при сближениях

Первая и самая очевидная причина появления «смещенного» графика заключается в природе метода RK45 и поведении гравитационной силы. Метод Рунге-Кутты 4/5 порядка является методом с фиксированным или адаптивным шагом общего назначения. Он прекрасно справляется с гладкими нелинейными функциями, но гравитационная сила содержит член $\|\mathbf{r}_j - \r_i\|^3$ в знаменателе.

Когда два тела в процессе орбитального движения сближаются на очень близкое расстояние, значение $\|\mathbf{r}_j - \r_i\|$ стремится к нулю. В этот момент:

  • Ускорение стремится к бесконечности.
  • Производные меняются с колоссальной скоростью.
  • Адаптивный алгоритм RK45 пытается уменьшить шаг интегрирования до микроскопических значений, чтобы сохранить заданную погрешность (tolerance).

Однако из-за ограничений чисел с плавающей запятой двойной точности (IEEE 754, float64) при уменьшении шага начинают доминировать ошибки округления. Солвер теряет точность вычислений, накапливает погрешность, и после прохождения точки сближения траектория тел «отскакивает» совсем не по той физической траектории, которую предсказывает аналитика. На графике это выглядит как резкое смещение фазы орбиты или параллельный сдвиг всей траектории движения центра масс.

Ловушка №2: Нарушение законов сохранения (Энергия и Импульс)

Это самый коварный архитектурный подвох. Общая теория относительности и классическая механика утверждают, что изолированная система трех тел должна сохранять:

  1. Полный импульс системы (центр масс должен оставаться в покое или двигаться равномерно и прямолинейно).
  2. Полный момент импульса (угловой моментум).
  3. Полную механическую энергию (сумма кинетической и потенциальной энергий должна быть константой).

Стандартные методы Рунге-Кутты (включая RK45) не являются консервативными. Это означает, что они не встроены с учетом законов сохранения на геометрическом уровне. На каждом малом шаге алгоритма допускается крошечная погрешность округления и аппроксимации. В диссипативных или кратковременных задачах этой погрешностью пренебрегают. Но в консервативных системах, таких как небесная механика, эти микроскопические ошибки накапливаются со временем.

Если вы построите график полной энергии системы $\mathcal{E}(t) = K(t) + U(t)$ для стандартного вызова solve_ivp(..., method='RK45'), вы заметите систематический дрейф энергии (energy drift) — она будет либо медленно расти, либо убывать. Физически это означает, что ваша виртуальная Вселенная либо получает энергию из ниоткуда, либо теряет ее. Именно этот дрейф энергии и заставляет орбиты «смещаться» на графике: радиусы орбит меняются, периоды вращения искажаются, и через несколько оборотов смоделированная система полностью рассыпается.

Как исправить смещение: Практические рекомендации и код

Теперь, когда мы диагностировали проблемы (ошибки вблизи близких сближений и дрейф энергии из-за неконсервативности RK45), давайте рассмотрим проверенные методы исправления ситуации.

1. Смягчение потенциала (Softening)

Если ваша задача допускает прохождение тел на очень близкие расстояния (или столкновения), стандартный закон Ньютона приводит к числовым взрывам. В астрофизическом моделировании для предотвращения этого используют параметр смягчения (softening length, $\epsilon$):

# Модифицированный расчет расстояния в симуляциях N-тел
dist_cubed = (np.linalg.norm(diff)2 + epsilon2)**1.5

Это убирает сингулярность при $\mathbf{r}_j = \mathbf{r}_i$, спасая RK45 от падения шага интегрирования в ноль и накопления катастрофических ошибок округления.

2. Ужесточение допусков солвера (Tolerances)

По умолчанию scipy.integrate.solve_ivp использует значения rtol=1e-3 и atol=1e-6. Для гравитационных задач эти значения катастрофически велики. Они подходят для грубых оценок, но для построения красивых и точных графиков задачи трех тел требования к точности должны быть на несколько порядков выше:

from scipy.integrate import solve_ivp

# Правильные допуски для прецизионного моделирования орбиты
sol = solve_ivp(
    fun=three_body_ode,
    t_span=(0, 50),
    y0=initial_state,
    method='RK45',
    args=(masses,),
    rtol=1e-10,
    atol=1e-12,
    t_eval=np.linspace(0, 50, 1000)
)

Установка rtol=1e-10 и atol=1e-12 заставит алгоритм RK45 делать гораздо более мелкие шаги, что существенно снизит уровень накопления ошибок, хотя и увеличит время выполнения скрипта.

3. Переход на симплектические интеграторы (Альтернатива RK45)

Если ваша главная цель — долгосрочное стабильное моделирование орбитальной динамики без смещения графиков и дрейфа энергии, метод RK45 вообще не лучший выбор. Для гамильтоновых систем созданы специальные симплектические интеграторы (например, алгоритм Верле — Velocity Verlet, методы Фреста или симплектические интеграторы высших порядков, такие как Yoshida 4th/6th order).

Симплектические методы геометрически сохраняют фазовый объем и структуру пространства движения, благодаря чему энергия системы колеблется около истинного значения, но не дрейфует со временем. График траектории при этом остается стабильным даже на миллионах шагов.

«Если вы пишете симуляцию Солнечной системы или сложного гравитационного балета трех тел, помните правило: RK45 хорош для точности на коротких интервалах, но симплектический интегратор Верле спасает физику на длинных дистанциях». Убедитесь, что вы контролируете параметры шага и не доверяете слепо дефолтным настройкам библиотек общего назначения.

Заключение

Появление «смещенного» или деградирующего графика при решении задачи трех тел с помощью солвера RK45 — это классическая техническая проблема, с которой сталкиваются все разработчики на стыке программирования и математического моделирования. Причины кроются в фундаментальных свойствах алгоритма Рунге-Кутты: он не сохраняет полную энергию системы автоматически, а при сближении тел сталкивается с математическими сингулярностями и ограничениями плавающей запятой float64.

Чтобы получить корректный график и физически достоверные результаты, необходимо:

  • Критически ужесточить параметры точности солвера (rtol до $10^{-10}$ и выше, atol до $10^{-12}$).
  • Использовать параметр смягчения (softening) при высокой вероятности тесных сближений масс.
  • Для долгосрочных симуляций осознанно переходить на симплектические методы интегрирования (например, Velocity Verlet), созданные специально для гамильтоновых систем.

Учет этих нюансов позволит вам избавиться от артефактов смещения и наслаждаться безупречной гармонией небесной механики в вашем коде.