質問編集履歴

6

アルゴリズムにミスがあったので更新しました

2018/10/24 13:57

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -6,7 +6,11 @@
6
6
 
7
7
 
8
8
 
9
+ **20181024 追記メトロポリスアルゴリズムにミスがあったため、全面的に書きなおしました。**
10
+
11
+
12
+
9
- **初期状態の形成**
13
+ **初期状態の形成**
10
14
 
11
15
 
12
16
 
@@ -30,8 +34,6 @@
30
34
 
31
35
 
32
36
 
33
- **メトロポリスアルゴリズム(10/24追記:(6)以下で一部間違っています。下に正しいものを載せました。)**
34
-
35
37
 
36
38
 
37
39
  (4)行列Sのある要素S(i,j)をランダムに抽出し、その要素だけを0~2πの乱数で変動させる。
@@ -40,55 +42,27 @@
40
42
 
41
43
  (5)変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
42
44
 
43
- Eal(i,j)=cos(Sal(i,j)-S(i+1,j))-cos(Sal(i,j)-S(i-1,j))-cos(Sal(i,j)-S(i,j+1))-cos(Sal(i,j)-S(i,j-1))
44
-
45
-
46
-
47
- つまり、変動させた(i,j)以外の要素はそのままにして、手順(3)で行ったものと同じ計算をするということです。
48
-
49
-
50
-
51
- ~~(6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較ます
52
-
53
- ~~
54
-
55
- ~~(6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]Sal[i,j]に更新し、(7)以降過程へ
56
-
57
- ~~
58
-
59
- ~~(6-2) E[i,j] < Eal[i,j] 場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算ます。そし確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
60
-
61
- ~~
62
-
63
- ~~(7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するうにします。
64
-
65
- ~~
66
-
67
- ~~(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum=np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
68
-
69
- 参考程度にですが、Ssum~64、Esum~-250 (Nx=8、Ny=8)  となれば成功と言えます。~~
70
-
71
- ~~
72
-
73
- **10/24追記 (6)以下のメトロポリスアルゴリズム(文献を再度確認しなおしました。)**
74
-
75
- (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
76
-
77
-
78
-
79
- (6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
80
-
81
-
82
-
83
- (6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
84
-
85
-
86
-
87
- (7)(4)~(6-2)を1000回繰り返す。(ここでは(8)でしているような計算は行わず、Sの更新だけを行います。)
88
-
89
-
90
-
91
- (8)(7)の後(4)~(6-2)を10000回繰り返す。しかし、この時1回ごとにS_sum = np.sum(S)、E_sum = np.sun(E)を計算し、その和を求める。つまりi回目の更新で得られたS_sumをS_sum(i)とすると、S_end(10000)=S_sum(1)+S_sum(2)+S_sum(3)...+S_sum(10000)を求めるということです。E_end(10000)についても同様です。
45
+
46
+
47
+ つまり、変動させた(i,j)以外の要素はそのままにして、手順(3)で行ったものと同じ計算をするということです。この時、行列Sの要素S[i,j]の変動に伴い、行列Eの要素E[i,j],E[i+1,j],E[i-1,j],E[i,j+1],E[1,j-1]も変動することに注意してください。
48
+
49
+
50
+
51
+ (6)変動前の行列Eおよび変動後の行列Eの行列式をそれぞれE、Ealとします。
52
+
53
+ (6-1)E > Eal の場合、要素S[i,j]→Sal[i,j]への更新に成功したとし、(7)以降の過程へ
54
+
55
+
56
+
57
+ (6-2) E < Eal の場合、deltaE=Eal-Eを計算します。そして確率exp(-deltaE/0.001)でもともとの要素S[i,j]Sal[i,j]への更新成功したとします。ここで(6)~(6-2)においてS[i,j]の更新に成功た場合は行列E[i,j],E[i+1,j],E[i-1,j],E[i,j+1],E[1,j-1]要素も当然変化します
58
+
59
+
60
+
61
+ (7)(4)~(6-2)を1回試行として、1000回繰り返す。
62
+
63
+
64
+
65
+ (8)(7)の後(4)~(6-2)を10000回繰り返す。しかし、の時1回ごにS_sum = np.sum(S)、E_sum = np.sun(E)を計算し、足し合わせいく。つまりi回目の更新で得られたS_sumS_sum(i)とすると、S_end(10000)=S_sum(1)+S_sum(2)+S_sum(3)...+S_sum(10000)を求めるといことです。E_end(10000)ついても同様です。
92
66
 
93
67
 
94
68
 

5

アルゴリズムにミスがあったので、訂正を行いました。

2018/10/24 13:57

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -1,6 +1,6 @@
1
1
  ### 仕様
2
2
 
3
- コードにもコメント文で書いていますが、以下に具体的にきたいと思います。
3
+ コードにもコメント文で書いていますが、以下に具体的にきたいと思います。
4
4
 
5
5
  メトロポリス法という物理学の中の統計力学におけるシミュレーションをを行う際の計算方法です。
6
6
 
@@ -30,7 +30,7 @@
30
30
 
31
31
 
32
32
 
33
- **メトロポリスアルゴリズム**
33
+ **メトロポリスアルゴリズム(10/24追記:(6)以下で一部間違っています。下に正しいものを載せました。)**
34
34
 
35
35
 
36
36
 
@@ -48,11 +48,35 @@
48
48
 
49
49
 
50
50
 
51
+ ~~(6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
52
+
53
+ ~~
54
+
55
+ ~~(6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
56
+
57
+ ~~
58
+
59
+ ~~(6-2) E[i,j] < Eal[i,j] の場合、deltaE[i,j]=Eal[i,j]-E[i,j]を計算します。そして確率exp(-deltaE[i,j]/0.001)でもともとの要素S[i,j]をSal[i,j]に更新します。
60
+
61
+ ~~
62
+
63
+ ~~(7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
64
+
65
+ ~~
66
+
67
+ ~~(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum=np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
68
+
69
+ 参考程度にですが、Ssum~64、Esum~-250 (Nx=8、Ny=8)  となれば成功と言えます。~~
70
+
71
+ ~~
72
+
73
+ **10/24追記 (6)以下のメトロポリスアルゴリズム(文献を再度確認しなおしました。)**
74
+
51
75
  (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
52
76
 
53
77
 
54
78
 
55
- (6-1E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
79
+ (6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
56
80
 
57
81
 
58
82
 
@@ -60,13 +84,25 @@
60
84
 
61
85
 
62
86
 
63
- (7)(4)~(6)を10000回繰り返すこおおそすべて要素S[i,j]がスピン更新を経験するようにします。
64
-
65
-
66
-
67
- (8)(7)の操作の行列S用いて行列Eを計算それぞれ行列の行列式の値Ssumnp.sum(S)、Esum=np.sum(E)をとめることきれば計算終了です。
68
-
69
- 参考程度にですが、Ssum~64、Esum~-250 (Nx=8、Ny=8)  となれば成功と言えます。
87
+ (7)(4)~(6-2)を1000回繰り返す。(は(8)でしているうな計算は行わず、Sの更新だけ行います。)
88
+
89
+
90
+
91
+ (8)(7)の後(4)~(6-2)10000回繰り返す。しかし、こ時1回ごとにS_sum = np.sum(S)、E_sum = np.sun(E)を計算し、その和を求める。つまりi回目の更新で得られたS_sumをS_sum(i)すると、S_end(10000)=S_sum(1)+S_sum(2)+S_sum(3)...+S_sum(10000)を求めるということです。E_end(10000)についても同様です。
92
+
93
+
94
+
95
+ (9)(8)で求めた値を用いて、S_ave = S_end(10000)/10000,E_ave = E_end(10000)/10000 を出力できればOKです。
96
+
97
+
98
+
99
+
100
+
101
+
102
+
103
+
104
+
105
+
70
106
 
71
107
 
72
108
 

4

仕様の項目に追記を行いました

2018/10/24 08:40

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -64,9 +64,9 @@
64
64
 
65
65
 
66
66
 
67
- (8)(7)の操作後の行列式Send=np.sum(S)及びEend=np.sum(E)を出力できればこのシミュレーションは終了です。
67
+ (8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列の値Ssum=np.sum(S)Esum=np.sum(E)をもとめることができれば計算終了です。
68
+
68
-
69
+ 参考程度にですが、Ssum~64、Esum~-250 (Nx=8、Ny=8)  となれば成功と言えます。
69
-
70
70
 
71
71
 
72
72
 

3

追記2を更新しました

2018/10/19 22:34

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -319,3 +319,65 @@
319
319
 
320
320
 
321
321
  ```
322
+
323
+
324
+
325
+ ###追記2(#更新部分について、教えていただいた内容をもとに条件分岐を訂正しました)
326
+
327
+ ```python
328
+
329
+ def Metropolis(Nx=8, Ny=8, iterations=10000):
330
+
331
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
332
+
333
+ E = update(S)
334
+
335
+
336
+
337
+ indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
338
+
339
+ for k in range(iterations): # モンテカルロステップ数を 10000 とする
340
+
341
+ # ランダムに要素を選択する。
342
+
343
+ index = tuple(indices[np.random.randint(len(indices))])
344
+
345
+
346
+
347
+ # 乱数を生成する。
348
+
349
+ E_tmp = E.copy()
350
+
351
+ E_tmp[index] = np.random.uniform(0, 2 * np.pi)
352
+
353
+
354
+
355
+ # 計算
356
+
357
+ E_tmp = update(E_tmp)
358
+
359
+
360
+
361
+ # 更新
362
+
363
+ if E[index] > E_tmp[index]:
364
+
365
+ E[index] = E_tmp[index]
366
+
367
+
368
+
369
+ else:
370
+
371
+ delta = E_tmp[index] - E[index]
372
+
373
+ if np.random.random() < np.exp(-delta / 0.001):
374
+
375
+ E[index] = E_tmp[index]
376
+
377
+
378
+
379
+
380
+
381
+ return S, E
382
+
383
+ ```

2

追記を追加しました。

2018/10/19 11:27

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
File without changes
test CHANGED
@@ -235,3 +235,87 @@
235
235
  https://teratail.com/questions/152409
236
236
 
237
237
  https://teratail.com/questions/152518
238
+
239
+
240
+
241
+ ###追記(elseの分岐以降でエラーが出ています)
242
+
243
+
244
+
245
+
246
+
247
+ ```python
248
+
249
+ def Metropolis(Nx=8, Ny=8, iterations=10000):
250
+
251
+ S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数
252
+
253
+ E = update(S)
254
+
255
+
256
+
257
+ indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
258
+
259
+ for k in range(iterations): # モンテカルロステップ数を 10000 とする
260
+
261
+ # ランダムに要素を選択する。
262
+
263
+ index = tuple(indices[np.random.randint(len(indices))])
264
+
265
+
266
+
267
+ # 乱数を生成する。
268
+
269
+ E_tmp = E.copy()
270
+
271
+ E_tmp[index] = np.random.uniform(0, 2 * np.pi)
272
+
273
+
274
+
275
+ # 計算
276
+
277
+ E_tmp = update(E_tmp)
278
+
279
+
280
+
281
+ # 更新
282
+
283
+ if E[index] > E_tmp[index]:
284
+
285
+ continue
286
+
287
+ else:
288
+
289
+ delta = E_tmp[index] - E[index]
290
+
291
+ if random() < np.exp(-delta / 0.001):#変更部分です
292
+
293
+ continue
294
+
295
+
296
+
297
+ return S, E
298
+
299
+
300
+
301
+ ```
302
+
303
+
304
+
305
+ ```output
306
+
307
+ Traceback (most recent call last):
308
+
309
+ File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 44, in <module>
310
+
311
+ S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
312
+
313
+ File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 37, in Metropolis
314
+
315
+ if random() < np.exp(-delta / 0.001):
316
+
317
+ TypeError: 'module' object is not callable
318
+
319
+
320
+
321
+ ```

1

質問の題にミスがあったので修正しました。

2018/10/19 10:15

投稿

astromelt0416
astromelt0416

スコア15

test CHANGED
@@ -1 +1 @@
1
- python NameErrorが解決できない
1
+ python TypeErrorが解決できない
test CHANGED
File without changes