
Percolation Simulation

This is a simulation of aspects of mathematical percolation theory.

For other percolation simulations, see Category:Percolation Simulations, or:
1D finite grid simulation
Mean run density 
2D finite grid simulations
Site percolation | Bond percolation | Mean cluster density

Let 



c


{\displaystyle c}

 be a 2D boolean square matrix of 



n
×
n


{\displaystyle n\times n}

 values of either 1 or 0 where the
probability of any value being 1 is 



p


{\displaystyle p}

, (and of 0 is therefore 



1
−
p


{\displaystyle 1-p}

).
We define a cluster of 1's as being a group of 1's connected vertically or
horizontally (i.e., using the Von Neumann neighborhood rule) and bounded by either 



0


{\displaystyle 0}

 or by the limits of the matrix.
Let the number of such clusters in such a randomly constructed matrix be 




C

n




{\displaystyle C_{n}}

.

Percolation theory states that 



K
(
p
)


{\displaystyle K(p)}

 (the mean cluster density) will satisfy 



K
(
p
)
=

C

n



/


n

2




{\displaystyle K(p)=C_{n}/n^{2}}

 as 



n


{\displaystyle n}

 tends to infinity. For 



p
=
0.5


{\displaystyle p=0.5}

, 



K
(
p
)


{\displaystyle K(p)}

 is found numerically to approximate 



0.065770


{\displaystyle 0.065770}

...

Task
Show the effect of varying 



n


{\displaystyle n}

 on the accuracy of simulated 



K
(
p
)


{\displaystyle K(p)}

 for 



p
=
0.5


{\displaystyle p=0.5}

 and
for values of 



n


{\displaystyle n}

 up to at least 



1000


{\displaystyle 1000}

.
Any calculation of 




C

n




{\displaystyle C_{n}}

 for finite 



n


{\displaystyle n}

 is subject to randomness, so an approximation should be
computed as the average of 



t


{\displaystyle t}

 runs, where 



t


{\displaystyle t}

 ≥ 



5


{\displaystyle 5}

.

For extra credit, graphically show clusters in a 



15
×
15


{\displaystyle 15\times 15}

, 



p
=
0.5


{\displaystyle p=0.5}

 grid.

Show your output here.

See also
 s-Cluster on Wolfram mathworld.
from __future__ import division
from random import random
import string
from math import fsum
 
n_range, p, t = (2**n2 for n2 in range(4, 14, 2)), 0.5, 5
N = M = 15
 
NOT_CLUSTERED = 1   # filled but not clustered cell
cell2char = ' #' + string.ascii_letters
 
def newgrid(n, p):
    return [[int(random() < p) for x in range(n)] for y in range(n)]
 
def pgrid(cell):
    for n in range(N):
        print( '%i)  ' % (n % 10) 
               + ' '.join(cell2char[cell[n][m]] for m in range(M)))
 
 
def cluster_density(n, p):
    cc = clustercount(newgrid(n, p))
    return cc / n / n
 
def clustercount(cell):
    walk_index = 1
    for n in range(N):
        for m in range(M):
            if cell[n][m] == NOT_CLUSTERED:
                walk_index += 1
                walk_maze(m, n, cell, walk_index)
    return walk_index - 1
 
 
def walk_maze(m, n, cell, indx):
    # fill cell 
    cell[n][m] = indx
    # down
    if n < N - 1 and cell[n+1][m] == NOT_CLUSTERED:
        walk_maze(m, n+1, cell, indx)
    # right
    if m < M - 1 and cell[n][m + 1] == NOT_CLUSTERED:
        walk_maze(m+1, n, cell, indx)
    # left
    if m and cell[n][m - 1] == NOT_CLUSTERED:
        walk_maze(m-1, n, cell, indx)
    # up
    if n and cell[n-1][m] == NOT_CLUSTERED:
        walk_maze(m, n-1, cell, indx)
 
 
if __name__ == '__main__':
    cell = newgrid(n=N, p=0.5)
    print('Found %i clusters in this %i by %i grid\n' 
          % (clustercount(cell), N, N))
    pgrid(cell)
    print('')
 
    for n in n_range:
        N = M = n
        sim = fsum(cluster_density(n, p) for i in range(t)) / t
        print('t=%3i p=%4.2f n=%5i sim=%7.5f'
              % (t, p, n, sim))