
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 std.stdio, std.algorithm, std.random, std.math, std.array,
       std.range, std.ascii;
 
alias Cell = ubyte;
alias Grid = Cell[][];
enum Cell notClustered = 1; // Filled cell, but not in a cluster.
 
Grid initialize(Grid grid, in double prob, ref Xorshift rng) nothrow {
    foreach (row; grid)
        foreach (ref cell; row)
            cell = Cell(rng.uniform01 < prob);
    return grid;
}
 
void show(in Grid grid) {
    immutable static cell2char = " #" ~ letters;
    writeln('+', "-".replicate(grid.length), '+');
    foreach (row; grid) {
        write('|');
        row.map!(c => c < cell2char.length ? cell2char[c] : '@').write;
        writeln('|');
    }
    writeln('+', "-".replicate(grid.length), '+');
}
 
size_t countClusters(bool justCount=false)(Grid grid)
pure nothrow @safe @nogc {
    immutable side = grid.length;
    static if (justCount)
        enum Cell clusterID = 2;
    else
        Cell clusterID = 1;
 
    void walk(in size_t r, in size_t c) nothrow @safe @nogc {
        grid[r][c] = clusterID; // Fill grid.
 
        if (r < side - 1 && grid[r + 1][c] == notClustered) // Down.
            walk(r + 1, c);
        if (c < side - 1 && grid[r][c + 1] == notClustered) // Right.
            walk(r, c + 1);
        if (c > 0 && grid[r][c - 1] == notClustered) // Left.
            walk(r, c - 1);
        if (r > 0 && grid[r - 1][c] == notClustered) // Up.
            walk(r - 1, c);
    }
 
    size_t nClusters = 0;
 
    foreach (immutable r; 0 .. side)
        foreach (immutable c; 0 .. side)
            if (grid[r][c] == notClustered) {
                static if (!justCount)
                    clusterID++;
                nClusters++;
                walk(r, c);
            }
    return nClusters;
}
 
double clusterDensity(Grid grid, in double prob, ref Xorshift rng) {
    return grid.initialize(prob, rng).countClusters!true /
           double(grid.length ^^ 2);
}
 
void showDemo(in size_t side, in double prob, ref Xorshift rng) {
    auto grid = new Grid(side, side);
    grid.initialize(prob, rng);
    writefln("Found %d clusters in this %d by %d grid:\n",
             grid.countClusters, side, side);
    grid.show;
}
 
void main() {
    immutable prob = 0.5;
    immutable nIters = 5;
    auto rng = Xorshift(unpredictableSeed);
 
    showDemo(15, prob, rng);
    writeln;
    foreach (immutable i; iota(4, 14, 2)) {
        immutable side = 2 ^^ i;
        auto grid = new Grid(side, side);
        immutable density = nIters
                            .iota
                            .map!(_ => grid.clusterDensity(prob, rng))
                            .sum / nIters;
        writefln("n_iters=%3d, p=%4.2f, n=%5d, sim=%7.8f",
                 nIters, prob, side, density);
    }
}