3

量子ウォークという、いわゆるランダムウォークの量子版という物を勉強しています。
全確率が保存するように、以下の周期境界条件をつかって、時間t=0で(0,0)にだけ確率1をもってブルーの2n*2nの正方格子上に時間毎に1移動していくようにコードを書いたのですが、t>2秒以降の時間発展で確率が保存されず、どんどん減っていってしまい原因が全く分からず、かれこれ1カ月経過しようてしており、お力を借りたく質問させて頂きました。(バグを誘発している書き方・・・?)

図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

と全体の確率が1になるはずなのですが、減っていってしまいわかりません。
どなた様か、ご指導お願いいただけないでしょうか..宜しくお願い致します。

7
  • 「elif i == (12,12):」の行に違和感を覚えました。 ここは12のような定数ではなく、nの関数になるのではないでしょうか?
    – Fumu 7
    2018年7月30日 5:51
  • phi_map = next_phi_mapphi_map = next_phi_map.copy() に変更してみて下さい。
    – user39889
    2018年7月30日 5:53
  • python の場合、参照値が代入されるという事をすっかり忘れていました。。。回答の方は忘れて下さい(失礼な言辞はお赦し下さい、削除しました)。
    – user39889
    2018年7月30日 5:58
  • @Fumu7 elif i =(12,12)はちょっと別の案件でわざと書いてますが、書いてあっても問題はなく、答えは同じになりました。
    – 0018
    2018年7月30日 6:26
  • 3
    @Rei おそらく、Python に詳しい方が正確な解説を書いてくださるかと思いますが、今回の件では phi_map と next_phi_map が同一のインスタンス(numpy の array オブジェクト)を参照していますので、next_phi_map[i] に計算結果を代入すると、phi_map[i] も同じ値になってしまうのでした(逆でも同様です)。
    – user39889
    2018年7月30日 6:37

1 件の回答 1

4

コメントで既に解決していると思いますが、回答として書いておきます。

next_phi_map に計算結果を溜めていって、最後にそれを phi_map に代入して、次回の計算で使う、というコードだと思うのですが、

phi_map = next_phi_map

これでは、next_phi_mapphi_map が同じオブジェクトになってしまいます。

phi_map = np.copy(next_phi_map)

であるべきだと思います。

また、条件によって添字を変化させていますが、この状態では見通しが悪く、間違いに気づきにくいです。違う場所が添え字だけなら、添字の計算と行列計算は分離させた方が良いでしょう。

デバッグした結果のコードとその出力を、以下に貼り付けます。

コード:

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:
        pass
    else:
        for i in itertools.product(x_list,y_list):
            x = i[0]
            y = i[1]

            x1 = (x - 1 + 2 * n + 1) % (2 * n + 1)
            x2 = (x + 1            ) % (2 * n + 1)
            y1 = (y - 1 + 2 * n + 1) % (2 * n + 1)
            y2 = (y + 1            ) % (2 * n + 1)

            next_phi_map[i] = np.array([np.dot(P, phi_map[x1,y]) + np.dot(Q, phi_map[x2,y]) + np.dot(R, phi_map[x,y2]) + np.dot(S, phi_map[x,y1])])
            p_map[i] = np.real(np.inner(next_phi_map[i], np.conj(next_phi_map[i])))
        phi_map = np.copy(next_phi_map)
    print(t,np.real(p_map),p_map.sum())

結果:

luna:~ % python3 test.py
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.25   0.     0.0625 0.     0.     0.0625 0.    ]
 [0.     0.125  0.     0.     0.     0.     0.125 ]
 [0.0625 0.     0.     0.     0.     0.     0.    ]
 [0.     0.     0.     0.     0.     0.     0.    ]
 [0.     0.     0.     0.     0.     0.     0.    ]
 [0.0625 0.     0.     0.     0.     0.     0.    ]
 [0.     0.125  0.     0.     0.     0.     0.125 ]] 1.0
luna:~ % 
1
  • 最後のx1 = (x - 1 + 2 * n + 1) % (2 * n + 1) x2 = (x + 1 ) % (2 * n + 1) y1 = (y - 1 + 2 * n + 1) % (2 * n + 1) y2 = (y + 1 ) % (2 * n + 1)の部分は凄い目から鱗です。剰余を使って周期境界条件を表現する考えは、僕には全く思いつかなかったです。。。本当にお時間を割いて頂き、ありがとうございます。頑張ります。
    – 0018
    2018年7月31日 3:21

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

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