2
$\begingroup$

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?

$\endgroup$

2 Answers 2

6
$\begingroup$

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]

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Dec 12, 2016 at 17:11
  • $\begingroup$ @tsuresuregusa It's in the docs for With, under "Scope," but I remember it from "A Good Trick to Know" in Introduction to Dynamic. $\endgroup$
    – Michael E2
    Dec 12, 2016 at 20:28
  • $\begingroup$ @tsuresuregusa You can also find it used without explanation in the docs for WhenEvent under "Applications" > "BouncingBall", the 20 bouncing balls example. $\endgroup$
    – Michael E2
    Dec 12, 2016 at 23:30
5
$\begingroup$

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]
$\endgroup$
1
  • $\begingroup$ thanks a lot! this is an awesome use of vectors, very neat form of writing vforce $\endgroup$ Dec 13, 2016 at 3:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.