2次元格子上のとある計算がおかしい

以下のコードで、確率が計算されません。
大分修正などを試みたのですが、なぜだかわかりませんでした。
【設定】_________________________________________________________________

2m(今回はm=2)の格子上を中心(2,2)スタートで移動していき、移動するときの式は以下で与えられます。
phi(x,y,t+1) = Pphi(x+1,y,t) + Qphi(x-1,y,t) + Rphi(x,y+1,t) + S*phi(x,y-1,t)
で与えられます。

以上の【設定】をコードで書いたものを以下に示します。
コードは汚いかもしれません、ご了承ください。

python
1import numpy as np
2import matplotlib.pyplot as plt
3import math
4from mpl_toolkits.mplot3d import Axes3D
5
6
7#環境設定
8n=2  #tの範囲
9m=2  #偶数(ｘ、ｙの範囲）
10
11P = [[-1/2, 1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0],[0,0,0,0]] 
12Q = [[0,0,0,0],[1/2, -1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0]]
13R = [[0,0,0,0],[0,0,0,0],[1/2, 1/2, -1/2, 1/2],[0,0,0,0]]
14S = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[1/2, 1/2, 1/2, -1/2]]
15
16t_list =[]
17
18#初期設定
19phi_map = np.zeros((2*m+1, 2*m+1,4)) #np.zeros((行,列,[]の中身の数))
20phi_map[m,m,0]= 1
21
22p_map=np.zeros([2*m+1,2*m+1])
23
24for i in range(0,2*m+1):
25    p = np.dot(phi_map[i,i], np.conj(phi_map[i,i]))
26    p_map[i,i]=p
27#print(p_map)
28
29#以下は時間発展確率計算
30for t in range(0,n+1):
31    t_list.append(t)
32    if t == 0:
33        phi_map
34        p_map
35    else:
36        next_phi_map = np.zeros((2*m+1,2*m+1, 4))　#エラーを起こさないように、以下の計算を収納するために用意
37        for x in range(0,2*m+1):
38            if x == 0:
39                for y in range(0,2*m+1):
40                    if y == 0:　
41                        next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(R, phi_map[x,y+1])])#next_phi_mapのx行y列に計算されたもの[？,？,？,？,]を代入
42                    elif y == 2*m:
43                        next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(S, phi_map[x,y-1])])
44                    else:
45                        next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(S, phi_map[x,y-1]) + np.inner(R, phi_map[x,y+1])])
46            elif x == 2*m:
47                for y in range(0,2*m+1):
48                    if y == 0:
49                        next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1])])
50                    elif y == 2*m:
51                        next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1])])
52                    else:
53                        next_phi_map[x,y] = np.array([np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1]) + np.inner(R, phi_map[x,y+1])])
54            else:
55                if y == 0:
56                    next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1])])
57                elif y == 2*m:
58                    next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(S, phi_map[x,y-1])])
59                else:
60                    next_phi_map[x,y] = np.array([np.inner(P, phi_map[x+1,y]) + np.inner(Q, phi_map[x-1,y]) + np.inner(R, phi_map[x,y+1]) + np.inner(S, phi_map[x,y-1])])
61                    p_map[x,y] = np.dot(next_phi_map[x,y], np.conj(next_phi_map[x,y])) #確率計算
62        phi_map = next_phi_map
63
64    print(t,p_map)

出力結果は

0 [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
1 [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
2 [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

となり、t=1以降の求めるべき確率の値

0 [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]
1 [[0. 0. 0. 0. 0.]
 [0. 0. 1/4. 0. 0.]
 [0. 1/4. 0. 1/4. 0.]
 [0. 0. 1/4. 0. 0.]
 [0. 0. 0. 0. 0.]]
2 [[0. 0. 1/16. 0. 0.]
 [0. 2/16. 0. 2/16. 0.]
 [1/16. 0. 1/4. 0. 1/16.]
 [0. 2/16. 0. 2/16. 0.]
 [0. 0. 1/16. 0. 0.]]

と、でません。
私的に、
next_phi_map = np.zeros((2m+1,2m+1, 4))
の部分が、計算されずにそのまま処理されていると考えています。
何かしらのエラーが起きているということなのでしょうか？
宜しくお願い致します。

行動規範の内容に同意します

回答1件

ベストアンサー

一箇所 for y in range(0,2*m+1): が抜けてました。

python
1for t in range(0,n+1):
2    # ...
3    if t == 0:
4        # ....
5    else:
6        # ....
7        for x in range(0,2*m+1):
8            if x == 0:
9                for y in range(0,2*m+1):
10                    # ...
11            elif x == 2*m:
12                for y in range(0,2*m+1):
13                    # ...
14            else:
15                if y == 0: # <-- ここに「for y in range(0,2*m+1):」が抜けております！

上記修正後、正常に動いてるっぽいです。
出力:

python
10 [[0. 0. 0. 0. 0.]
2 [0. 0. 0. 0. 0.]
3 [0. 0. 1. 0. 0.]
4 [0. 0. 0. 0. 0.]
5 [0. 0. 0. 0. 0.]]
61 [[0.   0.   0.   0.   0.  ]
7 [0.   0.   0.25 0.   0.  ]
8 [0.   0.25 0.   0.25 0.  ]
9 [0.   0.   0.25 0.   0.  ]
10 [0.   0.   0.   0.   0.  ]]
112 [[0.    0.    0.    0.    0.   ]
12 [0.    0.125 0.    0.125 0.   ]
13 [0.    0.    0.25  0.    0.   ]
14 [0.    0.125 0.    0.125 0.   ]
15 [0.    0.    0.    0.    0.   ]]