Rain droplets falling on a table

There's a nice approach in Coverage by Randomly Deployed Wireless Sensor Networks by Wan and Yi that allows us to get an expansion in $\mu:=R/r$ for the expected number of raindrops required.

The basic idea is to focus on the intersections among the circles. Since any uncovered area must be bounded by arcs bounded by intersections that aren't covered by other disks, the entire table is covered if and only if there is at least one intersection and there is no uncovered intersection. Here intersections include intersections between raindrop circles and the table circle.

First consider the probability of there being at least one intersection, uncovered or not, between a raindrop circle and the table circle. The centre of the raindrop must lie within $r$ of the table circle, and the probability for that is $1-\pi(R-r)^2/(\pi R^2)=2\mu^{-1}\left(1+O\left(\mu^{-1}\right)\right)$. Since we need at least $\mu^2$ drops to cover the table, the probability that none of them intersects the table circle goes to zero. Thus, asymptotically the probability that the table is covered becomes the probability that there are no uncovered intersections.

Since the probability for a raindrop circle to intersect the table circle is $2\mu^{-1}\left(1+O\left(\mu^{-1}\right)\right)$, the expected number of intersections between a raindrop circle and the table circle is $2n\mu^{-1}\left(1+O\left(\mu^{-1}\right)\right)$, where $n$ is the number of drops. Each of these can be covered by any of the remaining $n-1$ disks, with probability $\mu^{-2}$ each (neglecting border effects), so the expected number of uncovered intersections of this type is $2n\mu^{-1}\left(1-\mu^{-2}\right)^{n-1}\left(1+O\left(\mu^{-1}\right)\right)$. For two raindrop circles to intersect, their centres have to be at most $2r$ apart, and (again neglecting border effects) the probability for this is $\pi(2r)^2/\pi R^2=4\mu^{-2}$. Thus the expected number of intersections among $n$ raindrop circles is $4n(n-1)\mu^{-2}$ (since there are $n(n-1)/2$ pairs and each intersecting pair contributes $2$ intersections). Again, each of these can be covered by any of the remaining $n-2$ disks, with probability $\mu^{-2}$ each, so the expected number of uncovered intersections of this type is $I(n)=4n(n-1)\mu^{-2}\left(1-\mu^{-2}\right)^{n-2}$. Since $n\ge\mu^2$, this dominates the contribution from intersections with the table circle by a factor of $\mu$, so we can neglect intersections with the table circle if we only want the leading terms of the expansion; this is the same level of approximation as the one applied in taking the area fraction covered by a raindrop to be $\mu^{-2}$ and neglecting border effects.

To that approximation, $I(n)$ is an upper bound for the probability that at least one intersection, and thus some part of the table, is uncovered. At $n=\mu^2$, this is approximately $(4/\mathrm e)\mu^2$, and thus greater than $1$ for sufficiently large $\mu$. Since $n\ge\mu^2$, we can safely replace the upper bound by $1$ up to $\mu^2$, and then also up to the point $n_0$ where $I(n)$ falls below $1$ again.

Since $n\ge\mu^2$, we can approximate $n-1$ and $n-2$ by $n$ and obtain $n_0$ as the greater of the two real solutions of the resulting equation $4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n=1$:

$$ n_0\approx\frac{2W_{-1}\left(\mu\log\left(1-\mu^{-2}\right)/4\right)}{\log\left(1-\mu^{-2}\right)}\;, $$

where $W_{-1}$ is the lower branch of the Lambert W function.

Now the expected number of drops required to cover the table is

$$ \begin{align} \langle N\rangle &= \sum_{k=0}^\infty kP(N=k) \\ &= \sum_{k=0}^\infty P(N\gt k) \\ &\lesssim \sum_{k=0}^\infty\min(1,I(k)) \\ &\le n_0+\sum_{k=\lfloor n_0\rfloor}^\infty I(k)\;. \end{align} $$

We can approximate the remaining sum by an integral:

$$ \begin{align} \sum_{k=\lfloor n_0\rfloor}^\infty I(k) &\approx \int_{n_0}^\infty \left.4\mu^{-2}\frac{\partial^2}{\partial q^2}q^k\right|_{q=1-\mu^{-2}}\,\mathrm dk \\ &= -\left.4\mu^{-2}\frac{\partial^2}{\partial q^2}\frac{q^{n_0}}{\log q}\right|_{q=1-\mu^{-2}} \\ &= -\left.4\mu^{-2}n_0(n_0-1)\frac{q^{n_0-2}}{\log q}\right|_{q=1-\mu^{-2}}\left(1+O\left(\mu^{-2}\right)\right) \\ &= -\frac1{\log\left(1-\mu^{-2}\right)}\left(1+O\left(\mu^{-2}\right)\right) \\ &= \mu^2+O(1)\;. \end{align} $$

Now we can use the asymptotic expansion of $W_{-1}$,

$$ W_{-1}(x)=\log(-x)-\log(-\log(-x))+\frac{\log(-\log(-x))}{\log(-x)}+O\left(\left(\frac{\log(-\log(-x))}{\log(-x)}\right)^2\right)\;, $$

to get an expansion for $\langle N\rangle$. If we want only the terms up to $O\left(\mu^2\right)$ we can replace $\log\left(1-\mu^{-2}\right)$ by $-\mu^{-2}$ in both instances to obtain

$$ \begin{align} n_0 &\approx -2\mu^2W_{-1}\left(-\frac1{4\mu}\right) \\ &= 2\mu^2\left(\log(4\mu)+\log(\log(4\mu))+\frac{\log\log(4\mu)}{\log(4\mu)}+O\left(\left(\frac{\log\log\mu}{\log\mu}\right)^2\right)\right)\;. \end{align} $$

Then adding the term from the remaining sum yields

$$ \langle N\rangle\lesssim2\mu^2\left(\log(4\mu)+\log(\log(4\mu))+\frac12+\frac{\log\log(4\mu)}{\log(4\mu)}+O\left(\left(\frac{\log\log\mu}{\log\mu}\right)^2\right)\right)\;, $$

where the $\lesssim$ and the big $O$ should be taken with a grain of salt, since I analyzed neither the error from replacing $n-1$ and $n-2$ by $n$, nor the one from approximating the sum by an integral, nor whether the upper bound from the first moment method is asymptotic. Perhaps someone more rigorous-minded than I can fill in some of the gaps.

I did some exact (i.e. non-discretized) simulations; the results are in good agreement with the above expansion. The code is here. I ran $1000$ trials each for $1/20\le\mu^{-1}\le1$ in steps of $1/100$.

Here's a log-log plot of $\langle N\rangle$ over $\mu$:

(You can get the full resolution by displaying the image in another tab or window.) The red crosses are the data points, and the green, blue, magenta and cyan curves respectively include successive terms of the above expansion. The fact that successive terms get increasingly close to the data seems to indicate that the unanalyzed errors are at most $O(\mu^2)$.

The estimated standard errors for the data points from the simulation are too small to display. The values and errors at different $\mu$ are strongly correlated, since the code works by successively adding raindrop centres and lowering a threshold $r$ at which drops at these centres would cover the table; but the trial runs as a whole are independent, so the error estimates are valid.

A few words about the code: It recursively subdivides the square containing the unit circle; for each square, it computes three thresholds for $r$: the least $r$ for which the square would be partially covered by exactly one drop, the least $r$ for which the square would be partially covered by more than one drop, and the least $r$ for which the entire square would be wet. These thresholds divide the $r$ values into four ranges; in the first two, the square is known not to be fully covered, in the last one it is known to be fully covered, and in between (if there is more than one drop partially covering it) the recursion has to descend to the next level to decide. Raindrops fall one after the other, and after each drop the statistics for the average and the variance estimate are collected for the radii up to which the table is fully covered. Of course there are some complications, such as dealing with squares on the border of the table.

Disclaimer: This sort of code is difficult to fully debug – I did some testing with known configurations and I hope that I caught all the bugs, and the agreement with the analytic expansion is reassuring, but there's no guarantee that some cases aren't perhaps being treated incorrectly.

[Update:]

Since mjqxxxx elegantly proved that the expected number $\langle N\rangle$ of drops is bounded from below by $2\mu^2\log\mu(1+o(1))$, I thought it's worthwhile to make at least the part of my argument fully rigorous that shows that $\langle N\rangle$ is also bounded from above by $2\mu^2\log\mu(1+o(1))$, since together these imply $\langle N\rangle\sim2\mu^2\log\mu$.

So consider again the intersections among circles. Since we will be bounding the probability of the table not being fully covered by its maximal value $1$ for $n\lt\mu^2$, we can disregard the possibility of there being no intersections at all, since this is impossible for $n\ge\mu^2$. Thus we only need to consider the probability of there being an uncovered intersection. As shown above, this is asympotically dominated by the probability of there being an uncovered intersection between raindrop circles, which is $4n(n-1)\mu^{-2}\left(1-\mu^{-2}\right)^{n-2}(1+o(1))$ (where the $o(1)$ term accounts for border effects in the area fractions $\mu^{-2}$ and $4\mu^{-2}$). We can bound this from above by

$$ 4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n\left(1-\mu^{-2}\right)^{-2}(1+o(1))=4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n(1+o(1))\;. $$

From the solution for $4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n=1$ above in terms of $W_{-1}$ and its asymptotic expansion, we know that we can choose an integer $n_0=2\mu^2\log\mu(1+o(1))$ such that this bound is $1+o(1)$. Then we can bound the probability of the table not being fully covered by $1$ for $n\le n_0$ and by the above bound for $n\gt n_0$. This gives a contribution $n_0$ to the expected number of disks from $n\le n_0$, which has the right asymptotic behaviour, so we only need to show that the contribution from $n\gt n_0$ is $o\left(\mu^2\log\mu\right)$. To that end, take the logarithm, bound it linearly using $\log n\le\log n_0+(n-n_0)/n_0$, and sum the resulting geometric series:

$$ \begin{align} &\log\left(4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n(1+o(1))\right) \\ =& \log\left(4\mu^{-2}\right)+2\log n+n\log\left(1-\mu^{-2}\right)+o(1) \\ \le& \log\left(4\mu^{-2}\right)+2\log n_0+2(n-n_0)/n_0+n\log\left(1-\mu^{-2}\right)+o(1) \\ =& (n-n_0)\left(2/n_0+\log\left(1-\mu^{-2}\right)\right)+o(1)\;, \end{align} $$

$$ \begin{align} &\sum_{n=n_0}^\infty\left(4n^2\mu^{-2}\left(1-\mu^{-2}\right)^n\left(1+o(1)\right)\right) \\ \le& \sum_{n=n_0}^\infty\exp\left((n-n_0)\left(2/n_0+\log\left(1-\mu^{-2}\right)\right)+o(1)\right) \\ =& \frac1{1-\exp\left(2/n_0+\log\left(1-\mu^{-2}\right)\right)}(1+o(1)) \\ =& \frac1{1-\left(1-\mu^{-2}\right)\exp\left(2/n_0\right)}(1+o(1)) \\ =& \frac1{1-\left(1-\mu^{-2}\right)\left(1+o\left(\mu^{-2}\right)\right)}(1+o(1)) \\ =& \mu^2(1+o(1))\;. \end{align} $$

Thus $\langle N\rangle$ is bounded from above by $2\mu^2\log\mu(1+o(1))$, and together with mjqxxxx' result it follows that $\langle N\rangle\sim2\mu^2\log\mu$.

P.S.: Note that neither mjqxxxx' argument nor mine depends on the shapes of the drops or the table (as regards the leading term, not the corrections derived above). Thus we can conclude generally that the expected number of objects with area a fraction $\lambda$ of the table area that are required to cover the entire table is asymptotic to $\lambda\log\lambda$.

Allow me to make some discretization approximations to a circle. Define a circle of radius $r$ at location $\bf x$ as the set of all points $\bf y$ on a lattice in $\mathbb{Z}_d$ where $\|{\bf x}-{\bf y}\|_2<r$. This isn't to bad when $r \gg 1$, shown below is a 100 radius "circle":

The approximation gets worse when $r<10$, so we will limit all droplets to this size. To code up a simple Monte Carlo example we make a square area $A$ whose length is $2R + 2r$, and put a circular table in the center of radius $R$ with a positive mask. Choose points in $A$ and keep them if they are in within the radius $R$. For each selected point, subtract a circular mask of radius $r$. Continue until the circular table is empty and count the drops needed. An example of this in motion looks like:

A time-series of the log of the area remaining for an individual run:

I ran 10 trials for $r \in 10..R-1$ and plotted the mean with the std dev as the error bar.

I'm at a loss to describe the function - the best fit I could come up with (which really is not fit at all) was $c \exp{(-17.9 (r/R)^{-1.6})}$. At the very least, I can say that it drops faster than an exponential.

For posterity and completeness, the code I wrote is below. Note that you can generalize to $d$ dimensions!

from numpy import *

def circle(radius, dimension):
    A = zeros((2*radius+1,)*dimension)
    for idx in ndindex(A.shape):
        coord = array(idx) - radius # recenter
        r     = dot(coord, coord)
        A[idx] = r
    B = zeros(A.shape,dtype=bool)
    B[A<radius**2] = True
    return B

def spot_on_table(dimension):
    table_r_squared = table_r ** 2
    r2 = table_r_squared+1
    max_dim = 2*table_r
    while r2 > table_r_squared:
        pt = random.randint(-max_dim,max_dim,dimension)
        r2  = dot(pt,pt)
    return pt

def raindrop(table, droplet, loc, add=False):
    r = droplet.shape[0]/2
    center = array(table.shape[0])/2
    idx = [slice(center-r+x, center+r+1+x) for x in loc]
    idx = tuple(idx)

    if add: table[idx] += droplet
    else  : table[idx] *= ~droplet

def raindrop_question(dimension, droplet_r, table_r):
    random.seed()
    buff = droplet_r

    # Create the table
    table = zeros((2*table_r+1+2*buff,)*dimension)
    raindrop(table,circle(table_r,dimension), zeros(dimension), add=True)

    # Let it rain!
    counter = 0
    while table.any():
        raindrop(table,circle(droplet_r,dimension), spot_on_table(dimension))
        counter += 1

    return (dimension, droplet_r, table_r), counter

def CB(sol):
    print sol

import multiprocessing, itertools

dimension = 2
table_r   = 100
droplet_r = arange(10,table_r)

P = multiprocessing.Pool(16)

for r in droplet_r:
    for n in xrange(10):
        sol = P.apply_async(raindrop_question, (dimension, r, table_r),
                            callback=CB)
P.close()
P.join()

For any covering of the target region with shapes of diameter $\le r$, if one droplet center lies within each shape, this is sufficient to cover the entire region with droplets. This can be used to calculate an upper bound on the expected time to cover the table. Note that this upper bound will be tighter when the number of shapes is smaller; i.e., we want the area of each shape to be as large as possible. Suppose the covering uses regular hexagons of diameter $r$: these have area $3r^2\sqrt{3}/8$, so there will be about $$ N_{\rm hex} = \frac{8\pi}{3\sqrt{3}}\left(\frac{R}{r}\right)^2 $$ of them. The coupon collector's problem then states that the expected time to hit all the hexagons is $E[T_{\rm hex}] \sim N_{\rm hex}\ln N_{\rm hex}$.

On the other hand, for any packing of the target region with disks of radius $r$, it is necessary for at least one droplet center to fall within each disk in order for the target region to be covered (in particular, the disk centers cannot be covered otherwise). This can be used to provide a lower bound on the expected time. Here we again use a hexagonal packing, which covers a fraction $\pi/\sqrt{12}$ of the table's surface, but now the hexagon area is $2r^2\sqrt{3}$. The number of disks is then $$ N_{\rm disc}=\frac{\pi}{2\sqrt{3}}\left(\frac{R}{r}\right)^2=\frac{3}{16}N_{\rm hex}. $$ The coupon collector's problem in this case states that the expected time to hit all the disks is $E[T_{\rm disc}] \sim \frac{\sqrt{12}}{\pi}N_{\rm disc}\ln N_{\rm disc}=\left(\frac{R}{r}\right)^{2}\ln N_{\rm disc}$. Note the extra factor of $\sqrt{12}/\pi$, which comes from the fact that each droplet only "collects a coupon" with probability $\pi/\sqrt{12}$.

Putting these bounds together, we have that the expected time to cover the table satisfies $$ 2 \mu^2 \ln \mu + O(\mu^2) = E[T_{\rm disc}] \le E[T] \le E[T_{\rm hex}] = \frac{16\pi}{3\sqrt{3}} \mu^2 \ln \mu + O(\mu^2) $$ for $\mu \equiv (R/r) \gg 1$. In particular, the expected time is proven to be $\Theta(\mu^2 \ln \mu)$, with a multiplicative factor bounded between $2$ and about $9.68$.