Forward Euler N-Body Simulation в Python не возвращает правильные результаты

Question

Forward Euler N-Body Simulation

Эта программа берет позиции и скорости планет солнечной системы в конкретную дату из эфемерид, а затем моделирует движения планет, используя метод прямой эйлеровой интеграции.

Фундаментальная физика, лежащая в основе этого моделирования, - это закон Ньютона о всеобщей гравитации.

Для каждой итерации интеграции кумулятивно учитывается влияние каждой j-й планеты на i-ю планету (где j!= I).

#Christopher Iliffe Sprague, Solar System Simulation
#This simulation takes position and velocity valuse of various planets at a particular point in time, and then runs a forward Euler approximation to simulate actual planetary motion.
#This position and velocity values are supplied by NASA's Jet Propulsion Laboratory via the ephemeris module 'jplephem'

#Modules (Obtaining the ephemeris data requires Numpy)
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

# Ephemeris (positions and velocities of each body)
import de423; from jplephem import Ephemeris; eph = Ephemeris(de423)

#Initial Planetary Values
t0 = 2457180.500000 #Julian date, CE 2015 June 7 00:00:00.0 UT
Planets=['sun','mercury','venus','earthmoon','mars'] #Planet names (for use in jplephem)
Masses=[1.989E30,328.5E21,4.867E24,5.972E24,7.34767309E22,639E21,1.898E27,568E24,86.8E24,102E24,0.0131E24] #Mass of each body [kg]
G=6.67384E-11 #Gravitational Constant [m^3kg^-1s^-2]
x0=[]; y0=[]; z0=[]; vx0=[]; vy0=[]; vz0=[]; #Lists of initial values for each body
for i in range(len(Planets)):
    position, velocity = eph.position_and_velocity(Planets[i], t0) #[km], [km/day]
    x,y,z = position*(1000./1.) #[m]
    x0.append(x[0]); y0.append(y[0]); z0.append(z[0]); 
    vx,vy,vz = velocity*(1./24.)*(1./60.)*(1./60.)*(1000./1.) #[m/s]
    vx0.append(vx[0]); vy0.append(vy[0]); vz0.append(vz[0]);

#Create Arrays
N=100; nl=range(0,N); tll=[0]*len(nl); trl=[0]*len(nl); trl[0]=t0
X=[]; Y=[]; Z=[]; VX=[]; VY=[]; VZ=[]
X.append(x0); Y.append(y0); Z.append(z0); VX.append(vx0); VY.append(vy0); VZ.append(vz0)

#Equations of Motion
def dvxdt(xi,yi,zi,xj,yj,zj,mj): return (-G*mj*(xi-xj))/(((xi-xj)**2.+(yi-yj)**2.+(zi-zj)**2.)**(3./2.)) #x-velocity [m/s^2]
def dvydt(xi,yi,zi,xj,yj,zj,mj): return (-G*mj*(yi-yj))/(((xi-xj)**2.+(yi-yj)**2.+(zi-zj)**2.)**(3./2.)) #y-velocity [m/s^2]
def dvzdt(xi,yi,zi,xj,yj,zj,mj): return (-G*mj*(zi-zj))/(((xi-xj)**2.+(yi-yj)**2.+(zi-zj)**2.)**(3./2.)) #z-velocity [m/s^2]
def dxdt(vx): return vx #[m/s] #Redundant
def dydt(vy): return vy #[m/s] #Redundant
def dzdt(vz): return vz #[m/s] #Redundant

#Numerical Integration
dt=10 #Time step [s]
for n in nl[:-1]: #For each iteration
    tll[n+1]=tll[n]+dt #Increment time
    xn=[]; yn=[]; zn=[]; vxn=[]; vyn=[]; vzn=[] #Lists for each time iteration
    for i in range(len(Planets)): #For each planet
        xi=X[n][i]+dt*VX[n][i]; yi=Y[n][i]+dt*VY[n][i]; zi=Z[n][i]+dt*VZ[n][i]
        vxi=0; vyi=0; vzi=0 #Initial velocity values for each planet in each time step set to zero to accomidate the summation
        for j in range(len(Planets)): #For each other planet j in the context of planet i
            if i!=j: #Only consider the planets other than planet i
                #Add the effects of each other planet
                vxi+=VX[n][i]+dt*dvxdt(X[n][i],Y[n][i],Z[n][i],X[n][j],Y[n][j],Z[n][j],Masses[j])
                vyi+=VY[n][i]+dt*dvydt(X[n][i],Y[n][i],Z[n][i],X[n][j],Y[n][j],Z[n][j],Masses[j])
                vzi+=VZ[n][i]+dt*dvzdt(X[n][i],Y[n][i],Z[n][i],X[n][j],Y[n][j],Z[n][j],Masses[j])
        #Fill planetary values for the nth iteration
        xn.append(xi), yn.append(yi); zn.append(zi); vxn.append(vxi); vyn.append(vyi); vzn.append(vzi)
    #Append the nth iteration's values to the master list
    X.append(xn); Y.append(yn); Z.append(zn); VX.append(vxn); VY.append(vyn); VZ.append(vzn)
    #For example X[1][1] would return the Mercury's X-position in the second time step

Результаты

for n in X:
    print n

 [504036180.92439699, -881333815.26759899, -91468863977.716003, -36868767754.599091, 43716515029.563995]
 [504036226.3350125, -880944433.41161776, -91468682880.980255, -36868483835.137222, 43716286270.288544]
 [504036407.97747171, -879386905.93280208, -91467958493.070251, -36867348157.147964, 43715371233.140121]
 [504037134.5473057, -873156795.96266425, -91465060940.463226, -36862805445.049141, 43711711084.499817]
 [504040040.82663882, -848236356.02729917, -91453470729.06813, -36844634596.51207, 43697070489.891968]
 [504051665.94396847, -748554596.23127103, -91407109882.520767, -36771951202.222023, 43638508111.413956]
 [504098166.4132843, -349827556.99357504, -91221666495.364441, -36481217624.92012, 43404258597.455299]
 [504284168.29054493, 1245080600.0068531, -90479892945.772797, -35318283315.571083, 42467260541.574135]
 [505028175.79958457, 7624713228.0424356, -87512798746.441879, -30666546078.03463, 38719268318.003204]
 [508004205.83574033, 33143243740.155304, -75644421948.162399, -12059597127.752968, 23727299423.674278]
 [519908325.98036081, 135217365788.32568, -28170914754.127251, 62368198673.491524, -36240576153.682373]
 [567524806.55884087, 543513853980.01276, 161723114022.72266, 360079381878.51501, -276112078463.13312]
 [757990728.87276161, 2176699806746.2324, 921299229130.20862, 1550924114698.4771, -1235598087700.906]
 [1519854418.1284447, 8709443617811.0674, 3959603689560.1353, 6314303045978.2676, -5073542124651.9482]
 [4567309175.1511765, 34840418862070.402, 16112821531279.84, 25367818771097.43, -20425318272456.109]
 [16757128203.242104, 139364319839107.75, 64725692898158.656, 101581881671574.06, -81832422863672.75]
 [65516404315.605812, 557459923747257.12, 259177178365673.94, 406438133273480.62, -327460841228539.37]
 [260553508765.06064, 2229842339379854.5, 1036983120235735.0, 1625863139681107.0, -1309974514688005.7]
 [1040701926562.88, 8919372001910244.0, 4148206887715979.5, 6503563165311612.0, -5240029208525871.0]
 [4161295597754.1572, 35677490652031804.0, 16593101957636958.0, 26014363267833632.0, -20960247983877332.0]
 [16643670282519.266, 1.4270996525251805e+17, 66372682237320872.0, 1.0405756367792171e+17, -83841123085283184.0]
 [66573169021579.703, 5.7083986365446298e+17, 2.6549100335605651e+17, 4.1623036531827405e+17, -3.3536462349090656e+17]
 [266291163977821.44, 2.2833594572622428e+18, 1.061964287830999e+18, 1.6649215718796833e+18, -1.3414586251134001e+18]
 [1065163143802788.5, 9.1334378316933622e+18, 4.2478574257307694e+18, 6.6596863981253202e+18, -5.3658346316033741e+18]
 [4260651063102656.5, 3.6533751329417839e+19, 1.699142997732985e+19, 2.6638745703107871e+19, -2.1463338657563271e+19]
 [17042602740302128.0, 1.4613500532031575e+20, 6.7965720183726178e+19, 1.0655498292303806e+20, -8.585335476140286e+19]
 [68170409449100016.0, 5.8454002128390737e+20, 2.7186288100931148e+20, 4.2621993180275881e+20, -3.4341341917676123e+20]
 [2.7268163628429158e+17, 2.338160085138274e+21, 1.0874515243116528e+21, 1.7048797273216419e+21, -1.3736536768381945e+21]
 [1.0907265436250578e+18, 9.3526403405557405e+21, 4.3498060975210176e+21, 6.8195189093971746e+21, -5.4946147074839276e+21]
 [4.3629061729881226e+18, 3.7410561362225604e+22, 1.7399224390358477e+22, 2.7278075637699302e+22, -2.1978458830066862e+22]
 [1.7451624690440382e+19, 1.4964224544890507e+23, 6.9596897561708314e+22, 1.0911230255090782e+23, -8.7913835320398597e+22]
 [6.9806498760249418e+19, 5.9856898179562289e+23, 2.7838759024710767e+23, 4.3644921020374188e+23, -3.5165534128172552e+23]
 [2.7922599503948556e+20, 2.3942759271824942e+24, 1.113550360988705e+24, 1.7457968408150781e+24, -1.4066213651270333e+24]
 [1.1169039801564302e+21, 9.5771037087299791e+24, 4.4542014439550949e+24, 6.983187363260423e+24, -5.6264854605082643e+24]
 [4.4676159206242087e+21, 3.8308414834919921e+25, 1.7816805775820654e+25, 2.7932749453041804e+25, -2.250594184203319e+25]
 [1.7870463682495323e+22, 1.5323365933967968e+26, 7.1267223103282893e+25, 1.1173099781216732e+26, -9.0023767368132899e+25]
 [7.1481854729979781e+22, 6.1293463735871873e+26, 2.8506889241313184e+26, 4.4692399124866941e+26, -3.6009506947253173e+26]
 [2.8592741891991758e+23, 2.4517385494348749e+27, 1.1402755696525277e+27, 1.7876959649946776e+27, -1.4403802778901269e+27]
 [1.1437096756796688e+24, 9.8069541977394997e+27, 4.5611022786101106e+27, 7.1507838599787106e+27, -5.7615211115605078e+27]
 [4.5748387027186738e+24, 3.9227816790957999e+28, 1.8244409114440442e+28, 2.8603135439914842e+28, -2.3046084446242031e+28]
 [1.8299354810874695e+25, 1.56911267163832e+29, 7.297763645776177e+28, 1.1441254175965937e+29, -9.2184337784968124e+28]
 [7.3197419243498781e+25, 6.2764506865532798e+29, 2.9191054583104708e+29, 4.5765016703863748e+29, -3.687373511398725e+29]
 [2.9278967697399512e+26, 2.5105802746213119e+30, 1.1676421833241883e+30, 1.8306006681545499e+30, -1.47494940455949e+30]
 [1.1711587078959805e+27, 1.0042321098485248e+31, 4.6705687332967533e+30, 7.3224026726181996e+30, -5.8997976182379599e+30]
 [4.684634831583922e+27, 4.0169284393940991e+31, 1.8682274933187013e+31, 2.9289610690472798e+31, -2.359919047295184e+31]
 [1.8738539326335688e+28, 1.6067713757576396e+32, 7.4729099732748052e+31, 1.1715844276189119e+32, -9.4396761891807359e+31]
 [7.4954157305342751e+28, 6.4270855030305585e+32, 2.9891639893099221e+32, 4.6863377104756478e+32, -3.7758704756722944e+32]
 [2.9981662922137101e+29, 2.5708342012122234e+33, 1.1956655957239688e+33, 1.8745350841902591e+33, -1.5103481902689177e+33]
 [1.199266516885484e+30, 1.0283336804848894e+34, 4.7826623828958754e+33, 7.4981403367610364e+33, -6.041392761075671e+33]
 [4.7970660675419361e+30, 4.1133347219395575e+34, 1.9130649531583501e+34, 2.9992561347044146e+34, -2.4165571044302684e+34]
 [1.9188264270167744e+31, 1.645333888775823e+35, 7.6522598126334006e+34, 1.1997024538817658e+35, -9.6662284177210736e+34]
 [7.6753057080670977e+31, 6.5813355551032919e+35, 3.0609039250533602e+35, 4.7988098155270633e+35, -3.8664913670884294e+35]
 [3.0701222832268391e+32, 2.6325342220413168e+36, 1.2243615700213441e+36, 1.9195239262108253e+36, -1.5465965468353718e+36]
 [1.2280489132907356e+33, 1.0530136888165267e+37, 4.8974462800853764e+36, 7.6780957048433013e+36, -6.1863861873414871e+36]
 [4.9121956531629425e+33, 4.2120547552661068e+37, 1.9589785120341505e+37, 3.0712382819373205e+37, -2.4745544749365948e+37]
 [1.964878261265177e+34, 1.6848219021064427e+38, 7.8359140481366022e+37, 1.2284953127749282e+38, -9.8982178997463793e+37]
 [7.8595130450607081e+34, 6.7392876084257709e+38, 3.1343656192546409e+38, 4.9139812510997128e+38, -3.9592871598985517e+38]
 [3.1438052180242832e+35, 2.6957150433703084e+39, 1.2537462477018563e+39, 1.9655925004398851e+39, -1.5837148639594207e+39]
 [1.2575220872097133e+36, 1.0782860173481234e+40, 5.0149849908074254e+39, 7.8623700017595405e+39, -6.3348594558376828e+39]
 [5.0300883488388532e+36, 4.3131440693924934e+40, 2.0059939963229702e+40, 3.1449480007038162e+40, -2.5339437823350731e+40]
 [2.0120353395355413e+37, 1.7252576277569974e+41, 8.0239759852918806e+40, 1.2579792002815265e+41, -1.0135775129340292e+41]
 [8.0481413581421651e+37, 6.9010305110279894e+41, 3.2095903941167523e+41, 5.0319168011261059e+41, -4.054310051736117e+41]
 [3.219256543256866e+38, 2.7604122044111958e+42, 1.2838361576467009e+42, 2.0127667204504424e+42, -1.6217240206944468e+42]
 [1.2877026173027464e+39, 1.1041648817644783e+43, 5.1353446305868036e+42, 8.0510668818017695e+42, -6.4868960827777872e+42]
 [5.1508104692109856e+39, 4.4166595270579132e+43, 2.0541378522347214e+43, 3.2204267527207078e+43, -2.5947584331111149e+43]
 [2.0603241876843943e+40, 1.7666638108231653e+44, 8.2165514089388858e+43, 1.2881707010882831e+44, -1.0379033732444459e+44]
 [8.241296750737577e+40, 7.0666552432926612e+44, 3.2866205635755543e+44, 5.1526828043531325e+44, -4.1516134929777838e+44]
 [3.2965187002950308e+41, 2.8266620973170645e+45, 1.3146482254302217e+45, 2.061073121741253e+45, -1.6606453971911135e+45]
 [1.3186074801180123e+42, 1.1306648389268258e+46, 5.2585929017208869e+45, 8.244292486965012e+45, -6.642581588764454e+45]
 [5.2744299204720493e+42, 4.5226593557073032e+46, 2.1034371606883548e+46, 3.2977169947860048e+46, -2.6570326355057816e+46]
 [2.1097719681888197e+43, 1.8090637422829213e+47, 8.413748642753419e+46, 1.3190867979144019e+47, -1.0628130542023126e+47]
 [8.4390878727552789e+43, 7.2362549691316851e+47, 3.3654994571013676e+47, 5.2763471916576077e+47, -4.2512522168092506e+47]
 [3.3756351491021116e+44, 2.894501987652674e+48, 1.346199782840547e+48, 2.1105388766630431e+48, -1.7005008867237002e+48]
 [1.3502540596408446e+45, 1.1578007950610696e+49, 5.3847991313621882e+48, 8.4421555066521722e+48, -6.8020035468948009e+48]
 [5.4010162385633785e+45, 4.6312031802442784e+49, 2.1539196525448753e+49, 3.3768622026608689e+49, -2.7208014187579204e+49]
 [2.1604064954253514e+46, 1.8524812720977114e+50, 8.6156786101795011e+49, 1.3507448810643476e+50, -1.0883205675031682e+50]
 [8.6416259817014056e+46, 7.4099250883908455e+50, 3.4462714440718004e+50, 5.4029795242573902e+50, -4.3532822700126726e+50]
 [3.4566503926805622e+47, 2.9639700353563382e+51, 1.3785085776287202e+51, 2.1611918097029561e+51, -1.741312908005069e+51]
 [1.3826601570722249e+48, 1.1855880141425353e+52, 5.5140343105148807e+51, 8.6447672388118244e+51, -6.9652516320202762e+51]
 [5.5306406282888996e+48, 4.7423520565701411e+52, 2.2056137242059523e+52, 3.4579068955247297e+52, -2.7861006528081105e+52]
 [2.2122562513155598e+49, 1.8969408226280564e+53, 8.8224548968238091e+52, 1.3831627582098919e+53, -1.1144402611232442e+53]
 [8.8490250052622393e+49, 7.5877632905122258e+53, 3.5289819587295236e+53, 5.5326510328395676e+53, -4.4577610444929767e+53]
 [3.5396100021048957e+50, 3.0351053162048903e+54, 1.4115927834918095e+54, 2.213060413135827e+54, -1.7831044177971907e+54]
 [1.4158440008419583e+51, 1.2140421264819561e+55, 5.6463711339672378e+54, 8.8522416525433082e+54, -7.1324176711887628e+54]
 [5.6633760033678332e+51, 4.8561685059278245e+55, 2.2585484535868951e+55, 3.5408966610173233e+55, -2.8529670684755051e+55]
 [2.2653504013471333e+52, 1.9424674023711298e+56, 9.0341938143475805e+55, 1.4163586644069293e+56, -1.141186827390202e+56]
 [9.0614016053885331e+52, 7.7698696094845192e+56, 3.6136775257390322e+56, 5.6654346576277172e+56, -4.5647473095608082e+56]
 [3.6245606421554132e+53, 3.1079478437938077e+57, 1.4454710102956129e+57, 2.2661738630510869e+57, -1.8258989238243233e+57]
 [1.4498242568621653e+54, 1.2431791375175231e+58, 5.7818840411824515e+57, 9.0646954522043476e+57, -7.3035956952972931e+57]
 [5.7992970274486612e+54, 4.9727165500700923e+58, 2.3127536164729806e+58, 3.625878180881739e+58, -2.9214382781189172e+58]
 [2.3197188109794645e+55, 1.9890866200280369e+59, 9.2510144658919225e+58, 1.4503512723526956e+59, -1.1685753112475669e+59]
 [9.2788752439178578e+55, 7.9563464801121477e+59, 3.700405786356769e+59, 5.8014050894107824e+59, -4.6743012449902676e+59]
 [3.7115500975671431e+56, 3.1825385920448591e+60, 1.4801623145427076e+60, 2.320562035764313e+60, -1.869720497996107e+60]
 [1.4846200390268573e+57, 1.2730154368179436e+61, 5.9206492581708304e+60, 9.2822481430572519e+60, -7.4788819919844281e+60]
 [5.938480156107429e+57, 5.0920617472717745e+61, 2.3682597032683321e+61, 3.7128992572229008e+61, -2.9915527967937713e+61]
 [2.3753920624429716e+58, 2.0368246989087098e+62, 9.4730388130733286e+61, 1.4851597028891603e+62, -1.1966211187175085e+62]
 [9.5015682497718864e+58, 8.1472987956348392e+62, 3.7892155252293314e+62, 5.9406388115566412e+62, -4.786484474870034e+62]
 [3.8006272999087546e+59, 3.2589195182539357e+63, 1.5156862100917326e+63, 2.3762555246226565e+63, -1.9145937899480136e+63]
 [1.5202509199635018e+60, 1.3035678073015743e+64, 6.0627448403669303e+63, 9.5050220984906259e+63, -7.6583751597920544e+63]
 [6.0810036798540073e+60, 5.2142712292062971e+64, 2.4250979361467721e+64, 3.8020088393962504e+64, -3.0633500639168218e+64]

проблема

Как можно видеть, X-положение каждой планеты становится очень большим для каждой итерации, настолько большим, что на самом деле, если программа выполняется слишком много итераций, числа становятся настолько большими, что их невозможно представить больше. Очевидно, есть проблема, которую я не смог выявить.

Я чувствую, что либо не вижу простой синтаксической ошибки, либо неверна семантика построения моих списков.

Для ясности: для каждого окончательного списка (X, Y, Z, VX, VY, VZ) первый индекс идентифицирует итерацию, а второй индекс идентифицирует планету.

Таким образом, X[1][1] даст X-положение Меркурия на второй итерации.

Пожалуйста, помогите мне определить, что я делаю неправильно здесь.

Answer 1

Давайте посмотрим на расчет для одного компонента скорости:

vxi=0 ...
for j in range(len(Planets)):
    if i!=j:
        vxi+=VX[n][i]+...

Вы добавляете начальную скорость один раз для каждой другой планеты, а не один раз. Это заставляет скорость каждой планеты расти экспоненциально со временем.