优化零重力2d空间中粒子的重力计算

ztr*_*pez 5 python optimization physics numpy

我已经在python中创建了一个小的粒子可视化.我在没有重力的2D空间中对粒子的运动进行了解释.因为每个粒子基于粒子质量和距离吸引所有其他粒子.

我在pygame中做了一个视觉效果,一切都按计划(用caluclation)工作,但是我需要极大地优化计算.今天,该系统可以以相邻的帧速率计算大约100-150个粒子.我把所有的计算都放在一个单独的线程中,它给了我更多但不是我想要的东西.

我看着scipy和numpy,但因为我不是科学家或数学家,所以我只是感到困惑.看起来我在正确的轨道上,但我不知道怎么做.

我需要在循环中计算所有粒子上的所有吸引力.因为我需要找到是否有任何碰撞,我必须重新做同样的事情.

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

Numpy能够使用数组计算数组,但是我还没有找到任何计算数组中所有项目的内容,其中包含来自相同/另一个数组的所有项目.有吗?如果是这样,我可以创建几个数组并计算得更快,并且必须有一个函数将索引从2个数组中获取,其值匹配(Collitiondetect iow)

这是今天的吸引力/ collsion计算:

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)
Run Code Online (Sandbox Code Playgroud)

Eri*_*got 5

您可以首先尝试使用复数:相关的引力和动力学公式在这种形式中非常简单,并且也可以非常快(因为NumPy可以在内部进行计算,而不是单独处理x和y坐标).例如,z和z'处的两个粒子之间的力是简单的:

(z-z')/abs(z-z')**3
Run Code Online (Sandbox Code Playgroud)

对于所有z/z'对,NumPy可以非常快速地计算这样的量.例如,所有zz'值的矩阵简单地从1D Z坐标数组获得,因为Z-Z[:, numpy.newaxis](对角项[z = z']在计算时需要特别小心1/abs(z-z')**3:它们应该设置为零).

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

在任何情况下,钻研NumPy都非常有用,特别是如果你计划进行科学计算,因为NumPy非常快.


Mau*_*uro 1

(这可能应该放在评论中,但我没有这样做所需的声誉)

我不明白你是如何进行时间步进的。你有

P.speedx -= acceleration * xc
P.speedy -= acceleration * yc
Run Code Online (Sandbox Code Playgroud)

但要在时间 t+delta_t 获得新的速度,你会这样做

P.speedx -= acceleration * xc * delta_t
P.speedy -= acceleration * yc * delta_t
Run Code Online (Sandbox Code Playgroud)

然后像这样更新位置:

P.x = P.x + P.speedx * delta_t
P.y = P.y + P.speedy * delta_t
Run Code Online (Sandbox Code Playgroud)

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

另外,您是否看过维基百科,其中描述了一些加快计算速度的方法。

(根据迈克的评论进行编辑)