我目前正在尝试模拟盒子周围弹跳的许多粒子.
我考虑了@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 … props @kalhartt的一些知识),但我仍然不知道它在哪里搞砸了.我尝试更改实际更新的顺序,其中成对距离最后一次或position = velocity最后一次但无济于事.我添加了* 0.9以使整个过程更快消失,并尝试使用4来确保2 *半径(= 2)不太严格…但是似乎没有任何效果.
任何和所有帮助将不胜感激.
解决方法:
仅有两种错字妨碍您前进.首先是范围(len(position)/ 2)中的i：仅迭代一半以上的粒子.这就是为什么一半的粒子保留在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>一世.首先,您的代码实际上通过在range(i)中使用j来遍历较低的三角形区域. i,在这种情况下并不特别重要,因为不会重复i,j对.但是Numpy有两个内置函数可以替代,np.triu可以将对角线以下的所有值都设置为0,而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)):
...
第二个问题(如您所述)是速度只是被切换和缩放而不是被反映.我们真正想要做的是沿连接它们的轴求和否定每个粒子速度的分量.请注意,要做到这一点,我们将需要将它们连接到位置[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
...
由于这篇文章已经足够长了,我对代码的修改here is a gist.
标签：scipy,python-2-7,physics,python,numpy
