Random+sequential+deposition

toc =Introduction= During the first lectures, we have seen how in Monte-Carlo algorithms one defines a dynamics (through transition probabilities) which enables to sample a given distribution //π//(//C//) on a set of configuration {//C//}.
 * Class Session 09: Random sequential deposition **

One may consider a reversed point of view, where now (//e.g.// physically motivated) transition rates are given, and one is interested in finding/characterizing the unknown steady state distribution. In this exercice, we focus on the problem of random sequential deposition, defined for hard disks on a one-dimensional segment (see Class Session 02).



=Model= We consider the problem of hard "spheres" in one dimension, that we will call //clothe pins//, on a segment of length //L//. Pins of width 2//σ// are permanently deposited at random positions in the box provided they do not overlap with pins placed earlier. The deposition consists in trying to deposit a new pin with its center uniformly distributed in the interval [//σ,L//–//σ//]. This process goes on until no pins can be deposited anymore.

Time //t// is defined as follows: the first deposition takes place at //t// = 1, and //t// is increased by 1 at each (successfull or not) tentative deposition. One is interested for instance in the last deposition time, and in the statistical properties of the final configuration. One could also study the density as a function of time, and the dynamical approach towards the fully packed state.
 * Is this problem equivalent to the one-dimensional hard sphere problem of Class Session 02? Why (Why not)?

=Naive algorithm= The following (naive) algorithm implements the sequential deposition up to time tmax: code format="python" import math,random,pylab random.seed(3456789732) N = 1; L = 10; sigma = 0.1; tmax = 80; def dist(a,b): return math.fabs(b - a) def allowed(positions,newpos): dists=[ dist(pos,newpos) for pos in positions ] return min(dists) > 2*sigma positions=[ random.uniform(sigma,L-sigma) ] lastt = 1 for t in range(tmax): pos = random.uniform(sigma,L-sigma); if allowed(positions,pos): positions.append(pos) lastt = t       N += 1 print 'pin number ',N, ', deposition time: ',t print '\n last deposition time: ',lastt, ' < tmax =',tmax height = .33 * L square = pylab.Rectangle((sigma,0), L-2*sigma, height, fc='b') pylab.gca.add_patch(square) for x in positions: whiterec = pylab.Rectangle((x-2*sigma,0), 4*sigma,height, fc='w',ec='w') pylab.gca.add_patch(whiterec) for x in positions: redrec = pylab.Rectangle((x-sigma,0), 2*sigma,height, fc='r') pylab.gca.add_patch(redrec) pylab.axis('scaled') pylab.axis([0,L,0,height]) pylab.title('Clothe pins are in red, remaining available space for is in blue.') pylab.show code In this algorithm:
 * Is the rejection probability increasing with time?
 * Do you know when to stop the deposition?

=Faster-than-the-clock algorithm=

The idea
Consider a given configuration of the system.
 * Write the rejection probability //λ// in terms of the accessible length and of the total deposition length.
 * Starting from a given configuration, we define //t// accept as the number of time steps one has to wait before the deposition of a new pin (occurring in //t// accept +1) is accepted. Determine the probability distribution of //t// accept.
 * How would you sample this probability?

Implementation of the faster than the clock algorithm
code format="python" import math,random,pylab random.seed(34567897) N = 1; L = 10; sigma = .1; intervals = [(sigma,L-sigma)]  # list of the (ak,bk) dists = [b-a for (a,b) in intervals] positions = [] time = 0 def tower(data,Upsilon): sums = [0] for d in data: sums.append(d+sums[-1]) sums.pop(0) for k in range(len(sums)): if Upsilon < sums[k]: break return k while len(intervals) > 0: length = sum(dists) # choice of the time jump lambd = 1 - length/(L-2*sigma); # rejection rate deltat = 1 if (lambd > 0.): deltat += int(math.log(random.random)/math.log(lambd)) time += deltat N += 1 # evolution Upsilon = random.uniform(0,length) k = tower(dists,Upsilon) d2 = sum(dists[p] for p in range(k+1))-Upsilon (a,b) = intervals.pop(k) dists.pop(k) newpos = b - d2   positions.append(newpos) d1 = newpos - a;   if (d1 > 2*sigma): intervals.append((a,newpos-2*sigma));dists.append(newpos-2*sigma-a) if (d2 > 2*sigma): intervals.append((newpos+2*sigma,b));dists.append(b-newpos-2*sigma)

print 'total pin number: ',N,' ultimate deposition time: ',time height = .33*L square = pylab.Rectangle((sigma,0), L-2*sigma, height, fc='b') pylab.gca.add_patch(square) for x in positions: whiterec=pylab.Rectangle((x-2*sigma,0), 4*sigma,height, fc='w',ec='w') pylab.gca.add_patch(whiterec) for x in positions: redrec=pylab.Rectangle((x-sigma,0), 2*sigma,height, fc='r') pylab.gca.add_patch(redrec) pylab.axis('scaled') pylab.axis([0,L,0,height]) pylab.show code =References=
 * Y. Pomeau, //Some asymptotic estimates in the random parking problem//, J. Phys. A: Math. Gen. 13 L193 (1980)
 * V. Privman, J.-S. Wang, and P. Nielaba, //Continuum limit in random sequential adsorption//, Phys. Rev. B 43 3366 (1991)

url}?f=print| [Print this page]