2

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.

enter image description here

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

Community
  • 1
tsuresuregusa
  • 2,661
  • 2
  • 13
  • 31

2 Answers2

6

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]

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 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. – tsuresuregusa Dec 12 '16 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. – Michael E2 Dec 12 '16 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. – Michael E2 Dec 12 '16 at 23:30
5

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]

enter image description here

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]
Community
  • 1
george2079
  • 38,913
  • 1
  • 43
  • 110