編集履歴

質問編集履歴

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

2018/10/24 13:57

投稿

astromelt0416

スコア15

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -2,8 +2,10 @@
 コードにもコメント文で書いていますが、以下に具体的に書きたいと思います。
 メトロポリス法という物理学の中の統計力学におけるシミュレーションをを行う際の計算方法です。
- **初期状態の形成**
+**20181024 追記メトロポリスアルゴリズムにミスがあったため、全面的に書きなおしました。**
+**初期状態の形成**
 （1）Nx*Nyの行列を生成し、各成分を(i,j)であらわす。
 （2）行列の各要素は０～2*np.piの乱数を取るものとする。この値をS[i,j]とする。
@@ -14,37 +16,22 @@
 ～～周期的境界条件～～
 例えばNx=8,Ny=8の場合、S[0,1]=S[8,1]、S[1,0]=S[1,8]というように、行列の各端で値を折り返します。
-**メトロポリスアルゴリズム（10/24追記：(6)以下で一部間違っています。下に正しいものを載せました。）**
 （4）行列Sのある要素S（i,j)をランダムに抽出し、その要素だけを0～2πの乱数で変動させる。
 （5）変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
-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))
-つまり、変動させた(i,j)以外の要素はそのままにして、手順(3)で行ったものと同じ計算をするということです。
+つまり、変動させた(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]も変動することに注意してください。
-~~(6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
+（6）変動前の行列Eおよび変動後の行列Eの行列式をそれぞれE、Ealとします。
-~~
-~~(6-1）E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
+(6-1)E > Eal の場合、要素S[i,j]→Sal[i,j]への更新に成功したとし、(7)以降の過程へ。
-~~
-~~(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]に更新します。
-~~
-~~(7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
-~~
-~~(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
-参考程度にですが、Ssum～64、Esum～-250　(Nx=8、Ny=8) 　となれば成功と言えます。~~
-~~
-**10/24追記　(6)以下のメトロポリスアルゴリズム（文献を再度確認しなおしました。）**
-(6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
-(6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
+(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]の要素も当然変化します。
-(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]に更新します。
+(7)(4)~(6-2)を1回の試行として、1000回繰り返す。
-(7)(4)~(6-2)を1000回繰り返す。(ここでは(8)でしているような計算は行わず、Sの更新だけを行います。)
+(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)についても同様です。
-(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)についても同様です。
 (9)(8)で求めた値を用いて、S_ave = S_end(10000)/10000,E_ave = E_end(10000)/10000 を出力できればOKです。

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

2018/10/24 13:57

投稿

astromelt0416

スコア15

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -1,5 +1,5 @@
 ### 仕様
-コードにもコメント文で書いていますが、以下に具体的に描きたいと思います。
+コードにもコメント文で書いていますが、以下に具体的に書きたいと思います。
 メトロポリス法という物理学の中の統計力学におけるシミュレーションをを行う際の計算方法です。
  **初期状態の形成**
@@ -14,7 +14,7 @@
 ～～周期的境界条件～～
 例えばNx=8,Ny=8の場合、S[0,1]=S[8,1]、S[1,0]=S[1,8]というように、行列の各端で値を折り返します。
-**メトロポリスアルゴリズム**
+**メトロポリスアルゴリズム（10/24追記：(6)以下で一部間違っています。下に正しいものを載せました。）**
 （4）行列Sのある要素S（i,j)をランダムに抽出し、その要素だけを0～2πの乱数で変動させる。
@@ -23,17 +23,35 @@
 つまり、変動させた(i,j)以外の要素はそのままにして、手順(3)で行ったものと同じ計算をするということです。
+~~(6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
+~~
+~~(6-1）E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
+~~
+~~(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]に更新します。
+~~
+~~(7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
+~~
+~~(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
+参考程度にですが、Ssum～64、Esum～-250　(Nx=8、Ny=8) 　となれば成功と言えます。~~
+~~
+**10/24追記　(6)以下のメトロポリスアルゴリズム（文献を再度確認しなおしました。）**
 (6)ここで(3)で計算したE[i,j]と(5)で計算したEal[i,j]とを比較します。
-(6-1）E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
+(6-1)E[i,j] > Eal[i,j] の場合、もともとの要素S[i,j]をSal[i,j]に更新し、(7)以降の過程へ。
 (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]に更新します。
-(7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
+(7)(4)~(6-2)を1000回繰り返す。(ここでは(8)でしているような計算は行わず、Sの更新だけを行います。)
-(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
-参考程度にですが、Ssum～64、Esum～-250　(Nx=8、Ny=8) 　となれば成功と言えます。
+(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)についても同様です。
+(9)(8)で求めた値を用いて、S_ave = S_end(10000)/10000,E_ave = E_end(10000)/10000 を出力できればOKです。
 ### 実際に書いてみたコード
 ```python

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

2018/10/24 08:40

投稿

astromelt0416

スコア15

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -31,9 +31,9 @@
 (7)(4)~(6)を10000回繰り返すことでおおよそすべての要素S[i,j]がスピン更新を経験するようにします。
-(8)(7)の操作後の行列式Send=np.sum(S)及びEend=np.sum(E)を出力できればこのシミュレーションは終了です。
+(8)(7)の操作後の行列Sを用いて行列Eを計算する。それぞれの行列の行列式の値Ssum＝np.sum(S)、Esum=np.sum(E)をもとめることができれば計算終了です。
+参考程度にですが、Ssum～64、Esum～-250　(Nx=8、Ny=8) 　となれば成功と言えます。
 ### 実際に書いてみたコード
 ```python

追記２を更新しました

2018/10/19 22:34

投稿

astromelt0416

スコア15

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -158,4 +158,35 @@
     if random() < np.exp(-delta / 0.001):
 TypeError: 'module' object is not callable
+```
+###追記２（#更新部分について、教えていただいた内容をもとに条件分岐を訂正しました）
+```python
+def Metropolis(Nx=8, Ny=8, iterations=10000):
+    S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は０～２πの乱数
+    E = update(S)
+    indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
+    for k in range(iterations): # モンテカルロステップ数を 10000 とする
+        # ランダムに要素を選択する。
+        index = tuple(indices[np.random.randint(len(indices))])
+        # 乱数を生成する。
+        E_tmp = E.copy()
+        E_tmp[index] = np.random.uniform(0, 2 * np.pi)
+        # 計算
+        E_tmp = update(E_tmp)
+        # 更新
+        if E[index] > E_tmp[index]:
+            E[index] = E_tmp[index]
+        else:
+            delta = E_tmp[index] - E[index]
+            if np.random.random() < np.exp(-delta / 0.001):
+                E[index] = E_tmp[index]
+    return S, E
 ```

追記を追加しました。

2018/10/19 11:27

投稿

astromelt0416

スコア15

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -116,4 +116,46 @@
 ### 参考
 以前にも以下のような質問をし、本プログラム作成の際にお世話になりました。
 https://teratail.com/questions/152409
-https://teratail.com/questions/152518
+https://teratail.com/questions/152518
+###追記（elseの分岐以降でエラーが出ています）
+```python
+def Metropolis(Nx=8, Ny=8, iterations=10000):
+    S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は０～２πの乱数
+    E = update(S)
+    indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2)
+    for k in range(iterations): # モンテカルロステップ数を 10000 とする
+        # ランダムに要素を選択する。
+        index = tuple(indices[np.random.randint(len(indices))])
+        # 乱数を生成する。
+        E_tmp = E.copy()
+        E_tmp[index] = np.random.uniform(0, 2 * np.pi)
+        # 計算
+        E_tmp = update(E_tmp)
+        # 更新
+        if E[index] > E_tmp[index]:
+            continue
+        else:
+            delta = E_tmp[index] - E[index]
+            if random() < np.exp(-delta / 0.001):#変更部分です
+                continue
+    return S, E
+```
+↓
+```output
+Traceback (most recent call last):
+  File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 44, in <module>
+    S, E = Metropolis(Nx=8, Ny=8, iterations=10000)
+  File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 37, in Metropolis
+    if random() < np.exp(-delta / 0.001):
+TypeError: 'module' object is not callable
+```

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

2018/10/19 10:15

投稿

astromelt0416

スコア15

title CHANGED Viewed

	@@ -1,1 +1,1 @@
1	- python ~~NameError~~が解決できない
1	+ python TypeErrorが解決できない

body CHANGED Viewed

File without changes