当您使用数字网格并希望获得非常好的性能时,您应该考虑使用Numpy。它使用起来非常简单,让您可以考虑使用网格而不是网格上的循环来进行操作。性能来自这样一个事实,即这些操作随后使用优化的 SSE 代码在整个网格上运行。
例如,这里有一些使用我编写的代码的 numpy,它对通过弹簧连接的带电粒子进行蛮力数值模拟。此代码在 31 毫秒内计算具有 100 个节点和 99 条边的 3d 系统的时间步长。这比我能想到的最好的纯 Python 代码快 10 倍以上。
from numpy import array, sqrt, float32, newaxis
def evolve(points, velocities, edges, timestep=0.01, charge=0.1, mass=1., edgelen=0.5, dampen=0.95):
"""Evolve a n body system of electrostatically repulsive nodes connected by
springs by one timestep."""
velocities *= dampen
# calculate matrix of distance vectors between all points and their lengths squared
dists = array([[p2 - p1 for p2 in points] for p1 in points])
l_2 = (dists*dists).sum(axis=2)
# make the diagonal 1's to avoid division by zero
for i in xrange(points.shape[0]):
l_2[i,i] = 1
l_2_inv = 1/l_2
l_3_inv = l_2_inv*sqrt(l_2_inv)
# repulsive force: distance vectors divided by length cubed, summed and multiplied by scale
scale = timestep*charge*charge/mass
velocities -= scale*(l_3_inv[:,:,newaxis].repeat(points.shape[1], axis=2)*dists).sum(axis=1)
# calculate spring contributions for each point
for idx, (point, outedges) in enumerate(izip(points, edges)):
edgevecs = point - points.take(outedges, axis=0)
edgevec_lens = sqrt((edgevecs*edgevecs).sum(axis=1))
scale = timestep/mass
velocities[idx] += (edgevecs*((((edgelen*scale)/edgevec_lens - scale))[:,newaxis].repeat(points.shape[1],axis=1))).sum(axis=0)
# move points to new positions
points += velocities*timestep