確率保存の工夫のしかた
素朴な疑問で、以下の通れない道があった場合の対処法で良い案がある方がいらっしゃいましたら、ぜひご教授いただければ幸いなのですが、、、。
図の3×3の格子上の中央が通れない場合を考慮して、全体の確率を1に保存させたいのですが、上手くいきません(外側には周期境界条件を適用しています、それがダメなのかもしれませんが)。
import itertools
import numpy as np
import matplotlib.pyplot as plt
import math
#試行回数
itr = 5
#格子の大きさ
N =3
#場合分けのために2進数表示
stage =np.array([[1,1,1],
[1,0,1],
[1,1,1]])
#格子点に対応した状態をセット
psi = np.zeros((N,N,4),dtype="float")
psi[0,0]=[1,0,0,0]
#各場所での確率収納
p_map = np.zeros((3,3),dtype="float")
#右、左、上、下に移動するときの重み(これは定義されています)
C = np.array([[-1,1,1,1],[1,-1,1,1],[1,1,-1,1],[1,1,1,-1]])/2
right_m= np.zeros((4,4));right_m[0,:] = C[0,:]
down_m = np.zeros((4,4));down_m[1,:] = C[1,:]
left_m = np.zeros((4,4));left_m[2,:] = C[2,:]
up_m = np.zeros((4,4));up_m[3,:] = C[3,:]
#x,y軸
x_list = [i for i in range(N)]
y_list = [i for i in range(N)]
#メイン計算
for t in range(itr):
if t == 0:
pass
else:
#t+1秒後のpsiを収納して、毎回初期化
next_psi = np.zeros((N,N,4),dtype="float")
for k in itertools.product(x_list,y_list):
x = k[0]
y = k[1]
####boundary condition####
x0 = (x+1) % N
x1 = (x-1 + N) % N
y0 = (y+1) % N
y1 = (y-1 + N) % N
if stage[x,y] == 0:
continue
elif stage[x0,y]==0 and stage[x,y0]==1 and stage[x,y1]==1 and stage[x1,y]==1:
next_psi[x,y] = np.copy(np.array(np.dot(right_m,psi[x1,y]) + np.dot(up_m,psi[x,y1])+np.dot(down_m,psi[x,y0])))
elif stage[x0,y]==1 and stage[x,y0]==0 and stage[x,y1]==1 and stage[x1,y]==1:
next_psi[x,y] = np.copy(np.array(np.dot(right_m,psi[x1,y]) + np.dot(up_m,psi[x,y1])+np.dot(left_m,psi[x0,y])))
elif stage[x0,y]==1 and stage[x,y0]==1 and stage[x,y1]==0 and stage[x1,y]==1:
next_psi[x,y] = np.copy(np.array(np.dot(right_m,psi[x1,y]) + np.dot(left_m,psi[x0,y])+np.dot(down_m,psi[x,y0])))
elif stage[x0,y]==1 and stage[x,y0]==1 and stage[x,y1]==1 and stage[x1,y]==0:
next_psi[x,y] = np.copy(np.array(np.dot(left_m,psi[x0,y]) + np.dot(up_m,psi[x,y1])+np.dot(down_m,psi[x,y0])))
else:
next_psi[x,y] = np.copy( np.array( np.dot(right_m,psi[x1,y]) + np.dot(up_m,psi[x,y1]) + np.dot(left_m,psi[x0,y]) + np.dot(down_m,psi[x,y0])))
psi = np.copy(next_psi)
#############
for j in itertools.product(x_list,y_list):
p_map[j] = np.inner(psi[j],psi[j])
print(t,p_map.sum())
結果が
4 0.7265625
となってしまい、1に保存しません。
この通れない道の上手い対処の仕方がわかる方がいらっしゃいましたら、ぜひご教授下さい。宜しくお願い致します。