[Python] Why are lists performing faster than numpy arrays in my n-body integrator?

I made a simple n-body integrator using Euler's method in python. I didn't know about numpy arrays, so I treated each vector as a list of x, y, and z variables. I even made my own magnitude() function from scratch to square and add the values together. I was suggested to use numpy arrays: it would make my code look cleaner, I wouldn't constantly be dealing with three variables, and it would be faster. It was indeed all of those things, except it was slower. 2x slower.

Here's the code for the one with lists:

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)
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)



And here's the code for the one using 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)
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)



My results for one simulated day using lists was:

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]


and for 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]


As you can see, same results (although less sig figs...?), but much longer time.

A few notes: dt

• 1 second for both simulations
• All the code is the same, other than use of numpy arrays instead of lists.
• Instead of using my magnitude() function for the new version with numpy, I used linalg.norm. I expected this to be faster, maybe it isn't?
• All the arrays are 3 dimensional (x,y,z), maybe smaller arrays don't have performance benefits?
• All the numpy calculations I am doing are array = array + something, array = array * scalar, array = array / scalar, or linalg.norm(array)

I plan to try multiprocessing and other improvements to speed up simulation time in the future, so I want to make choices that are best for the speed of my code. A quick google search says numpy should be faster, but my results don't match that.

Continue reading...