1

我有一个适用于数千个粒子的 N 体模拟脚本。它输出粒子位置的最终 2D 投影，我想围绕给定轴以不同的旋转角度多次运行它，以便能够从不同的角度看到最终的 2D 结果。为此，我在开头添加了一些代码：

# Euler-Rodrigues formula definition for arbitrary 3D rotation
def rotation_matrix(axis,angle):
    axis = axis/math.sqrt(dot(axis,axis))
    a = math.cos(angle/2)
    b,c,d = -axis*math.sin(angle/2)
    return array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                  [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                  [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

接着：

# 3D rotation
angle = float(sys.argv[7])*math.pi/180.0
axis = array([0,0,1])
xAr,yAr,zAr=[],[],[]           #arrays for rotated particles
for xi,yi,zi in zip(xA,yA,zA): #cycle to rotate every particle
    ci = array([xi,yi,zi])
    cr = dot(rotation_matrix(axis,angle),ci)
    xAr.append(cr[0])
    yAr.append(cr[1])
    zAr.append(cr[2])
xA,yA,zA = array(xAr), array(yAr), array(zAr)

脚本的核心部分在此之后开始。总之，最初的脚本是这样做的：

  • 读取模型 > 运行模拟 > 输出 2D 投影

我添加了我的部分，所以现在：

  • 读取模型 >旋转 3D 中的每个粒子> 运行模拟 > 输出 2D 投影

但我发现轮换过程需要太多时间。Python中有没有办法或不同的方法来加速这个？（如果有帮助，我只希望最终的 2D 投影沿给定轴旋转）。

4

2 回答 2

2

您需要对操作进行矢量化，即在整个粒子集上运行它，而不是在它们上循环。假设xA等是一维数组，那就是：

particles = np.vstack([xA, xB, xC])
rot = rotation_matrix(axis, angle)
rotated = np.dot(rot, particles)
xA, yA, zA = rotated

这应该会给你几个数量级的加速。在更一般的说明中，您在循环的每次迭代中构造相同的旋转矩阵。那很浪费。

于 2014-02-04T15:16:14.797 回答
1

把rotation_matrix(axis,angle) 放在紧循环之外怎么样？

于 2014-02-04T15:19:29.397 回答