17
$\begingroup$

I want to make a bouncing bubbles animation, looking like this. Here is a start:

DynamicModule[{pts, dt},
 pts = RandomReal[{-1, 1}, 2];
 dt = {0.01, 0.01};
 Print[Dynamic[
   Graphics[Disk[pts, 0.1], Frame -> True, PlotRange -> 1.1]]];
 While[True,
  pts += dt;
  If[Abs[pts[[1]]] >= 1, dt[[1]] = -dt[[1]]];
  If[Abs[pts[[2]]] >= 1, dt[[2]] = -dt[[2]]];
  Pause[0.001];
  ]
 ]

enter image description here

My questions are:

  • How to make the movement of the disk look smoother?
  • How to program collisions between multiple disks?
  • Is it possible to create a 3D version?
$\endgroup$
  • 1
    $\begingroup$ related: particle bouncing around in an ellipse? and at least closely related: Particles bouncing in a 3D box $\endgroup$ – Kuba Aug 27 '16 at 10:59
  • $\begingroup$ if you want to simulate multiple object collisions you need to make a contact detection algorithm. An easy way is to make $n^2$ comparisons and define a force if interaction between the particles, then integrate that with NDSolve as done here, If you change the interaction force for a short range interaction you can get what you ask. If you want some procedural code as in the example you link you should check DEM simulations or MD simulation, depending where are you coming from. $\endgroup$ – tsuresuregusa Aug 27 '16 at 14:16
  • 1
    $\begingroup$ @tsuresuregusa The "contact detection algorithm" in this case is very simple since all the objects are circular. Nearest[objects, pos] will find all collisions with less than $n^2$ comparisons. I frankly don't think that NDSolve is any good for this problem, and Feyre's code doesn't convince me otherwise. It just forces you to formulate the problem in a form that isn't natural. I've never used DEM or MD but it seems outrageously unproportional to the problem of recreating the animation in OP's link (which uses Euler integration.) $\endgroup$ – C. E. Aug 27 '16 at 16:26
  • $\begingroup$ There is a demo ParticleMotionSimulationUsingAPrioriCollisionDetection/ $\endgroup$ – Nasser Aug 27 '16 at 21:20
27
$\begingroup$

This is my port of the Processing code that you referenced. It doesn't try to optimize, so I didn't try it either, for example I didn't use Nearest to find collisions even though that would be fast since it uses a quadtree. I found that I didn't need those optimizations to recreate the animation on the previously linked to website.

About your questions:

1) To make the animation look smooth you should use RunScheduledTask to make sure the updates happen uniformly. I use

RunScheduledTask[balls = step[balls], 0.02]

which means that I get an FPS of 50 frames per second. Of course I have to make sure that balls = step[balls] can be executed within 0.02 seconds; step[balls] // AbsoluteTiming tells me that one step takes 0.002 seconds to execute, meaning it's fast enough to keep up with the desired FPS.

2) This is what the code is all about.

3) It should be trivial to adapt the code for the three dimensional case, this is left as an exercise.

Code:

numBalls = 12;
spring = 0.05;
gravity = 0.03;
friction = -0.9;
width = 640;
height = 360;

balls = ball @@@ Transpose[{
     RandomReal[{0, width}, numBalls],
     RandomReal[{0, height}, numBalls],
     ConstantArray[0, numBalls],
     ConstantArray[0, numBalls],
     RandomReal[{30, 70}, numBalls]
     }];

ballCollide[ball[x1_, y1_, vx1_, vy1_, d1_], 
  ball[x2_, y2_, vx2_, vy2_, d2_]] := Module[{targetX, targetY}, If[
   Norm[{x2 - x1, y2 - y1}] < (d1/2 + d2/2),
   {targetX, targetY} = AngleVector[{x1, y1},{d1/2 + d2/2, ArcTan[x2 - x1, y2 - y1]}];
   ball[x1, y1, vx1 - (targetX - x2) spring, vy1 - (targetY - y2) spring, d1],
   ball[x1, y1, vx1, vy1, d1]
   ]]

boundaryCollide[ball[x_, y_, vx_, vy_, d_]] := Which[
  x + d/2 > width, ball[width - d/2, y, vx friction, vy, d],
  x - d/2 < 0, ball[d/2, y, vx friction, vy, d],
  y + d/2 > height, ball[x, height - d/2, vx, vy friction, d],
  y - d/2 < 0, ball[x, d/2, vx, vy friction, d],
  True, ball[x, y, vx, vy, d]
  ]

move[ball[x_, y_, vx_, vy_, d_]] := ball[x + vx, y + vy - gravity, vx, vy - gravity, d]

step[balls_] := Module[{new},
  new = Fold[ballCollide, #, Complement[balls, {#}]] & /@ balls;
  new = move /@ new;
  boundaryCollide /@ new
  ]

disk[ball[x_, y_, vx_, vy_, d_]] := Disk[{x, y}, d/2]

draw[balls_] := Graphics[{
   White, Opacity[0.8],
   disk /@ balls
   },
  PlotRange -> {{0, width}, {0, height}},
  Background -> Black,
  ImageSize -> width
  ]

To run it, first create a dynamic cell using

Dynamic@draw[balls]

then schedule updates at a uniform pace:

RunScheduledTask[balls = step[balls], 0.02]

To stop it, simply remove the scheduled task:

RemoveScheduledTask[ScheduledTasks[]]

What it looks like:

$\endgroup$
  • 2
    $\begingroup$ Big fun watching this, until they settle at bottom line packed together without any motion. $\endgroup$ – Narasimham Aug 27 '16 at 19:48
23
$\begingroup$

data generates n balls, here: 10

Note that it might be wise to make the box larger, if different WhenEvent[]s happen at the same time (both box and ball collision), one can overwrite the other.

data = Table[{RandomReal[{-0.85, 0.85}, {3}], 
RandomReal[{-1, 1}, {3}], 
RGBColor[RandomReal[{-1, 1}], RandomReal[{-1, 1}], 
 RandomReal[{-1, 1}]]}, {i, 10}];

I've adapted this from my n-body simulation:

Generating lists for use in the NDSolve[]

n = Length[data];
positions = Transpose[Table[data[[i, 1]], {i, n}]];
velocities = Transpose[Table[data[[i, 2]], {i, n}]];
(*Setting up stuff for NDSolve[]*)
{xt, yt, 
   zt} = {ToExpression[Table["x" <> ToString[i], {i, n}]], 
   ToExpression[Table["y" <> ToString[i], {i, n}]], 
   ToExpression[Table["z" <> ToString[i], {i, n}]]};
{xm, ym, zm} = {Through[xt[t]], Through[yt[t]], Through[zt[t]]};
{xz, yz, zz} = {Through[xt[0]], Through[yt[0]], Through[zt[0]]};
yt = ToExpression[Table["y" <> ToString[i], {i, n}]];
zt = ToExpression[Table["z" <> ToString[i], {i, n}]];
rm = Flatten[
    Table[If[i != j, 
      Sqrt[(xm[[j]] - xm[[i]])^2 + (ym[[j]] - ym[[i]])^2 + (zm[[j]] - 
           zm[[i]])^2]], {i, n}, {j, n}]] /. Null -> Sequence[];

No friction or gravity:

xf = Thread[D[D[xm, t], t] == ConstantArray[0, n]];
yf = Thread[D[D[ym, t], t] == ConstantArray[0, n]];
zf = Thread[D[D[zm, t], t] == ConstantArray[0, n]];

Friction and gravity:

xf = Thread[D[D[xm, t], t] == -0.01 D[xm, t]];
yf = Thread[D[D[ym, t], t] == -0.01 D[ym, t]];
zf = Thread[D[D[zm, t], t] == -0.01 D[zm, t] - ConstantArray[0.5, n]];

The final equations for the differential equation:

pos = {Thread[xz == positions[[1]]], Thread[yz == positions[[2]]], 
   Thread[zz == positions[[3]]]};
vel = {Thread[D[xm, t] == velocities[[1]]] /. t -> 0, 
   Thread[D[ym, t] == velocities[[2]]] /. t -> 0, 
   Thread[D[zm, t] == velocities[[3]]] /. t -> 0};
col = Table[data[[i, 3]], {i, n}];
we = {};

A list of WhenEvent[]s, note that the ball collisions aren't exactly right, I'm too tired right now to properly add elastic collision equations, feel free to edit the code. The first set is the ball collisions, the second is the boundary box collisions. In these WhenEvent[] collisions we can reduce total energy by little amounts per collision, as in the animation.

Table[AppendTo[we, 
   WhenEvent[
    Sqrt[(xm[[i]] - xm[[j]])^2 + (ym[[i]] - ym[[j]])^2 + (zm[[i]] - 
           zm[[j]])^2] <= 0.2 // Evaluate, {
      D[xm[[i]], t] -> 0.5 (D[xm[[i]], t] + D[xm[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[xm[[i]], t] + D[xm[[j]], t]), 
      D[xm[[i]], t] -> 0.5 (D[ym[[i]], t] + D[ym[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[ym[[i]], t] + D[ym[[j]], t]), 
      D[xm[[i]], t] -> 0.5 (D[zm[[i]], t] + D[zm[[j]], t]), 
      D[xm[[j]], t] -> 0.5 (D[zm[[i]], t] + D[ym[[j]], t])} // 
     Evaluate]], {i, n - 1}, {j, i + 1, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[xm[[i]]] >= 0.9 // Evaluate, 
    D[xm[[i]], t] -> -D[xm[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[ym[[i]]] >= 0.9 // Evaluate, 
    D[ym[[i]], t] -> -D[ym[[i]], t] // Evaluate]], {i, n}];
Table[AppendTo[we, 
   WhenEvent[Abs[zm[[i]]] >= 0.9 // Evaluate, 
    D[zm[[i]], t] -> -D[zm[[i]], t] // Evaluate]], {i, n}];
s = NDSolve[Flatten[{xf, yf, zf,
     vel, pos, we}], Flatten[{xt, yt, zt}], {t, 0, 50}, 
   MaxSteps -> ∞];

Manipulate[
 Graphics3D[
  Transpose[{col, 
    Thread[Sphere[
      Evaluate[
       Flatten[Thread[Transpose[{xm, ym, zm} /. t -> a] /. s], 1]], 
      0.1]]}], PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}], {a, 0, 50, 
  0.1}]

enter image description here

$\endgroup$
  • $\begingroup$ I was looking for how to make work WhenEvent with lists and could not find the answer, this is great. Would it be possible to make WhenEvent dynamically look only for particles that are nearby instead for all the possible combinations? $\endgroup$ – tsuresuregusa Aug 27 '16 at 15:28
  • $\begingroup$ @tsuresuregusa I almost asked a question about WhenEvent[] lists, it should really be in the documentation. The nearby particle thing, as far as I'm aware it's not possible, it'd be nice if there were a WhenCondition[] which could operate as a large If[] block for WhenEvent[]s. $\endgroup$ – Feyre Aug 27 '16 at 16:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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