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

我考慮了@kalhartt的建議,這是改進(jìn)的代碼,用于初始化框內(nèi)的粒子：

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)

初始化的位置如下所示：

初始化粒子后,我想在每個(gè)時(shí)間步更新它們.用于更新的代碼立即遵循先前的代碼,如下所示：

# 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

更新它并使其運(yùn)行完成后,我得到以下結(jié)果：

這比以前要好得多,但是仍然有一些補(bǔ)丁之間的距離太近-例如：

太近了.我認(rèn)為更新是有效的…感謝@kalhartt,我的代碼更好,更快了(而且我學(xué)到了有關(guān)numpy … props @kalhartt的一些知識(shí)),但我仍然不知道它在哪里搞砸了.我嘗試更改實(shí)際更新的順序,其中成對(duì)距離最后一次或position = velocity最后一次但無(wú)濟(jì)于事.我添加了* 0.9以使整個(gè)過(guò)程更快消失,并嘗試使用4來(lái)確保2 *半徑(= 2)不太嚴(yán)格…但是似乎沒(méi)有任何效果.

任何和所有幫助將不勝感激.

解決方法:

僅有兩種錯(cuò)字妨礙您前進(jìn).首先是范圍(len(position)/ 2)中的i：僅迭代一半以上的粒子.這就是為什么一半的粒子保留在x邊界中的原因(如果您觀察較大的迭代,則更清晰).第二,第二個(gè)y條件應(yīng)為最小(我假設(shè))position [i] [1]< 0.下面的塊為我綁定了粒子(我沒(méi)有使用碰撞代碼進(jìn)行測(cè)試,因此可能存在問(wèn)題).

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]

順便說(shuō)一句,嘗試盡可能利用numpy消除循環(huán).它更快,更高效,并且在我看來(lái)更具可讀性.例如,force_init看起來(lái)像這樣：

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

最后說(shuō)明,用position [xmax,0] = limit_x之類的東西將位置硬剪切到邊界框可能是一個(gè)好主意. position [xmin,0] =0.在某些情況下,速度較小,框外的粒子將被反射,但在下一次迭代中不會(huì)進(jìn)入框內(nèi).因此它將永遠(yuǎn)位于被永久反射的盒子外面.

編輯：碰撞

碰撞檢測(cè)是一個(gè)困難得多的問(wèn)題,但讓我們看看我們能做什么.讓我們看一下您當(dāng)前的實(shí)現(xiàn).

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

總體而言,這是一個(gè)非常好的方法,cdist將有效地計(jì)算距離

在點(diǎn)集之間,您會(huì)發(fā)現(xiàn)哪些點(diǎn)與pair_d = pair_dist< = 4發(fā)生碰撞.

嵌套的for循環(huán)是第一個(gè)問(wèn)題.我們需要遍歷pair_d的True值,其中j>一世.首先,您的代碼實(shí)際上通過(guò)在range(i)中使用j來(lái)遍歷較低的三角形區(qū)域. i,在這種情況下并不特別重要,因?yàn)椴粫?huì)重復(fù)i,j對(duì).但是Numpy有兩個(gè)內(nèi)置函數(shù)可以替代,np.triu可以將對(duì)角線以下的所有值都設(shè)置為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]:

...

相當(dāng)于

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)):

...

第二個(gè)問(wèn)題(如您所述)是速度只是被切換和縮放而不是被反映.我們真正想要做的是沿連接它們的軸求和否定每個(gè)粒子速度的分量.請(qǐng)注意,要做到這一點(diǎn),我們將需要將它們連接到位置[j]-position [i]的向量和連接它們的向量的長(zhǎng)度(我們已經(jīng)計(jì)算出).因此,不幸的是,部分cdist計(jì)算被重復(fù)了.讓我們退出使用cdist并自己做.這里的目標(biāo)是制作兩個(gè)數(shù)組diff和norm,其中diff [i] [j]是從粒子i指向j的向量(因此diff是3D數(shù)組),而norm [i] [j]是粒子i之間的距離和j.我們可以這樣用numpy做到這一點(diǎn)：

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

...

由于這篇文章已經(jīng)足夠長(zhǎng)了,我對(duì)代碼的修改here is a gist.

標(biāo)簽：scipy,python-2-7,physics,python,numpy

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。