
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.
import random, sequtils, strformat
 
const
  N = 15
  T = 5
  P = 0.5
  NotClustered = 1
  Cell2Char = " #abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
  NRange = [4, 64, 256, 1024, 4096]
 
type Grid = seq[seq[int]]
 
 
proc newGrid(n: Positive; p: float): Grid =
  result = newSeqWith(n, newSeq[int](n))
  for row in result.mitems:
    for cell in row.mitems:
      if rand(1.0) < p: cell = 1
 
 
func walkMaze(grid: var Grid; m, n, idx: int) =
  grid[n][m] = idx
  if n < grid.high and grid[n + 1][m] == NotClustered:
    grid.walkMaze(m, n + 1, idx)
  if m < grid[0].high and grid[n][m + 1] == NotClustered:
    grid.walkMaze(m + 1, n, idx)
  if m > 0 and grid[n][m - 1] == NotClustered:
    grid.walkMaze(m - 1, n, idx)
  if n > 0 and grid[n - 1][m] == NotClustered:
    grid.walkMaze(m, n - 1, idx)
 
 
func clusterCount(grid: var Grid): int =
  var walkIndex = 1
  for n in 0..grid.high:
    for m in 0..grid[0].high:
      if grid[n][m] == NotClustered:
        inc walkIndex
        grid.walkMaze(m, n, walkIndex)
  result = walkIndex - 1
 
 
proc clusterDensity(n: int; p: float): float =
  var grid = newGrid(n, p)
  result = grid.clusterCount() / (n * n)
 
 
proc print(grid: Grid) =
  for n, row in grid:
    stdout.write n mod 10, ") "
    for cell in row:
      stdout.write ' ', Cell2Char[cell]
    stdout.write '\n'
 
 
when isMainModule:
 
  randomize()
 
  var grid = newGrid(N, 0.5)
  echo &"Found {grid.clusterCount()} clusters in this {N} by {N} grid\n"
  grid.print()
  echo ""
 
  for n in NRange:
    var sum = 0.0
    for _ in 1..T:
      sum += clusterDensity(n, P)
    let sim = sum / T
    echo &"t = {T}  p = {P:4.2f}  n = {n:4}  sim = {sim:7.5f}"