
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.
#include <stdio.h>
#include <stdlib.h>
 
int *map, w, ww;
 
void make_map(double p)
{
	int i, thresh = RAND_MAX * p;
	i = ww = w * w;
 
	map = realloc(map, i * sizeof(int));
	while (i--) map[i] = -(rand() < thresh);
}
 
char alpha[] = "+.ABCDEFGHIJKLMNOPQRSTUVWXYZ"
		"abcdefghijklmnopqrstuvwxyz";
#define ALEN ((int)(sizeof(alpha) - 3))
 
void show_cluster(void)
{
	int i, j, *s = map;
 
	for (i = 0; i < w; i++) {
		for (j = 0; j < w; j++, s++)
			printf(" %c", *s < ALEN ? alpha[1 + *s] : '?');
		putchar('\n');
	}
}
 
void recur(int x, int v) {
	if (x >= 0 && x < ww && map[x] == -1) {
		map[x] = v;
		recur(x - w, v);
		recur(x - 1, v);
		recur(x + 1, v);
		recur(x + w, v);
	}
}
 
int count_clusters(void)
{
	int i, cls;
 
	for (cls = i = 0; i < ww; i++) {
		if (-1 != map[i]) continue;
		recur(i, ++cls);
	}
 
	return cls;
}
 
double tests(int n, double p)
{
	int i;
	double k;
 
	for (k = i = 0; i < n; i++) {
		make_map(p);
		k += (double)count_clusters() / ww;
	}
	return k / n;
}
 
int main(void)
{
	w = 15;
	make_map(.5);
	printf("width=15, p=0.5, %d clusters:\n", count_clusters());
	show_cluster();
 
	printf("\np=0.5, iter=5:\n");
	for (w = 1<<2; w <= 1<<14; w<<=2)
		printf("%5d %9.6f\n", w, tests(5, .5));
 
	free(map);
	return 0;
}