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

Question

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

我在pygame中制作了一个可视化程序，一切都按计划进行（通过计算），但是我需要从外部优化计算。今天，这个系统可以计算出大约100-150个粒子的帧率。我把所有的计算放在一个单独的线程，给了我更多，但不是我想要的。在

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

我需要计算所有粒子对一个循环的所有吸引力。 既然我需要找出是否有碰撞，我就得再做一遍。在

写这种代码让我心碎。。。。在

Numpy可以用array计算数组，但是我还没有找到用来自同一个/另一个数组的所有项来计算数组中所有项的内容。有吗？ 如果是这样的话，我可以创建和耦合数组并计算得更快，并且必须有一个函数将索引从2个数组中获取，它们的值匹配（Collitiondetect low）

今天的吸引力/吸引力计算如下：

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

我以前也做过这方面的工作，我在过去看到的加速碰撞计算的方法之一就是实际存储一个附近粒子的列表。在

基本上，在你的重力计算中，你可以做一些类似的事情：

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)。大多数八叉树实现都包含一个“开放参数”，它控制您将要进行的速度与精度的权衡

Answer 2

你可以先尝试使用复数：在这种形式下，相关的引力和动力学公式非常简单，而且速度也相当快（因为NumPy可以在内部进行计算，而不是单独处理x和y坐标）。例如，z和z'处两个粒子之间的力很简单：

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

对于所有的z/z'对，NumPy可以很快地计算出这样一个量。例如，所有z-z'值的矩阵只需从坐标的1D数组Z获得，即Z-Z[:, numpy.newaxis]（对角线项[z=z']在计算1/abs(z-z')**3时确实需要特别注意：它们应该设置为零）。在

至于时间演化，你当然可以使用SciPy的快速微分方程例程：它们比逐步欧拉积分精确得多。在

在任何情况下，钻研NumPy将是非常有用的，特别是如果你打算做科学计算，因为NumPy是非常快的。在