
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 



v


{\displaystyle v}

 be a vector of 



n


{\displaystyle n}

 values of either 1 or 0 where the probability of any
value being 1 is 



p


{\displaystyle p}

; the probability of a value being 0 is therefore 



1
−
p


{\displaystyle 1-p}

.
Define a run of 1s as being a group of consecutive 1s in the vector bounded
either by the limits of the vector or by a 0. Let the number of such runs in a given
vector of length 



n


{\displaystyle n}

 be 




R

n




{\displaystyle R_{n}}

.

For example, the following vector has 




R

10


=
3


{\displaystyle R_{10}=3}



[1 1 0 0 0 1 0 1 1 1]
 ^^^       ^   ^^^^^

Percolation theory states that





K
(
p
)
=

lim

n
→
∞



R

n



/

n
=
p
(
1
−
p
)


{\displaystyle K(p)=\lim _{n\to \infty }R_{n}/n=p(1-p)}


Task
Any calculation of 




R

n



/

n


{\displaystyle R_{n}/n}

 for finite 



n


{\displaystyle n}

 is subject to randomness so should be
computed as the average of 



t


{\displaystyle t}

 runs, where 



t
≥
100


{\displaystyle t\geq 100}

.

For values of 



p


{\displaystyle p}

 of 0.1, 0.3, 0.5, 0.7, and 0.9, show the effect of varying 



n


{\displaystyle n}


on the accuracy of simulated 



K
(
p
)


{\displaystyle K(p)}

.

Show your output here.

See also
 s-Run on Wolfram mathworld.
#include <algorithm>
#include <random>
#include <vector>
#include <iostream>
#include <numeric>
#include <iomanip>
using VecIt = std::vector<int>::const_iterator ;
 
//creates vector of length n, based on probability p for 1
std::vector<int> createVector( int n, double p ) {
   std::vector<int> result( n ) ;
   std::random_device rd ;
   std::mt19937 gen( rd( ) ) ;
   std::uniform_real_distribution<> dis( 0 , 1 ) ;
   for ( int i = 0 ; i < n ; i++ ) {
      double number = dis( gen ) ;
      if ( number <= p ) 
	 result[ i ] = 1 ;
      else 
	 result[ i ] = 0 ;
   }
   return result ;
}
 
//find number of 1 runs in the vector
int find_Runs( const std::vector<int> & numberVector ) {
   int runs = 0 ;
   VecIt found = numberVector.begin( ) ;
   while ( ( found = std::find( found , numberVector.end( ) , 1 ) ) 
	 != numberVector.end( ) ) {
      runs++ ;
      while ( found != numberVector.end( ) && ( *found == 1 ) ) 
	 std::advance( found , 1 ) ;
      if ( found == numberVector.end( ) ) 
	 break ;
   }
   return runs ;
}
 
int main( ) {
   std::cout << "t = 100\n" ;
   std::vector<double> p_values { 0.1 , 0.3 , 0.5 , 0.7 , 0.9 } ;
   for ( double p : p_values ) { 
      std::cout << "p = " << p << " , K(p) = " << p * ( 1 - p ) << std::endl ;
      for ( int n = 10 ; n < 100000 ; n *= 10 ) {
	 std::vector<double> runsFound ;
	 for ( int i = 0 ; i < 100 ; i++ ) {
	    std::vector<int> ones_and_zeroes = createVector( n , p ) ;
	    runsFound.push_back( find_Runs( ones_and_zeroes ) / static_cast<double>( n ) ) ;
	 }
	 double average = std::accumulate( runsFound.begin( ) , runsFound.end( ) , 0.0 ) / runsFound.size( ) ;
	 std::cout << "  R(" << std::setw( 6 ) << std::right << n << ", p) = " << average << std::endl ;
      }
   }
   return 0 ;
}