2
$\begingroup$

This is an N-Body simulation for the Sun and the following planetary bodies: Mercury, Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune and Pluto.

Initial Parameters

Ecc = {0.20563069, 0.00677323, 0.01671022, 0.09341233,0.04839266, 
0.05415060, 0.04716771, 0.00858587, 0.24880766};(*eccentricity of bodies*)

a = {0.38709893, 0.7233319899999999, 1.00000011, 1.52366231, 
5.2033630099999995, 9.537070319999998, 19.191263929999998, 
30.06896348, 39.48168677}; (*semi major axis of bodies*)

r = a (1 - Ecc^2)/(1 + Ecc Cos[\[Psi]]); (*orbital position*)
rx={0, 0.3075, 0.718433, 0.98329, 1.38133, 4.95156, 9.02063, 18.2861, 29.8108, 29.6583} (*x component of position*)
ry={0, 0., 0., 0., 0., 0., 0., 0., 0., 0.}
v = {0, 0.03406085426835039`, 0.020363076269733636`, 
0.017491554631468408`, 0.015304294697344465`, 
0.007915195286690359`, 0.005880353628887382`, 
0.004116410730170449`, 0.0031640275881454545`, 
0.0035297581940090896`};(*initial velocity*)
T = {0, 88.0, 224.7, 365.2, 687.0, 4331, 10747, 30589, 59800, 90560};(* period of orbit in days*)

Solving the differential equations

eq = {Table[
x[i]''[t] == 
 Sum[If[j == i, 
   0, (-\[Mu][[j]] (x[i][t] - 
        x[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - 
          y[j][t])^2)^(3/2)], {j, 10}], {i, 10}], 
 Table[y[i]''[t] == 
 Sum[If[j == i, 
   0, (-\[Mu][[j]] (y[i][t] - 
        y[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - 
          y[j][t])^2)^(3/2)], {j, 10}], {i, 10}]};

var = Join[Table[x[i], {i, 10}], Table[y[i], {i, 10}]];

orb = NDSolve[{eq, Table[x[i][0] == rx[[i]], {i, 10}], 
Table[y[i][0] == 0, {i, 10}], Table[x[i]'[0] == 0, {i, 10}], 
Table[y[i]'[0] == v[[i]], {i, 10}]}, var, {t, 0, 90600}];

Plotting the bodies

plot2D = Show[
Table[ParametricPlot[
 Evaluate[{x[i][t], y[i][t]} /. orb], {t, 0, 30000}(*,
 PlotStyle\[Rule]None*), PlotRange -> 5], {i, 10, 10}]];

Animating the bodies

Animate[Show[plot2D, 
Graphics[Table[{Red, PointSize[0.02], 
 Point[{x[i][t], y[i][t]} /. orb]}, {i, 1, 10}]]], {t, 0, 90000}, 
 AnimationRate -> 50, AnimationRunning -> False]

The Problem The simulation appears to work however, at around t=30000, the sun begins to drift which drags the rest of the bodies with it. Please see below for(first image) t=1000 and (second image)t=30000

As a result of this drift, pluto does not reach its aphelion of 49.3 AU.

Im aware that with N-body simulations, integration errors occur over time but could be an error in the code that might be causing this??

t=1000 t=30000

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ You did not set the initial parameters of the movement of the Sun. You must either set the movement of the Sun relative to the barycenter, or exclude the Sun from the simulation. $\endgroup$
    – Alex Trounev
    Oct 23, 2019 at 16:55
  • $\begingroup$ Assuming that the barycenter for the sun can be calculated in terms of the jupiter. Its position is then approximately 0.00496246 AU from coordinate origin. From this, is it then possible to calculate its respective velocity components for its orbit around the barycenter? $\endgroup$
    – isaac5122
    Oct 24, 2019 at 9:44
  • $\begingroup$ You can’t figure it out roughly. The barycenter is associated with ICRS, and ICRS is the current standard celestial reference system adopted by the International Astronomical Union (IAU). $\endgroup$
    – Alex Trounev
    Oct 24, 2019 at 11:34
  • $\begingroup$ I think its possible to approximate the barycenter if you sum: (mass of body 1*distance of body 1 from center of sun+mass of body 2*distanceof body 2 from center of sun)/(mass 1 + mass 2) for the Sun jupiter barycenter this is 0.00496246 AU. if we do this for all major planets(and pluto) then the barycenter is 0.00963163 AU $\endgroup$
    – isaac5122
    Oct 24, 2019 at 14:00

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.