Q&A
前回の質問に関連しています。
以下のシュレピンスキーカーペット(stage2とする)を...
以下の式で時間発展させます（詳細は除く）(初期状態のphiは[0,0,0,0]or[1,0,0,0])。
phi(x,y,t) = Pphi(x-1,y,t-1) + Qphi(x+1,y,t-1) + Rphi(x,y+1,t-1) + Sphi(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
. .
. .
...として計算していきました。
それが以下のコードになります。
python
1import numpy as np 2import matplotlib.pyplot as plt 3import math 4from mpl_toolkits.mplot3d import Axes3D 5from matplotlib import cm 6import matplotlib.colors as colors 7####### 8stage1=np.array([[1,1,1],[1,0,1],[1,1,1]],dtype="float") 9b=np.zeros([3,3]) 10####### 11c=np.vstack([stage1,stage1,stage1]) 12e=np.hstack([stage1,b,stage1]) 13####### 14stage2=np.hstack([c,e.T,c]) 15####### 16h = np.zeros([9,9]) 17f = np.vstack([stage2,stage2,stage2]) 18g = np.hstack([stage2,h,stage2]) 19stage3=np.hstack([f,g.T,f]) 20####### 21H = np.zeros([27,27]) 22J = np.vstack([stage3,stage3,stage3]) 23K = np.hstack([stage3,H,stage3]) 24stage4 = np.hstack([J,K.T,J]) 25####### 26n = 2 27m = 10 28###### 29x_list=[i for i in range(3*n+3)] 30y_list=[i for i in range(3*n+3)] 31t_list=[i for i in range(m+1)] 32####### 33P = [[-1/2, 1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0],[0,0,0,0]] 34Q = [[0,0,0,0],[1/2, -1/2, 1/2, 1/2],[0,0,0,0],[0,0,0,0]] 35R = [[0,0,0,0],[0,0,0,0],[1/2, 1/2, -1/2, 1/2],[0,0,0,0]] 36S = [[0,0,0,0],[0,0,0,0],[0,0,0,0],[1/2, 1/2, 1/2, -1/2]] 37####### 38p_map = np.zeros([3*n+3,3*n+3]) 39phi_map = np.zeros((3*n+3,3*n+3,4),dtype="complex") 40next_phi_map = np.zeros((3*n+3,3*n+3, 4),dtype="complex") 41phi_map[0,0]=np.array([1,0,0,0]) 42p_map[0,0]=np.real(np.inner(phi_map[0,0],np.conj(phi_map[0,0]))) 43########## 44for t in range(0,m): 45 if t == 0: 46 p_map 47 phi_map 48 else: 49 for x in range(0,3*n+3): 50 if x == 0: 51 for y in range(0,3*n+3): 52 if y == 0: 53 next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(Q,phi_map[x+1,y])]) 54 elif y == 3*n+2: 55 next_phi_map[x,y]=np.array([np.dot(S,phi_map[x,y-1])+np.dot(Q,phi_map[x+1,y])]) 56 else: 57 if stage2[x+1,y] == 0: 58 next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])]) 59 else: 60 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])]) 61 elif x == 3*n+2: 62 for y in range(0,3*n+3): 63 if y == 0: 64 next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(P,phi_map[x-1,y])]) 65 elif y == 3*n+2: 66 next_phi_map[x,y]=np.array([np.dot(S,phi_map[x,y-1])+np.dot(P,phi_map[x-1,y])]) 67 else: 68 if stage2[x-1,y]==0: 69 next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])]) 70 else: 71 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])]) 72 else: 73 for y in range(0,3*n+3): 74 if y == 0: 75 if stage2[x,y+1] == 0: 76 next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])]) 77 else: 78 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])]) 79 elif y == 3*n+2: 80 if stage2[x,y-1] == 0: 81 next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])]) 82 else: 83 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])]) 84 else: 85 if stage2[x,y-1] == 0: 86 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])]) 87 elif stage2[x,y+1] == 0: 88 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])]) 89 elif stage2[x-1,y] == 0: 90 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])]) 91 elif stage2[x+1,y] == 0: 92 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])]) 93 elif stage2[x,y-1] == 0 and stage3[x,y+1] == 0: 94 next_phi_map[x,y]=np.array([np.dot(P,phi_map[x-1,y])+np.dot(Q,phi_map[x+1,y])]) 95 elif stage2[x-1,y] == 0 and stage3[x+1,y] == 0: 96 next_phi_map[x,y]=np.array([np.dot(R,phi_map[x,y+1])+np.dot(S,phi_map[x,y-1])]) 97 elif stage2[x,y] == 0 : 98 next_phi_map[x,y] =np.array([0,0,0,0]) 99 p_map[x,y]=np.real(np.inner(next_phi_map[x,y], np.conj(next_phi_map[x,y]))) 100 phi_map=next_phi_map 101############## 102fig = plt.figure() 103ax = Axes3D(fig, rect=(0.1,0.1,0.8,0.8)) 104X,Y = np.meshgrid(x_list, y_list) 105ax.set_xlabel("x",labelpad=10,fontsize=24) 106ax.set_ylabel("y",labelpad=20,fontsize=24) 107ax.set_zlabel("$|\psi|^2$",labelpad=10,fontsize=24) 108ax.set_xlim(3*n+3,0) 109ax.set_ylim(0,3*n+3) 110ax.set_zlim(0,1) 111############# 112mask= p_map > 0.0 113offset = p_map[mask].ravel() + np.abs(p_map[mask].min()) 114fracs = offset.astype(float)/offset.max() 115norm = colors.Normalize(fracs.min(), fracs.max()) 116clrs = cm.cool(norm(fracs)) 117########### 118ax.bar3d(X[mask].ravel(), Y[mask].ravel(), p_map[mask].ravel() ,0.1, 0.1, -p_map[mask].ravel(),color =clrs) 119ax.w_xaxis.set_pane_color((0, 0, 0, 0)) 120ax.w_yaxis.set_pane_color((0, 0, 0, 0)) 121ax.w_zaxis.set_pane_color((0, 0, 0, 1)) 122ax.grid(color="white") 123ax.grid(False) 124############ 125plt.show() 126
結果
見るからに計算されておらず、結構考えたのですが、頭がパンクしそうで助けが欲しくなり、質問させて頂きました。
ちなみになのですが、前回の質問のやり方では、やはりstageが上がるにつれて面倒くさくなるので、やめました。
その時の結果は一応のっけて置きます(stage2,t=10)。
最後に、求めるべき確率の正確な答えはわかりませんが、計算されていないという点で上記の結果が間違っていると踏んだので質問させて頂きました。
参照論文
追記】
いろいろ直してみて以下のようになりました(stage4です)。
メイン計算がまだ余計に場合分けしている可能性が出てきたので、見つけた方、お手数ですがご指摘お願いします。自分でも改良します。
自己解決
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。