模拟求逾渗理论的阈值 | 寒山烟雨
现在的位置: 首页 > 小试身手 > 正文

模拟求逾渗理论的阈值

2014年08月12日 小试身手 ⁄ 共 1948字 ⁄ 字号 评论 3 条 ⁄ 阅读 4,031 views 次

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

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

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

2,依次走遍每一行。

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

4,依次走完每一行。

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

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

2,对撒的点进行分组

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

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

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

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

 Fortran | 
 
 copy code |
?

001
program main
002
    implicit none
003
    integer,parameter   :: s=4,e=256
004
    integer :: i,N,j,k
005
    real(8)    :: cal,avg
006
    open(10,file='data.dat')
007
    open(100,file='pc.dat')
008
    k=0
009
    do i=2,7
010
        N=2**i
011
        avg=0.0d0
012
        do j=1,10
013
            avg=avg+cal(N)
014
            write(100,*) k,cal(N)
015
            k=k+1
016
            print*,N,j
017
        enddo
018
        write(10,*)N,avg/10.d0 
019
    enddo
020
end program main
021
 
022
real(8) function cal(l)
023
    integer :: l,switch,i,j
024
    integer :: Space(l,l),Group(l,l)
025
    real(8)    :: pc,dpc
026
    pc=0.5d0
027
    dpc=0.005d0
028
    !赋初值
029
    do i=1,l
030
        do j=1,l
031
            Space(i,j)=0
032
            Group(i,j)=0
033
        enddo
034
    enddo
035
    call D(Space,Group,l,pc)
036
 
037
    do while(switch(Group,l)==0)  !判断是不是贯通,不是的话pc+dpc
038
        pc=pc+dpc
039
        call D(Space,Group,l,pc)
040
    enddo
041
    !cal=1.0*sum(SPace)/l**2.0
042
    cal=pc
043
end
044
 
045
subroutine D(Space,Group,l,pc)
046
    integer ::l,i,j,g,g1,g2,nn,xi,yi
047
    integer ::Space(l,l),Group(l,l)
048
    real(8)    ::pc,x,y ,r
049
    call RANDOM_SEED()
050
 
051
    nn=int(1.0*l**2.0*pc)+1
052
     i=0
053
     do
054
        call random_number(x)
055
        xi=int(l*x)+1
056
        call random_number(y)
057
        yi=int(l*y)+1
058
        if (space(xi,yi)==0)then
059
            Space(xi,yi)=1
060
        endif
061
       if (sum(space)>nn) exit
062
   enddo
063
 
064
    g=0
065
    !第一次分组
066
    do i=1,l
067
        do j=1,l
068
            if (i==1 .and. j==1 .and. space(i,j)==1) then
069
                g=g+1
070
                group(i,j)=1
071
            endif
072
 
073
            if (i==1 .and. j/=1 .and. space(i,j)==1) then
074
                if (space(i,j-1)==1) then
075
                    group(i,j)=group(i,j-1)
076
                else
077
                    g=g+1
078
                    group(i,j)=g
079
                endif
080
            endif
081
 
082
            if (i/=1 .and. space(i,j)==1) then
083
                if (j==1) then
084
                    if (group(i-1,j)==0) then
085
                        g=g+1
086
                        group(i,j)=g
087
                    else
088
                        group(i,j)=group(i,j)
089
                    endif
090
                else
091
                    if (group(i-1,j)==0 .and. group(i,j-1)==0)then
092
                        g=g+1
093
                        group(i,j)=g
094
                    else
095
                        if (group(i-1,j)*group(i,j-1)/=0) then
096
                            group(i,j)=min(group(i-1,j),group(i,j-1))
097
 
098
                        else
099
                            group(i,j)=group(i-1,j)+group(i,j-1)
100
                        endif
101
                    endif
102
                endif
103
            endif
104
        enddo
105
    enddo
106
    !第二次修正
107
    do i=1,l
108
        do j=1,l
109
            if (i==l .and. j==l) exit
110
            if (i==l .and. group(i,j)/=0) then
111
                if (group(i,j+1)/=0) then
112
                    if (group(i,j)<group(i,j+1)) then
113
                        g1=group(i,j)
114
                        g2=group(i,j+1)
115
                        call Change(Group,l,g1,g2)
116
                    else
117
                        g1=group(i,j+1)
118
                        g2=group(i,j)
119
                        call Change(Group,l,g1,g2)
120
                    endif
121
                endif
122
            endif
123
 
124
            if (j==l .and. group(i,j)/=0) then
125
                if (group(i+1,j)/=0) then
126
                    if (group(i,j)<group(i+1,j))then
127
                        g1=group(i,j)
128
                        g2=group(i+1,j)
129
                        call Change(Group,l,g1,g2)
130
                    else
131
                        g1=group(i+1,j)
132
                        g2=group(i,j)
133
                        call Change(Group,l,g1,g2)
134
                    endif
135
                endif
136
            endif
137
 
138
            if (i/=l .and. j/=l .and. group(i,j)/=0) then
139
                if (group(i+1,j)*group(i,j+1)/=0)then
140
                    g1=min(group(i,j+1),group(i,j))
141
                    g1=min(g1,group(i+1,j))
142
                    g2=max(group(i,j+1),group(i,j))
143
                    g2=max(g2,group(i+1,j))
144
                    call Change(Group,l,g1,g2)
145
                else
146
                    if (group(i+1,j)/=0) then
147
                        if (group(i,j)<group(i+1,j)) then
148
                            g1=group(i,j)
149
                            g2=group(i+1,j)
150
                            call Change(Group,l,g1,g2)
151
                        else
152
                            g1=group(i+1,j)
153
                            g2=group(i,j)
154
                            call Change(Group,l,g1,g2)
155
                        endif
156
                    endif
157
                    if (group(i,j+1)/=0) then
158
                        if (group(i,j)<group(i,j+1))then
159
                            g1=group(i,j)
160
                            g2=group(i,j+1)
161
                            call Change(Group,l,g1,g2)
162
                        else
163
                            g1=group(i,j+1)
164
                            g2=group(i,j)
165
                            call Change(Group,l,g1,g2)
166
                        endif
167
                    endif
168
                endif
169
            endif
170
        enddo
171
    enddo
172
end
173
 
174
subroutine Change(Group,l,g1,g2)
175
    integer :: l,group(l,l),g1,g2
176
    integer :: i,j
177
    do i=1,l
178
        do j=1,l
179
            if (group(i,j)==g2) group(i,j)=g1
180
        enddo
181
    enddo
182
end
183
 
184
integer function switch(Group,l)
185
    integer :: l
186
    integer :: Group(l,l),flag,i
187
    flag=0
188
    do i=1,l
189
        if (Group(1,i)/=0) then
190
            do j=1,l
191
                if (Group(1,i)==Group(l,j)) then
192
                    flag=1
193
                    exit
194
                endif
195
            enddo
196
        endif
197
        if (flag==1)    exit
198
    enddo
199
    switch=flag
200
end

 

1

目前有 3 条留言    访客:2 条, 博主:1 条

  1. 玻璃钢电缆支架 2015年04月28日 下午2:39  @回复  Δ-49楼 回复
    TheWorld Browser TheWorld Browser Windows Windows

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

  2. Zhengheng Li 2014年08月12日 下午10:42  @回复  Δ-48楼 回复
    Firefox Firefox Windows Windows

    代码200行!


    • 管理员
      hsyyf 2014年08月12日 下午10:56  @回复  ∇地下1层 回复
      Firefox Firefox Windows Windows

      这个程序几乎是我所写过的算法最复杂的一个了。。。

给我留言

留言无头像?


×