|
Two dimensional percolation
Authors: Nikolay Aleksandrovich Khokhlov, Kyaw Aunghein and Lin Ko Ko. July 16, 2017. Komsomolsk-na-Amure. http://physmath.tech.
Abstract: This is a first percolation package realised as a functional one. We consider a two-dimensional percolation with square cells and square grid.
We use three colors: white - dielectric cell, red - conductor cell, blue - conducting cluster (from the bottom to the top).
|
>
|
restart:with(plottools): with(plots):
|
We define a uniform distribution probability function:
|
>
|
randomize(): rnd:= rand(0.0 .. 1.0):
|
This function generates a random square percolation pattern. N is a grid size and p is the probability of being a conducter cell.
|
>
|
MY_C_initialize:=proc(N,p) local i,j,r,MY_C: MY_C:=[seq([seq(white,j=1..N)],i=1..N)]:for i from 1 to N do
for j from 1 to N do
r:=rnd():
if r < p then MY_C[i,j]:=red: else MY_C[i,j]:=white: fi:
od: od: return(MY_C): end proc:
|
This function assigns numbers to cells acorrding to their colors
|
>
|
numberOFcolor:=proc(MY_C) local MY_N,N,i,j: N:=nops(op(1,MY_C)); MY_N:=MY_C: for i from 1 to N do
for j from 1 to N do
if MY_C[i,j]=red then MY_N[i,j]:=1: else MY_N[i,j]:=0: fi:
if MY_C[i,j]=blue then MY_N[i,j]:=5: fi:
od: od: return(MY_N):end proc:
|
This function assigns colors to cells acorrding to their numbers
|
>
|
colorOFnumber:=proc(MY_N) local MY_C,i,j,N: MY_C:=MY_N: N:=nops(op(1,MY_C)); for i from 1 to N do
for j from 1 to N do
if MY_N[i,j]=1 then MY_C[i,j]:=red: else MY_C[i,j]:=white: fi:
if MY_N[i,j]=5 then MY_C[i,j]:=blue: fi:
od: od: return(MY_C):end proc:
|
This is a main function that searches for conducting cluster. The algorithm is as follows. First: we assign number 5 to all the botton row cells that are conducters. Second: we go from left to right and from bottom to the top and assign a number 5 to all the cells that are conducters and have a neighbor numbered by 5. You can see that condition check is different for the inner part and for the outer colomns
|
>
|
Change_color:=proc() global MY_N,MY_C: local N,NN,NN2,i,j,MY_N_OLD,MY_F: N:=nops(op(1,MY_C)); MY_N_OLD:=MY_N:
# We compute sum of all numbers. It will be needed to check growth of the conductor cluster
NN:=add(add(MY_N[i,j],i=1..N),j=1..N):
# All bottom conductor cells are presumed to be potential roots of the conducter cluster
for i from 1 to N do
if MY_N_OLD[i,1]=1 then MY_N[i,1]:=5: fi:
od:
# This is a double circle where we go from left to right and from bottom to top
for j from 2 to N-1 do
for i from 2 to N-1 do
# This is a condition check for the inner part of the grid
if (MY_N_OLD[i,j-1]+ MY_N_OLD[i,j+1]+MY_N_OLD[i-1,j]+MY_N_OLD[i+1,j] >= 5 and MY_N_OLD[i,j]=1)
then MY_N[i,j]:=5: fi:
od:
# Next two lines are condition checks for the right and left colomns and for the top row
if (MY_N_OLD[1,j-1]+ MY_N_OLD[1,j+1]+MY_N_OLD[2,j] >= 5 and MY_N_OLD[1,j]=1) then MY_N[1,j]:=5: fi:
if (MY_N_OLD[N,j-1]+ MY_N_OLD[N,j+1]+MY_N_OLD[N-1,j] >= 5 and MY_N_OLD[N,j]=1) then MY_N[N,j]:=5: fi:
od:
# This is the end of the double cicle
# This circle checks condition for the top row
for i from 2 to N-1 do
if (MY_N_OLD[i-1,N]+ MY_N_OLD[i+1,N]+MY_N_OLD[i,N-1] >= 5 and MY_N_OLD[i,N]=1) then MY_N[i,N]:=5: fi: od:
# Next two lines are condition checks for the two top corners
if (MY_N_OLD[2,N]+ MY_N_OLD[1,N-1] >= 5 and MY_N_OLD[1,N]=1) then MY_N[1,N]:=5: fi:
if (MY_N_OLD[N,N-1]+ MY_N_OLD[N-1,N] >= 5 and MY_N_OLD[N,N]=1) then MY_N[N,N]:=5: fi:
# Now we check growth of the conductor cluster. If it grows then we proceed by recursive call of the present cluster.
# If the growth stops then we check if there are conductor cluster cells in the top row. If so we print "Yes", otherwise we print "NO".
# At the end we draw a picture of our grid with white cells being dielectric, red ones being conducters, and blue ones being part of the conductor cluster.
NN2:=add(add(MY_N[i,j],i=1..N),j=1..N): if NN2>NN then Change_color() else for i from 1 to N do
if MY_N_OLD[i,N]=5 then print(YES): MY_F:=colorOFnumber(MY_N): return(MY_F) fi:
od: print(NO): MY_F:=colorOFnumber(MY_N): return(MY_F) fi:
end proc:
|
Now we check our algorithm. N is the size of our square grid. p is the probability of being a conducter cell. x1 and y1 are list of cells coordinates
|
>
|
N:=17: p:=0.57: x1:=[seq(i-1,i=1..N+1)]:y1:=[seq(i-1,i=1..N+1)]:
|
We perform a cicle of stochatic check if number p is the threshold probability.
The Ncheck is a number of random checks, you should increase it (as well as N) with a precaution - the computation may be long for big numbers.
You must know that threshold probability is defined only in the N-> infinity limit. So you must increase the size of the grid to find the threshold.
The line that draws the initial grid is commented, this line draws grid without the blue colored conductor cluster.
|
>
|
Ncheck:=2: for i from 1 to Ncheck do MY_C:=MY_C_initialize(N,p): MY_N:=numberOFcolor(MY_C):
#Sqs:=seq(seq(rectangle([x1[i],y1[j]],[x1[i+1],y1[j+1]],color=MY_C[i,j]),i=1..N),j=1..N): print(display(Sqs)):
MY_F:=Change_color(MY_N,MY_C): Sqs:=seq(seq(rectangle([x1[i],y1[j]],[x1[i+1],y1[j+1]],color=MY_F[i,j]),i=1..N),j=1..N): print(display(Sqs)); od:
|
|
