# N body simulation

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[{$$CellContextx[0],$$CellContexty[0]}, Part[{{{\$CellContextx -> 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.

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

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.0002959122082834964543838926328576059910.122911634152604,
8.88769244878024109007084412458.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

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.

• 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? Oct 23 '19 at 6:15
• @Luke4737 See update to my answer. Oct 23 '19 at 12:54
• 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? Oct 23 '19 at 13:38
• @Luke4737 Open a new theme for a simulation of 30,000 days. Oct 23 '19 at 13:44
• I think a symplectic method for NDSolve will do better to keep energy conservation. A built-in reference: "SymplecticPartitionedRungeKutta" Method for NDSolve Nov 30 '19 at 5:56