Почему в моем интеграторе n-body списки работают быстрее, чем массивы numpy? ⇐ Python

[Anonymous]

Я создал простой интегратор из n тел, используя метод Эйлера в Python. Я не знал о массивах numpy, поэтому рассматривал каждый вектор как список переменных x, y и z. Я даже создал свою собственную функцию value() с нуля до квадрата и сложил значения вместе. Мне предложили использовать массивы 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() #

Подробнее здесь: [url]https://stackoverflow.com/questions/79004875/why-are-lists-performing-faster-than-numpy-arrays-in-my-n-body-integrator[/url]