## Percolation threshold estimation

### August 30, 2013

The first assignment for Algorithms 1 is an estimation for the percolation threshold in an n-by-n array of sites; the system is said to percolate whenever there exists an uninterrupted series of connections between the top “edge” and the bottom one.

We need a fast lookup for the neighborhood of any given site; the neighbors() method checks whether a site is in the middle or on the sides or corners of the grid:

def neighbors(ii, n):
"returns list of neighbors of ii in n-by-n grid"
l = []
if mod(ii, n) > 0:
l.append(ii-1)
if mod(ii-n+1, n) > 0 or ii == 0:
l.append(ii + 1)
if ii > n-1:
l.append(ii - n)
if ii < n*(n-1):
l.append(ii + n)
return l

I use two lists: L is a pointer array as in the previous post and O contains the state tags (‘#’ for closed, ‘.’ for open site). There are two extra sites, serving as roots for the top, resp. bottom sites.

def percolation_frac(L, O, sites):
for ii, isite in enumerate(sites):
O[isite] = '.'
if isite < n:
L[isite] = ridxTop
elif n*(n-1) < isite <= N:
L[isite] = ridxBot
neighb = neighbors(isite, n)
for inb in neighb:
if O[inb]=='.':
L = wqupc2(L, isite, inb)
if L[-1]==L[-2]:
""" percolation:
top and bottom virtual sites
belong to same component
"""
pf = float(ii) / N
return pf

As there is no analytical means of predicting when the percolation phase transition is going to occur, we just run the experiment in Monte Carlo mode (many times with randomized initialization), and average the per-run percolation threshold estimate.
Theory tells us that in a square lattice, this value lies around 0.5927 .
My implementation bumps into a few false positives, possibly due to the order in which the union operation is performed, thereby skewing the average more toward 0.6something. Need to look into this.