
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

Given an 



M
×
N


{\displaystyle M\times N}

 rectangular array of cells numbered 




c
e
l
l

[
0..
M
−
1
,
0..
N
−
1
]


{\displaystyle \mathrm {cell} [0..M-1,0..N-1]}

assume 



M


{\displaystyle M}

 is horizontal and 



N


{\displaystyle N}

 is downwards.

Assume that the probability of any cell being filled is a constant 



p


{\displaystyle p}

 where

 



0.0
≤
p
≤
1.0


{\displaystyle 0.0\leq p\leq 1.0}


The task

Simulate creating the array of cells with probability 



p


{\displaystyle p}

 and then
testing if there is a route through adjacent filled cells from any on row 



0


{\displaystyle 0}

 to any on row 



N


{\displaystyle N}

, i.e. testing for site percolation.

Given 



p


{\displaystyle p}

 repeat the percolation 



t


{\displaystyle t}

 times to estimate the proportion of times that the fluid can percolate to the bottom for any given 



p


{\displaystyle p}

.

Show how the probability of percolating through the random grid changes with 



p


{\displaystyle p}

 going from 



0.0


{\displaystyle 0.0}

 to 



1.0


{\displaystyle 1.0}

 in 



0.1


{\displaystyle 0.1}

 increments and with the number of repetitions to estimate the fraction at any given 



p


{\displaystyle p}

 as 



t
>=
100


{\displaystyle t>=100}

.

Use an 



M
=
15
,
N
=
15


{\displaystyle M=15,N=15}

 grid of cells for all cases.

Optionally depict a percolation through a cell grid graphically.

Show all output on this page.

#lang racket
(require racket/require (only-in racket/fixnum for*/fxvector))
(require (filtered-in (lambda (name) (regexp-replace #rx"unsafe-" name ""))
                      racket/unsafe/ops))
 
(define cell-empty   0)
(define cell-filled  1)
(define cell-wall    2)
(define cell-visited 3)
(define cell-exit    4)
 
(define ((percol->generator p)) (if (< (random) p) cell-filled cell-empty))
 
(define t (make-parameter 1000))
 
(define ((make-percol-grid M N) p)
  (define p->10 (percol->generator p))
  (define M+1 (fx+ 1 M))
  (define M+2 (fx+ 2 M))
  (for*/fxvector
   #:length (fx* N M+2)
   ((n (in-range N)) (m (in-range M+2)))
   (cond
     [(fx= 0 m) cell-wall]
     [(fx= m M+1) cell-wall]
     [else (p->10)])))
 
(define (cell->str c) (substring " #|+*" c (fx+ 1 c)))
 
(define ((draw-percol-grid M N) g)
  (define M+2 (fx+ M 2))
  (for ((row N))
    (for ((col (in-range M+2)))
      (define idx (fx+ (fx* M+2 row) col))
      (printf "~a" (cell->str (fxvector-ref g idx))))
    (newline)))
 
(define ((percolate-percol-grid?! M N) g)
  (define M+2 (fx+ M 2))
  (define N-1 (fx- N 1))
  (define max-idx (fx* N M+2))
  (define (inner-percolate g idx)
    (define row (fxquotient idx M+2))    
    (cond
      ((fx< idx 0) #f)
      ((fx>= idx max-idx) #f)
      ((fx= N-1 row) (fxvector-set! g idx cell-exit) #t)
      ((fx= cell-filled (fxvector-ref g idx))
       (fxvector-set! g idx cell-visited)
       (or 
        ; gravity first (thanks Mr Newton)
        (inner-percolate g (fx+ idx M+2))
        ; stick-to-the-left
        (inner-percolate g (fx- idx 1))
        (inner-percolate g (fx+ idx 1))
        ; go uphill only if we have to!
        (inner-percolate g (fx- idx M+2))))
      (else #f)))
  (for/first ((m (in-range 1 M+2)) #:when (inner-percolate g m)) g))
 
(define make-15x15-grid (make-percol-grid 15 15))
(define draw-15x15-grid (draw-percol-grid 15 15))
(define perc-15x15-grid?! (percolate-percol-grid?! 15 15))
 
(define (display-sample-percolation p)
  (printf "Percolation sample: p=~a~%" p)
  (for*/first
      ((i (in-naturals))
       (g (in-value (make-15x15-grid 0.6)))
       #:when (perc-15x15-grid?! g))
    (draw-15x15-grid g))
  (newline))
 
(display-sample-percolation 0.4)
 
(for ((p (sequence-map (curry * 1/10) (in-range 0 (add1 10)))))
  (define n-percolated-grids
    (for/sum
     ((i (in-range (t))) #:when (perc-15x15-grid?! (make-15x15-grid p))) 1))
  (define proportion-percolated (/ n-percolated-grids (t)))
  (printf "p=~a\t->\t~a~%" p (real->decimal-string proportion-percolated 4)))