周期境界条件を適用した量子ウォークの確率が保存されない。

Question

問題点、欲しいもの

量子ウォークという、いわゆるランダムウォークの量子版という物を勉強しています。
全確率が保存するように、以下の周期境界条件をつかって、時間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は44の行列で、ψは[a,b,c,d]のような4状態で表すとします。
このψから確率pは**p(t,x,y) = ψ(t,x,y).ψ(t,x,y) **と算出できます。(*は複素共役)。
以上を踏まえて、以下にコードと結果をリポジっトします。

コード内容

import numpy as np
import matplotlib.pyplot as plt
import math
import itertools

n = 3
itr = 3
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: #x座標が0のときで場合分け
                if i[1]==0: #y =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:#x=2nの時で場合分け
                if i[1] == 0:　#y=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])])
            else: #xがそれ以外のときのyで場合分け
                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

とｔ＞２で既に確率が保存されていない状況が生まれてしまい、コードを見直しても原因がわかりません。
怪しいと思ったこと何でも言ってください。
よろしくお願い致します。

Answer 1

t = 1 でループを回った直後を考えてみてください。
phi_map = next_phi_map の直後なので、
phi_map と next_phi_map は同じものを指していますよね。

phi_map は時刻 t の、next_phi_map は時刻 t+1 のときの ψ の値を保持していて欲しいのですから、これは変です。
この状態で next_phi_map を書きかえているので、ψ_{t+1} だけを変更しているつもりが、
「過去の値」である ψ_{t} (phi_map) も書きかえられてしまっています。
……というのがうまくいかなかった理由です。

対策としては、ループ内で next_phi_map を毎回新たに作成すればよいです。
分かってしまえばよくある凡ミスですが、コードが入り組んでいて発見が難しかったのだと思います。
抽象化できる処理はできるだけ抽象化するのがオススメです。

// ちなみに、ψ は psi です。phi は φ のことですね。

# coding: utf-8

import numpy as np

# 一辺の長さは 2n+1
n = 3

# 繰り返し数
itr = 3

# 初期位置
(x0, y0) = (0, 0)

# =============================================
# ψ と p の初期化
# =============================================

# 遷移の重みを表す行列
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]] #↑

# ψ_{t}(x,y) の初期化
phi_map = np.zeros((2*n+1, 2*n+1,4), dtype="complex")

# 時刻 t における ψ_{t} の確率分布 p_{t}(x, y)
p_map = np.zeros([2*n+1,2*n+1])

# 初期条件の設定
phi_map[x0,y0]= np.array([1,0,0,0])
p_map[x0,y0] = 1.0

# =============================================
# 便利関数
# =============================================

# (x, y) where 0 <= x <= 2n, 0 <= y <= 2n を満たす (x,y) の組を生成するジェネレータを返す
def Lattice(width = 2*n+1, height = 2*n+1): return ( (x, y) for x in range(width) for y in range(height) )

# v ≡ x mod (2n+1) を満たす最小の非負整数 x を返す
def mod(v): return v % (2*n + 1)

# ψ の第 n 項を計算する
def T(v, x, y): return np.dot(v, phi_map[mod(x), mod(y)])
def T1(x, y): return T(P, x, y)
def T2(x, y): return T(Q, x, y)
def T3(x, y): return T(R, x, y)
def T4(x, y): return T(S, x, y)

# ψ_{t} を使って座標 (x, y) の ψ_{t+1} を計算する
def calc_phi(x, y):
    return np.array([ T1(x-1, y) + T2(x+1, y) + T3(x, y+1) + T4(x, y-1) ])

# ψ_{t} から p_{t} を計算する
def calc_p(phi):
    return np.real(np.vdot(phi, phi))

# 時刻 t とそのときの確率分布を表示する
def show(t, p):
    print(t)
    print(np.real(p), p.sum())

# =============================================
# main
# =============================================
np.set_printoptions(linewidth=256)

show(0, p_map)
for t in range(itr):
    # v こいつは毎回初期化が必要
    # ψ_{t+1}(x,y) の初期化
    next_phi_map = np.zeros((2*n+1, 2*n+1, 4), dtype="complex")
    for (x,y) in Lattice():
        next_phi = calc_phi(x, y)
        next_phi_map[x,y] = next_phi
        p_map[x,y] = calc_p(next_phi)
    
    phi_map = next_phi_map
    show(t, p_map)

Comments

[Fallout_18]

初めに、懇切丁寧なご指摘、お時間を割いて頂き、誠に感謝いたします。
ただ、まだピンとこないのですが、
僕としては、t+1のnext_phi_mapが計算し終わって、t+2に行くときにphi_mapにt+1を代入してあげたら次のループでt+2のnext_phi_mapが来ると想定しています（つまり書き換えられて良いと思ったのですが。。。）
その事を考え、
phi_map = np.copy(next_phi_map)としてみたら確率が保存されたのですが、つまりは、next_phi_mapのデータをphi_mapに移行するときに、完全に移行されていなかったということだと感じました（浅いコピー、深いコピーが関係している（？））。そのやりかたでも良いのでしょうか（自分では良いと思います）？

ただ、どちらにせよ、上記のpute様のコードの書き方は非常に為になり、とても勉強になります。貪欲に吸収させてい頂きます！！

[Fallout_18]

それと、psiとphiはいつもどっちだっけとなってしまいます。申し訳ないです。。

[Fallout_18]

next_phi_mapの初期化に関して、理解できました！

[pute]

よかったです。「=」が何をするのか、というのは割とつまづきどころかと思います。

ちなみに、過去の質問の Sierpinski carpet も今回と同じではないでしょうか。
https://teratail.com/questions/136122
確認はしていませんが、もしまだ未解決でしたら試してみてください。
それでは。

Answer 2

こんにちわ。

もうBA出ていますが、補足させてください。

コードの最後のあたり、

next_phi_mapをphi_mapにコピーしようとしている箇所ですが、

    phi_map = next_phi_map

これだと値をコピーしたことにならないので、

次のようにすればよいです。

    phi_map = copy.copy(next_phi_map)

修正後のコードを実行させたところ、以下の結果になりました。

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
3 [[0.       0.078125 0.       0.015625 0.015625 0.       0.078125]
 [0.015625 0.       0.078125 0.       0.       0.078125 0.      ]
 [0.       0.078125 0.       0.       0.       0.       0.078125]
 [0.015625 0.       0.       0.       0.       0.       0.      ]
 [0.015625 0.       0.       0.       0.       0.       0.      ]
 [0.       0.015625 0.       0.       0.       0.       0.015625]
 [0.390625 0.       0.015625 0.       0.       0.015625 0.      ]] 1.0

なお、修正コードはプログラムの頭に import copy を追加しています。

Comments

[Fallout_18]

参考にさせて頂きます！
ちなみになんですが、phi_map =np.copy(next_phi_map)でも解決しました！

[srsnsts]

それでもいいと思います。