
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.
$fill = 'x';
$D{$_} = $i++ for qw<DeadEnd Up Right Down Left>;
 
sub deq { defined $_[0] && $_[0] eq $_[1] }
 
sub perctest {
    my($grid) = @_;
    generate($grid);
    my $block = 1;
    for my $y (0..$grid-1) {
        for my $x (0..$grid-1) {
            fill($x, $y, $block++) if $perc[$y][$x] eq $fill
        }
    }
    ($block - 1) / $grid**2;
}
 
sub generate {
    my($grid) = @_;
    for my $y (0..$grid-1) {
        for my $x (0..$grid-1) {
            $perc[$y][$x] = rand() < .5 ? '.' : $fill;
        }
    }
}
 
sub fill {
    my($x, $y, $block) = @_;
    $perc[$y][$x] = $block;
    my @stack;
    while (1) {
        if (my $dir = direction( $x, $y )) {
            push @stack, [$x, $y];
            ($x,$y) = move($dir, $x, $y, $block)
        } else {
            return unless @stack;
            ($x,$y) = @{pop @stack};
        }
    }
}
 
sub direction {
    my($x, $y) = @_;
    return $D{Down}  if deq($perc[$y+1][$x  ], $fill);
    return $D{Left}  if deq($perc[$y  ][$x-1], $fill);
    return $D{Right} if deq($perc[$y  ][$x+1], $fill);
    return $D{Up}    if deq($perc[$y-1][$x  ], $fill);
    return $D{DeadEnd};
}
 
sub move {
    my($dir,$x,$y,$block) = @_;
    $perc[--$y][   $x] = $block if $dir == $D{Up};
    $perc[++$y][   $x] = $block if $dir == $D{Down};
    $perc[  $y][ --$x] = $block if $dir == $D{Left};
    $perc[  $y][ ++$x] = $block if $dir == $D{Right};
    ($x, $y)
}
 
my $K = perctest(15);
for my $row (@perc) {
    printf "%3s", $_ for @$row;
    print "\n";
}
printf  "𝘱 = 0.5, 𝘕 = 15, 𝘒 = %.4f\n\n", $K;
 
$trials = 5;
for $N (10, 30, 100, 300, 1000) {
    my $total = 0;
    $total += perctest($N) for 1..$trials;
    printf "𝘱 = 0.5, trials = $trials, 𝘕 = %4d, 𝘒 = %.4f\n", $N, $total / $trials;
}