#Forward Euler N-Body Simulation This program takes the positions and velocities of the solar system's planets at a specific date from an ephemeris then simulates the planets' movements using a forward euler integration method.
The fundamental physics behind this simulation is Newton's Law of Universal Gravitation.
For each iteration of the integration, the effects of each jth planet on the ith planet (where j !=i) are cumulatively accounted for.
# 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
#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 one can see, each planet's X-position is becoming very larger for every iteration, so large in fact that if the program is run for too many iterations the numbers become so large that they cannot be represented anymore. Obviously there is a problem, which I have not been able to identify.
I feel I am either not seeing a simple syntax error or that the semantics of the construction of my lists are incorrect.
To be clear: for each final list (X, Y, Z, VX, VY, VZ), the first index identifies the iteration, and the second index identifies the planet.
So X[1][1] would give the X-position of Mercury at the second iteration.
Please help me identify what I am doing incorrectly here.
Let's look at the calculation for one component of velocity:
You're adding in the initial velocity once for every other planet, rather than just once. This causes each planet's velocity to grow exponentially with time.