
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.
my @perc;
my $fill = 'x';
 
enum Direction <DeadEnd Up Right Down Left>;
 
my $𝘒 = perctest(15);
.fmt("%-2s").say for @perc;
say "𝘱 = 0.5, 𝘕 = 15, 𝘒 = $𝘒\n";
 
my $trials = 5;
for 10, 30, 100, 300, 1000 -> $𝘕 {
    my $𝘒 = ( [+] perctest($𝘕) xx $trials ) / $trials;
    say "𝘱 = 0.5, trials = $trials, 𝘕 = $𝘕, 𝘒 = $𝘒";
}
 
sub infix:<deq> ( $a, $b ) { $a.defined && ($a eq $b) }
 
sub perctest ( $grid ) {
    generate $grid;
    my $block = 1;
    for ^$grid X ^$grid -> ($y, $x) {
        fill( [$x, $y], $block++ ) if @perc[$y; $x] eq $fill
    }
    ($block - 1) / $grid²;
}
 
sub generate ( $grid ) {
    @perc = ();
    @perc.push: [ ( rand < .5 ?? '.' !! $fill ) xx $grid ] for ^$grid;
}
 
sub fill ( @cur, $block ) {
    @perc[@cur[1]; @cur[0]] = $block;
    my @stack;
    my $current = @cur;
 
    loop {
        if my $dir = direction( $current ) {
            @stack.push: $current;
            $current = move $dir, $current, $block
        }
        else {
            return unless @stack;
            $current = @stack.pop
        }
    }
 
    sub direction( [$x, $y] ) {
        ( Down  if @perc[$y + 1][$x] deq $fill ) ||
        ( Left  if @perc[$y][$x - 1] deq $fill ) ||
        ( Right if @perc[$y][$x + 1] deq $fill ) ||
        ( Up    if @perc[$y - 1][$x] deq $fill ) ||
        DeadEnd
    }
 
    sub move ( $dir, @cur, $block ) {
        my ( $x, $y ) = @cur;
        given $dir {
            when Up    { @perc[--$y; $x] = $block }
            when Down  { @perc[++$y; $x] = $block }
            when Left  { @perc[$y; --$x] = $block }
            when Right { @perc[$y; ++$x] = $block }
        }
        [$x, $y]
    }
}
 