Оптимизация расчета гравитации для частиц в невесомости 2d пространстве

Question

Я создал небольшую визуализацию частиц в python.Я рассчитываю движение частиц в двухмерном пространстве с невесомостью.Поскольку каждая частица притягивает все остальные частицы в зависимости от их массы и расстояния.

Я сделал визуализацию в Pygame, и все работает по плану (с расчетом), однако мне нужно очень быстро оптимизировать вычисления.Сегодня система может рассчитывать около 100-150 частиц в соседней частоте кадров.Я поместил все вычисления в отдельную ветку, которая дала мне немного больше, но не почти то, что я хочу.

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

Мне нужно вычислить все притяжение всех частиц, которые у меня есть, к петле в петле.И так как мне нужно выяснить, столкнулись ли они, я должен сделать то же самое снова и снова.

У меня разбито сердце, чтобы писать такой код ....

У Нампи естьВозможность рассчитать массив с массивом, однако я не нашел, что рассчитать все элементы в массиве со всеми элементами из того же / другого массива.Есть один?Если это так, я мог бы создать и пару массивов и вычислить намного быстрее, и должна быть функция для получения индекса из двух массивов, где их значения совпадают (Collitiondetect iow)

Вот сегодняшнее вычисление притяжения / коллизии:

class Particle:
    def __init__(self):
        self.x = random.randint(10,790)
        self.y = random.randint(10,590)
        self.speedx = 0.0
        self.speedy = 0.0
        self.mass = 4

#Attraction    
for p in Particles:
    for p2 in Particles:
        if p != p2:
            xdiff = P.x - P2.x
            ydiff = P.y - P2.y
            dist = math.sqrt((xdiff**2)+(ydiff**2))
            force = 0.125*(p.mass*p2.mass)/(dist**2)
            acceleration = force / p.mass
            xc = xdiff/dist
            yc = ydiff/dist
            P.speedx -= acceleration * xc
            P.speedy -= acceleration * yc
for p in Particles:
    p.x += p.speedx
    p.y += p.speedy

#Collision
for P in Particles:
   for P2 in Particles:
        if p != P2:
            Distance = math.sqrt(  ((p.x-P2.x)**2)  +  ((p.y-P2.y)**2)  )
            if Distance < (p.radius+P2.radius):
                p.speedx = ((p.mass*p.speedx)+(P2.mass*P2.speedx))/(p.mass+P2.mass)
                p.speedy = ((p.mass*p.speedy)+(P2.mass*P2.speedy))/(p.mass+P2.mass)
                p.x = ((p.mass*p.x)+(P2.mass*P2.x))/(p.mass+P2.mass)
                p.y = ((p.mass*p.y)+(P2.mass*P2.y))/(p.mass+P2.mass)
                p.mass += P2.mass
                p.radius = math.sqrt(p.mass)
                Particles.remove(P2)

Answer 1

Сначала вы можете попробовать поработать с комплексными числами : соответствующие формулы гравитации и динамики очень просты в этом формализме, а также могут быть довольно быстрыми (поскольку NumPy может выполнять вычисления внутри вас, а не вы обрабатывать отдельно координаты x и y). Например, сила между двумя частицами в z и z 'просто:

(z-z')/abs(z-z')**3

NumPy может очень быстро рассчитать такое количество для всех пар z / z '. Например, матрица всех значений zz 'просто получается из массива 1D Z координат как Z-Z[:, numpy.newaxis] (диагональные члены [z = z'] действительно требуют особой осторожности при вычислении 1/abs(z-z')**3: они должны быть установлен на ноль).

Что касается эволюции во времени, вы, безусловно, можете использовать Подпрограммы SciPy для быстрых дифференциальных уравнений : они намного более точны, чем пошаговая интеграция Эйлера.

В любом случае углубление в NumPy было бы очень полезно, особенно если вы планируете проводить научные вычисления, поскольку NumPy очень быстр.

Answer 2

Я работал над этим ранее, и одна из вещей, которые я видел в прошлом для ускорения расчетов столкновений, - это сохранение списка соседних частиц.

По сути, идея заключается в том, что внутри вашего расчета гравитации вы делаете что-то вроде:

for (int i = 0; i < n; i++)
{
    for (int j = i + 1; j < n; j++)
    {
        DoGravity(Particle[i], Particle[j]);
        if (IsClose(Particle[i], Particle[j]))
        {
            Particle[i].AddNeighbor(Particle[j]);
            Particle[j].AddNeighbor(Particle[i]);
        }
    }
}

Затем вы просто пропускаете все частицы и обнаруживает столкновения для каждого включения по очереди.Обычно это что-то вроде O(n) в лучшем случае, но в худшем случае оно может легко ухудшиться до O(n^2).

Другая альтернатива - попытаться поместить свои частицы внутрь Octree ,Создание одного это что-то вроде O(n), тогда вы можете запросить его, чтобы увидеть, есть ли что-нибудь рядом друг с другом.В этот момент вы просто обнаружите столкновения на парах.Делать это - O(n log n) Я считаю.

Не только это, но вы также можете использовать Octree для ускорения расчета гравитации.Вместо поведения O(n^2) оно также падает до O(n log n).В большинстве реализаций Octree есть «параметр открытия», который управляет компромиссом скорости и точности, который вы будете совершать.Таким образом, Octrees, как правило, менее точны, чем прямые парные вычисления, и сложны для кодирования, но они также делают возможным крупномасштабное моделирование.

Если вы используете Octree таким образом, вы будете делать то, что известно как Barnes-Hut Simulation .

Примечание. Поскольку вы работаете в 2D, 2D-аналог Octree известен как Quadtree.См. Следующую статью в Википедии для получения дополнительной информации: http://en.wikipedia.org/wiki/Quadtree

Answer 3

, чтобы сделать быстрый расчет, вам нужно хранить x, y, speedx, speedy, m в массивах numpy.Например:

import numpy as np

p = np.array([
    (0,0),
    (1,0),
    (0,1),
    (1,1),
    (2,2),
], dtype = np.float)

p - массив 5x2, в котором хранится положение частиц x, y.Чтобы рассчитать расстояние между каждой парой, вы можете использовать:

print np.sqrt(np.sum((p[:, np.newaxis] - p[np.newaxis, :])**2, axis=-1))

на выходе:

[[ 0.          1.          1.          1.41421356  2.82842712]
 [ 1.          0.          1.41421356  1.          2.23606798]
 [ 1.          1.41421356  0.          1.          2.23606798]
 [ 1.41421356  1.          1.          0.          1.41421356]
 [ 2.82842712  2.23606798  2.23606798  1.41421356  0.        ]]

или вы можете использовать cdist от scipy:

from scipy.spatial.distance import cdist
print cdist(p, p)

Answer 4

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

C ++ модуль для вычисления x и yкомпоненты гравитационного притяжения на двухмерной плоскости:

#include <Python.h>
#include <math.h>

double _acceleration(double &Vxa, double &Vya, double &Vxb, double &Vyb, double xa, double ya, double xb, double yb, double massa, double massb)
{
   double xdiff = xa - xb;
   double ydiff = ya - yb;
   double distance = sqrt(xdiff*xdiff + ydiff*ydiff) * pow(10, 5);

   if (distance <= 0)
      distance = 1;

   double force = (6.674 * pow(10, -11))*(massa*massb)/(distance*distance);

   double acca = force / massa;
   double accb = force / massb;
   double xcomponent = xdiff/distance;
   double ycomponent = ydiff/distance;

   Vxa -= acca * xcomponent;
   Vya -= acca * ycomponent;
   Vxb += accb * xcomponent;
   Vyb += accb * ycomponent;

   return distance;
}

static PyObject* gforces(PyObject* self, PyObject* args)
{
   double Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb, distance;

   if (!PyArg_ParseTuple(args, "dddddddddd", &Vxa, &Vya, &Vxb, &Vyb, &xa, &ya, &xb, &yb, &massa, &massb))
      return NULL;

   distance = _acceleration(Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb);

   return Py_BuildValue("ddddd", Vxa, Vya, Vxb, Vyb, distance);
}

static PyMethodDef GForcesMethods[] = {
{"gforces", gforces, METH_VARARGS, "Calculate the x and y acceleration of two masses and the distance between the two."},
{NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC
initgforces(void)
{
(void) Py_InitModule("gforces", GForcesMethods);
}

Если вы скомпилируете это как файл pyd, вы должны получить объект python, который вы можете импортировать.Вы должны правильно настроить все параметры компилятора и компоновщика.Я использую dev-C ++ и мои параметры компилятора установлены на -shared -o gforces.pyd, а компоновщик установлен на -lpython27 (убедитесь, что вы используете ту же версию, что вы установили), и добавьте путь к каталогу python для включаемых файлов и библиотек.tabs.

Объект принимает аргументы (p1.speedx, p1.speedy, p2.speedx, p2.speedy, p1.x, p1.y, p2.x, p2.y, p1.mass,p2.mass) и возвращает новые p1.speedx, p1.speedy, p2.speedx, p2.speedy и расстояние между p1 p2.

Используя вышеупомянутый модуль, я также попытался вырезатьнесколько шагов для обнаружения столкновения путем сравнения возвращенного расстояния с суммой радиусов частиц как таковой:

def updateForces(self):         #part of a handler class for the particles
    prevDone = []
    for i in self.planetGroup:  #planetGroup is a sprite group from pygame
        prevDone.append(i.ID)
        for j in self.planetGroup:
            if not (j.ID in prevDone):               #my particles have unique IDs
                distance = i.calcGForce( j )         #calcGForce just calls the above  
                if distance <= i.radius + j.radius:  #object and assigns the returned 
                    #collision handling goes here    #values for the x and y speed
                                                     #components to the particles

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

Answer 5

(Это может быть в комментарии, но у меня нет необходимой репутации, чтобы сделать это)

Я не понимаю, как вы делаете шаг времени.У вас есть

P.speedx -= acceleration * xc
P.speedy -= acceleration * yc

, но чтобы получить новую скорость в момент времени t + delta_t, вы должны сделать

P.speedx -= acceleration * xc * delta_t
P.speedy -= acceleration * yc * delta_t

, а затем обновить позицию следующим образом:

P.x = P.x + P.speedx * delta_t
P.y = P.y + P.speedy * delta_t

Тогда к твоей скорости относятся.Может быть, лучше хранить информацию о частицах не в классе, а в массивах?Но я не думаю, что вы можете избежать циклов.

Кроме того, вы смотрели на wikipedia , там описаны некоторые методы для ускорения расчетов.

(отредактированоблагодаря комментарию Майка)