for ループを用いて自動で計算を行いたい

閲覧ありがとうございます。

以下のコードでは初めに変数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件

ベストアンサー

python
1   eq()
2   Monte()

上記をループすれば良いです。以下のような感じ

python
1for t in range(0.1,2.5,0.05):
2 T = t
3 eq()
4 Monte8)

投稿2018/12/07 02:24

t_obara

総合スコア5488

eq()とMonte()の呼び出しをforかwhileで囲み、Tの足し込みを行う。
0.05は2進数で表現できないので、10から250まで5きざみで変化させ、100で割って使う方が誤差が小さくなる。整数を整数で割ると整数になることにも注意。

投稿2018/12/07 02:21

Q71

総合スコア995

あなたの回答

tips

プレビュー

行動規範の内容に同意します

質問の解決につながる回答をしましょう。サンプルコードなど、より具体的な説明があると質問者の理解の助けになります。また、読む側のことを考えた、分かりやすい文章を心がけましょう。

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

ただいまの回答率
85.48%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

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

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

for ループを用いて自動で計算を行いたい

関連した質問