Q&A
ごめんなさい、わかりそうでわからないです。
もう少し詳しくお願いします。
x軸上にいる物体が左右方向に動いた時の確率分布のコードが完成しました。。。？
【前提条件】__________________________________________________
ｘ座標、ｘにいるときの状態（ここではphai）は関数phai(t,x)とすると、P,Qをある定数として
phai(t,x) = Pphai(t,x+1)+Qphai(t,x-1)と表します。
(物理学的背景は除く)
あとは初期設定の数値を変更するだけで動くと思ったのですが、、問題が発生しています。
下記のx=nのとき、phai=[0,1]に変更すると、下図のグラフがちょうどx=nに線対称な、右側に山なりになるはずなのですが、実際に動かしても設定を変える前と変わりません。
一方で、θ（theta)の値を変えると、しっかりとネットでいわれているような確率分布になってくれるので、自分が書いたコードが根本的におかしいとは考えられず、解決できずにいます。
以下に、現在の進行状況を挙げます。
python
1import numpy as np 2import matplotlib.pyplot as plt 3import math 4 5#初期設定 6n=150 7theta= 3*(math.pi)/12 #θの値 8 9P = [[np.cos(theta),np.sin(theta)],[0,0]] #Shift Operator 10Q = [[0,0],[np.sin(theta),-np.cos(theta)]] #Shift Operator 11x_list=[]#x座標 12t_list=[]#time 13p_list=[]#probability 14s_list=[]#state 15 16for x in range(0,2*n+1): 17 if x == n: 18 phai = [1,0] #←このphaiの値を[0,1]にしたら、変わるはずなのですが、変わりません。 19 else: 20 phai = [0,0] 21 p = np.dot(phai,phai) 22 23 x_list.append(x) 24 s_list.append(phai) 25 p_list.append(p) 26 27 28 29for t in range(0,3*n): #以下の文は時間ごとにそれぞれの位置での確率を計算 30 t_list.append(t) 31 if t ==0: 32 s_list 33 p_list 34 else: 35 for x in range(0,2*n+1): 36 if x == 0: 37 s_list[0] = np.inner(P, s_list[1]) 38 if x == 2*n: 39 s_list[2*n] = np.inner(Q, s_list[2*n-1]) 40 else: 41 s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1]) 42 p_list[x] = np.dot(s_list[x],s_list[x]) 43 44 45 46 print(x_list,p_list) 47 #plt.bar(x_list,p_list,width=0.1,color='black') 48 plt.xlabel("x") 49 plt.ylabel("probability") 50 plt.ylim([0,0.1]) 51 plt.xlim([-n,3*n]) 52 plt.plot(x_list,p_list,color="red",linewidth=1.0) 53 plt.pause(0.001) 54 plt.cla() 55
と行うと、以下のような確率の挙動(今回は一部を貼り付け)がわかりました。
θ＝π/4
θ=5π/12
となりました。
ご指摘の程、宜しくお願いします。
python初心者ですので、お手柔らかにお願いします。。
回答3件
0
s_listの更新のタイミングが悪いことと、
確保してあるs_listのサイズが足りないことが問題だと思います。
端の場所がおかしいのも問題ですね…
wikipedia先生によると対称になるのは[1,i]/sqrt(2)の時ですが…
それが正しければ、ノルムの計算方法もまずいです。
python
1import numpy as np 2import matplotlib.pyplot as plt 3 4n = 100 5H = np.array([[1, 1],[1, -1]]) / np.sqrt(2) 6P = np.zeros((2, 2)); P[0, :] = H[0, :] 7Q = np.zeros((2, 2)); Q[1, :] = H[1, :] 8 9pos = np.zeros((2*n+1, 2), dtype=np.complex) 10pos[n, :] = np.array([1, 1j], dtype=np.complex) / np.sqrt(2) 11x = np.arange(-n, n+1) 12 13history = np.zeros((n, 2*n+1, 2), dtype=np.complex) 14history[0, :, :] = pos[:, :] 15prob = np.zeros((n, 2*n+1)) 16prob[0, :] = np.sum(np.conj(history[0, :, :]) * history[0, :, :], axis=-1)[:] 17for t in range(1, n): 18 history[t, :-1, :] += np.inner(history[t-1, 1:, :], P.T) 19 history[t, 1:, :] += np.inner(history[t-1, :-1, :], Q.T) 20 prob[t, :] = np.sum(np.conj(history[t, :, :]) * history[t, :, :], axis=-1)[:] 21 22history_ = np.zeros((n, 2*n+1)) 23history_[0, n] = 1. 24for t in range(1, n): 25 history_[t, :-1] += history_[t-1, 1:]/2. 26 history_[t, 1:] += history_[t-1, :-1]/2. 27 28def show(x, y, y_=None, nonzeros=True): 29 if nonzeros: 30 s = 1 31 ss = 2 32 else: 33 s = 0 34 ss = 1 35 plt.xlabel("x") 36 plt.ylabel("probability") 37 plt.ylim((0, np.max(y)*1.1)) 38 plt.xlim((np.min(x), np.max(x))) 39 plt.plot(x[s::ss], y[s::ss], color="red", linewidth=1.0) 40 if y_ is not None: 41 plt.plot(x[s::ss], y_[s::ss], color="blue", linewidth=1.0) 42 plt.grid() 43 plt.savefig('q.png') 44 plt.show() 45 46show(x, prob[-1], y_=history_[-1])
参考：
https://arxiv.org/pdf/quant-ph/0303081.pdf
Google画像: Quantum Random Walk
投稿2018/05/07 22:20
編集2018/05/08 12:31
総合スコア8562
0
ベストアンサー
【勘違いでした】
x == n のときに phai の値を変えて、x == n + 1 のときに元に戻してることが原因に見えました。
追記① 一点気づいたバグ（線対称にならない件とは関係ないです）
python
1if x == 0: 2 s_list[0] = np.inner(P, s_list[1]) 3if x == 2*n: 4 s_list[2*n] = np.inner(Q, s_list[2*n-1]) 5else: 6 s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1]) 7 p_list[x] = np.dot(s_list[x],s_list[x])
elifじゃなくてifになってるので、x==0のときにelse節が実行されます。
このとき s_list[x-1] == s_list[-1] == s_listの末尾となるので、
x==300における状態がx==0の状態に影響を及ぼしています。
追記②
python
1for x in range(0,2*n+1): 2 # 中略 3 s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])
↑の部分ですが、s_list[x-1]は既に更新済みなので、「左側の既に更新済みの状態を使って現在位置の状態を更新する」という操作を左から右に向かって行っていることになり、たぶんこれが主原因です。
次のようにnext_s_listってのを新設して対応してみるとー、、、
python
1for t in range(0,3*n): 2 t_list.append(t) 3 if t ==0: 4 s_list 5 p_list 6 else: 7 next_s_list = [0] * len(s_list) # next state list 8 for x in range(0,2*n+1): 9 if x == 0: 10 next_s_list[0] = np.inner(P, s_list[1]) 11 elif x == 2*n: # debug 12 next_s_list[2*n] = np.inner(Q, s_list[2*n-1]) 13 else: 14 next_s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1]) 15 p_list[x] = np.dot(next_s_list[x],next_s_list[x]) 16 s_list = next_s_list
これでわりとそれっぽくなってる気がします（右側が山なり）。
（コード変更によって phai = [1,0] のときのグラフも結構形が変わります。）
違ってたらごめんなさい。
投稿2018/05/07 15:26
編集2018/05/07 17:51
総合スコア2413
考えたんですけど、わかりません。。。泣
やはり、コードがおかしいです。。。
ごめんなさい僕が勘違いしてたようです。本当にすみません。
いえ、大丈夫です！
もし閃いたら、その時、ぜひお力を貸してください！
ありがとうございます！
一度、set0gut1さんの通りにやってみて、試してみます！！
つまり、phai=[1,0]の時の計算状態が残っていて、
それを、知らずにphai=[0,1]を計算しようとしているという解釈でよろしいですか？
tが1増える度にxは両隣（x-1とx+1）からのみ影響を受けるのが正しいようですが、質問に貼られた実装だと一瞬で右側全体に影響が及んでます。
少し考えてみます。。。
0
set0gut1さん、mkgreiさん、コメントやアドバイスありがとうございました！
おかげで理解し、次に進むことができます。
また、わからなくなったらお力添えお願い致します笑
mkgreiさんの、ノルム計算が怪しいという指摘に気づき、複素数にも対応できるようにコードを少し改良しました。おかげで完成できたと思います。
結果を以下にrepositoryしておきます。
何かおかしな点がございましたら、お手数ですが助言の程、よろしくお願い致します。
とても、有意義な質問、時間になりました！！
ありがとうございます。
python
1import numpy as np 2import matplotlib.pyplot as plt 3import math 4from mpl_toolkits.mplot3d import Axes3D 5import matplotlib.animation as animation 6 7#環境設定 8n=400 9theta = 3*(math.pi)/12 10 11P = [[np.cos(theta),np.sin(theta)],[0,0]] 12Q = [[0,0],[np.sin(theta),-np.cos(theta)]] 13x_list=[]#xline 14t_list=[]#time 15p_list=[]#probability 16s_list=[]#state 17a = 1/math.sqrt(2) 18b = 1j/math.sqrt(2) 19 20for x in range(0,2*n+1): 21 if x == n: 22 phai = [a ,b] 23 else: 24 phai = [0,0] 25 p = np.dot(phai,np.conj(phai)) 26 27 x_list.append(x) 28 s_list.append(phai) 29 p_list.append(p) 30 31 32for t in range(0,500): 33 t_list.append(t) 34 if t ==0: 35 s_list 36 p_list 37 else: 38 next_s_list = [0]*len(s_list) #s_listと同じ要素の数ですべて0を用意（初期化） 39 for x in range(0,2*n+1): 40 if x == 0: 41 next_s_list[0] = np.inner(P, s_list[1]) 42 43 elif x == 2*n: 44 next_s_list[2*n] = np.inner(Q, s_list[2*n-1]) 45 46 else: 47 next_s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1]) 48 p_list[x] = np.dot(next_s_list[x],np.conj(next_s_list[x])) 49 s_list = next_s_list 50 51 52 53 print(p_list) 54 plt.xlabel("x") 55 plt.ylabel("probability") 56 plt.ylim([0,0.1]) 57 plt.xlim([-n,3*n]) 58 plt.plot(x_list,np.real(p_list),color="red",linewidth=1.0,) 59 plt.pause(0.01) 60 plt.cla()
投稿2018/05/08 07:52
編集2018/05/08 07:54総合スコア124
とりあえず前進したようで何よりです。
コードとグラフが対応していないので、なんとも言えないです。
グラフの方は少し込み入っているので、詳しくはわかりません。
一つ飛びにゼロになるはずなので、一個飛ばしにプロットすると見通しが少し良くなります。
コードの方ですが、
xの数が2*n+1ならば、境界条件を設けない限りtがnより大きくなるのは問題です。
また、使用する予定のない変数があったりと、少し読みにくくなっています。
if文の分岐を再度検討してみてください。
ほとんど意味がないものになっています。
このままだと、今後手を加えていくとバグを誘発しやすいです。
お気をつけください。
---
ファイは正しくphiと書いたほうがいいです。
たまにうるさく突っ込まれます。
ありがとうございます！！
頑張って、勉強します！
また、改良できたら、掲載しようと思います
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2018/05/08 06:35
2018/05/08 06:59