我试图计算由于重力引起的加速度,对于一个三空间的n体问题(我使用辛欧拉)。在
我有每个时间步的位置和速度向量,并使用下面的(工作)代码来计算加速度并更新速度和位置。注意加速度是3维空间中的向量,而不仅仅是量值。在
我想知道是否有一种更有效的方法来用numpy来计算,以避免循环。在
def accelerations(positions, masses):
'''Params:
- positions: numpy array of size (n,3)
- masses: numpy array of size (n,)
Returns:
- accelerations: numpy of size (n,3), the acceleration vectors in 3-space
'''
n_bodies = len(masses)
accelerations = numpy.zeros([n_bodies,3]) # n_bodies * (x,y,z)
# vectors from mass(i) to mass(j)
D = numpy.zeros([n_bodies,n_bodies,3]) # n_bodies * n_bodies * (x,y,z)
for i, j in itertools.product(range(n_bodies), range(n_bodies)):
D[i][j] = positions[j]-positions[i]
# Acceleration due to gravitational force between each pair of bodies
A = numpy.zeros((n_bodies, n_bodies,3))
for i, j in itertools.product(range(n_bodies), range(n_bodies)):
if numpy.linalg.norm(D[i][j]) > epsilon:
A[i][j] = gravitational_constant * masses[j] * D[i][j] \
/ numpy.linalg.norm(D[i][j])**3
# Calculate net acceleration of each body (vectors in 3-space)
accelerations = numpy.sum(A, axis=1) # sum of accel vectors for each body of shape (n_bodies,3)
return accelerations
我对你原帖的评论如下:
这是一个使用
blas
的优化版本。blas
对对称矩阵或厄米特矩阵上的线性代数有特殊的例程。它们使用专用的压缩存储,只保留上三角或下三角,而不保留(冗余)镜像条目。这样,blas不仅节省了一半的存储空间,而且节省了大约一半的失败。在我写了不少评论让它可读。在
示例运行;与OP相比,我的纯numpy向量化和@p Mende的
^{pr2}$我们可以看到
1)p Mende在矢量化方面略优于I
2)
blas
的速度是~5倍;请注意,我的blas不是很好;我怀疑使用优化的BLA,你可能会变得更好(不过,在更好的BLA上,numpy也会运行得更快)3)任何一个答案都比循环快得多
需要考虑的一些事项:
您只需要一半的距离;一旦您计算了
D[i][j]
,这与-D[j][i]
相同。在你可以做
df2 = df.apply(lambda x:gravitational_constant/x**3)
你可以为每一个产品生成一对数据体。您只需执行一次,然后每次调用它时都可以将其传递给
accelearations
。在然后
df.product(df2).product(mass_products).sum().div(masses)
给你加速度。在相关问题 更多 >
编程相关推荐