为什么我在重力模型中的粒子会相互排斥？答案

【问题标题】：Why do my particles in a gravity model repel with each other?为什么我在重力模型中的粒子会相互排斥？
【发布时间】：2021-04-13 10:00:15
【问题描述】：

所以我正在尝试建立一个简单的重力粒子相互作用模型（稍后我想让它变得更复杂，但现在不是）。

我正在处理 Processing，特别是 processing.py。我创建了 2 个文件，一个是主要的 setup 和 draw 函数。另一个是粒子类，我在其中定义粒子属性和重力定律。对于重力，我使用牛顿 F = Gm1m2 / r^2 。

这是主文件：

#MAIN --> simulation of a particles system influenced by gravity and temperatures
from particle import Particle

num_particles = 20
particle_array = []

show = False #To show the distance with particles that are interacting (LATER)

def setup():
    size(1000, 600)
    for p in range(num_particles):
        p = Particle()
        particle_array.append(p)
    
def draw():
    background(255)
    
    for i in range(len(particle_array)):
        particle_array[i].display()
        #I'm creating another array without the particle considered in the loop. So i can confront the first one with
        #every other particle except itself.
        part_array = [elem for elem in particle_array if elem != particle_array[i]]
        #check interaction with every particles
        for k in range(len(part_array)):
            particle_array[i].check_interaction(part_array[k], show)
    
    #different loops to first determine the acceleration without changing positions
    #and later on updating the positions
    for elem in particle_array: 
        elem.update()
        
    #fill(173, 164, 164)
    #rect(10, 10, 20, 20) #Show Interactions box

这是课程： https://pastebin.com/nJJE45Rs

那么发生了什么？当我开始模拟时，坐标发生了一些奇怪的事情。

我收到“ValueError：无法将浮点 NaN 转换为整数”。我试图从它的未来看，似乎一些坐标随便变成了 NaN 并且 int() 函数崩溃了，因为它不能进行从 NaN 到 int 的转换。所以很奇怪。（已解决，见评论）

此外，粒子的行为似乎非常奇怪，它看起来不像是重力。它们似乎是排斥而不是吸引。

对于这个模型的数学，我已经开始得到力的加速度模块。然后，假设加速度矢量指向另一个粒子（正在计算相互作用的粒子），我必须确定我将用于计算 Ax 和 Ay 的 x 轴角度，如下所示：

angle = acos((x2 - x1) / r) #Orientation of the acceleration vector
            print("Angle: " + str(angle) + "x1, y1, x2, y2: " + str(x1) +" "+ str(y1) +" "+ str(x2) +" "+ str(y2), "R: " + str(r))
            #print("New Angle: " + str(angle), "X1; Y1, X2, Y2: " + str(x1) +str(y1) + str(x2) + str(y2))
            ax = a * cos(angle)
            ay = a * sin(angle)
            self.a.add(ax, ay, 0) #update Acceleration vector

我使用过这样的 acos 函数，因为在一个矩形中：

A = Point (x1, y1)
B = Point (x2, y2)
AB * cos(angle) = (x2 - x1)

但我不知道，事情似乎失败了。

【问题讨论】：

我不是 python 爱好者，但我要做的第一件事是打印 x1、x2、y1、y2 并查看它们包含的内容。那里可能会有一些要理解的东西，即使只是他们很好。如果你这样做，请标记我并将信息添加到你的帖子中，这样我就可以知道你更新了它。
@laancelot 我已经这样做了。通常代码运行顺利，但有时，四个值中的 2 个，尤其是 x2 和 y2 突然变成了 NaN。即使我强制它们为 0 而不是 NaN。
好的，请检查您的代码中有关在崩溃方法中到达或影响到达 NaN 值的数字的以下可能性：dividing by zero、sqrt(num) where num < 0、number is not a float。您可以通过不同方式检查 NaN，除其他外（在伪代码中）：if num != num then return 可以让您跳过异常，如果没有其他方法。
好的，我找到了错误。这是一个数学错误。我将一个大于 1 的值传递给 asin 函数，这是不可能的。因此该函数返回一个 NaN 值，该值会导致下一个循环出现问题。这是由应用于 r 分配的 int() 引起的。 R 是近似的，这导致值大于一。但我认为物理学非常糟糕。粒子似乎相互排斥而不是相互吸引。啊！
恭喜您查明问题！提醒我，在计算 GPS 坐标与光照覆盖范围时，我遇到了类似不可能三角问题的应用程序崩溃。如果您想增加获得帮助的机会，我建议您针对当前问题发表一篇新文章。

标签： python processing simulation particles physics-engine

【解决方案1】：

我建议完全避免使用角度。我发现向量更容易使用。如果以后想知道角度，可以使用PVector.heading()方法输出向量与(0, 1)向量之间的角度。

如果我从头开始制作，我也会将位置转换为向量。由于这会影响您的其他方法，因此我将如何调整您的代码：

（免责声明：我不在 Python 中使用处理，所以我无法运行它来检查错误。如果遇到任何问题，请发表评论）

def check_interaction(self, p2, show_interactions):
    r_vec = PVector(p2.x - self.x, p2.y - self.y) #vector from self to other particle
    if r.mag() < 200:
        if show_interactions: #draw the line
            stroke(110)
            line(self.x, self.y, p2.x, p2.y)
        if 1 < r.magSq():
            grav_acc = G * p2.mass / r_vec.magSq() #acceleration due to the other particle
            r_vec.setMag(grav_acc) #Change magnitude of r_vec so it's now the acceleration vector (due to p2)
            self.a = self.a.add(r_vec) #update acceleration

如果你想转换

self.x = random(screen_width)
self.y = random(screen_height)

到

self.pos = PVector(random(screen_width), random(screen_height))

在def __init__(self)，你可以这样定义r_vec：

r_vec = p2.pos.sub(self.pos)

【讨论】：