第6章 モンテカルロ法

６．１　適用と乱数分布関数例

【リスト6-1】正規分布関数とその乱数発生頻度グラフの表示
from tkinter import *
import random, math
def randNomal():#正規分布乱数発生関数
    T=-6
    for i in range(12):T+=random.random()
    return T
M=100;N=1000000;X=[0 for i in range(M*2)]#以下，グラフ表示
for k in range(N):
    RR=randNomal();ID=int((RR+6)*M/10);X[ID]+=1
root=Tk(); root.title("正規分布")
canvas=Canvas(root, width=460,height=160, bg="white");canvas.pack()
for i in range(7):
    Y0=140-i*20
    canvas.create_line(40,Y0, 440,Y0,fill="blue",width=1)
    canvas.create_text(25,Y0,text="%4.2f" % (i/100),
                       font=("ＭＳ　ゴシック",8))
for i in range(0,M+1,10):
    X0=40+i*4
    canvas.create_line(X0,20, X0,140,fill="blue",width=1)
    canvas.create_text(X0-4,150,text="%4.1f" % (i/50-1),
                       font=("ＭＳ　ゴシック",8))
X0=40;Y0=140
for i in range(10,int(M+11)):
    X1=i*4; Y1=140-X[i]/N*M*20
    canvas.create_line(X0,Y0, X1,Y1,fill="red",width=2)
    X0=X1; Y0=Y1

【リスト6-2】ポアソン分布関数とその乱数発生頻度グラフの表示
from tkinter import *
from tkinter import messagebox
import random, math, time
def poissonDisp(M):#ポアソン分布関数
    MM=math.exp(M)*random.random(); R=0
    while MM>1: MM*=random.random();R+=1
    return R
M=100;N=1000000;X=[0 for i in range(M+1)]#以下，グラフ表示
for k in range(N): ID=int(poissonDisp(50)); X[ID]+=1
root=Tk(); root.title("ポアソン分布")
canvas=Canvas(root, width=460,height=160, bg="white");canvas.pack()
for i in range(7):
    Y0=140-i*20
    canvas.create_line(40,Y0, 440,Y0,fill="blue",width=1)
    canvas.create_text(25,Y0,text="%4.2f" % (i/100), 
　　　　font=("ＭＳ　ゴシック",8))
for i in range(0,M+1,10):
    X0=40+i*4
    canvas.create_line(X0,20, X0,140,fill="blue",width=1)
    canvas.create_text(X0-4,150,text="%4.0f" % i, 
　　　　font=("ＭＳ　ゴシック",8))
X0=40;Y0=140
for i in range(1,M+1):
    X1=40+i*4; Y1=140-X[i]/N*M*20
    canvas.create_line(X0,Y0, X1,Y1,fill="red",width=2)
    X0=X1; Y0=Y1from tkinter import *
import random, math
def randNomal():

６．２　パーコレーション問題

【リスト6-3】パーコレーション問題
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
N=100; M=100
A=array(N+2,M+2);D=array(N+2,M+2)
def setCell():
    for i in range(2*N*M):
        X=int(N*random.random()+1);Y=int(M*random.random()+1)
        if A[X][Y]==0 : break
    A[X][Y]=1
def direction(i, j):
    if   j > M                  : return True
    elif j < 1 or i < 1 or i > M: return False
    elif D[i][j]!=0             : return False
    else:
        D[i][j] = 1  #以前来たことがあることを示すフラグ
        if A[i][j]==0 : return False
        R = direction(i, j + 1)
        if R: D[i][j+1]=2; return R
        R = direction(i, j - 1)
        if R: D[i][j-1]=2;return R
        R = direction(i + 1, j)
        if R: D[i+1][j]=2;return R
        R=direction(i - 1, j)
        if R: D[i-1][j]=2;return R
        return R
def judge():
    for i in range(1,N+1):
        for j in range(M+1): D[i][j]=0
    for i in range(1, N+1):
        if A[i][1]!=0:
            if direction(i,1): D[i][1]=2;return True
    return False
def percolation(root):
    for i in range(N*M+1):
        if i % 100 ==0: root.title("粒子数=%d" % i);root.update()
        setCell()
        if judge():
            root.title("粒子数=%d" % i);root.update();return i
    return 0
def dsp(canvas):
    for i in range(1,N+1):
        ii=i*2+10
        for j in range(1, M+1):
            jj=j*2+10
            if A[i][j]!=0:
                CL="#7777FF"
                if D[i][j]==2: CL="#FF0000"
                canvas.create_rectangle(jj,ii,jj+2,ii+2,fill=CL,outline="")    
#■メイン処理
root=Tk(); root.title("パーコレーション(浸透)問題")
canvas=Canvas(root, width=220,height=230, bg="white");canvas.pack()
if percolation(root)!=0: dsp(canvas); root.update()

６．３　酔歩(random walk)問題

【リスト6-4】酔っ払った人が橋を渡れることができる確率
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
XB=0; YB=0
def dsp(canvas, X,Y):
    global XB,YB
    XX=X*50+10; YY=Y*50+10
    canvas.create_line(XB, YB, XX,YY,fill="red",width=2)
    XB=XX; YB=YY
def dspInit(canvas,X,Y): global XB,YB; XB=X*50+10; YB=Y*50+10
def UnitRandom(root, canvas):#　1回試行
    W = 4; L = 8; X = 0; Y = W / 2; K = 2
    canvas.delete(ALL); PI2=2*math.pi
    canvas.create_rectangle(10, 10, L*50+10,W*50+10,fill="#FFFF77",
                            outline="#7777FF",width=2)
    dspInit(canvas, X,Y)
    while X < L:
        TH = PI2 * random.random()
        X += 0.5 * math.cos(TH); Y += 0.5 * math.sin(TH)
        K+=1; dsp(canvas, X,Y)
        if X < 0 or Y < 0 or Y > W : return False
        root.update()
    return True
def RandomWalk(root, canvas):
    N=10000;M=0; DM=100 # DM:成功したときの終了回数(画面キャプチャ用)
    for i in range(1,N+1):
        if UnitRandom(root, canvas):
            M +=1
            root.title("試行回数 %d, 成功回数 %d(%8.4f)"%(i, M, (M/i)))
            if M==DM: return M/i
        else: root.title("試行回数 %d, 成功回数 %d(%8.4f)"%(i, M, (M/i)))
    return M / N
root=Tk()
canvas=Canvas(root, width=420,height=230, bg="white");canvas.pack()
RandomWalk(root,canvas)

６．４　崩壊の影響範囲

【リスト6-5】落石／土砂流出影響範囲シミュレーション
from tkinter import *
from tkinter import messagebox
import random, math, time
def array(N1, N2=0, N3=0):
    if N2 == 0 and N3==0: return [0 for i in range(N1)]
    elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
    else: return [[[0 for i in range(N3)] for j in range(N2)]
                  for k in range(N1)]
KX=80; KY=80; MX=KX+1; MY=KY+1
Myu    =0.12# 粘性抵抗(様々な障害物の影響を代表させるので大きくする)
Grad   =1   # 重力加速度
MyuRand=0.0 # 粘性抵抗の揺らぎ(0～0.25程度:粘性抵抗が負にならないように)
DT=1        # 時間刻み
Z=array(MX,MY)
MaxB=40; nBall = -1; Ball=array(MaxB, 5)
for j in range(KY):#■初期条件の設定
    for i in range(KX):
        if i>25: Z[i][j]=0
        else: Z[i][j]=10-i*20/50
def genBall(canvas,KX,KY, CS):#岩石／土砂粒子の生成
    global nBall,MaxB, Ball
    if nBall<(MaxB-1):
        nBall+=1; ii=CS+10; jj=int(KY*CS/2+20)
        Ball[nBall][0]=canvas.create_oval(ii,jj,ii+CS,jj+CS,fill="#FFFFFF")
        Ball[nBall][1]=1           #位置
        Ball[nBall][2]=int(KY/2)
        Ball[nBall][3]=0           #速度
        Ball[nBall][4]=0
def chkC(R):                       #■R/G/B値のチェック
    if R>255: return 255
    if R<0: return 0
    return int(R)       
def dspCell(canvas, Z,KX,KY,CS):   #■地形図をカラーマップで表示
    ZMax = Z[1][1]; ZMin = Z[1][1] #最大・最小値を求める
    for i in range(KX):
        ZMax=max(ZMax, max(Z[i]))
        ZMin=min(ZMin, min(Z[i]))
    ZDX= ZMax - ZMin; B=0
    for j in range(1,KY):
        jj=j*CS+20
        for i in range(1,KX):
            ii=i*CS+10
            V=int(((Z[i][j]-ZMin)/ZDX)*765)#表示色の決定
            if   V<256: G=int(55+200*V/255); R=0
            elif V<512: V-=256; R=int(V)   ; G=255       
            else      : V-=512; R=255      ; G=255-int(V)
            C="#%02X%02X%02X" % (chkC(R),chkC(G),B)
            canvas.create_rectangle(ii,jj,ii+CS,jj+CS,
                                    fill=C, outline="")
def randNomal():
    T=-6
    for i in range(12):T+=random.random()
    return T
def simulate(canvas,CS):
    global nBall, Ball,Z,Myu,minMyu,maxMyu
    X=Ball[nBall][1]; Y=Ball[nBall][2]
    i=int(X); j=int(Y)  #運動方程式計算
    DX=Z[i+1][j]-Z[i-1][j]; DY=Z[i][j+1]-Z[i][j-1]
    GRD2=DX*DX+DY*DY; F=math.sqrt(GRD2)/math.sqrt(1+GRD2)
    TH=math.atan2(DY,DX)+math.pi*randNomal()/6;
    Fx=-F*math.cos(TH);Fy=F*math.sin(TH)
    rMyu=Myu*(1+randNomal()*MyuRand)
    #if rMyu<minMyu: minMyu=rMyu #以下3行チェック用
    #if rMyu>maxMyu: maxMyu=rMyu
    #print("%16.10f, %16.10f, %16.10f " % (rMyu, minMyu, maxMyu))
    Ball[nBall][3]=(Ball[nBall][3]+Fx*Grad*DT)*(1-rMyu)#速度計算
    Ball[nBall][4]=(Ball[nBall][4]+Fy*Grad*DT)*(1-rMyu)
    X=X+Ball[nBall][3]*DT; Ball[nBall][1]=X            #位置計算
    Y=Y+Ball[nBall][4]*DT; Ball[nBall][2]=Y
    ii=X*CS+10;jj=Y*CS+20
    #canvas.coords(Ball[nBall][0], ii,jj,ii+CS, jj+CS)
    Ball[nBall][0]=canvas.create_oval(ii,jj,ii+CS,jj+CS,fill="#FFFFFF")
    return Ball[nBall][3]*Ball[nBall][3]+Ball[nBall][4]*Ball[nBall][4]
#■メイン処理
root=Tk(); root.title("斜面上の落石シミュレーション")
canvas=Canvas(root, width=340,height=340, bg="white");canvas.pack()
N=0; minMyu=99999;maxMyu=-999999
S="粘性抵抗基準値 % 8.4f  揺らぎ %8.4f" % (Myu, MyuRand)
while True:
    dspCell(canvas, Z,KX,KY, 4)
    canvas.create_text(160,10, text=S,font=("ＭＳ　ゴシック",10), fill="blue")
    for k in range(MaxB):
        genBall(canvas,KX,KY,4); root.update()
        while True:
            epsi=simulate(canvas,4); root.update()
            if epsi<0.00001:break
            time.sleep(0.01)
    if not messagebox.askyesno('落石','繰り返しますか'): break