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

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

ただいまの
回答率

90.34%

  • Python

    9180questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • Python 3.x

    7387questions

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

sierpinski carpet 上を歩く計算がうまくいかない。

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 124

Fallout_18

score 67

前回の質問に関連しています。
以下のシュレピンスキーカーペット(stage2とする)を...
イメージ説明
以下の式で時間発展させます(詳細は除く)(初期状態のphiは[0,0,0,0]or[1,0,0,0])。


phi(x,y,t) = P*phi(x-1,y,t-1) + Q*phi(x+1,y,t-1) + R*phi(x,y+1,t-1) + S*phi(x,y-1,t-1)...① (状態計算)
p(x,y,t) = phi(x,y,t).np.conj(phi(x,y,t))....② (確率計算)


ただし、上のstage2を参照して、0のところは通れないとします。
なので
stage2[x,y]=0 のところは常に[0,0,0,0]
x,yを上の座標に変化させていった時に、
stage2[x+1,y]=0 ①式のQは作用させない
stage2[x-1,y]=0      P
.                    R
.                    S
.                    .
.                    .

...として計算していきました。
それが以下のコードになります。

import numpy as np
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.colors as colors
#######
stage1=np.array([[1,1,1],[1,0,1],[1,1,1]],dtype="float")
b=np.zeros([3,3])
#######
c=np.vstack([stage1,stage1,stage1])
e=np.hstack([stage1,b,stage1])
#######
stage2=np.hstack([c,e.T,c])
#######
h = np.zeros([9,9])
f = np.vstack([stage2,stage2,stage2])
g = np.hstack([stage2,h,stage2])
stage3=np.hstack([f,g.T,f])
#######
H = np.zeros([27,27])
J = np.vstack([stage3,stage3,stage3])
K = np.hstack([stage3,H,stage3])
stage4 = np.hstack([J,K.T,J])
#######
n = 2
m = 10
######
x_list=[i for i in range(3*n+3)]
y_list=[i for i in range(3*n+3)]
t_list=[i for i in range(m+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]] 
#######
p_map = np.zeros([3*n+3,3*n+3])
phi_map = np.zeros((3*n+3,3*n+3,4),dtype="complex")
next_phi_map = np.zeros((3*n+3,3*n+3, 4),dtype="complex")
phi_map[0,0]=np.array([1,0,0,0])
p_map[0,0]=np.real(np.inner(phi_map[0,0],np.conj(phi_map[0,0])))
##########
for t in range(0,m):
    if t == 0:
        p_map
        phi_map
    else:
        for x in range(0,3*n+3):
            if x == 0:
                for y in range(0,3*n+3):
                    if y == 0:
                        next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(Q,phi_map[x+1,y])])
                    elif y == 3*n+2:
                        next_phi_map[x,y]=np.array([np.dot(S,phi_map[x,y-1])+np.dot(Q,phi_map[x+1,y])])
                    else:
                        if stage2[x+1,y] == 0:
                            next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
                        else:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
            elif x == 3*n+2:
                for y in range(0,3*n+3):
                    if y == 0:
                        next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(P,phi_map[x-1,y])])
                    elif y == 3*n+2:
                        next_phi_map[x,y]=np.array([np.dot(S,phi_map[x,y-1])+np.dot(P,phi_map[x-1,y])])
                    else:
                        if stage2[x-1,y]==0:
                            next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
                        else:
                            next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
            else:
                for y in range(0,3*n+3):
                    if y == 0:
                        if stage2[x,y+1] == 0:
                            next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])])
                        else:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(R,phi_map[x,y+1])+np.dot(P,phi_map[x-1,y])])
                    elif y == 3*n+2:
                        if stage2[x,y-1] == 0:
                            next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])])
                        else:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(Q,phi_map[x,y-1])+np.dot(P,phi_map[x-1,y])])
                    else:
                        if stage2[x,y-1] == 0:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(R,phi_map[x,y+1])+np.dot(P,phi_map[x-1,y])])
                        elif stage2[x,y+1] == 0:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(S,phi_map[x,y-1])+np.dot(P,phi_map[x-1,y])])
                        elif stage2[x-1,y] == 0:
                            next_phi_map[x,y]=np.array([np.dot(Q,phi_map[x+1,y])+np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
                        elif stage2[x+1,y] == 0:
                            next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
                        elif stage2[x,y-1] == 0 and stage3[x,y+1] == 0:
                            next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])])
                        elif stage2[x-1,y] == 0 and stage3[x+1,y] == 0:
                            next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])])
                        elif stage2[x,y] == 0 :
                            next_phi_map[x,y] =np.array([0,0,0,0])
                    p_map[x,y]=np.real(np.inner(next_phi_map[x,y], np.conj(next_phi_map[x,y])))
        phi_map=next_phi_map
##############
fig = plt.figure()
ax = Axes3D(fig, rect=(0.1,0.1,0.8,0.8))
X,Y = np.meshgrid(x_list, y_list)
ax.set_xlabel("x",labelpad=10,fontsize=24)
ax.set_ylabel("y",labelpad=20,fontsize=24)
ax.set_zlabel("$|\psi|^2$",labelpad=10,fontsize=24)
ax.set_xlim(3*n+3,0)
ax.set_ylim(0,3*n+3)
ax.set_zlim(0,1)
#############
mask= p_map > 0.0
offset = p_map[mask].ravel() + np.abs(p_map[mask].min())
fracs = offset.astype(float)/offset.max()
norm = colors.Normalize(fracs.min(), fracs.max())
clrs = cm.cool(norm(fracs))
###########
ax.bar3d(X[mask].ravel(), Y[mask].ravel(), p_map[mask].ravel() ,0.1, 0.1, -p_map[mask].ravel(),color =clrs)
ax.w_xaxis.set_pane_color((0, 0, 0, 0))
ax.w_yaxis.set_pane_color((0, 0, 0, 0))
ax.w_zaxis.set_pane_color((0, 0, 0, 1))
ax.grid(color="white")
ax.grid(False)
############
plt.show()


結果
![結果](4ea506de3f4486e4ef45169354e63439.png)

見るからに計算されておらず、結構考えたのですが、頭がパンクしそうで助けが欲しくなり、質問させて頂きました。
ちなみになのですが、前回の質問のやり方では、やはりstageが上がるにつれて面倒くさくなるので、やめました。
その時の結果は一応のっけて置きます(stage2,t=10)。
イメージ説明
最後に、求めるべき確率の正確な答えはわかりませんが、計算されていないという点で上記の結果が間違っていると踏んだので質問させて頂きました。
参照論文


追記】
いろいろ直してみて以下のようになりました(stage4です)。
イメージ説明
メイン計算がまだ余計に場合分けしている可能性が出てきたので、見つけた方、お手数ですがご指摘お願いします。自分でも改良します。

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

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

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

    クリップを取り消します

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

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

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

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

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

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

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

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

    質問の評価を下げる

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

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

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

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

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

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

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

    詳細な説明はこちら

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

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

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

回答 1

check解決した方法

0

イメージ説明
上記の質問に対して、自分で結果が得られました。

投稿

  • 回答の評価を上げる

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

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

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

  • 回答の評価を下げる

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

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

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

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

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

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

  • Python

    9180questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • Python 3.x

    7387questions

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