# Mass distribution in NDSolve

I'm trying to model the solar system with NDSolve, basically a numerical n-body solution.

I model the initial positions of bodies from a data set thusly:

{"sun",
{-1.068000648301820*^9, -4.176802125684930*^8, 3.084467020687090*^7},
{9.305300847631915, -1.283176670344807*10,-1.631528028381386*^-1},
1.988544*^30, White}


Where the first triplet is the body's initial position with regard to the solar system's barycentre, the second triplet its initial velocity, and the fourth position number its mass.

The NDSolve, then evaluates the law of gravitation with respect to all tuples of masses in Cartesian component form.

This works fine for the major bodies. I'm currently getting roughly 3-4 accurate significant numbers after a setting the time of the model to 1 year.

My problem is with modelling the asteroid belt. Naturally, I can't model it as a single point particle. So I was wondering if there is a way to model it as a distribution with mass $3\times10^{21}$. I'm really only interested in finding a way to have the asteroid belt excert its gravitational influence on other bodies.

Here is a minimal working example, if I understand how to do it here, I can work out how to incorporate it in the actual model.

Supposing I started with a very simple model of one planet (sun+earth):

G = 6.6740831*10^-11;
{m1,m2} = {1.98855*^30, 5.97237*^24};
a = {x1''[
t] == (G m2) (-x1[t] +
x2[t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2),
x2''[
t] == (G m1) (-x2[t] +
x1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2),
y1''[
t] == (G m2) (-y1[t] +
y2[t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2),
y2''[
t] == (G m1) (-y2[t] +
y1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2),
z1''[
t] == (G m2) (-z1[t] +
z2[t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2),
z2''[
t] == (G m1) (-z2[t] +
z1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2)};
vel = {{Derivative[x1] == 0,
Derivative[x2] == 0}, {Derivative[y1] == -0.0894406,
Derivative[y2] == 29780}, {Derivative[z1] == 0,
Derivative[z2] == 0}};
pos = {{x1 == 0, x2 == 152100000000}, {y1 == 0,
y2 == 0}, {z1 == 0, z2 == 0}};
var = {x1, x2, y1, y2, z1, z2};
tmax = 3600*24*365.24;
s = NDSolve[Flatten[{a, vel, pos}], var, {t, 0, tmax}];
ParametricPlot3D[
Flatten[Thread[{{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}} /.s]],
{t, 0, tmax},
PlotRange -> {{-2*^11, 2*^11}, {-2*^11, 2*^11}, {-0.5*^11, 0.5*^11}}] This doesn't make a complete orbit in one year, but obviously the accuracy here is low. Where a is the acceleration vectors ,vel the velocities, and pos the positions.

• While it doesn't answer your mass question, you really should be in vector form here: X1[t_] = {x1[t], y1[t], z1[t]};X2[t_] = {x2[t], y2[t], z2[t]}; eqns = {X1''[t] == (g m2 (X2[t] - X1[t]))/Norm[(X2[t] - X1[t])]^3, X2''[t] == (g m1 (X1[t] - X2[t]))/Norm[(X2[t] - X1[t])]^3}; s = NDSolve[{eqns, vel, pos}, {X1[t], X2[t]}, {t, 0, tmax}]; – KraZug Jul 19 '16 at 14:47
• You could also label the vectors with an index, e.g. X[t_] or X["sun"][t_] etc, and that lets you iterate over those. – KraZug Jul 19 '16 at 14:49
• In a computation giving only 3-4 significant digits over a year, I strongly suspect the asteroid belt objects have no detectable effect on the major bodies of the solar system. – m_goldberg Jul 19 '16 at 14:52
• I would recommend working in something like N-body units. Keeping around all those factors of ten can be a problem, precision-wise, I think. – march Jul 19 '16 at 15:22
• Well, scaling quantities in a reasonable set of units is a good skill to have in any situation, really. You want as many of the parameters to be on the order of 1 as possible. – march Jul 19 '16 at 16:30

Because the Asteroid belt's mass is negligible (only 4% of the Moon's mass), it is not very efficient to include it in your simulations. However for doing so I have two methods for approximation:

• Only consider the top massive objects in asteroid belt (i.e Ceres, Vesta, etc), which include more than 50 percent of the belt's total mass
• Use the gravitational field of a thin torus to consider the effects of the Asteroid belt. you can use different approximation formulas for planets outside and inside the asteroid belt, which can be found here

UPDATE

I've tried to use an approximation for the gravitational field of a torus shaped asteroid belt for the Sun and Earth's orbit. If I have implemented the code correctly, then it shows how little can be the effect of asteroid belt on solar system. I have estimated the belt's Radius to be about 2.5 AU.

G = 6.6740831*10^-11;
AU =2.5*152100000000;
F[r_, R_] := (Block[{x = r/R}, (3.14159 x^2 +
3.5342917352885173 x^4 + 3.6815538909255383 x^6)/(2 Pi R^2)]);
{m1, m2, m3} = {1.98855*^30, 5.97237*^24, 3.*10^21};
a = {x1''[t] == (G m2) (-x1[t] +
x2[t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}], AU] (
Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t], z3[t]}][]),
x2''[t] == (G m1) (-x2[t] +
x1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}],
AU] (Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t],
z3[t]}][]),
x3''[t] == -G m1 F[
EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}],
AU] Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t],
z3[t]}][] -
G m2 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}],
AU] (Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t],
z3[t]}][])
, y1''[t] == (G m2) (-y1[t] +
y2[t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}], AU] (
Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t], z3[t]}][]),
y2''[t] == (G m1) (-y2[t] +
y1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}],
AU] (Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t],
z3[t]}][]),
y3''[t] == -G m1 F[
EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}],
AU] Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t],
z3[t]}][] -
G m2 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}], AU] (
Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t], z3[t]}][]),
z1''[t] == (G m2) (-z1[t] +
z2[
t])/((-x1[t] + x2[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}], AU] (
Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t], z3[t]}][]),
z2''[t] == (G m1) (-z2[t] +
z1[t])/((-x2[t] + x1[t])^2 + (-y1[t] + y2[t])^2 + (-z1[t] +
z2[t])^2)^(3/2) +
G m3 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}],
AU] (Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t],
z3[t]}][])
, z3''[t] == -G m1 F[
EuclideanDistance[{x1[t], y1[t], z1[t]}, {x3[t], y3[t],
z3[t]}],
AU] (Normalize[{x1[t], y1[t], z1[t]} - {x3[t], y3[t],
z3[t]}][]) -
G m2 F[EuclideanDistance[{x2[t], y2[t], z2[t]}, {x3[t], y3[t],
z3[t]}], AU] (
Normalize[{x2[t], y2[t], z2[t]} - {x3[t], y3[t], z3[t]}][])};
vel = {{Derivative[x1] == 0, Derivative[x2] == 0,
Derivative[x3] == 0}, {Derivative[y1] == -0.0894406,
Derivative[y2] == 29780,
Derivative[y3] == 0}, {Derivative[z1] == 0,
Derivative[z2] == 0, Derivative[z3] == 0}};
pos = {{x1 == 0, x2 == 152100000000, x3 == 0.1}, {y1 == 0,
y2 == 0, y3 == 0}, {z1 == 0, z2 == 0, z3 == 0}};
var = {x1, x2, x3, y1, y2, y3, z1, z2, z3};
tmax = 3600*24*365.24;
s = NDSolve[Flatten[{a, vel, pos}], var, {t, 0, tmax},
MaxStepFraction -> 1/1000000];
ParametricPlot3D[
Flatten[Thread[{{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}, {x3[t],
y3[t], z3[t]}} /. s]], {t, 0, tmax},
PlotRange -> {{-2*^11, 2*^11}, {-2*^11, 2*^11}, {-0.5*^11, 0.5*^11}}] Note: Depending on whether the objects are inside or outside the torus, you should use different approximation for the torus' gravitational field, which can be found in the link mentioned above.

In:= Simplify[
Series[1 x Integrate[(Cos[t] - x)/(1 - 2 x Cos[t] + (x)^2)^1.5, {t,
0, 2 Pi}], {x, 0, 6}], x < 1]

Out= SeriesData[x, 0, {3.1415926535897927, 0., \
3.5342917352885173, 0., 3.6815538909255383}, 2, 7, 1]

In:= Simplify[
Series[1/x Integrate[(Cos[t] - 1/x)/(1 -
2/x Cos[t] + (1/x)^2)^1.5, {t, 0, 2 Pi}], {x, 0, 6}], 0 < x < 1]

Out= SeriesData[x, 0, {-6.283185307179585, 0., \
-4.712388980384691, 0., -4.417864669110656}, 1, 7, 1]

In:= Plot[{3.1415926535897927 x^2 + 3.5342917352885173 x^4 +
3.6815538909255383 x^6,
x NIntegrate[(Cos[t] - x)/(1 - 2 x Cos[t] + x^2)^1.5, {t, 0,
2 Pi}], -6.283185307179585 /x - 4.712388980384691 (1/x)^3 -
4.417864669110656 (1/x)^5}, {x, 0, 2}] • First of all, I certainly appreciate the effort gone into this. If I make a VectorPlot3D[] of the acceleration for {fx,fy,fz} = G m3 F[EuclideanDistance[{x1, y1, z1}, {x3, y3, z3}], AU] (Normalize[{x1, y1, z1} - {x3, y3, z3}][[{1,2,3}]]) /. {x3 -> 0, y3 -> 0, z3 -> 0}; VectorPlot3D[{fx, fy, fz}, {x1, -a AU, a AU}, {y1, -a AU, a AU}, {z1, -a/2 AU, a/ AU}, BoxRatios -> {1, 1, 0.5}] For any a` gives the same result, of radiating outward horizontally, which isn't an expected result. – Feyre Jul 20 '16 at 10:58
• Please Notice that the Approximation that I've made was for an object inside the asteroid belt. As you can see you should use another form of approximation for the objects outside the torus: – Siav Josep Jul 20 '16 at 12:17
• Alright, I can't verify now, will accept once I have worked this through. – Feyre Jul 20 '16 at 12:59