我目前正在嘗試模擬盒子周圍彈跳的許多粒子.

我考慮了@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

總結

以上是生活随笔為你收集整理的python 粒子动画_python-盒子中有很多粒子-物理模拟的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。