
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.
{-# language FlexibleContexts #-}
import Data.List
import Data.Maybe
import System.Random
import Control.Monad.State
import Text.Printf
import Data.Set (Set)
import qualified Data.Set  as S
 
type Matrix = [[Bool]]
type Cell = (Int, Int)
type Cluster = Set (Int, Int)
 
clusters :: Matrix -> [Cluster]
clusters m = unfoldr findCuster cells
  where
    cells = S.fromList [ (i,j) | (r, i) <- zip m [0..]
                               , (x, j) <- zip r [0..], x]
 
    findCuster s = do
      (p, ps) <- S.minView s
      return $ runState (expand p) ps
 
    expand p = do
      ns <- state $ extract (neigbours p)
      xs <- mapM expand $ S.elems ns
      return $ S.insert p $ mconcat xs
 
    extract s1 s2 = (s2 `S.intersection` s1, s2 S.\\ s1)
    neigbours (i,j) = S.fromList [(i-1,j),(i+1,j),(i,j-1),(i,j+1)]
    n = length m
 
showClusters :: Matrix -> String
showClusters m = unlines [ unwords [ mark (i,j)
                                   | j <- [0..n-1] ]
                         | i <- [0..n-1] ]
  where
    cls = clusters m
    n = length m
    mark c = maybe "." snd $ find (S.member c . fst) $ zip cls syms
    syms = sequence [['a'..'z'] ++ ['A'..'Z']]
------------------------------------------------------------
 
randomMatrices :: Int -> StdGen -> [Matrix]
randomMatrices n = clipBy n . clipBy n . randoms
  where
    clipBy n = unfoldr (Just . splitAt n)
 
randomMatrix n = head . randomMatrices n
 
tests :: Int -> StdGen -> [Int]
tests n = map (length . clusters) . randomMatrices n
 
task :: Int -> StdGen -> (Int, Double)
task n g = (n, result) 
  where
    result = mean $ take 10 $ map density $ tests n g
    density c = fromIntegral c / fromIntegral n**2
    mean lst = sum lst / genericLength lst
 
main = newStdGen >>= mapM_ (uncurry (printf "%d\t%.5f\n")) . res
  where
    res = mapM task [10,50,100,500]