0

素朴な疑問で、以下の通れない道があった場合の対処法で良い案がある方がいらっしゃいましたら、ぜひご教授いただければ幸いなのですが、、、。
格子上
図の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に保存しません。
この通れない道の上手い対処の仕方がわかる方がいらっしゃいましたら、ぜひご教授下さい。宜しくお願い致します。

5
  • 全体の確率って、何に対して、何の確率を求めているのですか? コードから意図を推測するのは難しいです。
    – Fumu 7
    Commented 2018年10月9日 8:15
  • 説明足らずで申し訳ないです。この場合ですと、毎秒常にこの格子上の全ての確率の和が1になっていることです
    – 0018
    Commented 2018年10月9日 23:25
  • 例えば下に通れない格子があった場合のCは、どういう値になりますか。 Commented 2018年10月9日 23:47
  • 僕が知っている限りですと、c=[[-1,2,2],[2,-1,2],[2,2,-1]]/3であるのですが、今は各格子点に4状態が置いてあるので計算が合わないです。
    – 0018
    Commented 2018年10月9日 23:50
  • じゃあ、3状態で良いじゃないかといわれればそうなのですが、これよりも複雑なステージになった場合を想定して、今は4状態の場合で考えています。(すいません勝手ながら、これから用事で返信が今日はできません)
    – 0018
    Commented 2018年10月9日 23:53

1 件の回答 1

0

すいません。自己解決しました!!!
各方向にそのまま重みを足せばよかっただけでした。

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。