閲覧ありがとうございます。
以下のコードでは初めに変数Tを指定することでHという値を出力するプログラムになっています。
python
1# import library 2import numpy as np 3import matplotlib.pyplot as plt 4from random import random,randrange 5 6np.random.seed(0) 7 8L = 64 9Nx = int(L) 10Ny = int(L) 11T = 0.1 12ESTEP = 40000 13STEP = 10000 14 15 16#S = np.random.uniform(0,2*np.pi,(L,L)) 17S = np.zeros((Nx,Ny),dtype = float) 18 19def eq(): 20 for x in range(ESTEP): 21 22 23 i = randrange(Nx) 24 j = randrange(Ny) 25 26 S_al = np.random.uniform(0,2*np.pi) 27 28 E = -np.cos(S[i , j] - S[(i+1)%Nx , j])\ 29 -np.cos(S[i , j] - S[(i-1)%Nx , j])\ 30 -np.cos(S[i , j] - S[i , (j+1)%Ny])\ 31 -np.cos(S[i , j] - S[i , (j-1)%Ny]) 32 33 E_al = -np.cos(S_al - S[(i+1)%Nx , j])\ 34 -np.cos(S_al - S[(i-1)%Nx , j])\ 35 -np.cos(S_al - S[i , (j+1)%Ny])\ 36 -np.cos(S_al - S[i , (j-1)%Ny]) 37 38 deltaE = E_al -E 39 40 if deltaE < 0: 41 42 #スピンを更新 43 S[i,j] = S_al 44 #エネルギーを計算 45 for i in range(Nx): 46 for j in range(Ny): 47 e = -np.cos(S[i , j] - S[(i+1)%Nx , j])\ 48 -np.cos(S[i , j] - S[(i-1)%Nx , j])\ 49 -np.cos(S[i , j] - S[i , (j+1)%Ny])\ 50 -np.cos(S[i , j] - S[i , (j-1)%Ny]) 51 52 53 54 55 56 57 else: 58 if np.random.random() < np.exp(-deltaE/T): 59 60 #スピンを更新 61 S[i,j] = S_al 62 #エネルギーを計算 63 for i in range(Nx): 64 for j in range(Ny): 65 e = -np.cos(S[i , j] - S[(i+1)%Nx , j])\ 66 -np.cos(S[i , j] - S[(i-1)%Nx , j])\ 67 -np.cos(S[i , j] - S[i , (j+1)%Ny])\ 68 -np.cos(S[i , j] - S[i , (j-1)%Ny]) 69 70 71 #print(np.sum(esum)/L**2/2,x) 72 #print(np.sum(mag)/L**2,x) 73 74 75 76 return S 77 78 79def Monte(): 80 S = eq() 81 82 Y11 = [] 83 Y22 = [] 84 85 86 for y in range(STEP): 87 88 Y1 = [] 89 Y2 = [] 90 91 i = randrange(Nx) 92 j = randrange(Ny) 93 94 S_al = np.random.uniform(0,2*np.pi) 95 96 E = -np.cos(S[i , j] - S[(i+1)%Nx , j])\ 97 -np.cos(S[i , j] - S[(i-1)%Nx , j])\ 98 -np.cos(S[i , j] - S[i , (j+1)%Ny])\ 99 -np.cos(S[i , j] - S[i , (j-1)%Ny]) 100 101 E_al = -np.cos(S_al - S[(i+1)%Nx , j])\ 102 -np.cos(S_al - S[(i-1)%Nx , j])\ 103 -np.cos(S_al - S[i , (j+1)%Ny])\ 104 -np.cos(S_al - S[i , (j-1)%Ny]) 105 106 deltaE = E_al -E 107 108 if deltaE < 0: 109 S[i,j] = S_al 110 111 #helicity 112 for i in range(Nx): 113 for j in range(Ny): 114 y1 = np.cos(S[i,j]-S[(i+1)%Nx,j]) 115 y2 = np.sin(S[i,j]-S[(i+1)%Nx,j]) 116 Y1.append(y1) 117 Y2.append(y2) 118 119 120 Y11.append(np.sum(Y1)/L**2) 121 Y22.append(np.sum(Y2)/L**2) 122 123 124 else: 125 if np.random.random() < np.exp(-deltaE/T): 126 127 #スピンを更新 128 S[i,j] = S_al 129 130 #helicity 131 for i in range(Nx): 132 for j in range(Ny): 133 y1 = np.cos(S[i,j]-S[(i+1)%Nx,j]) 134 y2 = np.sin(S[i,j]-S[(i+1)%Nx,j]) 135 Y1.append(y1) 136 Y2.append(y2) 137 138 139 Y11.append(np.sum(Y1)/L**2) 140 Y22.append(np.sum(Y2)/L**2) 141 142 143 144 145 146 147 148 del Y1[:] 149 del Y2[:] 150 151 H = (np.sum(Y11)/len(Y11)) -(((np.sum(Y22))/(len(Y22)))**2)*(L**2/T) 152 print(H) 153 print(T,L) 154 155 156 157if __name__ == '__main__': 158 eq() 159 Monte() 160
Tの値は0.1~2.5まで0.05の幅で変化させ、その度にHを得るという作業をしているのですが、どうにかしてこの作業を自動化できないでしょうか。具体的には、forループを用いて、一度runさせるとT=0.1~2.5のすべてのTにおけるHをprintするというプログラムに書き換えたいと考えています。詳しい方、助言をいただけると助かります。
回答2件
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。