
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.
using Printf, Distributions
 
newgrid(p::Float64, r::Int, c::Int=r) = rand(Bernoulli(p), r, c)
 
function walkmaze!(grid::Matrix{Int}, r::Int, c::Int, indx::Int)
    NOT_CLUSTERED = 1 # const
    N, M = size(grid)
    dirs = [[1, 0], [-1, 0], [0, 1], [0, -1]]
    # fill cell
    grid[r, c] = indx
    # check for each direction
    for d in dirs
        rr, cc = (r, c) .+ d
        if checkbounds(Bool, grid, rr, cc) && grid[rr, cc] == NOT_CLUSTERED
            walkmaze!(grid, rr, cc, indx)
        end
    end
end
 
function clustercount!(grid::Matrix{Int})
    NOT_CLUSTERED = 1 # const
    walkind = 1
    for r in 1:size(grid, 1), c in 1:size(grid, 2)
        if grid[r, c] == NOT_CLUSTERED
            walkind += 1
            walkmaze!(grid, r, c, walkind)
        end
    end
    return walkind - 1
end
clusterdensity(p::Float64, n::Int) = clustercount!(newgrid(p, n)) / n ^ 2
 
function printgrid(G::Matrix{Int})
    LETTERS = vcat(' ', '#', 'A':'Z', 'a':'z')
    for r in 1:size(G, 1)
        println(r % 10, ") ", join(LETTERS[G[r, :] .+ 1], ' '))
    end
end
 
G = newgrid(0.5, 15)
@printf("Found %i clusters in this %i×%i grid\n\n", clustercount!(G), size(G, 1), size(G, 2))
printgrid(G)
println()
 
const nrange = 2 .^ (4:2:12)
const p = 0.5
const nrep = 5
for n in nrange
    sim = mean(clusterdensity(p, n) for _ in 1:nrep)
    @printf("nrep = %2i p = %.2f dim = %-13s sim = %.5f\n", nrep, p, "$n × $n", sim)
end