# How to WhenEvent for indexed variables

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?

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]


• thanks a lot! I was stuck with this for ages... where can I find more info in this "standard" problem? I read the documentation of WhenEvent once and again and could not figure out how to fix this. Commented Dec 12, 2016 at 17:11
• @tsuresuregusa It's in the docs for With, under "Scope," but I remember it from "A Good Trick to Know" in Introduction to Dynamic. Commented Dec 12, 2016 at 20:28
• @tsuresuregusa You can also find it used without explanation in the docs for WhenEvent under "Applications" > "BouncingBall", the 20 bouncing balls example. Commented Dec 12, 2016 at 23:30

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]


...

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]

• thanks a lot! this is an awesome use of vectors, very neat form of writing vforce Commented Dec 13, 2016 at 3:19