在我的粒子模拟中到底发生了什么？-解网

问：

我的碰撞球模拟工作正常，但粒子发疯了。

他们中的一些人相互粘在一起并旋转。其他人则被卡在边缘。

我想不出为什么......谁能帮我？谢谢

此代码使用 Taichi Gui。我打算将来使用 taichi 装饰器，但现在我只是用普通的 python 测试它

import taichi as ti
import numpy as np

ti.init(arch=ti.vulkan)

gui = ti.GUI('Circles Test', res=(600, 600))

numpoints = 20
centres = np.random.random((numpoints, 2))
velocities = np.random.uniform(-1,1, (numpoints,2))
radius = 20
guiradius = 0.04
damping = 0.9
timestep = 1/60

'''
@ti.kernel
def randompoints():
    for i, j in centres:
        centres[i,j] =ti.random()
'''
#randompoints()
print(centres[1])

def wallcollision(centres, velocities):
    for i in range(len(centres)):
        if centres[i][0] - guiradius <=0:
            velocities[i][0] = -damping*velocities[i][0]
        if centres[i][0] + guiradius >=1:
            velocities[i][0] = -damping*velocities[i][0]
        
        if centres[i][1] - guiradius <=0:
            velocities[i][1] = -damping*velocities[i][1]
        if centres[i][1] + guiradius >=1:
            velocities[i][1] = -damping*velocities[i][1]
    return velocities

def selfcollision(centres, velocities):
    for i in range(len(centres)):
        for j in range(i+1,len(centres)):
            if np.linalg.norm(centres[i]-centres[j]) < 2*guiradius:
                print("Collision!")
                vitemp = np.copy(velocities[i])
                vitemp = velocities[i] - (np.dot(velocities[i]-velocities[j], centres[i]-centres[j])*(centres[i]-centres[j])/(np.linalg.norm(centres[i]-centres[j])**2))
                velocities[j] = velocities[j] - (np.dot(velocities[j]-velocities[i], centres[j]-centres[i])*(centres[j]-centres[i])/(np.linalg.norm(centres[j]-centres[i])**2))
                velocities[i] = vitemp
                
                ##velocities[i] = -velocities[i]
    return velocities

def update(centres, velocities):
    centres += velocities*timestep
    velocities = wallcollision(centres, velocities)
    velocities = selfcollision(centres, velocities)

while gui.running:
    gui.circles(centres, radius, color=0xEEEEF0)
    gui.show()
    update(centres, velocities)

Python Numpy 物理

欢迎使用 Stack Overflow。请阅读如何提问。我们在这里不写“找到错误”的答案;我们需要一个具体的问题 - 这将是你理解和定位特定问题的最佳尝试，并在一个最小的可重复的例子中展示它。适合 Stack Overflow 的问题是，你已经弄清楚了代码中执行与预期不同的操作的特定部分（您应该具体期望某些内容），并且不明白为什么。

1赞 Jérôme Richard 11/10/2023

@KarlKnechtel 虽然我基本同意第一部分，但 OP 实际上确实提供了一个最小的可重复示例！它是完全可重复的，所以我很容易做到。OP 很难使它变小，因为（相同的）错误同时位于和中。您如何使其更具体，同时仍具有可重复性？其他新人提出的大多数问题通常甚至没有提供像这样的 MRE（当他们这样做时，他们不会像这里那样得到那么多的反对票）。wallcollisionselfcollision

答：

2赞 Jérôme Richard 11/10/2023 #1

问题在于，当检测到碰撞时，速度会被修改，但中心不会被修改，并且速度并不总是大到足以让粒子在下一帧中逃离碰撞区域。这意味着粒子可以在许多帧中保持碰撞状态。这导致速度分量在负值和正值之间振荡，使碰撞逃逸有时无法进行。这就是为什么粒子以某种方式粘合在一起的原因。同样的问题也发生在壁上：当粒子的中心超出壁时，速度分量会振荡，因此粒子被粘在壁上。这种问题在基础物理模拟中经常出现。一个简单的解决方法是简单地校正中心，使粒子不会在下一帧中碰撞。

以下是更正后的代码：

# My parameters:
# guiradius = 0.04
# radius = 600*guiradius # More accurate visually

def wallcollision(centres, velocities):
    for i in range(len(centres)):
        if centres[i][0] - guiradius <=0:
            velocities[i][0] = -damping*velocities[i][0]
            centres[i][0] = guiradius
        if centres[i][0] + guiradius >=1:
            velocities[i][0] = -damping*velocities[i][0]
            centres[i][0] = 1 - guiradius
        if centres[i][1] - guiradius <=0:
            velocities[i][1] = -damping*velocities[i][1]
            centres[i][1] = guiradius
        if centres[i][1] + guiradius >=1:
            velocities[i][1] = -damping*velocities[i][1]
            centres[i][1] = 1 - guiradius
    return velocities

def selfcollision(centres, velocities):
    for i in range(len(centres)):
        for j in range(i+1,len(centres)):
            distance = np.linalg.norm(centres[i]-centres[j])
            if distance < guiradius * 2:
                print("Collision!")
                vitemp = np.copy(velocities[i])
                vitemp = velocities[i] - (np.dot(velocities[i]-velocities[j], centres[i]-centres[j])*(centres[i]-centres[j])/distance**2)
                velocities[j] = velocities[j] - (np.dot(velocities[j]-velocities[i], centres[j]-centres[i])*(centres[j]-centres[i])/distance**2)
                velocities[i] = vitemp

                # Fix the centers
                overlapping_distance = guiradius * 2 - distance
                direction = (centres[j] - centres[i]) / distance
                epsilon = 1e-6 # Additional safe gap to avoid particles to collide again
                shift = overlapping_distance * direction * (0.5 + epsilon)
                centres[i] -= shift
                centres[j] += shift
    return velocities

请注意，在更真实的模拟中，力取决于碰撞表面（如果要模拟负责球排斥的电子，请使用非平凡的公式），因此碰撞不应长时间发生。

@ollied 例如，有很多很好的问题，人们在代码中看到数字不稳定问题，这些问题可以通过与数学相关的技巧来解决。这样的问题通常不适合 mathoverflow，因为问题是由于数字的离散化造成的。事实上，这是一个计算科学研究的领域。但不幸的是，许多程序员经常认为这是纯粹的数学问题。这不太适合 cs.stackexchange 和 mathoverflow，虽然 scicomp.stackexchange 可能，但 AFAIK 它仍然很新/实验性。

0赞 Jérôme Richard 11/11/2023

@ollied 可以在这里和那里看到这个问题的更具体的例子（后者可能会让您感兴趣）。这些问题显然不仅仅是编程问题，而且已经得到了社区的好评。

上一个：欧拉积分绘制水星轨道，仅绘制圆形轨道和不正确的周期

下一个：带有悬链线的curve_fit会无缘无故地返回巨大的帕斯

在我的粒子模拟中到底发生了什么？

What on earth is going in in my particle simulation?

评论

评论