質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.49%

  • Python 3.x

    6850questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

1次元量子ウォーク(ランダムウォークの量子版)が後少しでできそうですけど、何かがおかしい。その何かがわかりません!!

解決済

回答 3

投稿 編集

  • 評価
  • クリップ 1
  • VIEW 315

Fallout_18

score 65

x軸上にいる物体が左右方向に動いた時の確率分布のコードが完成しました。。。?
【前提条件】______________________________________________
x座標、xにいるときの状態(ここではphai)は関数phai(t,x)とすると、P,Qをある定数として
phai(t,x) = P*phai(t,x+1)+Q*phai(t,x-1)と表します。
(物理学的背景は除く)


あとは初期設定の数値を変更するだけで動くと思ったのですが、、問題が発生しています。
下記のx=nのとき、phai=[0,1]に変更すると、下図のグラフがちょうどx=nに線対称な、右側に山なりになるはずなのですが、実際に動かしても設定を変える前と変わりません。
一方で、θ(theta)の値を変えると、しっかりとネットでいわれているような確率分布になってくれるので、自分が書いたコードが根本的におかしいとは考えられず、解決できずにいます。
以下に、現在の進行状況を挙げます。

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

#初期設定
n=150
theta= 3*(math.pi)/12 #θの値

P = [[np.cos(theta),np.sin(theta)],[0,0]]  #Shift Operator
Q = [[0,0],[np.sin(theta),-np.cos(theta)]] #Shift Operator
x_list=[]#x座標
t_list=[]#time
p_list=[]#probability
s_list=[]#state

for x in range(0,2*n+1):
    if x == n:
        phai = [1,0] #←このphaiの値を[0,1]にしたら、変わるはずなのですが、変わりません。
    else:
        phai = [0,0]
    p = np.dot(phai,phai)

    x_list.append(x)
    s_list.append(phai)
    p_list.append(p)



for t in range(0,3*n):  #以下の文は時間ごとにそれぞれの位置での確率を計算
    t_list.append(t)
    if t ==0:
        s_list
        p_list
    else:
        for x in range(0,2*n+1):
            if x == 0:
                s_list[0] = np.inner(P, s_list[1])
            if x == 2*n:
                s_list[2*n] = np.inner(Q, s_list[2*n-1])
            else:
                s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])
                p_list[x] = np.dot(s_list[x],s_list[x])



    print(x_list,p_list)
    #plt.bar(x_list,p_list,width=0.1,color='black')
    plt.xlabel("x")
    plt.ylabel("probability")
    plt.ylim([0,0.1])
    plt.xlim([-n,3*n])
    plt.plot(x_list,p_list,color="red",linewidth=1.0)
    plt.pause(0.001)
    plt.cla()


と行うと、以下のような確率の挙動(今回は一部を貼り付け)がわかりました。
θ=π/4
確率分布(θ=π/4)
θ=5π/12
イメージ説明
となりました。
ご指摘の程、宜しくお願いします。
python初心者ですので、お手柔らかにお願いします。。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

checkベストアンサー

+3

【勘違いでした】
x == n のときに phai の値を変えて、x == n + 1 のときに元に戻してることが原因に見えました。


追記① 一点気づいたバグ(線対称にならない件とは関係ないです)

if x == 0:
    s_list[0] = np.inner(P, s_list[1])
if x == 2*n:
    s_list[2*n] = np.inner(Q, s_list[2*n-1])
else:
    s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])
    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の状態に影響を及ぼしています。


追記②

for x in range(0,2*n+1):
        # 中略
        s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])

↑の部分ですが、s_list[x-1]は既に更新済みなので、「左側の既に更新済みの状態を使って現在位置の状態を更新する」という操作を左から右に向かって行っていることになり、たぶんこれが主原因です。

次のようにnext_s_listってのを新設して対応してみるとー、、、

for t in range(0,3*n):
    t_list.append(t)
    if t ==0:
        s_list
        p_list
    else:
        next_s_list = [0] * len(s_list) # next state list
        for x in range(0,2*n+1):
            if x == 0:
                next_s_list[0] = np.inner(P, s_list[1])
            elif x == 2*n: # debug
                next_s_list[2*n] = np.inner(Q, s_list[2*n-1])
            else:
                next_s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])
            p_list[x] = np.dot(next_s_list[x],next_s_list[x])
        s_list = next_s_list

これでわりとそれっぽくなってる気がします(右側が山なり)。
(コード変更によって phai = [1,0] のときのグラフも結構形が変わります。)
違ってたらごめんなさい。

わりとそれっぽいグラフ

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/08 00:33

    ごめんなさい、わかりそうでわからないです。
    もう少し詳しくお願いします。

    キャンセル

  • 2018/05/08 00:48

    考えたんですけど、わかりません。。。泣
    やはり、コードがおかしいです。。。

    キャンセル

  • 2018/05/08 01:21

    ごめんなさい僕が勘違いしてたようです。本当にすみません。

    キャンセル

  • 2018/05/08 01:31

    いえ、大丈夫です!
    もし閃いたら、その時、ぜひお力を貸してください!

    キャンセル

  • 2018/05/08 15:35

    ありがとうございます!
    一度、set0gut1さんの通りにやってみて、試してみます!!

    キャンセル

  • 2018/05/08 15:38

    つまり、phai=[1,0]の時の計算状態が残っていて、
    それを、知らずにphai=[0,1]を計算しようとしているという解釈でよろしいですか?

    キャンセル

  • 2018/05/08 15:41

    tが1増える度にxは両隣(x-1とx+1)からのみ影響を受けるのが正しいようですが、質問に貼られた実装だと一瞬で右側全体に影響が及んでます。

    キャンセル

  • 2018/05/08 15:44

    少し考えてみます。。。

    キャンセル

+3

s_listの更新のタイミングが悪いことと、
確保してあるs_listのサイズが足りないことが問題だと思います。

端の場所がおかしいのも問題ですね…


wikipedia先生によると対称になるのは[1,i]/sqrt(2)の時ですが…
それが正しければ、ノルムの計算方法もまずいです。


import numpy as np
import matplotlib.pyplot as plt

n = 100
H = np.array([[1, 1],[1, -1]]) / np.sqrt(2)
P = np.zeros((2, 2)); P[0, :] = H[0, :]
Q = np.zeros((2, 2)); Q[1, :] = H[1, :]

pos = np.zeros((2*n+1, 2), dtype=np.complex)
pos[n, :] = np.array([1, 1j], dtype=np.complex) / np.sqrt(2)
x = np.arange(-n, n+1)

history = np.zeros((n, 2*n+1, 2), dtype=np.complex)
history[0, :, :] = pos[:, :]
prob = np.zeros((n, 2*n+1))
prob[0, :] = np.sum(np.conj(history[0, :, :]) * history[0, :, :], axis=-1)[:]
for t in range(1, n):
    history[t, :-1, :] += np.inner(history[t-1, 1:, :], P.T)
    history[t, 1:, :]  += np.inner(history[t-1, :-1, :], Q.T)
    prob[t, :] = np.sum(np.conj(history[t, :, :]) * history[t, :, :], axis=-1)[:]

history_ = np.zeros((n, 2*n+1))
history_[0, n] = 1.
for t in range(1, n):
    history_[t, :-1] += history_[t-1, 1:]/2.
    history_[t, 1:]  += history_[t-1, :-1]/2.

def show(x, y, y_=None, nonzeros=True):
    if nonzeros:
        s = 1
        ss = 2
    else:
        s = 0
        ss = 1
    plt.xlabel("x")
    plt.ylabel("probability")
    plt.ylim((0, np.max(y)*1.1))
    plt.xlim((np.min(x), np.max(x)))
    plt.plot(x[s::ss], y[s::ss], color="red", linewidth=1.0)
    if y_ is not None:
        plt.plot(x[s::ss], y_[s::ss], color="blue", linewidth=1.0)
    plt.grid()
    plt.savefig('q.png')
    plt.show()

show(x, prob[-1], y_=history_[-1])

参考:
https://arxiv.org/pdf/quant-ph/0303081.pdf
Google画像: Quantum Random Walk

QRW:赤、RW:青

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/08 15:35

    すみません、それだと下記のset0gut1さんのやり方でも、だめだということでしょうか?
    質問が当り障りで申し訳ないです。

    キャンセル

  • 2018/05/08 15:59

    少しわかった気がします

    キャンセル

+1

set0gut1さん、mkgreiさん、コメントやアドバイスありがとうございました!
おかげで理解し、次に進むことができます。
また、わからなくなったらお力添えお願い致します笑
mkgreiさんの、ノルム計算が怪しいという指摘に気づき、複素数にも対応できるようにコードを少し改良しました。おかげで完成できたと思います。
結果を以下にrepositoryしておきます。
何かおかしな点がございましたら、お手数ですが助言の程、よろしくお願い致します。
とても、有意義な質問、時間になりました!!
ありがとうございます。

import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

#環境設定
n=400
theta = 3*(math.pi)/12

P = [[np.cos(theta),np.sin(theta)],[0,0]]
Q = [[0,0],[np.sin(theta),-np.cos(theta)]]
x_list=[]#xline
t_list=[]#time
p_list=[]#probability
s_list=[]#state
a = 1/math.sqrt(2)
b = 1j/math.sqrt(2)

for x in range(0,2*n+1):
    if x == n:
        phai = [a ,b]
    else:
        phai = [0,0]
    p = np.dot(phai,np.conj(phai))

    x_list.append(x)
    s_list.append(phai)
    p_list.append(p)


for t in range(0,500):
    t_list.append(t)
    if t ==0:
        s_list
        p_list
    else:
        next_s_list = [0]*len(s_list)  #s_listと同じ要素の数ですべて0を用意(初期化)
        for x in range(0,2*n+1):
            if x == 0:
                next_s_list[0] = np.inner(P, s_list[1])

            elif x == 2*n:
                next_s_list[2*n] = np.inner(Q, s_list[2*n-1])

            else:
                next_s_list[x] = np.inner(P, s_list[x+1]) + np.inner(Q, s_list[x-1])
            p_list[x] = np.dot(next_s_list[x],np.conj(next_s_list[x]))
        s_list = next_s_list



    print(p_list)
    plt.xlabel("x")
    plt.ylabel("probability")
    plt.ylim([0,0.1])
    plt.xlim([-n,3*n])
    plt.plot(x_list,np.real(p_list),color="red",linewidth=1.0,)
    plt.pause(0.01)
    plt.cla()


イメージ説明

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/08 22:10 編集

    とりあえず前進したようで何よりです。

    コードとグラフが対応していないので、なんとも言えないです。
    グラフの方は少し込み入っているので、詳しくはわかりません。
    一つ飛びにゼロになるはずなので、一個飛ばしにプロットすると見通しが少し良くなります。

    コードの方ですが、
    xの数が2*n+1ならば、境界条件を設けない限りtがnより大きくなるのは問題です。
    また、使用する予定のない変数があったりと、少し読みにくくなっています。
    if文の分岐を再度検討してみてください。
    ほとんど意味がないものになっています。
    このままだと、今後手を加えていくとバグを誘発しやすいです。
    お気をつけください。

    ---

    ファイは正しくphiと書いたほうがいいです。
    たまにうるさく突っ込まれます。

    キャンセル

  • 2018/05/08 22:56

    ありがとうございます!!
    頑張って、勉強します!

    キャンセル

  • 2018/05/08 23:04

    また、改良できたら、掲載しようと思います

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.49%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る

  • Python 3.x

    6850questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

  • トップ
  • Python 3.xに関する質問
  • 1次元量子ウォーク(ランダムウォークの量子版)が後少しでできそうですけど、何かがおかしい。その何かがわかりません!!