全確率が保存されない原因がわからない
量子ウォークという、いわゆるランダムウォークの量子版という物を勉強しています。
全確率が保存するように、以下の周期境界条件をつかって、時間t=0で(0,0)にだけ確率１をもってブルーの2n*2nの正方格子上に時間毎に１移動していくようにコードを書いたのですが、t>2秒以降の時間発展で確率が保存されず、どんどん減っていってしまい原因が全く分からず、かれこれ1カ月経過しようてしており、お力を借りたく質問させて頂きました。（バグを誘発している書き方・・・？）
物理的背景を除いて、量子状態の時間発展の式は
ψ(t+1,x,y) = P*ψ(t,x-1,y) + Q*ψ(t,x+1,y) + R*ψ(t,x,y+1) + S*ψ(t,x,y-1)
と表せるとします。P~Rは4*4の行列で、ψは[a,b,c,d]のような4状態で表すとします。
このψから確率pはp(t,x,y) = ψ*(t,x,y).ψ(t,x,y) と算出できます。(*は複素共役)。
以上を踏まえて、以下にコードと結果をリポジっトします。
import numpy as np
import itertools
n = 3
N = 4*n**2
itr = 2
x_list=[i for i in range(0,2*n+1)]
y_list=[i for i in range(0,2*n+1)]
#####
P = [[-1/2, 1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0],[0,0,0,0]] #→
Q = [[0,0,0,0],[1/2, -1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0]] #←
R = [[0,0,0,0],[0,0,0,0],[1/2, 1/2, -1/2, 1/2],[0,0,0,0]] #↓
S = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[1/2, 1/2, 1/2, -1/2]] #↑
####準備
phi_map = np.zeros((2*n+1, 2*n+1,4),dtype="complex")
phi_map[0,0]= np.array([1,0,0,0])
#t+1用の状態を収納する
next_phi_map = np.zeros((2*n+1, 2*n+1, 4),dtype="complex")
p_map=np.zeros([2*n+1,2*n+1])
p_map[0,0] = np.real(np.inner(phi_map[0,0], np.conj(phi_map[0,0])))
###main
for t in range(0,itr+1):
if t == 0:
p_map
phi_map
else:
for i in itertools.product(x_list,y_list):
if i[0]==0:
if i[1]==0:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]+2*n,0]) + np.dot(Q, phi_map[i[0]+1,i[0]]) + np.dot(R, phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],i[1]+2*n])])
elif i[1] == 2*n:
next_phi_map[i] = np.array([np.dot(P, phi_map[2*n,i[1]]) + np.dot(Q, phi_map[i[0]+1, i[1]])
+ np.dot(R, phi_map[i[0], 0])
+ np.dot(S, phi_map[0,2*n-1])])
else:
next_phi_map[i] = np.array([np.dot(P, phi_map[2*n,i[1]]) + np.dot(Q, phi_map[i[0]+1,i[1]]) + np.dot(R, phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],i[1]-1])])
elif i[0] == 2*n:
if i[1] == 0:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]-1,i[1]]) + np.dot(Q, phi_map[0,i[1]]) + np.dot(R, phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],2*n])])
elif i[1] == 2*n:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]-1,i[1]]) + np.dot(Q, phi_map[0, i[1]]) + np.dot(R, phi_map[i[0], 0]) + np.dot(S, phi_map[i[0],i[1]-1])])
else:
next_phi_map[i] = np.array([np.dot(P, phi_map[2*n-1,i[1]]) + np.dot(Q, phi_map[0,i[1]]) + np.dot(R, phi_map[2*n, i[1]+1])+ np.dot(S, phi_map[2*n, i[1]-1])])
elif i == (12,12):
next_phi_map[i]= np.array([np.dot(P, phi_map[i[0]-1,i[1]]) +np.dot(Q,phi_map[i[0]+1,i[1]])+np.dot(R,phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],i[0]-1])])
else:
if i[1] == 0:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]-1,i[1]]) + np.dot(Q, phi_map[i[0]+1,i[1]]) + np.dot(R, phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],2*n])])
elif i[1] == 2*n:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]-1, i[1]]) + np.dot(Q, phi_map[i[0]+1,i[1]]) + np.dot(R, phi_map[i[0],0])
+ np.dot(S, phi_map[i[0],i[1]-1])])
else:
next_phi_map[i] = np.array([np.dot(P, phi_map[i[0]-1,i[1]]) + np.dot(Q, phi_map[i[0]+1,i[1]]) + np.dot(R, phi_map[i[0],i[1]+1])
+ np.dot(S, phi_map[i[0],i[1]-1])])
p_map[i] = np.real(np.inner(next_phi_map[i], np.conj(next_phi_map[i])))
phi_map = next_phi_map
print(t,np.real(p_map),p_map.sum())
結果は
0 [[1. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0.]] 1.0
1 [[0. 0.25 0. 0. 0. 0. 0.25]
[0.25 0. 0. 0. 0. 0. 0. ]
[0. 0. 0. 0. 0. 0. 0. ]
[0. 0. 0. 0. 0. 0. 0. ]
[0. 0. 0. 0. 0. 0. 0. ]
[0. 0. 0. 0. 0. 0. 0. ]
[0.25 0. 0. 0. 0. 0. 0. ]] 1.0
2 (省略)0.8534011840820312
3 (省略）0.7919882354326546
と全体の確率が１になるはずなのですが、減っていってしまいわかりません。
どなた様か、ご指導お願いいただけないでしょうか..宜しくお願い致します。