What is the statistical equilibrium for this simulation of happy bubbles?

Question

Happy Bubbles

I hope this is not too specific or practical, but I just made a simulation of sorts and seem to have hit quite close to an equilibrium (by accident). Now I am wondering if and how you would calculate the exact function or exact value for this situations statistical "sweetspot".

The Problem

I was playing around with D3.js and created a canvas with randomly appearing bubbles with a 5% chance every iteration. Also every iteration bubbles touching any border will disappear. This can happen because the bubbles move around every iteration.

The Parameters

• Starting with a single bubble
• Every bubble has a x- and a y-coordinate which represents the center of the bubble
• Bubble Size: 20px radius
• Iterationlength: 100 ms
• Canvas Size: 600px by 600px
• Bubble Speed: -5px to 5px on X-Axis and Y-Axis
  • (Math.random() * 10) - 5 (the random number is a floating point number between 0 (including) and 1 (excluding))
• Starting position per bubble: Math.random() * 600 for x and y coordinate (seperately)
• Bubble pops when it gets within 20px of any outer edge (again using the center as reference)

The Equilibrium

With these Parameters new bubbles and disappearing bubbles seem to cancel each other out. The number of bubbles ranges from 100 to 130 but has not breached those limits in more then an hour.

The Question

Do these Parameters ensure such an equilibirum or is it nothing special, that the number seems to be rather stable? How would one calculate this?

Caveats

One I can think of is that the bubble could move from 20.000001 to 15.000001 in one iteration. But I guess that is a small detail. Maybe this isn't even relevant.

Edit:

Here is the jsFiddle! (Be carefull with less then 100 ms iteration time, it might stress the browser)

Answer 1

This is a simple model to estimate $\mathcal{N}$, the expected number of bubbles after some sort of equilibrium is reached. In this model, the bubble will be randomly displaced at each iteration.

The first claim is $\mathcal{N}$ is proportional to $\mathcal{T}$, the expected life time of a bubble . After we start the system and wait long enough for the system to stabilize, there are two ways to estimate $\mathcal{N}$.

1. Observe the system for $N \gg 1 $ iterations. Sum the number of bubbles observed in these iterations and then divide the final result by $N$.

2. Follow all bubbles created and destroyed during these $N$ iterations. Sum over their life time and then divide by $N$.

In general, these two ways of counting will differ by a number which is of the order of $\mathcal{T}$. If $\mathcal{T}$ is finite. This two method will produce the same estimation of $\mathcal{N}$ as $N \to \infty$.

If on average, $\alpha$ bubble is created during each iteration, one will have

$$\mathcal{N} = \alpha \mathcal{T}$$

In this problem $\alpha = 0.05$. If the bubble is created within $20$px from the edge, the bubble will be destroyed in the next iteration. The same thing happens to the bubble at the end of its life. All this process will add $1$ to the expected life time $\mathcal{T}$.

For a probability $\left(\frac{560}{600}\right)^2$, the bubble will be created within the safe zone $[20,580] \times [20,580]$.
Let $L = 580-20 = 560$ and $t$ be the number of iteration since the bubble is created. Introduce random variables $X, Y, \eta_x(t), \eta_y(t)$ and we will update the position of the bubble by following recipe

$$\begin{cases} x(0) &= 20+X\\ y(0) &= 20+Y \end{cases} \quad\text{ and }\quad \begin{cases} x(t) &= x(t-1) + \eta_x(t)\\ y(t) &= y(t-1) + \eta_y(t)\\ \end{cases} \quad\text{ for }\quad t > 0 \tag{*1} $$

Subject to the constraint the bubble is created in the safe-zone, the variables $X, Y$ and $\eta_x(t), \eta_y(t)$ are uniform random variables taking values over $[0,L]$ and $[-5,5]$ respectively.

Let $T$ be the random variable corresponds to the life time of a bubble (excluding the last iteration discussed above). Introduce a function $$u(x,y) = {\bf E}\big[ T \big| (X,Y) = (x,y) \big]\quad\text{ for }\quad (x,y) \in [0,L]\times[0,L]$$ If $u(x,y)$ can be approximated by a continuous function, then the $\mathcal{T}$ we want will be given by a integral

$$\mathcal{T} = 1 + \left(\frac{560}{600}\right)^2 {\bf E}[T] \quad\text{ with }\quad {\bf E}[T] = \frac{1}{L^2} \int_0^L\int_0^L u(x,y) dx dy $$

We can view the random update in $(*1)$ as a Brownian motion. If equilibrium has been reached, we should have

$$u(x,y) = {\bf E}\big[ u(x + \eta_x, y + \eta_y) \big] + 1$$ This implies $u(x,y)$ should approximately satisfies a PDE

$$u(x,y) = u(x,y) + \frac{{\bf E}[\eta_x^2]}{2} \frac{\partial^2 u(x,y)}{\partial x^2} + \frac{{\bf E}[\eta_y^2]}{2} \frac{\partial^2 u(x,y)}{\partial y^2} + 1 $$ Let $\sigma^2 = {\bf E}[\eta_x^2] = \frac{5^2}{3}$, this has the form of a Possion equation $$\nabla^2 u(x,y) = -\frac{2}{\sigma^2}$$ and the appropriate condition for this problem is $u(x,y) = 0$ on the boundaries.

If I'm not mistaken, the solution of these boundary value problem is given by

$$u(x,y) = \frac{32L^2}{\sigma^2\pi^4} \sum_{n=1,\text{odd}}^\infty\sum_{m=1,\text{odd}}^{\infty} \frac{1}{n m (n^2+m^2)} \sin\left(\frac{\pi n x}{L}\right) \sin\left(\frac{\pi m y}{L}\right) $$ Integrate this, we get

$$\mathcal{T} = 1 + \left(\frac{560}{600}\right)^2\frac{128L^2}{\sigma^2\pi^6} K \quad\text{ where }\quad K = \sum_{n=1,\text{odd}}^\infty\sum_{m=1,\text{odd}}^{\infty} \frac{1}{n^2 m^2 (n^2+m^2)} \approx 0.52973 $$ Substitute this into the formula $\mathcal{N} = \alpha T$, we get $$\mathcal{N}_{600} \approx 115.65$$ which now falls within the observed range of OP's numerical result. If one repeat the same argument to other canvas size $s$, we have

$$\mathcal{N}_s = 1 + \left(1 - \frac{40}{s}\right)^2 \frac{384(s - 40)^2}{25\pi^6} K $$ Following is a table for $\mathcal{N}_s$ for other canvas sizes $s$:

$$ \begin{array}{r:r} s & \mathcal{N}_s\\ \hline 300 & 21.54\\ 400 & 44.47\\ 500 & 75.84\\ 600 & 115.65\\ 700 & 163.92\\ 800 & 220.64\\ 900 & 285.82 \end{array} $$

Answer 2

I would make sense that there is some kind of equilibrium. I've been performing some calculations but have not gotten far. Perhaps it is would be interesting to start with a number of bubbles higher than, say, 500, and see whether the number of bubbles converges back to the said range.

Edit: I wrote a script based on your description of the system and obtained very similar results (i.e. the number of bubbles hovered around $120$). However, I wrote another script to simulate the average number of iterations it took for a bubble to reach one of the borders. This number is somewhere around $2500$. In other words, the average bubble has about $\frac{1}{2500}$ chance of reaching a border in the next iteration. In equilibrium, $$n\cdot \frac{1}{\ell_{\mathrm{path}}} = p_{\mathrm{new}},$$ where $n$ is the number of bubbles currently existing, $\ell_{\mathrm{path}}$ is the average path length, and $p_{\mathrm{new}}$ is the probability that a new bubble appears in any given iteration. And as it turns out, $$120 \cdot \frac{1}{2500} = 0.048 \approx 0.05 = p_{\mathrm{new}}.$$ This is, however, not an analytic approach yet.

Answer 3

There will always be an equilibrium. You have a fixed rate of bubble productions, but the popping rate will be proportional to how many bubbles are on screen at the time. Thus the number of bubbles will rise until the popping rate matches.

Answer 4

I fitted the average path lengths I obtained.

$x$ axis: Size of square canvas $\times 100$ pixels.

$y$ axis: Average number of iterations per bubble.