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

Question

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

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

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

我需要计算所有粒子的所有吸引力，我必须在一个循环中循环。由于我需要找出是否有任何碰撞，我必须重新做同样的事情。

写那种代码让我心碎....

Numpy 有能力用数组计算数组，但是我还没有找到任何东西来计算数组中的所有项目以及来自相同/另一个数组的所有项目。有吗？如果是这样的话，我可以创建和耦合数组并计算得更快，并且必须有一个函数可以从 2 个数组中获取索引，其中它们的值匹配(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 可以非常快速地计算出这样的数量。例如，所有 z-z' 值的矩阵只是从坐标的一维数组 Z 中作为 Z-Z[:, numpy.newaxis] 获得(对角线项 [z= z'] 在计算 1/abs(z-z')**3 时确实需要特别注意:它们应该设置为零)。

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

无论如何，深入研究 NumPy 都会非常有用，特别是如果您打算进行科学计算，因为 NumPy 非常快。