gravitational two body problem for the orbit of the sun and jupiter around their barycenter

Question

This is a simulation of the sun and jupiter orbiting the respective barycenter.

Initial conditions

m = {1, 0.0009546133303706552`};(*masses of sun and jupiter in solar masses*)
G = 0.00029589743849552926`;(*gravitational constant in relevant units*)
\[Mu] = G*m;(*standard gravitational parameters of sun and jupiter*)
rx = {-0.004962462459288476`, 4.951558433000493`};(*Initial position from barycenter at(0,0)*)
v = {-7.203*10^-6, 0.007915195286690359`};(*relative velocity*)
T = {4331, 4331};(*period*)

The center of mass was calculated by:

Where r_s=0 and r_j=4.951558433000493 AU and m_s=1 and m_j=0.0009546133303706552 Solar masses

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, 2}], {i, 2}], 
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, 2}], {i, 2}]};

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

Plotting the orbits

plot2D = Show[
Table[ParametricPlot[
 Evaluate[{x[i][t], y[i][t]} /. orb], {t, 0, 
  90000},(*PlotStyle\[Rule]None,*)PlotRange -> 6], {i, 2}]];
Animate[Show[plot2D, 
Graphics[Table[{Red, PointSize[0.02], 
 Point[{x[i][t], y[i][t]} /. orb]}, {i, 2}]]], {t, 30000, 1}, 
AnimationRate -> 50, AnimationRunning -> False]

The problem

Upon initial inspection, the orbits seem to stable around the barycenter for up to the given period.

After this, the bodies begin to drift upwards

This can be seen by setting the values of plot2d to the following.

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

This is the sun orbiting the barycenter during the given period

This is the orbiting the barycenter but drifting upwards after the given period.

What i think is wrong

• When i calculated the velocity of the sun i assumed that the sun would share the same orbtial period of jupiter which may be wrong
• When i calculated the barycenter, i assumed that the sun would be displaced in the -x direction.
• I may have calculated the intial positons wrong
• I was following this tutorial for initial positon barycenters: https://www.youtube.com/watch?v=4cv8IeeBMtc
• and this tutorial for calcualting the velocity: https://www.youtube.com/watch?v=Lp4u2L8HNPI

Why is the sun drifiting instead of orbting its barycenter after one period of 4331 days? Have i made an error calculating the suns velocity that would cause this?

What im trying to achieve is a barycentric orbit like in the image below

Answer 1

We cannot depict the orbits of the Sun and Jupiter in the same figure on the same scale, since the radius of the orbit of the Sun is about 0.001 the radius of the orbit of Jupiter. But we can show in one animation their synchronous movement around the barycenter.

m = {1, 0.0009546133303706552`};(*masses of sun and jupiter in solar \
masses*)G = 0.00029589743849552926`;(*gravitational constant in \
relevant units*)\[Mu] = 
 G*m;(*standard gravitational parameters of sun and jupiter*)rx = \
{-0.004962462459288476`, 
  4.951558433000493`};(*Initial position*)v = {-7.203*10^-6, 
  0.007915195286690359`};(*relative velocity*)T = {4331, 
  4331};(*period*)
dv = m.v/(m.{1, 1});
v = v - {dv, dv};
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, 2}], {i, 2}], 
   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, 2}], {i, 2}]};

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

The orbit of Jupiter and the Sun around a stationary barycenter (I reset the speed of the barycenter to 0)

plot2D = Table[
  ParametricPlot[Evaluate[{x[i][t], y[i][t]} /. orb], {t, 0, 30000}, 
   PlotRange -> All], {i, {2, 1}}]

Animation

StarData["Sun", "Diameter"]

ae = 149597870700; ds = 1.3914`4.*^6*10^3/ae

so = ParametricPlot[
  Evaluate[{x[1][t], y[1][t]} /. orb], {t, 0, 30000}, 
  PlotRange -> {{-.01, .01}, {-.01, .01}}]
frame = Table[{Show[so, 
     Graphics[{Red, Disk[First[{x[1][t], y[1][t]} /. orb], ds/2]}]], 
    Show[plot2D[[1]], 
     Graphics[{Red, PointSize[0.05], 
       Point[{x[2][t], y[2][t]} /. orb]}]]}, {t, 0, 4331, 43.31}];

ListAnimate[frame]