# Detecting a collision in n-body simulation with NDSolve

Consider this three-body system:

g = 10; (* gravitational constant *)
m = {100, 1, 1}; (* masses *)
μ0 = g*m; (* gravitational parameters *)
r0 = {{0, 0, 0}, {1, 0, 0}, {-1.1, 0, 0}}; (* initial location vectors *)
v0 = {{0, 0, 0}, {0, 30, 0}, {0, -30.1, 0}}; (* initial velocity vectors *)


Differential equations for velocities and positions:

jmax = Length[m];
v = Evaluate[#] & /@
Table[ToExpression["v" <> ToString[j]], {j, jmax}];
r = Evaluate[#] & /@
Table[ToExpression["r" <> ToString[j]], {j, jmax}];
eqns[μ_List] := Flatten[Table[{
v[[j]]'[t] ==
Sum[(-Normalize[r[[j]][t] - r[[i]][t]]*μ[[i]])/
EuclideanDistance[r[[j]][t], r[[i]][t]]^2, {i,
Delete[Range[jmax], j]}],
r[[j]]'[t] == v[[j]][t],
v[[j]] == v0[[j]],
r[[j]] == r0[[j]]
}, {j, jmax}]]
fns = Flatten[Table[{v[[j]][t], r[[j]][t]}, {j, jmax}]]


The numerical integration stops when two bodies collide:

μ = μ0;
sol = NDSolve[eqns[μ], fns, {t, 0, 10}][]
vel = (fns /. sol)[[1 ;; -2 ;; 2]];
pos = vel = (fns /. sol)[[2 ;; -1 ;; 2]];


NDSolve::ndsz: At t == 2.393483348087306, step size is effectively zero; singularity or stiff system suspected. >>

The solution:

tmax = 2.393483348087306;
ParametricPlot3D[pos, {t, tmax - .1, tmax},
PlotRange -> {All, All, {-.01, .01}}, ImageSize -> 800,
PlotPoints -> 500, MaxRecursion -> 15] I would like to detect collisions and modify the computation so that the masses (and in effect, gravitational parameters) and velocities of the collided bodies change: one body's new mass becomes the sum of both masses and the new velocity becomes the sum of both velocities; the other body's new mass and velocity become zero.

I have tried to do it with WhenEvent:

μ = μ0;
sol = NDSolve[
Append[eqns[μ],
WhenEvent[Length[$MessageList] > 0, Print[tmax = t]; dist = Table[EuclideanDistance[r[[j + 1]][t], r[[j]][t]], {j, Length[r] - 1}]; distDif = Table[Abs[dist[[n + 1]] - dist[[n]]], {n, Length[dist] - 1}]; p = Position[distDif, Min[distDif]][[1, 1]]; μ[[p]] += μ[[p + 1]]; μ[[p + 1]] = 0; v[[p]][t] += v[[p + 1]][t]; v[[p + 1]][t] = 0; "RestartIntegration"]], fns, {t, 0, 10}][] vel = (fns /. sol)[[1 ;; -2 ;; 2]]; pos = vel = (fns /. sol)[[2 ;; -1 ;; 2]];  However, this does not work. It seems that I do not catch the error properly, the program aborts before the message can be caught. How to catch this properly? • I think you just can't use the message that NDSolve issues as an event trigger. Presumably the message is only generated after the integration has stopped for that reason. What I'd do is to check for the actual minimal distance between all bodies and if that is smaller than their extents that means collision. This will also be faster as it avoids the very small step sizes which lead to the message... – Albert Retey Dec 12 '13 at 11:41 ## 1 Answer As mentioned in my comment I think it makes more sense if you try to detect the collision by investigating the current positions than examining $Messages. Here is how this can be done by introducing additional discrete variables for the masses:

g = 10;
m0 = {100, 1, 1};
r0 = {{0, 0, 0}, {1, 0, 0}, {-1.1, 0, 0}};
v0 = {{0, 0, 0}, {0, 30, 0}, {0, -30.1, 0}};

numbodies = Length[m0];

odesys = Table[{
v[j]'[t] == UnitStep[m[j][t] - 0.01]*Sum[
(-Normalize[r[j][t]-r[i][t]]*g*m[i][t])/EuclideanDistance[r[j][t], r[i][t]]^2,
{i, Delete[Range[numbodies], j]}
],
r[j]'[t] == v[j][t],
v[j] == v0[[j]],
r[j] == r0[[j]]
},
{j, numbodies}
];

initialmasses = Table[m[k] == m0[[k]], {k, numbodies}];

whenevnt =
With[{rmin =
Min[EuclideanDistance @@@
Subsets[r[#][t] & /@ Range[numbodies], {2}]],
numbodies = numbodies},
WhenEvent[rmin < 0.0001,
pids =
SortBy[Subsets[Range[numbodies], {2}],
EuclideanDistance[r[#[]][t], r[#[]][t]] &][];
Print["collision at t=", t, " bodies:", pids];
{m[pids[]][t], v[pids[]][t], m[pids[]][t],
v[pids[]][t]} -> {
Total[
m[#][t] & /@ pids], (m[#][t] & /@ pids).(v[#][t] & /@ pids)/
Total[m[#][t] & /@ pids],
0, {0, 0, 0}
}
]
];

depvars = Flatten[Table[{v[j], r[j]}, {j, numbodies}]];

sol = NDSolve[
Flatten[{odesys, initialmasses, whenevnt}],
depvars, {t, 0, 5},
DiscreteVariables -> Table[m[k], {k, numbodies}]
][];

vel = Array[v, {3}] /. sol;
pos = Array[r, {3}] /. sol;
tmax = (Head[r[t] /. sol]@"Domain")[[1, 2]]

Animate[
Show[
ParametricPlot3D[Evaluate[#@tau & /@ pos], {tau, t - 0.1, t},
PlotRange -> {3*{-1, 1}, 3*{-1, 1}, {-.01, .01}},
ImageSize -> 800, PlotPoints -> 500, MaxRecursion -> 15,
Boxed -> False, Axes -> False,
Background -> Black
],
Graphics3D[{White, PointSize[Large],
Point[pos[[#]][t] & /@ Range[numbodies]]}]
],
{t, 0.1, tmax},
AnimationRepetitions -> 1, AnimationRate -> 0.5
]


Except for the mentioned change in the WhenEvent I have also made the following changes:

I am using expressions like r instead of symbols like r1 as variables. It is one of the great things about NDSolve that this works and it of course makes the programmatic generation of variables much easier.

I have also changed the equations slightly: by setting one of the colliding masses to zero you want to make that effectively not interact anymore. That will fail because your equations implicitly were divided by m, which of course goes wrong when m=0. I have corrected that in a somewhat ad hoc way by using a UnitStep[m-eps], where you need to adjust eps so that it is smaller than any mass you want to use. Without that term, the zero-mass body would still see the forces of the others and that will make NDSolve stop because of "effectively zero stepsize" as before.

Of course for the general case you should add momentum, not velocities which I have also changed from your original code. For your particular case (nonrelativistic with equal masses colliding) it doesn't make a difference, though.

For the collision detection I have also used an arbitrary "small" number which in reality proably has to do something with the extent of the bodies and should be adjusted accordingly...

I also understand your units and the value of g to be arbitrary, and expect you are aware that you should use the gravitational constant (usually denoted G), not the free fall acceleration (usually g) in those equations...

• +1 Also there is a possibility of three-body collisions. – Silvia Dec 12 '13 at 16:31
• Thank you for such a detailed answer, I appreciate all suggestions and improvements to my code. However, in this solution the Euclidean distance between all bodies is calculated for every computation step in NDSolve. I wanted to avoid it and only do the calculation when NDSolve throws an error. I'd like to know if this can be achieved, if not I will accept this answer. – shrx Dec 12 '13 at 17:11
• @shrx: AFAIK it isn't possible to see the message in the WhenEvent, but I may well be wrong. Even if it would be possible I think waiting for the message to appear wouldn't be a very efficient solution as it might take a lot of very small steps before NDSolve stops. What you could try is to check the step size yourself and react to a small step size in the WhenEvent`, probably only then checking the distances as well. Maybe you should add it as a requirement to your question that you explicitly want to avoid the calculation of the euclidian distance. – Albert Retey Dec 13 '13 at 9:59
• @Silvia: and for n bodies there would be a possibility for an m-body collision with 2<=m<=n. Of course the code I've shown doesn't cover that, but I think it would just mean to provide more sophisticated code for the event-action to cover such cases... – Albert Retey Dec 13 '13 at 10:02
• I don't need to deal with more than 2 body collisions for now, this is good enough. – shrx Dec 13 '13 at 13:50