47105


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

[![図１](https://i.stack.imgur.com/uy8r2.png)](https://i.stack.imgur.com/uy8r2.png)

物理的背景を除いて、量子状態の時間発展の式は

[code]

    ψ(t+1,x,y) = P*ψ(t,x-1,y) + Q*ψ(t,x+1,y) + R*ψ(t,x,y+1) + S*ψ(t,x,y-1)
    
[/code]

* * *

と表せるとします。P~Rは4*4の行列で、ψは[a,b,c,d]のような4状態で表すとします。  
このψから確率pは`p(t,x,y) = ψ*(t,x,y).ψ(t,x,y)` と算出できます。(*は複素共役)。  
以上を踏まえて、以下にコードと結果をリポジっトします。

[code]

    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())
    
[/code]

結果は

[code]

    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
    
[/code]

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


CC BY-SA 4.0
2018-07-30T03:21:47.643
47091
2020-11-03T10:39:21.727
2020-11-03T10:39:21.727
29111
question
python python3 quantum-computing
全確率が保存されない原因がわからない


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

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

[code]

    phi_map = next_phi_map
    
[/code]

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

[code]

    phi_map = np.copy(next_phi_map)
    
[/code]

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

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

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

コード:

[code]

    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())
    
[/code]

結果:

[code]

    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:~ % 
    
[/code]


CC BY-SA 4.0
2018-07-30T14:37:05.133
47105
2018-07-30T14:37:05.133
5288
47091
answer