How to WhenEvent for indexed variables

Question

I'm trying to solve a many bodies problem with NDSolve and I want to integrate just untill one particle has moved beyond a certain point, say the circle of radius 10 in 2D.

ClearAll[force, r, numbodies, pos0]
numbodies = 10;
pos0 = .5 RandomReal[{-10, 10}, {numbodies, 2}];

force[j_] := 
  Sum[(Normalize[r[j][t] - r[i][t]])/
     EuclideanDistance[r[j][t], r[i][t]]^2 +
    1/EuclideanDistance[r[j][t], r[i][t]]^12, {i, 
    Delete[Range[numbodies], j]}];

odesys = Table[{r[j]'[t] == force[j], r[j][0] == pos0[[j]](*,
    WhenEvent[Norm[r[j][t]]>10//Evaluate,"StopIntegration"]*)}, {j, numbodies}];


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

tfin = 550;

sol = NDSolve[odesys, depvars, {t, 0, tfin}, AccuracyGoal -> 2][[1]];

pos = Array[r, {numbodies}] /. sol;

Show[ListLinePlot[
  Transpose@Table[#[t] & /@ pos, {t, 0, tfin, tfin/100}], 
  Axes -> False, AspectRatio -> 1, Background -> Hue[.5], 
  PlotStyle -> Table[Hue[.5 + i/(2 numbodies)], {i, numbodies}], 
  ImageSize -> Large], Graphics[Circle[{0, 0}, 10]]]

The problem is when I uncomment the WhenEvent, and depending if I add an Evaluate to it, as in this answer it will give either do nothing or break the integration. I guess my problem is that Evaluate should act only on j and not on t but I don't know how to acomplish that.

So, what is the correct way to stop integrating my system when one particle crosses the circle?

Answer 1

The problem is a standard one and a standard solution is to use With[{j = j},...] to inject the value of j into the held argument:

odesys = Table[
   With[{j = j},
    {r[j]'[t] == force[j], r[j][0] == pos0[[j]], 
     WhenEvent[Norm[r[j][t]] > 10, "StopIntegration"]}
    ], {j, numbodies}];

Then executing the OP's code and resetting tfin, we get:

tfin = Part[r[1]["Domain"] /. sol, 1, -1]
(*  27.3852  *)

Show[ListLinePlot[
  Transpose@Table[Through[pos[t]], {t, 0., tfin, tfin/100}], 
  Axes -> False, AspectRatio -> 1, Background -> Hue[.5], 
  PlotStyle -> Table[Hue[.5 + i/(2 numbodies)], {i, numbodies}], 
  ImageSize -> Large], Graphics[Circle[{0, 0}, 10]], PlotRange -> All]

Answer 2

this is done more cleanly in vector form:

numbodies = 10;
pos0 = .5 RandomReal[{-10, 10}, {numbodies, 2}];
vforce[r : {{_?NumberQ, _?NumberQ} ..}] := 
  Table[Sum[(Normalize[r[[j]] - r[[i]]])/
      EuclideanDistance[r[[j]], r[[i]]]^2 +
     1/EuclideanDistance[r[[j]], r[[i]]]^12,
    {i, Delete[Range[numbodies], j]}], {j, numbodies}];
maxr[r : {{_?NumberQ, _?NumberQ} ..}] := Max[Norm /@ r]
vsoln = NDSolveValue[{r'[t] == vforce[r[t]], r[0] == pos0 , 
    WhenEvent[maxr[r[t]] > 10, "StopIntegration"]}, 
   r[t], {t, 0, Infinity}];
tmax = (vsoln /. t -> "Domain")[[1, 2]];
Show[{ParametricPlot[vsoln, {t, 0, tmax}, Axes -> False], 
  Graphics[Circle[{0, 0}, 10]]}, PlotRange -> All]

see: https://mathematica.stackexchange.com/a/78644/2079

...

a bit cleaner still define vforce as:

fi[x_List, x_List] = {0, 0}
fi[x_List, y_List] := Normalize[x - y] 
      (1/EuclideanDistance[y, x]^2 + 1/EuclideanDistance[y, x]^12)
vforce[r : {{_?NumberQ, _?NumberQ} ..}] := Total /@ Outer[ fi, r, r, 1]