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};
   Graphics[Disk[pts, 0.1], Frame -> True, PlotRange -> 1.1]]];
  pts += dt;
  If[Abs[pts[[1]]] >= 1, dt[[1]] = -dt[[1]]];
  If[Abs[pts[[2]]] >= 1, dt[[2]] = -dt[[2]]];

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?
  • 1
    $\begingroup$ related: particle bouncing around in an ellipse? and at least closely related: Particles bouncing in a 3D box $\endgroup$
    – Kuba
    Aug 27, 2016 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$ Aug 27, 2016 at 14:16
  • 2
    $\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, 2016 at 16:26
  • $\begingroup$ There is a demo ParticleMotionSimulationUsingAPrioriCollisionDetection/ $\endgroup$
    – Nasser
    Aug 27, 2016 at 21:20

2 Answers 2


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.


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


then schedule updates at a uniform pace:

RunScheduledTask[balls = step[balls], 0.02]

To stop it, simply remove the scheduled task:


What it looks like:

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

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.

    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}];
   WhenEvent[Abs[xm[[i]]] >= 0.9 // Evaluate, 
    D[xm[[i]], t] -> -D[xm[[i]], t] // Evaluate]], {i, n}];
   WhenEvent[Abs[ym[[i]]] >= 0.9 // Evaluate, 
    D[ym[[i]], t] -> -D[ym[[i]], t] // Evaluate]], {i, n}];
   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 -> ∞];

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

enter image description here

  • $\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$ Aug 27, 2016 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, 2016 at 16:30
  • $\begingroup$ @Feyre This is a beautiful simulation. Would you happen to know how one could modify your even structure for particle collisions such that the particles during any collision and number of particles stick together (perfectly inelastically), such that mr'[i] + mr'[j] -> 2*m*r'[k]? $\endgroup$
    – Zach
    Jul 21, 2022 at 18:24

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.