3

#Forward Euler N-Body Simulation 该程序从星历中获取太阳系行星在特定日期的位置和速度,然后使用前向欧拉积分方法模拟行星的运动。

这种模拟背后的基本物理学是牛顿万有引力定律。

对于积分的每次迭代,每个 jth 行星对第 i 个行星(其中 j !=i)的影响被累积计算。

# 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

#X 中 n 的结果:打印 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 正如我们所看到的,每颗行星的 X 位置在每次迭代中都变得非常大,实际上如此之大,以至于如果程序运行太多次迭代,数字就会变得如此之大,以至于无法再表示它们。显然有一个问题,我无法确定。

我觉得我要么没有看到简单的语法错误,要么我的列表构造的语义不正确。

需要明确的是:对于每个最终列表(X、Y、Z、VX、VY、VZ),第一个索引标识迭代,第二个索引标识行星。

所以 X[1][1] 将在第二次迭代中给出 Mercury 的 X 位置。

请帮助我确定我在这里做错了什么。

4

1 回答 1

2

让我们看一下速度的一个分量的计算:

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

您为每个其他行星添加一次初始速度,而不仅仅是一次。这导致每个行星的速度随着时间呈指数增长。

于 2015-06-25T01:10:20.710 回答