Алгоритм для гравитационной задачи N тел

#алгоритм #физика #численные_методы #моделирование

К своему стыду очень слаб в численных методах.
Суть задачи в том, чтобы смоделировать движение тел с течением времени под воздействием
силы гравитации друг на друга в количестве >= 3.

В данный момент проект готов. Работает он следующим образом: за атомарную единицу
времени принимается N микросекунд и считается, что в течении этого времени все тела
двигались линейно, линейно же воздействуя друг на друга.

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

Предположительное решение: для каждой пары тел высчитывать своё значение атома времени
в зависимости от силы притяжения. Но у этого подхода есть ряд проблем, которым мне
не удаётся придумать решение.


Например тела А и В имеют такое влияние друг на друга, что их атом времени составляет
5 секунд. При этом атом тел А и С составляет 1 / 1000000 секунды. Во время течения
атома A-B из-за гравитационного маневра тело А начало стремительное движение в сторону
тела В. В этот момент надо менять атом времени между A и В, иначе они пролетят мимо
друг друга. Получается нужна реактивность: при изменении стейта тела необходимо менять
атомы времени всех его связей с другими телами. Я правильно рассуждаю?
Должна ли при регулировке атома времени учитываться скорость сближения объектов?
Ведь два объекта могут быть достаточно далеко друг от друга чтобы не оказывать практически
никакое влияние друг на друга, но двигаться на сближение с огромной скоростью. Атом
времени между ними может оказаться настолько большим, что они сначала пролетят мимо
друг друга, и только потом высчитается их гравитационное влияние.


Сломал голову придумывая алгоритм реализации. Может быть есть у кого-нибудь пример
решения такой задачи: с радостью бы посмотрел код. А может-быть я услышу волшебное
"почитай про метод Вальфгауна-Штауца" и уйду гуглить в правильном направлении.

P.S. Для того чтобы немного поднять интерес прикреплю GUI существующего проекта.
Практической ценности к вопросу не несет никакого:)

https://makarov-andrey.github.io/n-body-problem/
    

Ответы

Ответ 1

Задача трех тел, есть частный случай задачи Коши для обыкновенных дифференциальных
уравнений: дан вектор начального состояния и функция который вычисляет производную
это вектора по времени (т.н. "функция правой части"). 
(Гуглить "численное решение задачи Коши".)

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

Типичный подход для оценки шага по времени (кроме привлечения физических соображений)
такой:


из начального состояния делаем шаг длинной t запоминаем
результат.  
из начального состояния делаем два шага длинной t/2
запоминаем результат. 
если результаты (1) и (2)  отличаются меньше чем заданный пользователем параметр
"точность на шаге", то принимаем этот шаг если больше то возвращаемся к начальному
состоянию и пытаемся сделать шаг меньшего размера. 
Насколько изменить шаг при пересчете шага, и как выбрать начально приближение для
следующего шага - зависит от выбранного метода решения. Для метода Ньютона, который
ты реализовал, чаще всего, делят или умножают длину шага на 2.


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

Важное замечание: задача трех - задача с разбеганием траекторий. Т.е. невозможно
заранее сказать с какой точностью нужно было посчитать решение для момента t1, чтобы
на его основании можно было получить решение для t2 с заданной точностью, пока не посчитаем
t2. Т.е. невозможно заранее сказать какая "точность на шаге" нужна для того чтобы получить
решение в момент t с заданной точностью. Типичный подход: задачу Коши решают несколько
раз, решают несколько раз, увеличивая "точность на шаге", пока результаты решений в
конечной точке по времени не будут отличаться друг от друга меньше чем заданная точность
интегрального решения.

Если цель изучения численных методов, рекомендую взять готовую библиотеку решения
задачи Коши. 
Для C++/C/fortran рекомендую:
 SUNDIALS https://computation.llnl.gov/projects/sundials 
Несколько самых простых методов, также реализованы в boost https://www.boost.org/doc/libs/1_64_0/libs/numeric/odeint/doc/html/index.html
(она попроще в использовании, но и методы послабее.)

Ключевые слова для поиска библиотек для других языков: метод Ньютона, метод Эйлера,
метод Розенброка, Stiffl, DASKR. 

Если цель изучить численные методы и построение солверов ОДУ, рекомендую книжку Э.Хайрер,
С.Нёрсет Г.Ваннер "Решение обыкновенных дифференциальных уравнений". (Поскольку система
уравнений в задаче трех тел - жесткая, тебе потребуются методы из начала 2-го тома:
"Жесткие и дифференциально-алгебраические задачи",  но без первого тома его читать
бесполезно.)