
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.
package main
 
import (
    "fmt"
    "math/rand"
    "time"
)
 
var (
    n_range = []int{4, 64, 256, 1024, 4096}
    M       = 15
    N       = 15
)
 
const (
    p             = .5
    t             = 5
    NOT_CLUSTERED = 1
    cell2char     = " #abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"
)
 
func newgrid(n int, p float64) [][]int {
    g := make([][]int, n)
    for y := range g {
        gy := make([]int, n)
        for x := range gy {
            if rand.Float64() < p {
                gy[x] = 1
            }
        }
        g[y] = gy
    }
    return g
}
 
func pgrid(cell [][]int) {
    for n := 0; n < N; n++ {
        fmt.Print(n%10, ") ")
        for m := 0; m < M; m++ {
            fmt.Printf(" %c", cell2char[cell[n][m]])
        }
        fmt.Println()
    }
}
 
func cluster_density(n int, p float64) float64 {
    cc := clustercount(newgrid(n, p))
    return float64(cc) / float64(n) / float64(n)
}
 
func clustercount(cell [][]int) int {
    walk_index := 1
    for n := 0; n < N; n++ {
        for m := 0; m < M; m++ {
            if cell[n][m] == NOT_CLUSTERED {
                walk_index++
                walk_maze(m, n, cell, walk_index)
            }
        }
    }
    return walk_index - 1
}
 
func walk_maze(m, n int, cell [][]int, indx int) {
    cell[n][m] = indx
    if n < N-1 && cell[n+1][m] == NOT_CLUSTERED {
        walk_maze(m, n+1, cell, indx)
    }
    if m < M-1 && cell[n][m+1] == NOT_CLUSTERED {
        walk_maze(m+1, n, cell, indx)
    }
    if m > 0 && cell[n][m-1] == NOT_CLUSTERED {
        walk_maze(m-1, n, cell, indx)
    }
    if n > 0 && cell[n-1][m] == NOT_CLUSTERED {
        walk_maze(m, n-1, cell, indx)
    }
}
 
func main() {
    rand.Seed(time.Now().Unix())
    cell := newgrid(N, .5)
    fmt.Printf("Found %d clusters in this %d by %d grid\n\n",
        clustercount(cell), N, N)
    pgrid(cell)
    fmt.Println()
 
    for _, n := range n_range {
        M = n
        N = n
        sum := 0.
        for i := 0; i < t; i++ {
            sum += cluster_density(n, p)
        }
        sim := sum / float64(t)
        fmt.Printf("t=%3d p=%4.2f n=%5d sim=%7.5f\n", t, p, n, sim)
    }
}