Simulation of N-Body movements of direct euler

This program takes the positions and velocities of the planets in the solar system on a specific day from the ephemeris and then simulates the planetary motions using the direct Euler integration method.

The fundamental physics behind this simulation is Newton's law of universal gravity.

For each iteration of integration, the effects of each j-th planet on the i-th planet (where j! = I) are cumulatively taken into account.

#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 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 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 X-position in the second time step

      

        
        
        
      

    

results

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]

      

        
        
        
      

    

Problem

As you can see, each X-position of the planet becomes very large for each iteration, so large in fact that if the program is run for too many iterations, the numbers become so large that they can no longer be represented. Obviously there is a problem that I have not been able to identify.

I feel like I am either not seeing a simple syntax error, or that the semantics of my lists are wrong.

To be clear: for each final list (X, Y, Z, VX, VY, VZ), the first index identifies the iteration and the second index indicates the planet.

So X [1] [1] will give the X-position of Mercury in the second iteration.

Please help me determine what I am doing wrong here.