Гравитационная задача N тел - классическая задача небесной механики, то есть астрономии, изучающей и вычисляющей траектории небесных тел при помощи законов механики. Выглядит это так: у нас есть N тел в замкнутой системе (те без воздействия внешних сил) с массами и скоростями. Нужно узнать, где они будут находится спустя какое-то время и как они туда попадут (по каким траекториям). Из школьного курса физики мы знаем, что сила притяжения и скорости считаются по следующим формулам:
\[ \overline{r_{t+1, j}} = \overline{r_{t, j}} + \overline{v_{t, j}}dt \quad \overline{v_{t+1, j}} = \overline{v_{t, j}} + \overline{a_{t, j}}dt \quad \overline{F_{t, j}} = m_j\overline{a_{t, j}} \quad \overline{F_{t, j}} = \sum_{i=0, i \neq j}^{N} Gm_im_j\left(\frac{r_i-r_j}{|r_i-r_j|^3} \right) \] Где

\( \overline{r_{t, j}} - \) координаты j тела на t шаге

\( \overline{v_{t, j}} - \) скорость j тела на t шаге

\( \overline{a_{t, j}} - \) ускорение j тела на t шаге

\( m_j - \) масса j тела

\( \overline{F_{t, j}} - \) сила притяжения j тела на t шаге

\( G - \) гравитационная постоянная


Соответственно ускорение можно посчитать для каждого тела как \[ \sum_{i=0, i \neq j}^{N} Gm_i\left(\frac{r_i-r_j}{|r_i-r_j|^3} \right) \]
И вроде бы все хорошо и просто: получаем ускорение из масс и координат, домнажаем на \( dt \), получаем скорость, еще раз на \( dt \) и вот у нас уже новые координаты. Но проблема в том, что, по сути, это дифференциальное уравнение 2 порядка, имеющее вид \[ \frac{dr_j}{dt} = v_j, \qquad \frac{dv_j}{dt} = \sum_{i=0, i \neq j}^{N} Gm_i\left(\frac{r_i-r_j}{|r_i-r_j|^3} \right), \qquad \] не имеющего общего решения, когда количество тел больше 2. Мы не можем взять пару интегралов и сказать сразу, где будут тела в любой момент времени. Только считать изменения "на ходу" последовательно для всех значений времени. Но и этого мало, из-за того, что мы переходим от непрерывной модели поведения к дискретной, у нас появляется погрешность в зависимости от \( dt \). Чем оно меньше, соответственно, тем меньше погрешность.

Данный численный метод решения задачи - метод Эйлера, имеет довольно таки большую погрешность - \( O(dt^2) \) на каждом шаге и \( O(dt) \) в целом. Особенно сильно это проявляется при сильном сближении тел.
Куда лучше справляется метод Рунге-Кутты. Он немного сложнее в осознании и вычисление, зато дает куда более маленькую погрешность - \( O(dt^5) \) на каждом шаге и \( O(dt^4) \) в целом.
Выглядит он как: \[ \Large \overline{k}_{1_{v_{(t+1, j)}}} = \overline{a}(\overline{r}_{(t, j)}) \qquad \overline{k}_{2_{v_{(t+1, j)}}} = \overline{a}(\overline{r}_{(t, j)} + \overline{k}_{1_{r_{(t+1, j)}}} \frac{dt}{2}) \qquad \] \[ \Large \overline{k}_{3_{v_{(t+1, j)}}} = \overline{a}(\overline{r}_{(t, j)} + \overline{k}_{2_{r_{(t+1, j)}}} \frac{dt}{2}) \qquad \overline{k}_{4_{v_{(t+1, j)}}} = \overline{a}(\overline{r}_{(t, j)} + \overline{k}_{3_{r_{(t+1, j)}}} dt) \] \[ \Large \overline{k}_{1_{r_{(t+1, j)}}} = \overline{v}_{(t, j)} \qquad \overline{k}_{2_{r_{(t+1, j)}}} = \overline{v}_{(t, j)} \overline{k}_{1_{v_{(t+1, j)}}} \frac{dt}{2} \qquad \] \[ \Large \overline{k}_{3_{r_{(t+1, j)}}} = \overline{v}_{(t, j)} \overline{k}_{1_{v_{(t+1, j)}}} \frac{dt}{2} \qquad \overline{k}_{4_{r_{(t+1, j)}}} = \overline{v}_{(t, j)} \overline{k}_{1_{v_{(t+1, j)}}} dt \] И итоговые изменения коэффицентов: \[ \Large \overline{v}_{(t+1, j)} = \overline{v}_{(t, j)} + \frac{dt}{6} \left( \overline{k}_{1_{v_{(t+1, j)}}} + 2\overline{k}_{2_{v_{(t+1, j)}}} + 2\overline{k}_{3_{v_{(t+1, j)}}} + \overline{k}_{4_{v_{(t+1, j)}}} \right) \] \[ \Large \overline{r}_{(t+1, j)} = \overline{r}_{(t, j)} + \frac{dt}{6} \left( \overline{k}_{1_{r_{(t+1, j)}}} + 2\overline{k}_{2_{r_{(t+1, j)}}} + 2\overline{k}_{3_{r_{(t+1, j)}}} + \overline{k}_{4_{r_{(t+1, j)}}} \right) \]
Я запрограммировал оба метода. Слева решение через Эйлера, справа - Рунге-Кутта. Программа сразу задает 3 тела их скорости и массы, которые (в идеале) движутся по замкнутому символу бесконечности. Однако из-за неточных начальных условий (единственные, которые удалось найти), тела движутся не совсем по замкнутой траектории. Но близкой к ней =)
Подобные начальные значения и траектории называются хореографией. Существует достаточно большое найденных хореографий, но, к сожалению, найти и понять начальные условия весьма проблематично, а данной теме посвящают отдельные, весьма немаленькие главы в учебниках. Один из самых популярных примеров - сборник К. Симо "Современные проблемы хаоса и нелинейности"
Построенными более сложными хореографиями можно полюбоваться тут

Но это все были слова и цифры. Куда интереснее попробовать самому позапускать планеты друг в друга =)
Верхняя панель отвечает за левую рисовалку, нижняя за правую. Стрелочки указывают скорость. Массы, гравитационную постоянную и шаг \( dt \) можно изменять вручную.


Открыть на весь экран

Как мы можем наблюдать, метод Рунге-Кутты работает медленнее. Так же масса тел для хореографии нужна на порядок ниже. Однако если изменить ускорение хотя бы в несколько раз, то при расчете Эйлером тела практически сразу разлетаются, в то время как движение при помощи Рунге-Кутты дает хоть и погрешность, но не на столько большую что бы тела разлетались мгновенно. При ускорении в 5 раз траектории намного плавнее, а разлет и отклонение от восьмерок намного ниже. Как минимум тела не разлетаются во все стороны.

Спасибо за внимание.
Проект выполнил: студент 2 курса 176 группы ПМИ,
Тропин Федор