我的Python代码计算太慢！可以用numpy加速吗？怎么做？

Question

我想用Verlet算法来计算行星的运动。

问题是，我的代码运行一个完整的周期大约需要30分钟。这是基本的计算，我不知道为什么会这么慢。

我在这里看了一些资料，听说用numpy应该会更快，对吧？

但是我该怎么用它呢？另外，代码应该是计算完后回到起点并停止在0（或者2*pi*r_earth，如果是在计算行驶的距离），但它根本不停止，只是继续运行，有人能告诉我为什么我的限制没有起作用吗？

这是我的代码：

grav    =   6.673481*(10**-11)             # = 6.673481e-11  # a direct way
m_sun   =   1.989*(10**30)                 # = 1.989e+30
m_earth =   5.972*(10**24)                 # = 5.972e+24    
r_earth = 149.59787*(10**9)                # = 1.4959787e+11

def verlet_v( prev_v, prev_f, mass ):

    current_v = ( prev_v + ( prev_f / mass ) )
    print "New Velocity: %f" % current_v
    return current_v

def verlet_x( prev_x, current_v ):

    current_x = prev_x + current_v
    print "New Position: %f" % current_x
    return current_x

verlet_v( 20, 50, 3 )

v = 29.8*(10**3)                          # = 2.98e+04
x = 0
f = ( -grav * ( ( m_earth * m_sun ) / r_earth**2 ) )
v_history = []
x_history = []

while( abs(x) > ( -0.1 ) ):               # ABS( <_anything_> ) is ALWAYS > -0.1

    print "Mod(x): %f" % abs(x)
    print "Limit: %f"  % (0)
    v = verlet_v( v, f, m_earth )
    x = verlet_x( x, v )
    v_history.append(v)
    x_history.append(x)

print v_history
print x_history

Answer 1

一个建议是在你的循环之前定义“历史”数组，而不是每次循环都往里面添加数据。这样的话，内存会提前分配好，这样可以提高程序的运行效率。不过，你需要提前计算好 v_history 和 x_history 的大小。

比如，使用一维数组：

v_history = np.zeros((k,))
x_history = np.zeros((k,))

这里的 (k,) 是数组的形状。

然后你需要用一个索引值来把计算出来的值存储到这个数组里。

# outside loop
x=0

# inside loop
v_history[x] = v
x +=1

另外，你可能还想开始了解一下 广播 的概念。

Answer 2

看起来你陷入了一个无限循环。

abs(x) 是用来计算一个数字的绝对值的，绝对值总是正数。这意味着它的值永远不会低于 -0.1，所以你的循环永远不会结束。

你需要把 abs(x) 或者 -0.1 改成其他的值。