9
$\begingroup$
n = 4;
countTime = 5;
SeedRandom[5];
initPos = RandomPoint[Disk[], n];
data = NBodySimulation[
   Association["PairwisePotential" -> "Coulomb", "Region" -> Disk[], 
    "ExternalForce" -> (Quantity[-0.5 Normalize[
          QuantityMagnitude[#["Velocity"]]](*An extra damping*), "Newtons"] &)], <|
      "Mass" -> Quantity[1, "Kilograms"], 
      "Position" -> Quantity[#, "Meters"], 
      "Velocity" -> Quantity[{0, 0}, "Meters"/"Seconds"], 
      "Charge" -> Quantity[10^-5, "Coulombs"]|> & /@ initPos, 
   countTime];
colors = RandomColor[n];
Manipulate[
 Graphics[{Circle[], Red, PointSize[0.02], 
   Riffle[colors, Point /@ data[All, "Position", time]]}, 
  Axes -> True], {time, $MachineEpsilon, countTime, 
  Appearance -> "Labeled"}]

enter image description here

enter image description here

I noticed that because one of the balls is coming to rest, the other balls can no longer be simulated because the iteration spacing is 0 after t>4.09486. How can I solve this problem so that the simulation continues until all the balls are at rest?

I guess we can change certain system options by Block? But I don't know which system option is causing this error...

$\endgroup$
4
  • $\begingroup$ You can control the method of integration for the underlying NDSolve by providing the option Method. Because the systems are Hamiltonian, the default choice in NBodySimulation is "SymplecticPartitionedRungeKutta". You can try with Method -> "StiffnessSwitching" which goes past the stiffness point, but soon after that, unfortunately, the balls wander outside the circle ... Internally, reflections from the walls are performed via WhenEvents. Perhaps there is already an answer here at MMA.SE on how to properly set up WhenEvents for such cases ... $\endgroup$
    – Domen
    Feb 24, 2022 at 19:50
  • $\begingroup$ @Domen I note the method of ExplicitMidpoint ,StiffnessSwitching and Extrapolation all will avoid error result the balls wander outside the circle. I don't know if this is a bug or not in NBodySimulation $\endgroup$
    – yode
    Feb 24, 2022 at 20:24
  • $\begingroup$ @yode What force do you try to describe with option "ExternalForce" -> (Quantity[-0.5 Normalize[ QuantityMagnitude[#["Velocity"]]](*An extra damping*), "Newtons"] &)? $\endgroup$ Feb 25, 2022 at 5:07
  • $\begingroup$ @AlexTrounev If I don't add damping to the charge, those charges will never stop. This force I am describing is like friction, always in the opposite direction to the speed and constant in magnitude. $\endgroup$
    – yode
    Feb 25, 2022 at 5:41

3 Answers 3

6
$\begingroup$

Solution 1

Funny, I'm not sure if this is the only solution, but Method -> "ExplicitEuler" (we know this is a rather primary method for ODE solving) solves the problem. This seems to be the first time I found a problem can be resolved with ExplicitEuler!

I've also added StartingStepSize -> 10^-3 to speed up the calculation a bit. (The default step choosing is around 10^-6, which turns out to be unnecessarily small. ) The calculation takes about 6 seconds on my laptop, tested in v12.3.1.

Result:

enter image description here

ListPlot[#, PlotRange -> All] & /@ data[All, "Position"]

Mathematica graphics


Solution 2

Not as straightforward and efficient as Solution 1, but a possible work-around is to define a smoother Normalize:

norm = Compile[{{v, _Real, 1}}, 
  If[v == {0, 0}, {0., 0.}, 2/Pi ArcTan[10^4 Total[v^2]] v/Sqrt@Total[v^2]], 
  RuntimeOptions -> EvaluateSymbolically -> False]

Use norm instead of Normalize in the code produces a solution visually the same as that of Solution 1, it takes about 26 seconds on my laptop, though.


Remark

I attempted to spot the root of problem, without success. But NBodySimulation has probably set up the system properly, because even if the problem is set up with NDSolve and WhenEvent manually, the issue remains.

The following is the code. Definition of data is the same as that of OP, definition of reflect is taken from this post, showStatus is from this post:

reflect[vector_, normal_] = -(vector - 2 (vector - Projection[vector, normal]));

f[x_, y_] = x^2 + y^2 - 1;

eq = data["Equations"] /. {Subscript[\[FormalQ], i_] :> Subscript[x, i], 
    Subscript[\[FormalP], i_] :> Subscript[x, i]', \[FormalT] -> t};

event = WhenEvent[
     f @@ #[t] == 
      0, {#'[t] -> 
       reflect[#'[t], {Derivative[1, 0][f] @@ #[t], 
         Derivative[0, 1][f] @@ #[t]}]}] & /@ (Subscript[x, #] & /@ Range[4]);

var = (Subscript[x, #] & /@ Range[4])[t] // Through;

init = <|"Mass" -> Quantity[1, "Kilograms"], "Position" -> Quantity[#, "Meters"], 
     "Velocity" -> Quantity[{0, 0}, "Meters"/"Seconds"], 
     "Charge" -> Quantity[10^-5, "Coulombs"]|> & /@ initPos;  
ic = {var == Through@init["Position"], 
      D[var, t] == Through@init["Velocity"]} /. t -> 0;

showStatus[status_]:=LinkWrite[$ParentLink,
  SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],
                     ToString[status]]];
clearStatus[]:=showStatus[""];
clearStatus[]
jianshi[t_]:=EvaluationMonitor:>showStatus["t = "<>ToString[CForm[t]]]

sol = NDSolveValue[{eq, event, ic} /. 
      HoldPattern@Quantity[a__] :> QuantityMagnitude@Quantity@a /. 
     QuantityMagnitude -> Identity, Head /@ var, {t, 0, countTime}, jianshi[t], 
    MaxSteps -> Infinity(*, Method -> "ExplicitEuler", 
    StartingStepSize -> 10^-3*)]; // AbsoluteTiming

NDSolveValue stops at t==4.09, too.

$\endgroup$
4
  • 1
    $\begingroup$ I note a built-in function NDSolve`NBodySimulationDump`reflect is equal to your -reflect[vector, normal]. but what is jianshi[t] in your code? $\endgroup$
    – yode
    Feb 25, 2022 at 12:17
  • $\begingroup$ @yode Oops, it's something in my tool bag. Code added. $\endgroup$
    – xzczd
    Feb 25, 2022 at 12:28
  • $\begingroup$ BTW, that is a good name. :) And do you think this is a bug about NDSolveValue as your analyze? $\endgroup$
    – yode
    Feb 25, 2022 at 12:29
  • $\begingroup$ @yode Glad you like it. :D Not necessarily a bug, but at least a performance issue, I think. Something might be wrong with the adaptive step choosing of default ODE solver, but of course this is no more than a guess of mine. I suggest you reporting this to WRI. $\endgroup$
    – xzczd
    Feb 25, 2022 at 12:43
6
$\begingroup$

If we put countTime = 20 and remove Normalize then we have desired rest state

n = 4;
countTime = 20;
SeedRandom[5];
initPos = RandomPoint[Disk[], n];
data = NBodySimulation[
   Association["PairwisePotential" -> "Coulomb", "Region" -> Disk[], 
    "ExternalForce" -> (Quantity[-0.5 QuantityMagnitude[#[
           "Velocity"]](*An extra damping*), "Newtons"] &)], <|
      "Mass" -> Quantity[1, "Kilograms"], 
      "Position" -> Quantity[#, "Meters"], 
      "Velocity" -> Quantity[{0, 0}, "Meters"/"Seconds"], 
      "Charge" -> Quantity[10^-5, "Coulombs"]|> & /@ initPos, 
   countTime];
colors = RandomColor[n];

Manipulate[
 Graphics[{Circle[], Red, PointSize[0.02], 
   Riffle[colors, Point /@ data[All, "Position", time]]}, 
  Axes -> True], {time, $MachineEpsilon, countTime, 
  Appearance -> "Labeled"}]

Figure 1

We also can regularized force proposed by yode as follows

n = 4;
countTime = 8;
SeedRandom[5];
initPos = RandomPoint[Disk[], n];
data = NBodySimulation[
   Association["PairwisePotential" -> "Coulomb", "Region" -> Disk[], 
    "ExternalForce" -> (Quantity[-.5 (1 - 
           Exp[-10 Norm[QuantityMagnitude[#["Velocity"]]]]) Normalize[
          QuantityMagnitude[#["Velocity"]]](*An extra damping*), 
        "Newtons"] &)], <|"Mass" -> Quantity[1, "Kilograms"], 
      "Position" -> Quantity[#, "Meters"], 
      "Velocity" -> Quantity[{1. 10^-12, 0}, "Meters"/"Seconds"], 
      "Charge" -> Quantity[10^-5, "Coulombs"]|> & /@ initPos, 
   countTime];
colors = RandomColor[n];
Manipulate[
 Graphics[{Circle[], Red, PointSize[0.02], 
   Riffle[colors, Point /@ data[All, "Position", time]]}, 
  Axes -> True], {time, $MachineEpsilon, countTime, 
  Appearance -> "Labeled"}]

Figure 2

$\endgroup$
6
  • 2
    $\begingroup$ Hi, if we remove Normalize that means that the faster the ball is, the more damped it is, which is probably not a physical reality. The damping should be velocity independent. $\endgroup$
    – yode
    Feb 25, 2022 at 9:47
  • $\begingroup$ @yode Actually damping is velocity dependent, for instance, in the water or air the friction force is proportional to $-k\vec{v}$ at small v and $-k_1|v|\vec{v}$ at large v. $\endgroup$ Feb 25, 2022 at 9:59
  • $\begingroup$ Well, I mean in a disc that ignores air resistance as this post. :) $\endgroup$
    – yode
    Feb 25, 2022 at 10:02
  • $\begingroup$ @yode Do you try to describe dry friction with some minimal velocity to start? $\endgroup$ Feb 25, 2022 at 10:09
  • 1
    $\begingroup$ @yode I see, that for this specific unphysical force $−0.5\vec{v} /v$ there is unphysical result. :) But with regularized force $−0.5(1-e^{10 v})\vec{v} /v$ we have desired rest state - see update to my answer. $\endgroup$ Feb 25, 2022 at 23:24
5
$\begingroup$

This is because of your external force:

(Quantity[
    -0.5 Normalize[QuantityMagnitude[#["Velocity"]]](*An extra damping*)
, "Newtons"] &)

blowing up at small velocity, because of the normalizing factor going infinity $$\hat{v}=\left(\frac{1}{\sqrt{v\cdot v}}\right)v\;,\quad \lim_{v\to\vec{0}}\frac{1}{\sqrt{v\cdot v}}=\infty\;.$$ To avoid this, you could define the external force as a piecewise function:

(Piecewise[{
    {Quantity[-0.5 Normalize[QuantityMagnitude[#["Velocity"]]], "Newtons"], Norm[QuantityMagnitude[#["Velocity"]]] > 10^-1},
    {0, Norm[QuantityMagnitude[#["Velocity"]]] <= 10^-1}
}] &)

As a side note, specifying the units and mapping the quantities to their magnitude is not very necessary but has a significant performance impact. Compare, for example, with this:

n = 4;
countTime = 5;
SeedRandom[5];
initPos = SetPrecision[RandomPoint[Disk[], n], 10];
data = NBodySimulation[
    Association[
        "PairwisePotential" -> "Coulomb",
        "Region" -> Disk[],
        "ExternalForce" -> (Piecewise[{
            {-0.5 Normalize[#["Velocity"]], Norm[#["Velocity"]] > 10^-1},
            {0, Norm[#["Velocity"]] <= 10^-1}
        }] &)
    ],
    Association[
        "Mass" -> 1,
        "Position" -> #,
        "Velocity" -> {0, 0},
        "Charge" -> 10^-5
    ]& /@ initPos
, countTime];
colors = RandomColor[n];
Manipulate[
    Graphics[
        {Circle[], Red, PointSize[0.02], Riffle[colors, Point /@ data[All, "Position", time]]}
    , Axes -> True]
, {time, $MachineEpsilon, countTime, Appearance -> "Labeled"}]

enter image description here


Addendum. This is an explicit explanation as to why the Normalize part blows up in NDSolve. User @xzczd insists otherwise, further claiming that Normalize doesn't transform into anything, such as x/Abs[x]. While that's partially true, it's only when the input is already given as a number. When NDSolve computes a step, it first turns the terms containing variables/target function in the equation into their effective numerical form as much as possible, such as Normalize[{x'[t],y'[t]}] to {x'[t]/Norm[{x'[t],y'[t]}],y'[t]/Norm[{x'[t],y'[t]}]} and then to {x'[t]/Sqrt[Abs[x'[t]^2+y'[t]^2]],y'[t]/Sqrt[Abs[x'[t]^2+y'[t]^2]]}. One can easily check this:

Trace[NDSolve[{Normalize[{f[t], 0}] == {f'[t], 0}, f[0] == 1}, f, {t, 0, 1}]][[1, 1, 1]]

(* {Normalize[{f[t], 0}], {f[t]/Abs[f[t]], 0}} *)

Now, in the case of damping, if "Velocity" or numerical f'[t] gets small at some sufficiently large step (sufficient enough to kill off the precision) during NDSolve, then it's straightforward f'[t]/Abs[f'[t]] blows up:

q1 = 1`10;
q2 = 1`10 + 1/5*10^-9;
dq = q2 - q1;

Precision /@ {q1, q2, dq}
(* {10., 10., 0.} *)

dq/Norm[dq]
(* Error encountered; Indeterminate *)

One workaround is to define the function (damper) piecewise, for small Norm. This is precisely what I did above, which effectively solved the problem.

$\endgroup$
28
  • $\begingroup$ The solution works, but the explanation isn't quite right, actually Normalize[{0, 0}] evaluates to {0, 0}. $\endgroup$
    – xzczd
    Feb 25, 2022 at 9:36
  • $\begingroup$ @xzczd Not quite. 0 is specially returning 0. Otherwise Normalize[x] is equivalent to x/Abs[x]. $\endgroup$
    – Shin Kim
    Feb 25, 2022 at 9:39
  • $\begingroup$ That's what I meant. There's no blow up. $\endgroup$
    – xzczd
    Feb 25, 2022 at 9:47
  • $\begingroup$ @xzczd Normalize[...] never jumps to 0 when ... is a floating number. You're suggesting that Normalize[{0,0}]={0,0} but this doesn't have much point because {0,0} has infinite precision, whereas NDSolve handles floating numbers. A better argument would have been suggesting Normalize[{0., 0.}]={0.,0.}, but when NBodySimulation gets a result having overshooting precision (i.e., higher than the working precision), it returns singularity or stiff system suspected or, in other words, blowing up. $\endgroup$
    – Shin Kim
    Feb 25, 2022 at 10:15
  • $\begingroup$ But you claimed that "This is because of your external force blowing up". The reason why ndsz warning pops up in this case isn't immediatly clear, but I'm pretty sure it's not caused by blowing up of external force. Also, observing the solution up to t==4.09 suggests the solution doesn't blow up, either. (ndsz warning doesn't necessarily mean there's a blow-up, see e.g. mathematica.stackexchange.com/a/219656/1871 , you can find more by searching ndsz in this site. ) $\endgroup$
    – xzczd
    Feb 25, 2022 at 10:26

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.