Bouncy Bubbles animation

Question

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];
  ]
 ]

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?

Answer 1

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:

Answer 2

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}]