Редактирование: Оптимазация вычисления сил в задаче N тел с помощью алгоритмов древовидных структур

Оптимазация вычисления сил в задаче N тел с помощью древовидных алгоритмов
==Постановка задачи==
* Исследовать эволюцию плоского облака частиц, в гравитационном центральном поле;
* Выделить параметры, влияющие на количество образовавшихся устойчивых скоплений частиц - кластеров;
* Оптимизировать методы моделирования системы.

=Решение=
==Математическая модель==
Сформулируем задачу в терминах динамики частиц.
Пусть имеется N частиц. В общем случае их распределение в расчетной области определяет функция плотности распределения p(x,y). Из соображений физической реалистичности эта функция должна быть симметричной. Вся совокупность частиц образует так называемый протопланетный диск. Обращаю отдельное внимание, что формально, это может оказаться и не диском с геометрической точки зрения. Скорость задается как сумма двух компонент: скорости вращения и хаотической компоненты скорости. Первая отвечает за движение каждой частицы вокруг центрального тела, вторая – привносит в систему хаотичность. Системы без второй компоненты являются недостаточно физически реалистичными. В центре протопланетного диска (в начале координат) находится неподвижный источник гравитационного поля. Взаимодействие частиц между собой задаются потенциалом.
==Имплементация==
* Частицы – объекты класса с полями скорости, положения, массы и пр.
* Метод интегрирования – leapfrog:
<math>\dot{x}(t+dt) = \dot{x}(t)-x(t)dt</math><br />
<math>x(t+dt) = x(t)+\dot{x}(t+dt)dt</math><br />
* Начальное распределение частиц: равномерное внутри диска;
* Скорости вращения – первая космическая скорость, направлена нормально к радиус-вектору;
* Хаотические компоненты скорости выражаются в виде случайного угла поворота скоростей вращения.
* Потенциал взаимодействия: <math> f(r,\dot{r})= \gamma \frac{m^2}{a^2} \left[ \left( \frac{a}{r} \right)^{13} \left(1- \beta \frac{\dot{r}}{r} \right) - \left(\frac{a}{r}\right)^2 \right] </math>. Данный потенциал помог исследовать образование системы Луна-Земля  в результате гравитационного коллапса. Подробнее можно прочитать в статье: [http://tm.spbstu.ru/Origin_of_the_Moon._New_Concept]
==Оптимизация==
Задача оптимизации появляется из-за чересчур высокой асимптотической сложности обсчета взаимодействий каждой частицы с каждой. Пусть всего в системе N частиц, тогда асимптотическая сложность вычисление взаимодействий в системе - <math>O(n^2)</math>. <br/>
Первым шагом оптимизации был отказ от взаимодействий между всеми частицами. Простейший способ – разбить плоскость на одинаковые квадраты и считать взаимодействие с соседними квадратами. Проблема данного подхода в сложности работы с дальнодействующими потенциалами, так как слишком маленький размер сетки приводил бы к большой ошибке. А кластеры являются настолько плотными структурами, что при сильных неоднородностях в распределении частиц сложность алгоритма опять бы оказывалась асимптотически равна <math>O(n^2)</math>.<br/>

Более правильный способ построения сеток и вычисления на них предложили Josh Barnes и Piet Hut. Этот алгоритм получил название алгоритма Барнс-Хата [ссылка на статью]. Визуализацию можно посмотреть здесь [http://tm.spbstu.ru/Метод_Барнса-Хата]. Основная идея метода строится на том, что взаимодействие частицы с удаленной группой частиц можно заменить на взаимодействие этой частицы с центром масс частиц. Pic.1 Этот подход возможно реализовать благодаря использованию древовидных структур данных. Заявляется о снижении асимптотической сложности до <math>O(n \cdot \log(n))</math>. В статье говорится о том, что эта оценка верна при большом количестве частиц и близком к равномерному распределении. Образование планет подразумевает много плотных скоплений. Барнс-Хат, в отличии от предыдущего метода описывает дальнодействие с сильно меньшей ошибкой благодаря иерархическому подходу к построению сетки. <br/>
Проблема в неоптимальном расчете этого взаимодействия. Предположим, что кластеры состоят из m и l частиц и находятся далеко друг от друга. Самый точный алгоритм совершит <math>O(m \cdot l))</math>  операций, так как будет находить взаимодействие всех частиц со всеми. Второй описанный здесь алгоритм не позволит вычислить взаимодействие вовсе, из-за большого расстояния. Алгоритм Барнс-Хата на само вычисление силы потратит <math>O(m +  l))</math> операций, так как каждая частица в системе один раз будет использоваться в расчетах (взаимодействие ее с центром масс). Идеальный алгоритм в этом случае должен был бы потратить <math>O(1))</math>, так должно быть рассмотрено взаимодействие кластера с кластером, то есть нужно посчитать одно взаимодействие между центрами масс. Эта идея была реализована в следующей статье, там же приведено сравнение быстродействия. Pic.2 Картинки оттуда [http://mech.spbstu.ru/images/6/6c/GalimovSbornik-LeZakharovKrivtsov-ProblemyEvolutsii-2011-v05.pdf] <br/>
Таким образом, вместо того, чтобы решать задачу N-тел методом «вычислить взаимодействие между всевозможными парами частиц» предлагается решать ее следующим образом:
#Построить представление частиц в виде древовидной структуры.
#Используя древовидную структуру подсчитать приближенные значения сил взаимодействия. <br/>
Необходимо отметить, что в такой формулировке сложность заключается именно в первом пункте, так как алгоритм выполнения второго, независимо от построенного представления будет одинаковым.


==Древовидные структуры==
Древовидной структурой называются ненаправленные, связанные, ацикличные, простые графы, для которых выполняются следующие свойства:
#Каждая вершина содержит некоторое (возможное нулевое) количество точек и соединена с одним родителем и некоторым количеством потомков (возможно нулевым);
#В структуре существует одна вершина без родителя – она называется корнем дерева;
#Каждая точка принадлежит хотя бы одной вершине;
#Множество точек, принадлежащих данной вершине представимо в виде пересечения множеств точек, принадлежащих потомкам этой вершины.
Главным плюсом древовидных структур является возможность задания очень гибких правил построения в сочетании с высокой производительностью. Возвращаясь к самой примитивной сетке – ее можно представить в виде древовидной структуры с одним корнем и <math>(\frac{d}{h})^2 </math> потомков (где d-размер расчетного поля, а h – размер клетки). Если пронумеровать потомков от нуля до <math>(\frac{d}{h})^2\:-\:1 </math>, то каждая точка будет находится в единственной клетке с номером <math> y \: div \: h \cdot \frac{d}{h} + x \: div \: h </math>. Pic. о формировании дерева Сразу видно, что при переходе к древовидным структурам становятся очевидны ограничения, наложенные предложенной сеткой – d –должно быть кратно h – самое существенное из них. Можно предложить усовершенствование при некратности, даже можно придумать отображение с неквадратных полей, но все это уже сильно ставит под вопрос физику моделируемых процессов (второе например вводит анизотропию, которая в классической задаче N-тел не предусмотрена). <br/>
Алгоритм Барнса-Хата использует так называемое квадродерево, каждый узел в котором имеет ровно 4 потомка. Терминальные вершины (вершины, не имеющие потомков) в оригинальном методе Барнса-Хата имеют 1 точку. Более детальное описание алгоритма можно посмотреть здесь [http://arborjs.org/docs/barnes-hut]. <br/>
Разумеется можно использовать и другие деревья [https://habr.com/ru/post/341208/], в частности, крайне важная статья [http://www.cs.cmu.edu/~agray/nips-final.pdf]. <br/>
Мы пока будем говорить об усовершенствовании построения именно квадродеревьев. Напомню, что алгоритм Барнс-Хата работает тем хуже, чем менее равномерно распределение частиц. Таким образом, сконцентрируемся на таком построении деревьев, при котором эффективно считаются большие скопления частиц, но при этом при равномерном распределении деревья вырождались в классический метод Барнса-Хата. <br/>
Для этого добавим новые правила в слабое место Барнса-Хата - спуск по иерархии при невыполнении MAC. Для этого сначала введем два типа скоплений частиц: кластеры и твердые тела. 
Будем называть кластерами достаточно такое скопление частиц, чтобы создаваемое им воздействие можно было заменить потенциалом. Например, если удастся доказать, что частицы в некоторой вершине сформировали круг, то нет смысла считать взаимодействие окружающих частиц с каждой из частиц, входящих в эту вершину – достаточно одного вычисления. Нельзя не отметить, что такой подход – прямое развитие оригинальной идеи Барнса и Хата, только с менее жестким ограничением на расстояние. Напомню, что аппроксимация силы в оригинальном методе происходила лишь на достаточно большом расстоянии, что теряло смысл при больших неоднородностях в распределении частиц. 
По определению твердые тела – тела, взаимные расстояния между частицами которого не меняются. Распределение скоростей точек движущегося абсолютно твёрдого тела описывается формулой Эйлера. Твердые тела – по сути своей те же кластеры, только взаимодействия внутри твердого тела можно не считать. <br/>
Взаимодействие близко расположенных твердых тел так же можно считать используя модель V-model[http://tm.spbstu.ru/V-model] <br/>
Очевидно, что точное решение задачи определения скопления как кластера и тем более, как твердого тела асимптотически сложнее, чем прямой расчет всех взаимодействий, поэтому сейчас разрабатывается несколько теорий: 
#В рамках исследования образований планет можно считать, что кластеры (и твердые тела соответственно) имеют круглую форму, так из-за проблемы «метрового барьера» [https://www.researchgate.net/publication/264276355_Overcoming_the_Meter_Barrier_and_The_Formation_of_Systems_with_Tightly-packed_Inner_Planets_STIPs] планетообразование происходит в результате коллапса, то есть имеют стахостическую природу.
#Возможно достаточно точно определить наличие кластера или твердого тела в вершине на основе различных статистик.
Первая теория лишь определяет, какой именно потенциал нужно будет брать для аппроксимации взаимодействий. Возможно, ее придется совмещать со второй теорией. В поддержку второй теории выступает большое распространение методов, построенных на случайных лесах для решения задач о нахождении k ближайших соседей. <br/>
Выгода применения статистики можно продемонстрировать на примере - за O(k) можно найти среднюю плотность k частиц, но так же можно построить диаграмму распределений k <math>  \cdot </math> const случайных расстояний, что включит в себя и среднюю плотность, и информацию об отсутствии/наличии неоднородостей. <br/> 
==Алгорим подсчета сил==
Предположим, что древовидная структура уже построена. Тогда запустив следующая рекурсивная функция позволит вычислить силу, действующую на точку P со стороны вершины Node. (в рамках задачи N-тел данная функция будет запущена для всех частиц, кроме принадлежащих твердым телам и Node будет корнем):
#function CalGrav(P, Node)
#gravity = 0
#для каждой дочерней вершины N в Node
#Проверить, является ли N кластером
##Если да, то проверить, принадлежит ли P к кластеру.
###Если да, то проверить принадлежит ли P к N
####Если да, то провести внутрикластерную аппроксимацию
####Если нет, то провести аппроксимацию "кластер-кластер"
###Если нет, то провести аппроксимацию "частица-кластер"
##Если нет, то проверить, можно ли N считать материальной точкой для P
###Если да, то провести апроксимацию Барнса-Хата
###Если нет, то gravity = gravity + CalcGrav(P, N)
#Вернуть gravity