
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" for Random
import "/fmt" for Fmt
 
var rand = Random.new()
var RAND_MAX = 32767
 
var list = []
var w = 0
var ww = 0
 
var ALPHA = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
var ALEN  = ALPHA.count - 3
 
var makeList = Fn.new { |p|
    var thresh = (p * RAND_MAX).truncate
    ww = w * w
    var i = ww
    list = List.filled(i, 0)
    while (i != 0) {
        i = i - 1
        var r = rand.int(RAND_MAX+1)
        if (r < thresh) list[i] = -1
    }
}
 
var showCluster = Fn.new {
    var k = 0
    for (i in 0...w) {
        for (j in 0...w) {
            var s = list[k]
            k = k + 1
            var c = (s < ALEN) ? ALPHA[1 + s] : "?"
            System.write(" %(c)")
        }
        System.print()
    }
}
 
var recur // recursive
recur = Fn.new { |x, v|
    if (x >= 0 && x < ww && list[x] == -1) {
        list[x] = v
        recur.call(x - w, v)
        recur.call(x - 1, v)
        recur.call(x + 1, v)
        recur.call(x + w, v)
    }
}
 
var countClusters = Fn.new {
    var cls = 0
    for (i in 0...ww) {
        if (list[i] == -1) {
            cls = cls + 1
            recur.call(i, cls)
        }
    }
    return cls
}
 
var tests = Fn.new { |n, p|
    var k = 0
    for (i in 0...n) {
        makeList.call(p)
        k = k + countClusters.call() / ww
    }
    return k / n
}
 
w = 15
makeList.call(0.5)
var cls = countClusters.call()
System.print("width = 15, p = 0.5, %(cls) clusters:")
showCluster.call()
 
System.print("\np = 0.5, iter = 5:")
w = 1 << 2
while (w <= (1 << 13)) {
    var t = tests.call(5, 0.5)
    Fmt.print("$5d $9.6f", w, t)
    w = w << 1
}