質問編集履歴
6
アルゴリズムにミスがあったので更新しました
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
|
-
|
44
|
-
|
45
|
-
|
46
|
-
|
47
|
-
|
48
|
-
|
49
|
-
|
50
|
-
|
51
|
-
|
52
|
-
|
53
|
-
|
54
|
-
|
55
|
-
|
56
|
-
|
57
|
-
|
58
|
-
|
59
|
-
~
|
60
|
-
|
61
|
-
|
62
|
-
|
63
|
-
|
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_sumをS_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
アルゴリズムにミスがあったので、訂正を行いました。
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-1
|
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)を1000
|
64
|
-
|
65
|
-
|
66
|
-
|
67
|
-
(8)(7)の
|
68
|
-
|
69
|
-
|
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
仕様の項目に追記を行いました
test
CHANGED
File without changes
|
test
CHANGED
@@ -64,9 +64,9 @@
|
|
64
64
|
|
65
65
|
|
66
66
|
|
67
|
-
(8)(7)の操作後の行列式S
|
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を更新しました
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
追記を追加しました。
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
質問の題にミスがあったので修正しました。
test
CHANGED
@@ -1 +1 @@
|
|
1
|
-
python
|
1
|
+
python TypeErrorが解決できない
|
test
CHANGED
File without changes
|