3
$\begingroup$

I'm trying to solve a 3 body problem consisting of the earth moon and sun. I'm wanting to show the earth orbiting the sun and the moon orbiting the earth as a result of the gravitational forces.

Where Force due to gravity is given by:

$$\mathbf F_{ij}=\frac{Gm_i m_j}{\|\mathbf q_j-\mathbf q_i\|^2}\cdot\frac{(\mathbf q_j-\mathbf q_i)}{\|\mathbf q_j-\mathbf q_i\|}=\frac{Gm_i m_j(\mathbf q_j-\mathbf q_i)}{\|\mathbf q_j-\mathbf q_i\|^3}$$

And the acceleration of the bodies can be given as

$$\vec{a}_j=\sum_{i\ne j}^n G\frac{M_i}{\left|\vec{r}_i-\vec{r}_j\right|^3}\left(\vec{r}_i-\vec{r}_j\right)$$

INITIAL CONDITIONS

The initial conditions of the bodies are as follows:

Ecc = {0.01671022, 0.00549}(*eccentricity of earth and moon orbits*);
Inc = {0.00005, 5.145};
RAAN = {-11.260640, 0};
ArgPer = {102.94719, 0};
G = 0.00029589743849552926`;(*gravitational constant in AU*)
m = {1, 3.004*10^-6, 
  3.694*10^-8}; (*mass of sun, earth and moon in solar    masses*)
μ = Table[
G*m[[i]], {i, 1, 3}];(*standard gravitational parameter of *)
(*μ={8.8878*10^-10, 1.093*10^-11};*)(*standard gravitational \
constant\[Rule] Sunm, earth, moon*)    
ψ = 0;
a = {1.00000011, 0.99743};(*semi major axis in AU*);
r = a (1 - Ecc^2)/(1 + 
Ecc Cos[ψ])(* disatances for sun to earth and sun to moon 
in AU*)
rx = r Cos[ψ](* x component of distance for: sun to earth and \
sun to moon*);
ry = r Sin[ψ](* y component of distance forsun to earth and sun \
to moon*);
v = Table[
Sqrt[μ[[i]] (2/r[[i]] - 1/a[[i]])], {i, 
2}](*earth velocity due to sun , moon velocity due to earth units of AU/day*)
T = {365.2, 27}(*earth moon period in days*);

EQUATIONS OF MOTION

Nbody = Table[
  NDSolve[{(x''[t] + (
  G*m[[i]]*x[t])/((x[t])^2 + (y[t])^2)^(3/2) + (
  G*m[[j]]* (x[t]))/((x[t])^2 + (y[t])^2)^(3/2) == 
 0, (y^''[
  t] == -((G*m[[i]] (y[t]))/((x[t])^2 + (y[t])^2)^(3/2)) - (
  G*m[[j]] (y[t]))/((x[t])^2 + (y[t])^2)^(3/2), x[0] == rx[[i]], 
y[0] == ry[[i]], [x]'[0] == 0, 
[y]'[0] == v[[i]]}, {x, y}, {t, 0, 365.2}], {i, 
2}, {j, 2}]

I can plot the results using `

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

However I get this graph Along with the error message

ReplaceAll::reps: {{{{x->InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}],y->InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}]}}}[[2]]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

Coordinate ReplaceAll[{$CellContext`x[0], $CellContext`y[0]}, Part[{{{$CellContext`x -> InterpolatingFunction[{{0., 365.2}}, {5, 7, 2, {84}, {4}, 0, 0, 0, 0, Automatic, {}, {}, False}, {{0., 0.011199207293206922`, 0.022398414586413843`, 0.27915745480600035`, 0.5359 should be a pair of numbers, or a Scaled or Offset form.

What I think is wrong

From the graph, the orbit of the earth is correct. The orbit of the moon should be around the earth which obviously it's not!

From section of Code for Nbody;

The masses for the sun and the moon are correct However I think the respective vectors for each body are not.

What I think I need to do is find a way for apply the initial conditions to both bodies individually, which I don't think I've done using this code.

Would anyone be able to help in correcting this?

Thank you for your time.

plot with error

Improve this question
$\endgroup$
8
  • $\begingroup$ The variable planet is used in the Graphics[ ... ] expression, but planet is not defined in the code above. $\endgroup$
    – LouisB
    Oct 22, 2019 at 19:00
  • $\begingroup$ Thank you for your reply, i've corrected the error however, the original problem still remains. $\endgroup$
    – isaac5122
    Oct 22, 2019 at 19:12
  • $\begingroup$ (x^\[Prime]\[Prime])[t] is wrong, output only looks like Derivative[2][x][t]! That's why NBodyisn't evaluated! $\endgroup$
    – Ulrich Neumann
    Oct 22, 2019 at 19:35
  • $\begingroup$ Even with x''[t], y''[t], x'[0], y'[0], the problem remains, thank you for your reply, and your patience! $\endgroup$
    – isaac5122
    Oct 22, 2019 at 19:45
  • 4
    $\begingroup$ If you have access to M12, why don't you try using NBodySimulation instead? $\endgroup$
    – Carl Woll
    Oct 22, 2019 at 19:53

1 Answer 1

Reset to default
9
$\begingroup$

Solving the problem requires high accuracy in determining the initial data. All parameters can be determined using StarData["Sun", "..."], PlanetData["Earth", "..."], \ PlanetaryMoonData["Moon", "..."],For example, mass

m = {StarData["Sun", "Mass"], PlanetData["Earth", "Mass"], 
  PlanetaryMoonData["Moon", "Mass"]}

(*Out[]= {Quantity[1.988435*10^30, "Kilograms"], 
 Quantity[5.972*10^24, "Kilograms"], 
 Quantity[7.3459*10^22, "Kilograms"]}*)

Gravitational Constant Mass Product

 mG = {StarData["Sun", "GravitationalConstantMassProduct"], 
  PlanetData["Earth", "GravitationalConstantMassProduct"], 
  PlanetaryMoonData["Moon", "GravitationalConstantMassProduct"]};

 mG = mG // QuantityMagnitude

(*Out[]= {1.327124400*10^20, 3.98600436*10^14, 4.9028*10^12}*}

So I defined all the parameters

μ = {0.00029591220828349645438389263285760599`10.122911634152604, 
  8.8876924487802410900708441245`8.697447781980085*^-10, 
  1.0931894685338682`*^-11}; (*standard gravitational parameter of \
sun,earth and moon in AU*)
rx = {0, 0.9832672274945706, 
  0.9847793657772581}(*x and y components of initial position*);
ry = {0, 0, 0.0020731551812556903};
vy = {0, 0.017491119160350586, 0.017965913470733775}; vx = {0, 
  0, -0.0003463101361750346`}(*initial velocities of earth and moon*);
eq = {Table[
    x[i]''[t] == 
     Sum[If[j == i, 
       0, (-μ[[j]] (x[i][t] - 
          x[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - y[j][t])^2)^(
       3/2)], {j, 3}], {i, 3}], 
   Table[y[i]''[t] == 
     Sum[If[j == i, 
       0, (-μ[[j]] (y[i][t] - 
          y[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - y[j][t])^2)^(
       3/2)], {j, 3}], {i, 3}]};
var = Join[Table[x[i], {i, 3}], Table[y[i], {i, 3}]];

orb = NDSolve[{eq, 
   Table[x[i][0] == rx[[i]], {i, 3}], 
   Table[y[i][0] == ry[[i]], {i, 3}], 
   Table[x[i]'[0] == vx[[i]], {i, 3}], 
   Table[y[i]'[0] == vy[[i]], {i, 3}]}, var, {t, 0, 366}, 
  MaxStepSize -> 10^-3]

The orbit of the earth and moon

{ParametricPlot[
  Evaluate[{{var[[2]][t], var[[5]][t]}, {var[[3]][t], 
      var[[6]][t]}} /. First[orb]], {t, 0, 365}], 
 ParametricPlot[
  Evaluate[{var[[2]][t] - var[[3]][t], var[[5]][t] - var[[6]][t]} /. 
    First[orb]], {t, 0, 27.322}]}

There is still insufficient accuracy for the orbit of the moon Figure 1

Some details of calculating input parameters.

1.Calculate the solstice date 

PlanetData["Earth", "PeriapsisTime"]

(*Out[]= DateObject[{2020, 1, 5}, "Day", "Gregorian", -4.]*)

  1. Calculate the coordinates and speed of the earth and moon

    PlanetData["Earth", 
 EntityProperty["Planet", 
  "VelocityAroundSun", {"Date" -> DateObject[DateList[{2020, 1, 8}]]}]]
PlanetData["Earth", 
 EntityProperty["Planet", 
  "DistanceFromSun", {"Date" -> DateObject[DateList[{2020, 1, 8}]]}]]
(*Quantity[30.28511785356967`, ("Kilometers")/("Seconds")]
Quantity[0.9832672274945706`, "AstronomicalUnit"]*)


PlanetaryMoonData["Moon", 
 EntityProperty["PlanetaryMoon", 
  "DistanceFromSun", {"Date" -> DateObject[DateList[{2020, 1, 8}]]}]]
 (*Quantity[0.985403, "AstronomicalUnit"]*)
PlanetaryMoonData["Moon", 
 EntityProperty["PlanetaryMoon", 
  "SunElongation", {"Date" -> DateObject[DateList[{2020, 1, 8}]]}]]

(*Quantity[MixedMagnitude[{143, 49, 27.613}], 
 MixedUnit[{"AngularDegrees", "ArcMinutes", "ArcSeconds"}]]*)

 PlanetaryMoonData["Moon", "AverageOrbitVelocity"]

 (*Quantity[1.02, ("Kilometers")/("Seconds")]*)
 PlanetaryMoonData["Moon", 
 EntityProperty["PlanetaryMoon", 
  "DistanceFromEarth", {"Date" -> DateObject[DateList[{2020, 1, 8}]]}]]

(*Quantity[383873., "Kilometers"]*)

    We use geometry to calculate the coordinates and speed of the moon.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you for your answer Alex. May i ask how you derived the expressions for vx and vy for the moon? Did you use the expression for orbital velocity that I used? $\endgroup$
    – isaac5122
    Oct 23, 2019 at 6:15
  • $\begingroup$ @Luke4737 See update to my answer. $\endgroup$
    – Alex Trounev
    Oct 23, 2019 at 12:54
  • $\begingroup$ Thank you very much for your insight alex, ive been able to make an N-body simulation for the solar system. Theres an interesting quirk where all planets begin to drift with the sun after an approximate period of 30000 days. As a result of this drift, pluto never reaches its aphelion of 49.3 AU yet still completes a full orbit. If i can discover the error behind this i believe it'll work. would you like to see? $\endgroup$
    – isaac5122
    Oct 23, 2019 at 13:38
  • 1
    $\begingroup$ @Luke4737 Open a new theme for a simulation of 30,000 days. $\endgroup$
    – Alex Trounev
    Oct 23, 2019 at 13:44
  • 1
    $\begingroup$ I think a symplectic method for NDSolve will do better to keep energy conservation. A built-in reference: "SymplecticPartitionedRungeKutta" Method for NDSolve $\endgroup$
    – Silvia
    Nov 30, 2019 at 5:56

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.

Not the answer you're looking for? Browse other questions tagged or ask your own question.