逾渗理论作为一个经典的理论,在宏观、微观领域都有极其重要的应用,如重整化群、电阻网络、相变等等。二维格子是逾渗理论最基本的一种模型,其中逾渗阈值则是其中尤为重要的一个概念。

逾渗模拟最大的难点在于团簇分组,机器不比人,寻常递归分组过于繁琐,我这用的是最笨的二次分组法。

1,从左上角开始向右走,若有格点,向左和向上看最近临点,如果有一个点格点,分为改为改组组号,二者都有格点,选择较小者。

2,依次走遍每一行。

3,走完一遍之后,再从左上角走,这次看相下和相右的格点,遇到比自己组号小的,则把自己组所有的成员的组号变成小号。

4,依次走完每一行。

研究完分组的算法,再来看逾渗阈值的算法就简单多了。

1,在方格内以某个较小的概率值开始撒点

2,对撒的点进行分组

3,判断是否上下贯通,如果贯通,则求出格点密度(概率值、阈值),如果未贯通,增加概率值,重新散点,调回第二步

4,多次循环1~3步,求出平均值

5,扩大方格,循环1~4步,统计方格和阈值关系

这个题目的计算量超大,方格边长从2~256,统计10次求平均,大概需要半个小时左右。机器性能堪忧的同学谨慎运行。

program main
implicit none
integer,parameter :: s=4,e=256
integer :: i,N,j,k
real(8) :: cal,avg
open(10,file='data.dat')
open(100,file='pc.dat')
k=0
do i=2,7
N=2**i
avg=0.0d0
do j=1,10
avg=avg+cal(N)
write(100,*) k,cal(N)
k=k+1
print*,N,j
enddo
write(10,*)N,avg/10.d0
enddo
end program main

real(8) function cal(l)
integer :: l,switch,i,j
integer :: Space(l,l),Group(l,l)
real(8) :: pc,dpc
pc=0.5d0
dpc=0.005d0
!赋初值
do i=1,l
do j=1,l
Space(i,j)=0
Group(i,j)=0
enddo
enddo
call D(Space,Group,l,pc)

do while(switch(Group,l)==0) !判断是不是贯通,不是的话pc+dpc
pc=pc+dpc
call D(Space,Group,l,pc)
enddo
!cal=1.0*sum(SPace)/l**2.0
cal=pc
end

subroutine D(Space,Group,l,pc)
integer ::l,i,j,g,g1,g2,nn,xi,yi
integer ::Space(l,l),Group(l,l)
real(8) ::pc,x,y ,r
call RANDOM_SEED()

nn=int(1.0*l**2.0*pc)+1
i=0
do
call random_number(x)
xi=int(l*x)+1
call random_number(y)
yi=int(l*y)+1
if (space(xi,yi)==0)then
Space(xi,yi)=1
endif
if (sum(space)>nn) exit
enddo

g=0
!第一次分组
do i=1,l
do j=1,l
if (i==1 .and. j==1 .and. space(i,j)==1) then
g=g+1
group(i,j)=1
endif

if (i==1 .and. j/=1 .and. space(i,j)==1) then
if (space(i,j-1)==1) then
group(i,j)=group(i,j-1)
else
g=g+1
group(i,j)=g
endif
endif

if (i/=1 .and. space(i,j)==1) then
if (j==1) then
if (group(i-1,j)==0) then
g=g+1
group(i,j)=g
else
group(i,j)=group(i,j)
endif
else
if (group(i-1,j)==0 .and. group(i,j-1)==0)then
g=g+1
group(i,j)=g
else
if (group(i-1,j)*group(i,j-1)/=0) then
group(i,j)=min(group(i-1,j),group(i,j-1))

else
group(i,j)=group(i-1,j)+group(i,j-1)
endif
endif
endif
endif
enddo
enddo
!第二次修正
do i=1,l
do j=1,l
if (i==l .and. j==l) exit
if (i==l .and. group(i,j)/=0) then
if (group(i,j+1)/=0) then
if (group(i,j)

 

作者 hsyyf

《模拟求逾渗理论的阈值》有3条评论

回复 玻璃钢电缆支架 取消回复

您的电子邮箱地址不会被公开。 必填项已用*标注