仕様
コードにもコメント文で書いていますが、以下に具体的に書きたいと思います。
メトロポリス法という物理学の中の統計力学におけるシミュレーションをを行う際の計算方法です。
20181024 追記メトロポリスアルゴリズムにミスがあったため、全面的に書きなおしました。
初期状態の形成
(1)Nx*Nyの行列を生成し、各成分を(i,j)であらわす。
(2)行列の各要素は0~2*np.piの乱数を取るものとする。この値をS[i,j]とする。
(3)E[i,j]を以下のように定義し、すべてのS[i,j]に対して計算を行う。
E[i,j] =cos(S(i,j)-S(i+1,j))-cos(S(i,j)-S(i-1,j))-cos(S(i,j)-S(i,j+1))-cos(S(i,j)-S(i,j-1)))
しかし、このままでは行列の端では計算ができなくなる(例えばE[1.1]を計算するとき、S[0,1]やS[1,0]が定義されていないので計算ができない)ので、以下のような周期的境界条件を課します。
~~周期的境界条件~~
例えばNx=8,Ny=8の場合、S[0,1]=S[8,1]、S[1,0]=S[1,8]というように、行列の各端で値を折り返します。
(4)行列Sのある要素S(i,j)をランダムに抽出し、その要素だけを0~2πの乱数で変動させる。
(5)変動させたS(i,j)をSal(i,j)として、以下の計算を行い、Eal(i,j)を計算する。
つまり、変動させた(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)変動前の行列Eおよび変動後の行列Eの行列式をそれぞれE、Ealとします。
(6-1)E > Eal の場合、要素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]の要素も当然変化します。
(7)(4)~(6-2)を1回の試行として、1000回繰り返す。
(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
1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3 4 5import numpy as np 6from random import random,randrange 7 8Nx = int(input("Nx=")) 9Ny = int(input("Ny=")) 10 11np.set_printoptions(linewidth=200) 12 13def Initial(): #ランダムな初期状態の形成 14 S = np.random.uniform(0,2*np.pi,(Nx,Ny))#Nx*Nyの行列Sを形成。各要素は0~2πの乱数 15 print(S) 16 17 Einitial = np.empty_like(S)#Sと同じ形の行列Eを作る 18 S_pad = np.pad(S,(1,1),"wrap")#行列Sの両端で折り返す周期的境界条件 19 20 for i, j in np.dstack(np.mgrid[1:Nx+1, 1:Ny+1]).reshape(-1, 2):#Eの各要素(i,j)を計算する 21 Einitial[i - 1, j - 1] = -np.cos(S_pad[i, j] - S_pad[i - 1, j]) \ 22 - np.cos(S_pad[i, j] - S_pad[i + 1, j]) \ 23 - np.cos(S_pad[i, j] - S_pad[i, j + 1]) \ 24 - np.cos(S_pad[i, j] - S_pad[i, j - 1]) 25 26 Einitialsum = np.sum(Einitial) 27 print("Einitial=", Einitial) 28 print("Einitialsum=", Einitialsum) 29 30 31def Metropolis():#メトロポリスアルゴリズム 32 for k in range(10000):#モンテカルロステップ数を10000とする 33 i = randrange(1,Nx-1) 34 j = randrange(1,Ny-1) 35 #ランダムに抽出した要素S(i,j)の値をランダムで決定する 36 Sal = np.random.uniform(0,2*np.pi) 37 38 #変動させたS(i,j)を用いて計算。S(i,j)以外はもともとの値 39 Eal[i - 1, j - 1] = -np.cos(Sal[i, j] - S_pad[i - 1, j]) \ 40 - np.cos(Sal[i, j] - S_pad[i + 1, j]) \ 41 - np.cos(Sal[i, j] - S_pad[i, j + 1]) \ 42 - np.cos(Sal[i, j] - S_pad[i, j - 1]) 43 44 if E[i,j] > Eal[i,j]: 45 S[i,j] = Sal[i,j] 46 E[i,j] = Eal[i,j] 47 else: 48 deltaE[i,j] = Eal[i,j] - E[i,j] 49 if random() < np.exp(-deltaE[i,j]/0.001): 50 S[i,j] = Sal[i,j] 51 E[i,j] = Eal[i,j] 52 53 Send = np.sum(S) 54 Eend = np.sum(E) 55 56 print(Send) 57 print(Eend) 58 59 60 61if __name__ == '__main__': 62 Initial() 63 Metropolis()
↓
output
1Traceback (most recent call last): 2 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 63, in <module> 3 Metropolis() 4 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest9.py", line 42, in Metropolis 5 - np.cos(Sal[i, j] - S_pad[i, j - 1]) 6TypeError: 'float' object is not subscriptable
type Errorが出るのですが、現状解決方法がわかりません。また、メトロポリスアルゴリズムにおいても、Sal[i,j]やE[i,j]のように書いていますがこれで最後まで実行できるかわかりません。pythonに詳しい方、よろしくお願いします。
参考
以前にも以下のような質問をし、本プログラム作成の際にお世話になりました。
https://teratail.com/questions/152409
https://teratail.com/questions/152518
###追記(elseの分岐以降でエラーが出ています)
python
1def Metropolis(Nx=8, Ny=8, iterations=10000): 2 S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数 3 E = update(S) 4 5 indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2) 6 for k in range(iterations): # モンテカルロステップ数を 10000 とする 7 # ランダムに要素を選択する。 8 index = tuple(indices[np.random.randint(len(indices))]) 9 10 # 乱数を生成する。 11 E_tmp = E.copy() 12 E_tmp[index] = np.random.uniform(0, 2 * np.pi) 13 14 # 計算 15 E_tmp = update(E_tmp) 16 17 # 更新 18 if E[index] > E_tmp[index]: 19 continue 20 else: 21 delta = E_tmp[index] - E[index] 22 if random() < np.exp(-delta / 0.001):#変更部分です 23 continue 24 25 return S, E 26
↓
output
1Traceback (most recent call last): 2 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 44, in <module> 3 S, E = Metropolis(Nx=8, Ny=8, iterations=10000) 4 File "C:/Users/user2/Desktop/2DXYTEST/testtesttesttest13.py", line 37, in Metropolis 5 if random() < np.exp(-delta / 0.001): 6TypeError: 'module' object is not callable 7
###追記2(#更新部分について、教えていただいた内容をもとに条件分岐を訂正しました)
python
1def Metropolis(Nx=8, Ny=8, iterations=10000): 2 S = np.random.uniform(0, 2 * np.pi, (Nx, Ny)) # Nx * Ny の行列Sを形成。各要素は0~2πの乱数 3 E = update(S) 4 5 indices = np.dstack(np.mgrid[0:Nx, 0:Ny]).reshape(-1, 2) 6 for k in range(iterations): # モンテカルロステップ数を 10000 とする 7 # ランダムに要素を選択する。 8 index = tuple(indices[np.random.randint(len(indices))]) 9 10 # 乱数を生成する。 11 E_tmp = E.copy() 12 E_tmp[index] = np.random.uniform(0, 2 * np.pi) 13 14 # 計算 15 E_tmp = update(E_tmp) 16 17 # 更新 18 if E[index] > E_tmp[index]: 19 E[index] = E_tmp[index] 20 21 else: 22 delta = E_tmp[index] - E[index] 23 if np.random.random() < np.exp(-delta / 0.001): 24 E[index] = E_tmp[index] 25 26 27 return S, E
回答1件
あなたの回答
tips
プレビュー