为什么在我的 n 体积分器中列表的执行速度比 numpy 数组更快？

Question

我使用Python中的欧拉方法制作了一个简单的n体积分器。我不知道

numpy

数组，所以我将每个向量视为 x、y 和 z 变量的列表。我什至从头开始制作了自己的

magnitude()

函数，并将这些值加在一起。有人建议我使用

numpy

数组：它会让我的代码看起来更干净，我不会经常处理三个变量，而且会更快。确实是所有这些事情，只是速度较慢。慢 2 倍。

这是带有列表的代码：

print('Initializing...')

import math
import time
import calendar

dt=1

class Body:
  def __init__(self, pos, vel, acc, GM, color):
    self.pos = pos
    self.vel = vel
    self.acc = acc
    self.GM = GM
    self.color = color

sun=Body(
  [0,0,0],
  [0,0,0],
  [],
  1.32712440041279419e20, #1.32712440018e20 #1.3271244004193938e20
  'orange')  
mercury=Body(
  [-4.108411877039495E+07,2.997375954154480E+07,6.217890408222714E+06],
  [-3.865743010383652E+01,-3.733889075044869E+01,4.944436024774976E-01],
  [],
  2.2031868551E+13,
  'black')
venus=Body(
  [-1.069987422398024E+08,-1.145572515113905E+07,6.016588327139664E+06],
  [3.513460276994624E+00,-3.497755629371660E+01,-6.830913209445484E-01],
  [],
  3.24858592000E+14,
  'yellow')
earth=Body(
  [-2.481099325965390E+07,1.449948612736719E+08,-8.215203670851886E+03],
  [-2.984146365518679E+01,-5.126262286859617E+00,1.184224839788195E-03],
  [],
  3.98600435507E+14, #3.98600435436E+14
  'blue')
mars=Body(
  [-4.388577457389553E+07,-2.170849264747288E+08,-3.473007284527391E+06],
  [2.466191455174334E+01,-5.594344409323426E+00,4.240886274745919E-01],
  [],
  4.2828375816E+13, #4.2828375214E+13
  'red')
jupiter=Body( #Barycenter
  [5.225707808892267E+08,5.318268084719173E+08,-1.390073672466460E+07],
  [-9.479547147917305E+00,9.781812771839792E+00,1.714496767625100E-01],
  [],
  1.26712764100000E+17, #1.26686531900E+17
  'brown')
saturn=Body( #Barycenter
  [1.345793470545739E+09,-5.559292621384816E+08,-4.389272824857706E+07],
  [3.145491294741479E+00,8.918908096728853E+00,-2.803661611085468E-01],
  [],
  3.7940584841800E+16, #3.7931206234E+16
  'green')
uranus=Body( #Barycenter
  [1.835714287613103E+09,2.288891426581001E+09,-1.529866616820955E+07],
  [-5.371909870919789E+00,3.954382268177655E+00,8.420336496942427E-02],
  [],
  5.794556400000E+15, #5.793951256E+15
  'cyan')
neptune=Body( #Barycenter
  [4.464446478443181E+09,-2.679156449818206E+08,-9.736552658539964E+07],
  [2.818332382913437E-01,5.469932267577336E+00,-1.190017337507621E-01],
  [],
  6.836527100580E+15, #6.83509997E+15
  'darkblue')
moon=Body(
  [-2.517894579079584E+07,1.451613930689555E+08,1.696211394185573E+04],
  [-3.025123150480888E+01,-6.001889331002243E+00,-5.808094571415667E-02],
  [],
  4.902800118E+12, #4.902800066E+12
  'gray')
#Pos Vel Data: https://ssd.jpl.nasa.gov/horizons/app.html#/
#GM Data: https://ssd.jpl.nasa.gov/astro_par.html


bodies=[sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,moon]

#functions

def magnitude(x,y,z):
  return math.sqrt(x**2+y**2+z**2)


def grav_acc(body1,body2):
  r = magnitude(
    body2.pos[0] - body1.pos[0],
    body2.pos[1] - body1.pos[1],
    body2.pos[2] - body1.pos[2])
  a=body2.GM/math.pow(r*1000,2)/1000
  return (
    (body2.pos[0]-body1.pos[0])/r*a,
    (body2.pos[1]-body1.pos[1])/r*a,
    (body2.pos[2]-body1.pos[2])/r*a)


def total_acc(body_number):
    body=bodies[body_number]
    body.acc=[0,0,0]
    for i in range(len(bodies)):
      if i != body_number:
        acc_data=grav_acc(body,bodies[i])
        body.acc[0]+=acc_data[0]
        body.acc[1]+=acc_data[1]
        body.acc[2]+=acc_data[2]
        
      

def suvat(body_number):
    total_acc(body_number)
    body=bodies[body_number]
   
    body.pos[0]+=body.vel[0]*dt+0.5*body.acc[0]*math.pow(dt,2)
    body.pos[1]+=body.vel[1]*dt+0.5*body.acc[1]*math.pow(dt,2)
    body.pos[2]+=body.vel[2]*dt+0.5*body.acc[2]*math.pow(dt,2)
  
    body.vel[0]+=body.acc[0]*dt
    body.vel[1]+=body.acc[1]*dt
    body.vel[2]+=body.acc[2]*dt


def timestep():
  for i in range(len(bodies)):
    suvat(i)

print('Initialization Complete!\n\nInitial Data:\n')
print(earth.pos)
print(earth.vel)
print(earth.acc,'\n')


days=input('Days: ')
start_time=calendar.timegm(time.strptime("01 Jan 24", "%d %b %y"))

print('\nStart!\n')
start_time = time.time()
for x in range (0,int(float(days)*(24*60*60)/dt)):
  if x*dt/24/60/60 % 1 == 0:
    print("Day",int(x*dt/24/60/60))
    print("Earth:\n",earth.pos)
    print(earth.vel)
    print(earth.acc)
    print("\nSun:")
    print(sun.pos)
    print(sun.vel)
    print(sun.acc)
  else:
    timestep() #<-- the important thing

print(int(days),'\nComplete!\n\nEarth Data:\n')
end_time = time.time()
print("Total Time: ", end_time-start_time)
print(earth.pos)
print(earth.vel)
print(earth.acc)
print("\nSun:")
print(sun.pos)
print(sun.vel)
print(sun.acc)

这是使用 numpy 的代码：

print('Initializing...')

import math
import time
import calendar

import numpy as np


dt=1

class Body:
  def __init__(self, pos, vel, acc, GM, color):
    self.pos = pos
    self.vel = vel
    self.acc = acc
    self.GM = GM
    self.color = color

sun=Body(
  np.array([0,0,0]),
  np.array([0,0,0]),
  np.array([0,0,0]),
  1.32712440041279419e20, #1.32712440018e20 #1.3271244004193938e20
  'orange')  
mercury=Body(
  np.array([-4.108411877039495E+07,2.997375954154480E+07,6.217890408222714E+06]),
  np.array([-3.865743010383652E+01,-3.733889075044869E+01,4.944436024774976E-01]),
  np.array([0,0,0]),
  2.2031868551E+13,
  'black')
venus=Body(
  np.array([-1.069987422398024E+08,-1.145572515113905E+07,6.016588327139664E+06]),
  np.array([3.513460276994624E+00,-3.497755629371660E+01,-6.830913209445484E-01]),
  np.array([0,0,0]),
  3.24858592000E+14,
  'yellow')
earth=Body(
  np.array([-2.481099325965390E+07,1.449948612736719E+08,-8.215203670851886E+03]),
  np.array([-2.984146365518679E+01,-5.126262286859617E+00,1.184224839788195E-03]),
  np.array([0,0,0]),
  3.98600435507E+14, #3.98600435436E+14
  'blue')
mars=Body(
  np.array([-4.388577457389553E+07,-2.170849264747288E+08,-3.473007284527391E+06]),
  np.array([2.466191455174334E+01,-5.594344409323426E+00,4.240886274745919E-01]),
  np.array([0,0,0]),
  4.2828375816E+13, #4.2828375214E+13
  'red')
jupiter=Body( #Barycenter
  np.array([5.225707808892267E+08,5.318268084719173E+08,-1.390073672466460E+07]),
  np.array([-9.479547147917305E+00,9.781812771839792E+00,1.714496767625100E-01]),
  np.array([0,0,0]),
  1.26712764100000E+17, #1.26686531900E+17
  'brown')
saturn=Body( #Barycenter
  np.array([1.345793470545739E+09,-5.559292621384816E+08,-4.389272824857706E+07]),
  np.array([3.145491294741479E+00,8.918908096728853E+00,-2.803661611085468E-01]),
  np.array([0,0,0]),
  3.7940584841800E+16, #3.7931206234E+16
  'green')
uranus=Body( #Barycenter
  np.array([1.835714287613103E+09,2.288891426581001E+09,-1.529866616820955E+07]),
  np.array([-5.371909870919789E+00,3.954382268177655E+00,8.420336496942427E-02]),
  np.array([]),
  5.794556400000E+15, #5.793951256E+15
  'cyan')
neptune=Body( #Barycenter
  np.array([4.464446478443181E+09,-2.679156449818206E+08,-9.736552658539964E+07]),
  np.array([2.818332382913437E-01,5.469932267577336E+00,-1.190017337507621E-01]),
  np.array([0,0,0]),
  6.836527100580E+15, #6.83509997E+15
  'darkblue')
moon=Body(
  np.array([-2.517894579079584E+07,1.451613930689555E+08,1.696211394185573E+04]),
  np.array([-3.025123150480888E+01,-6.001889331002243E+00,-5.808094571415667E-02]),
  np.array([0,0,0]),
  4.902800118E+12, #4.902800066E+12
  'gray')
#Pos Vel Data: https://ssd.jpl.nasa.gov/horizons/app.html#/
#GM Data: https://ssd.jpl.nasa.gov/astro_par.html


bodies=[sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,moon]

#functions

def grav_acc(body1,body2):
  r = np.linalg.norm(body2.pos-body1.pos)
  a=body2.GM/math.pow(r*1000,2)/1000
  return ((body2.pos-body1.pos)/r*a)


def total_acc(body_number):
    body=bodies[body_number]
    body.acc=[0,0,0]
    for i in range(len(bodies)):
      if i != body_number:
        acc_data=grav_acc(body,bodies[i])
        body.acc = body.acc + acc_data
        
      

def suvat(body_number):
    total_acc(body_number)
    body=bodies[body_number]
   
    body.pos= body.pos + body.vel*dt+ 0.5 * body.acc* math.pow(dt,2)
    body.vel=body.vel + body.acc*dt


def timestep():
  for i in range(len(bodies)):
    suvat(i)

 
print('Initialization Complete!\n\nInitial Data:\n')
print(earth.pos)
print(earth.vel)
print(earth.acc,'\n')


days=input('Days: ')
start_time=calendar.timegm(time.strptime("01 Jan 24", "%d %b %y"))

print('\nStart!\n')
start_time = time.time()
for x in range (0,int(float(days)*(24*60*60)/dt)):
  if x*dt/24/60/60 % 1 == 0:
    print("Day",int(x*dt/24/60/60))
    print("Earth:\n",earth.pos)
    print(earth.vel)
    print(earth.acc)
    print("\nSun:")
    print(sun.pos)
    print(sun.vel)
    print(sun.acc)
  else:
    timestep() #<-- the important thing

print(int(days),'\nComplete!\n\nEarth Data:\n')
end_time = time.time()
print("Total Time: ", end_time-start_time)
print(earth.pos)
print(earth.vel)
print(earth.acc)
print("\nSun:")
print(sun.pos)
print(sun.vel)
print(sun.acc)

我使用列表模拟一天的结果是：

Complete!

Earth Data:

Total Time:  23.128202438354492
[-27385404.548457384, 144529451.5861996, -8105.186673260389]
[-29.74987226023974, -5.646917018209525, 0.0013563469814852206]
[1.112851770138038e-06, -6.019252799282894e-06, 1.7747241298097921e-09]

对于

numpy

：

Complete!

Earth Data:

Total Time:  52.90268087387085
[-2.73854045e+07  1.44529452e+08 -8.10518667e+03]
[-2.97498723e+01 -5.64691702e+00  1.35634698e-03]
[ 1.11285177e-06 -6.01925280e-06  1.77472413e-09]

正如你所看到的，相同的结果（虽然更少的数字......？），但时间更长。

一些注意事项： dt

• 两次模拟均需 1 秒
• 除了使用 numpy 数组而不是列表之外，所有代码都是相同的。
• 我没有将新版本的

  magnitude()

  函数与

  numpy

  一起使用，而是使用

  linalg.norm

  。我预计这会更快，也许事实并非如此？
• 所有数组都是 3 维（x，y，z），也许较小的数组没有性能优势？
• 我正在做的所有 numpy 计算都是

  array = array + something

  、

  array = array * scalar

  、

  array = array / scalar

  或

  linalg.norm(array)

我计划在未来尝试多处理和其他改进来加快模拟时间，所以我想做出最适合我的代码速度的选择。快速谷歌搜索表明 numpy 应该更快，但我的结果与此不匹配。

Answer 1

与列表相比，Numpy 应该更快，但您需要使用矢量化来实现。您可以在代码中进行一些优化：

1. 您可以重写 grav_acc 和total_acc 函数，因为它们占用了您的大部分代码时间（使用 https://kernprof.readthedocs.io/en/latest/ 来分析您的代码）。据我了解，您正在计算范数并在这里进行标准化。目前，每对都进行两次！你可以减少一半。
2. 最好将 numpy 与向量化一起使用，否则优势将很小。例如。使用 numpy，您不需要单独计算每个向量对的范数。您可以定义一组“pos”向量并一次性计算所有对。

希望这有帮助。