How can I simulate repulsion between multiple point charges (ball bearings) in Mathematica?

Question

I'm trying to write a program in Mathematica that will simulate the way charged ball bearings spread out when they are charged with like charges (they repel each other). My program so far keeps the ball bearings from moving off the screen, and counts the number of times they hit the side of the box. I have the ball bearings moving randomly around the box so far, but I need to know how to make them repel each other.

Here's my code so far:

Manipulate[
 (*If the number of points has been reduced, discard  points*)
 If[ballcount < Length[contents], 
   contents = Take[contents, ballcount]];

 (*If the number of points has been increased, generate some random points*)
 If[ballcount > Length[contents], 
  contents = 
   Join[contents, 
    Table[{RandomReal[{-size, size}, {2}], {Cos[#], Sin[#]} &[
       RandomReal[{0, 2 \[Pi]}]]}, {ballcount - Length[contents]}]]];

 Grid[{{Graphics[{PointSize[0.02],

  (*Draw the container*)
  Line[size {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}], 
  Blend[{Blue, Red}, charge/0.3],
  Point[

   (*Start the main dynamic actions*)
   Dynamic[

    (*Reset the collision counter*)
    collision = 0;

    (*Check for mouse interaction and add points if there has been one*)
    Refresh[
     If[pt =!= lastpt, If[ballcount =!= 50, ballcount++]; 
      AppendTo[
       contents, {pt, {Cos[#], Sin[#]} &[
         RandomReal[{0, 2 \[Pi]}]]}]; lastpt = pt], 
     TrackedSymbols -> {pt}];

    (*Update the position of the points using their velocity values*)
    contents = Map[{#[[1]] + #[[2]] charge, #[[2]]} &, contents];

    (*Check for and fix points that have exceeded the box in Y
      direction, incrementing the collision counter for each one*)
    contents = Map[
      If[Abs[#[[1, 2]]] > size, 
        collision++; {{#[[1, 1]], 
          2 size Sign[#[[1, 2]]] - #[[1, 2]]}, {1, -1} #[[
           2]]}, #] &,
      contents];


    (*Check for and fix points that have exceeded the box in X 
      direction, incrementing the collision counter for each one*)
    contents = Map[
      If[Abs[#[[1, 1]]] > size, 
        collision++; {{2 size Sign[#[[1, 1]]] - #[[1, 1]], #[[1, 
           2]]}, {-1, 1} #[[2]]}, #] &,
      contents];

    hits = Take[PadLeft[Append[hits, collision/size], 200], 200];
    Map[First, contents]]]},
 PlotRange -> {{-1.01, 1.01}, {-1.01, 1.01}}, 
 ImageSize -> {250, 250}],

(*Show the hits*)
Dynamic@Show
  [
   ListPlot
   [
   Take[MovingAverage[hits, smooth], -100
    ]
   ,
   Joined -> True, ImageSize -> {250, 250}, AspectRatio -> 1, 
   PlotLabel -> "number of hits", AxesLabel -> {"time", "hits"}, 
   PlotRange -> {0, Max[Max[hits], 1]}], Graphics[]
  ]
}}
  ]
 ,
 {{pt, {0, 1}}, {-1, -1}, {1, 1}, Locator, Appearance -> None},
 {{ballcount, 5, "number of ball bearings"}, 1, 50, 1},
 {{charge, 0.05, "charge"}, 0.002, 0.3},
 {smooth, 1, ControlType -> None, Appearance -> None},
 {size, 1, ControlType -> None, Appearance -> None},
 {hits, {{}}, ControlType -> None},
 {contents, {{}}, ControlType -> None},
 {lastpt, {{0, 0}}, ControlType -> None}
 ]

Pack

You'd probably want NDSolve to handle the equations of motion. Could approximate a near field by using Nearest at each step and only using particles within a certain distance of a give particle for contributing to it's force field. Maybe not important if you only have a few dozen particles though. – Drowsy

Answer 1

What you need for your simulation is a "collision detection algorithm". The field of those algorithms is widespread since it is as old as computer games (Pong) are and it is impossible to give a complete answer here.

Your simulation as it is now is very basic because you advance your charged balls every time step which makes them "jump" from position to position. If the movement is as simple as it is with the constant velocity and zero acceleration, you know the exact equation of the movement and could calculate all positions by simply putting the time into the equations. When a ball bounces off the wall, it gets a new equation.

With this, you could predict, when two balls will collide. You simply solve for two balls, whether they have at the same time the same position. This is called A Priori detection. When you take your simulation as it is now, you would have to check at every timestep, whether or not two balls are so close together, that they may collide.

The problem there is, that your simulation speed is not infinitely high and the faster your balls are, the bigger the jumps in your simulation. It is then not unlikely, that two balls over-jump each other and you miss a collision.

With this in mind, you could start by reading the Wikipedia article to that topic, to get an overview. Next, you could read some scientific articles about it or check out, how the cracks do it. The Chipmunk physics engine for instance is an amazing 2d-physics engine. To ensure that such stuff works, I'm pretty sure they had to put a lot thoughts in their collision detection.

Dantzler

Answer 2

I can't do the detection part since that is a project on its own. But it seems now the interaction is faster, you had few Dynamics which were not needed. (if you see the side of the cell busy all the time, this always mean a Dynamic is busy where it should not be).

I also added a STOP/START button. I could not understand everything you were doing, but enough to make the changes I made. You are also using AppendTo everything. You should try to allocate contents before hand, and use Part[] to access it, would be much faster, since you seem to know the maximum points allowed?

I like to spread the code out more, this helps me see the logic more.

Here is a screen shot, and the updated version code is below. Hope you find it faster.

Please see code below, in update (1)

Update (1)

(*updated version 12/30/11 9:40 AM*)
Manipulate[(*If the number of points has been reduced,discard points*)
\


 Module[{tbl, rand, npt, ballsToAdd},

  If[running,
   (
    tick += $MachineEpsilon;
    If[ballcount < Length[contents], 
     contents = Take[contents, ballcount]];

    (*If the number of points has been increased,
    generate some random points*)

    If[ballcount > Length[contents],
     (
      ballsToAdd = ballcount - Length[contents];
      tbl = 
       Table[{RandomReal[{-size, size}, {2}], {Cos[#], Sin[#]} &[
          RandomReal[{0, 2 \[Pi]}]]}, {ballsToAdd}];
      contents = Join[contents, tbl]
      )
     ];

    image = Grid[{
       {LocatorPane[Dynamic[pt], Graphics[{

           PointSize[0.02],(*Draw the container*)
           Line[size {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}],
           Blend[{Blue, Red}, charge/0.3],

           Point[(*Start the main dynamic actions*)

            (*Reset the collision counter*)
            collision = 0;

            (*Check for mouse interaction and add points if there has \
been one*)
            If[EuclideanDistance[pt, lastpt] > 0.001, (*adjust*)
             (
              If[ballcount < MAXPOINTS,
               ballcount++
               ];

              rand = RandomReal[{0, 2 \[Pi]}];
              npt = {Cos[rand], Sin[rand]};
              AppendTo[contents, {pt, npt}  ];
              lastpt = pt
              )
             ];

            (*Update the position of the points using their velocity \
values*)

            contents = 
             Map[{#[[1]] + #[[2]] charge, #[[2]]} &, contents];

            (*Check for and fix points that have exceeded the box in \
Y direction,incrementing the collision counter for each one*)

            contents = Map[
              If[Abs[#[[1, 2]]] > size,
                (
                 collision++;
                 {{#[[1, 1]], 
                   2 size Sign[#[[1, 2]]] - #[[1, 2]]}, {1, -1} #[[2]]}
                 ),
                (
                 #
                 )
                ] &, contents
              ];

            (*Check for and fix points that have exceeded the box in \
X direction,
            incrementing the collision counter for each one*)


            contents = 
             Map[If[Abs[#[[1, 1]]] > size, 
                collision++; {{2 size Sign[#[[1, 1]]] - #[[1, 1]], #[[
                   1, 2]]}, {-1, 1} #[[2]]}, #] &, contents
              ];


            hits = Take[PadLeft[Append[hits, collision/size], 200], 
              200];
            Map[First, contents]
            ]
           },
          PlotRange -> {{-1.01, 1.01}, {-1.01, 1.01}}, 
          ImageSize -> {250, 250}
          ], Appearance -> None
         ],(*Show the hits*)

        Show[ListPlot[Take[MovingAverage[hits, smooth], -100],
          Joined -> True, ImageSize -> {250, 250}, AspectRatio -> 1,
          PlotLabel -> "number of hits", AxesLabel -> {"time", "hits"},
          PlotRange -> {0, Max[Max[hits], 1]}
          ]
         ]
        }
       }
      ]
    )
   ];

  image
  ],

 {{MAXPOINTS, 50}, None},
 {pt, {{0, 1}}, None},
 {{ballcount, 5, "number of ball bearings"}, 1, MAXPOINTS, 1, 
  Appearance -> "Labeled", ImageSize -> Small},
 {{charge, 0.05, "charge"}, 0.002, 0.3, Appearance -> "Labeled", 
  ImageSize -> Small},
 Row[{Button["START", {running = True; tick += $MachineEpsilon}], 
   Button["STOP", running = False]}],
 {{tick, 0}, None},
 {smooth, 1, None},
 {size, 1, None},
 {hits, {{}}, None},
 {{contents, {}}, None},
 {lastpt, {{0, 0}}, None},
 {{collision, 0}, None},
 {image, None},
 {{running, True}, None},
 TrackedSymbols ->     { tick},
 ContinuousAction -> False,
 SynchronousUpdating -> True

 ]

Carving

The Point construction in your code produces an error on my version of Mathematica. I don't know what causes the error but placing the code for the main dynamics outside the Point statement (i.e. doing something like <code for updating points>;Point[Map[First, contents]] instead of Point[<code for updating points>;Map[First, contents]]) seems to solve it. – Encephalo

@Heike, thanks for letting know. But this is strange, as it runs OK on my V 8.04, it must be when I copied the code and pasted it to SO, I broke something. I remember moving a little a little after the copy to make it look better. May be I broke it. What I'll do now, is paste what I have, which is unchanged from earlier, and this time, not touch anything after I paste it, and if you could try the new one, and if that works for you, then it was the pasting error, and then I can delete the earlier version. Thanks – Carving

btw, the problem with pasting M code from a notebook to SO and then copy it back to M, is that it gets clobbered up and reformatted differently from the way it looked in the notebook. I think I should use CODE cell next time I do a copy from/to SO/notebook, since a CODE cell in the notebook should not have this reformatting problem. – Carving

Yes, I just verified, it was my silly editing the code in SO after I copied it. I added an extra "," by mistake. Silly me. The code in above update (1) is the correct one. I will now delete the earlier copy. (next time I will learn not to touch code in SO buffer after I paste it). Thanks – Carving