优化零重力二维空间中粒子的引力计算

Question

我用 python 创建了一个小的粒子可视化。 我正在计算零重力二维空间中粒子的运动。 每个粒子都会根据粒子质量和距离吸引所有其他粒子。

我在 pygame 中做了一个可视化，一切都按计划进行（通过计算），但是我需要极大地优化计算。如今，该系统可以在相当好的帧速率下计算大约 100-150 个粒子。我将所有计算放在一个单独的线程中，这给了我更多的计算，但远不是我想要的。

我看过 scipy 和 numpy 但由于我不是科学家或数学大师，我只是感到困惑。看起来我走在正确的轨道上，但我不知道如何做。

我需要计算所有粒子对循环的所有吸引力。 因为我需要找出是否有碰撞，所以我必须再次做同样的事情。

写出这种代码让我心碎……

Numpy 能够用数组计算数组，但是我还没有找到任何方法来计算数组中的所有项目和同一个/另一个数组中的所有项目。有吗？ 如果是这样，我可以创建几个数组并计算得更快，并且必须有一个函数可以从两个数组中获取索引，其中它们的值匹配（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

对于所有 z/z' 对，NumPy 可以非常快速地计算出这样的量。例如，所有 zz' 值的矩阵只需从坐标为 ZZ[:, numpy.newaxis] 的一维数组 Z 获得（对角项 [z= z'] 在计算 1/abs(zz')**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)。

不仅如此，您还可以使用八叉树来加速重力计算。它不再是 O(n^2) 行为，而是下降到 O(n log n) 。大多数八叉树实现都包含一个“开放参数”，用于控制您将进行的速度与准确性权衡。因此，八叉树往往不如直接的成对计算准确，并且编码复杂，但它们也使大规模模拟成为可能。

如果您以这种方式使用八叉树，您将执行所谓的 Barnes-小屋模拟。

注意：由于您在 2D 环境中工作，因此 八叉树 的 2D 模拟称为四叉树。有关详细信息，请参阅以下维基百科文章：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.        ]]

或者您可以使用 scipy 中的 cdist：

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

Answer 4

不确定这是否会对您有很大帮助，但这是我一直致力于解决同一问题的解决方案的一部分。我没有注意到这样做在性能上有巨大的提升，仍然开始在 200 个粒子左右停止，但也许它会给你一些想法。

用于计算 2d 平面上引力的 x 和 y 分量的 C++ 模块：

#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 目录路径添加到包含文件和库中选项卡。

该对象接受参数（ 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

然后是您的速度问题。也许将粒子信息存储在 numpy 数组中而不是类中会更好？但我认为你无法避免循环。

另外，您是否看过wikipedia，那里描述了一些加快计算速度的方法。

（根据迈克的评论进行编辑）