盒中的许多粒子 - 物理模拟
Many particles in box - physics simulation
我目前正在尝试模拟盒子中弹跳的许多粒子。
我已经考虑了@kalhartt 的建议，这是用于初始化盒子内粒子的改进代码：
import numpy as np
import scipy.spatial.distance as d
import matplotlib.pyplot as plt
# 2D container parameters
# Actual container is 50x50 but chose 49x49 to account for particle radius.
limit_x = 20
limit_y = 20
#Number and radius of particles
number_of_particles = 350
radius = 1
def force_init(n):
# equivalent to np.array(list(range(number_of_particles)))
count = np.linspace(0, number_of_particles-1, number_of_particles)
x = (count + 2) % (limit_x-1) + radius
y = (count + 2) / (limit_x-1) + radius
return np.column_stack((x, y))
position = force_init(number_of_particles)
velocity = np.random.randn(number_of_particles, 2)
初始化的位置如下所示：
初始化粒子后，我想在每个时间步更新它们。更新的代码紧跟之前的代码，如下：
# Updating
while np.amax(abs(velocity)) > 0.01:
# Assume that velocity slowly dying out
position += velocity
velocity *= 0.995
#Get pair-wise distance matrix
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
#If pdist [i,j] is <=4 then the particles are too close and so treat as collision
for i in range(len(pair_d)):
for j in range(i):
# Only looking at upper triangular matrix (not inc. diagonal)
if pair_d[i,j] ==True:
# If two particles are too close then swap velocities
# It's a bad hack but it'll work for now.
vel_1 = velocity[j][:]
velocity[j] = velocity[i][:]*0.9
velocity[i] = vel_1*0.9
# Masks for particles beyond the boundary
xmax = position[:, 0] > limit_x
xmin = position[:, 0] < 0
ymax = position[:, 1] > limit_y
ymin = position[:, 1] < 0
# flip velocity and assume that it looses 10% of energy
velocity[xmax | xmin, 0] *= -0.9
velocity[ymax | ymin, 1] *= -0.9
# Force maximum positions of being +/- 2*radius from edge
position[xmax, 0] = limit_x-2*radius
position[xmin, 0] = 2*radius
position[ymax, 0] = limit_y-2*radius
position[ymin, 0] = 2*radius
更新并让它 运行 完成后，我得到了这个结果：
这比以前好多了，但仍然有一些补丁靠得太近了——例如：
靠得太近了。我认为更新有效......感谢@kalhartt，我的代码变得更好更快（我学到了一些关于 numpy 的东西......道具@kalhartt）但我仍然不知道它在哪里搞砸了。我试过更改实际更新的顺序，将成对距离放在最后或 position +=velocity 最后，但无济于事。我添加了 *0.9 以使整个事物更快地消失，我尝试使用 4 来确保 2*radius (=2) 不是太严格的标准......但似乎没有任何效果。
如有任何帮助，我们将不胜感激。
只有两个错别字妨碍了您。首先 for i in range(len(positions)/2): 只迭代一半以上的粒子。这就是为什么一半的粒子停留在 x 边界内（如果你观察大的迭代它会更清楚）。其次，第二个 y 条件应该是最小值（我假设）position[i][1] < 0。以下块用于为我绑定粒子（我没有使用碰撞代码进行测试，因此那里可能存在问题）。
for i in range(len(position)):
if position[i][0] > limit_x or position[i][0] < 0:
velocity[i][0] = -velocity[i][0]
if position[i][1] > limit_y or position[i][1] < 0:
velocity[i][1] = -velocity[i][1]
顺便说一句，尽可能利用 numpy 来消除循环。它更快、更高效，而且在我看来更具可读性。例如 force_init 看起来像这样：
def force_init(n):
# equivalent to np.array(list(range(number_of_particles)))
count = np.linspace(0, number_of_particles-1, number_of_particles)
x = (count * 2) % limit_x + radius
y = (count * 2) / limit_x + radius
return np.column_stack((x, y))
你的边界条件应该是这样的：
while np.amax(abs(velocity)) > 0.01:
position += velocity
velocity *= 0.995
# Masks for particles beyond the boundary
xmax = position[:, 0] > limit_x
xmin = position[:, 0] < 0
ymax = position[:, 1] > limit_y
ymin = position[:, 1] < 0
# flip velocity
velocity[xmax | xmin, 0] *= -1
velocity[ymax | ymin, 1] *= -1
最后一点，用 position[xmax, 0] = limit_x; position[xmin, 0] = 0 之类的东西将位置硬剪裁到边界框可能是个好主意。在某些情况下，速度可能很小，盒子外的粒子会被反射，但在下一次迭代中不会进入盒子内。所以它只会坐在盒子外面被永远反射。
编辑：碰撞
碰撞检测是一个更难的问题，但让我们看看我们能做些什么。让我们看看您当前的实现。
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
for i in range(len(pair_d)):
for j in range(i):
# Only looking at upper triangular matrix (not inc. diagonal)
if pair_d[i,j] ==True:
# If two particles are too close then swap velocities
# It's a bad hack but it'll work for now.
vel_1 = velocity[j][:]
velocity[j] = velocity[i][:]*0.9
velocity[i] = vel_1*0.9
总的来说这是一个很好的方法，cdist 将有效地计算距离
在点集之间，你会发现哪些点与 pair_d = pair_dist<=4.
相撞
嵌套的for循环是第一个问题。我们需要迭代 pair_d 的 True 值，其中 j > i。首先，您的代码实际上通过使用 for j in range(i) 遍历下三角区域，因此 j < i 在这种情况下并不特别重要，因为 i,j 对不重复。然而 Numpy 有两个我们可以使用的内置函数，np.triu lets us set all values below a diagonal to 0 and np.nonzero 会给我们矩阵中非零元素的索引。所以这个：
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
for i in range(len(pair_d)):
for j in range(i+1, len(pair_d)):
if pair_d[i, j]:
...
相当于
pair_dist = d.cdist(position, position)
pair_d = np.triu(pair_dist<=4, k=1) # k=1 to exclude the diagonal
for i, j in zip(*np.nonzero(pair_d)):
...
第二个问题（如您所述）是速度只是切换和缩放而不是反射。我们真正想要做的是否定和缩放每个粒子速度沿连接它们的轴的分量。请注意，为此我们需要连接它们的向量 position[j] - position[i] 和连接它们的向量的长度（我们已经计算过了）。所以不幸的是 cdist 计算的一部分被重复了。让我们停止使用 cdist 并自己动手吧。这里的目标是制作两个数组 diff 和 norm，其中 diff[i][j] 是从粒子 i 指向 j 的向量（因此 diff 是一个 3D 数组）和 norm[i][j] 是粒子 i 和 j 之间的距离。我们可以像这样用 numpy 做到这一点：
nop = number_of_particles
# Give pos a 3rd index so we can use np.repeat below
# equivalent to `pos3d = np.array([ position ])
pos3d = position.reshape(1, nop, 2)
# 3D arras with a repeated index so we can form combinations
# diff_i[i][j] = position[i] (for all j)
# diff_j[i][j] = position[j] (for all i)
diff_i = np.repeat(pos3d, nop, axis=1).reshape(nop, nop, 2)
diff_j = np.repeat(pos3d, nop, axis=0)
# diff[i][j] = vector pointing from position[i] to position[j]
diff = diff_j - diff_i
# norm[i][j] = sqrt( diff[i][j]**2 )
norm = np.linalg.norm(diff, axis=2)
# check for collisions and take the region above the diagonal
collided = np.triu(norm < radius, k=1)
for i, j in zip(*np.nonzero(collided)):
# unit vector from i to j
unit = diff[i][j] / norm[i][j]
# flip velocity
velocity[i] -= 1.9 * np.dot(unit, velocity[i]) * unit
velocity[j] -= 1.9 * np.dot(unit, velocity[j]) * unit
# push particle j to be radius units from i
# This isn't particularly effective when 3+ points are close together
position[j] += (radius - norm[i][j]) * unit
...
因为这个 post 已经足够长了，here is a gist 的代码加上我的修改。
我目前正在尝试模拟盒子中弹跳的许多粒子。
我已经考虑了@kalhartt 的建议，这是用于初始化盒子内粒子的改进代码：
import numpy as np
import scipy.spatial.distance as d
import matplotlib.pyplot as plt
# 2D container parameters
# Actual container is 50x50 but chose 49x49 to account for particle radius.
limit_x = 20
limit_y = 20
#Number and radius of particles
number_of_particles = 350
radius = 1
def force_init(n):
# equivalent to np.array(list(range(number_of_particles)))
count = np.linspace(0, number_of_particles-1, number_of_particles)
x = (count + 2) % (limit_x-1) + radius
y = (count + 2) / (limit_x-1) + radius
return np.column_stack((x, y))
position = force_init(number_of_particles)
velocity = np.random.randn(number_of_particles, 2)
初始化的位置如下所示：
初始化粒子后，我想在每个时间步更新它们。更新的代码紧跟之前的代码，如下：
# Updating
while np.amax(abs(velocity)) > 0.01:
# Assume that velocity slowly dying out
position += velocity
velocity *= 0.995
#Get pair-wise distance matrix
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
#If pdist [i,j] is <=4 then the particles are too close and so treat as collision
for i in range(len(pair_d)):
for j in range(i):
# Only looking at upper triangular matrix (not inc. diagonal)
if pair_d[i,j] ==True:
# If two particles are too close then swap velocities
# It's a bad hack but it'll work for now.
vel_1 = velocity[j][:]
velocity[j] = velocity[i][:]*0.9
velocity[i] = vel_1*0.9
# Masks for particles beyond the boundary
xmax = position[:, 0] > limit_x
xmin = position[:, 0] < 0
ymax = position[:, 1] > limit_y
ymin = position[:, 1] < 0
# flip velocity and assume that it looses 10% of energy
velocity[xmax | xmin, 0] *= -0.9
velocity[ymax | ymin, 1] *= -0.9
# Force maximum positions of being +/- 2*radius from edge
position[xmax, 0] = limit_x-2*radius
position[xmin, 0] = 2*radius
position[ymax, 0] = limit_y-2*radius
position[ymin, 0] = 2*radius
更新并让它 运行 完成后，我得到了这个结果：
这比以前好多了，但仍然有一些补丁靠得太近了——例如：
靠得太近了。我认为更新有效......感谢@kalhartt，我的代码变得更好更快（我学到了一些关于 numpy 的东西......道具@kalhartt）但我仍然不知道它在哪里搞砸了。我试过更改实际更新的顺序，将成对距离放在最后或 position +=velocity 最后，但无济于事。我添加了 *0.9 以使整个事物更快地消失，我尝试使用 4 来确保 2*radius (=2) 不是太严格的标准......但似乎没有任何效果。
如有任何帮助，我们将不胜感激。
只有两个错别字妨碍了您。首先 for i in range(len(positions)/2): 只迭代一半以上的粒子。这就是为什么一半的粒子停留在 x 边界内（如果你观察大的迭代它会更清楚）。其次，第二个 y 条件应该是最小值（我假设）position[i][1] < 0。以下块用于为我绑定粒子（我没有使用碰撞代码进行测试，因此那里可能存在问题）。
for i in range(len(position)):
if position[i][0] > limit_x or position[i][0] < 0:
velocity[i][0] = -velocity[i][0]
if position[i][1] > limit_y or position[i][1] < 0:
velocity[i][1] = -velocity[i][1]
顺便说一句，尽可能利用 numpy 来消除循环。它更快、更高效，而且在我看来更具可读性。例如 force_init 看起来像这样：
def force_init(n):
# equivalent to np.array(list(range(number_of_particles)))
count = np.linspace(0, number_of_particles-1, number_of_particles)
x = (count * 2) % limit_x + radius
y = (count * 2) / limit_x + radius
return np.column_stack((x, y))
你的边界条件应该是这样的：
while np.amax(abs(velocity)) > 0.01:
position += velocity
velocity *= 0.995
# Masks for particles beyond the boundary
xmax = position[:, 0] > limit_x
xmin = position[:, 0] < 0
ymax = position[:, 1] > limit_y
ymin = position[:, 1] < 0
# flip velocity
velocity[xmax | xmin, 0] *= -1
velocity[ymax | ymin, 1] *= -1
最后一点，用 position[xmax, 0] = limit_x; position[xmin, 0] = 0 之类的东西将位置硬剪裁到边界框可能是个好主意。在某些情况下，速度可能很小，盒子外的粒子会被反射，但在下一次迭代中不会进入盒子内。所以它只会坐在盒子外面被永远反射。
编辑：碰撞
碰撞检测是一个更难的问题，但让我们看看我们能做些什么。让我们看看您当前的实现。
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
for i in range(len(pair_d)):
for j in range(i):
# Only looking at upper triangular matrix (not inc. diagonal)
if pair_d[i,j] ==True:
# If two particles are too close then swap velocities
# It's a bad hack but it'll work for now.
vel_1 = velocity[j][:]
velocity[j] = velocity[i][:]*0.9
velocity[i] = vel_1*0.9
总的来说这是一个很好的方法，cdist 将有效地计算距离
在点集之间，你会发现哪些点与 pair_d = pair_dist<=4.
嵌套的for循环是第一个问题。我们需要迭代 pair_d 的 True 值，其中 j > i。首先，您的代码实际上通过使用 for j in range(i) 遍历下三角区域，因此 j < i 在这种情况下并不特别重要，因为 i,j 对不重复。然而 Numpy 有两个我们可以使用的内置函数，np.triu lets us set all values below a diagonal to 0 and np.nonzero 会给我们矩阵中非零元素的索引。所以这个：
pair_dist = d.cdist(position, position)
pair_d = pair_dist<=4
for i in range(len(pair_d)):
for j in range(i+1, len(pair_d)):
if pair_d[i, j]:
...
相当于
pair_dist = d.cdist(position, position)
pair_d = np.triu(pair_dist<=4, k=1) # k=1 to exclude the diagonal
for i, j in zip(*np.nonzero(pair_d)):
...
第二个问题（如您所述）是速度只是切换和缩放而不是反射。我们真正想要做的是否定和缩放每个粒子速度沿连接它们的轴的分量。请注意，为此我们需要连接它们的向量 position[j] - position[i] 和连接它们的向量的长度（我们已经计算过了）。所以不幸的是 cdist 计算的一部分被重复了。让我们停止使用 cdist 并自己动手吧。这里的目标是制作两个数组 diff 和 norm，其中 diff[i][j] 是从粒子 i 指向 j 的向量（因此 diff 是一个 3D 数组）和 norm[i][j] 是粒子 i 和 j 之间的距离。我们可以像这样用 numpy 做到这一点：
nop = number_of_particles
# Give pos a 3rd index so we can use np.repeat below
# equivalent to `pos3d = np.array([ position ])
pos3d = position.reshape(1, nop, 2)
# 3D arras with a repeated index so we can form combinations
# diff_i[i][j] = position[i] (for all j)
# diff_j[i][j] = position[j] (for all i)
diff_i = np.repeat(pos3d, nop, axis=1).reshape(nop, nop, 2)
diff_j = np.repeat(pos3d, nop, axis=0)
# diff[i][j] = vector pointing from position[i] to position[j]
diff = diff_j - diff_i
# norm[i][j] = sqrt( diff[i][j]**2 )
norm = np.linalg.norm(diff, axis=2)
# check for collisions and take the region above the diagonal
collided = np.triu(norm < radius, k=1)
for i, j in zip(*np.nonzero(collided)):
# unit vector from i to j
unit = diff[i][j] / norm[i][j]
# flip velocity
velocity[i] -= 1.9 * np.dot(unit, velocity[i]) * unit
velocity[j] -= 1.9 * np.dot(unit, velocity[j]) * unit
# push particle j to be radius units from i
# This isn't particularly effective when 3+ points are close together
position[j] += (radius - norm[i][j]) * unit
...
因为这个 post 已经足够长了，here is a gist 的代码加上我的修改。