蒙特卡洛模擬XY伊辛模型（python）

?

?

? 前言故事

? 世界上最早的通用電子計算機之一----ENIAC在發明后即被用于曼哈頓計劃，烏拉姆敏銳地意識到在計算機的幫助下，可通過重復數百次模擬過程的方式來對概率變量進行統計估計。馮諾依曼立即認識到這個想法的重要性并給予支持。1947年，烏拉姆提出這種統計方法并應用于計算裂變的連鎖反應。由于烏拉姆常說他的叔叔又在蒙特卡洛賭場輸錢了，因此他的同事Nicolas Metropolis 戲稱該方法為“蒙特卡洛”，不料卻流傳開去。

算法說明

參考文檔：維基百科

1，把握這個算法，首先要體會到隨機的意思。這是一個通過隨機初始化，然后進行優化的一種算法。

2，直接上算法流程，之后再進行解釋。

一，初始化：選擇任意一個x0，假設我們的選擇服從某種概率分布，比如高斯分布或者平均分布。

二，for i from 0 to N? ? ?開始循環迭代

? ? ? ? ? ? ? ? ? ?（1）Generate? 產生一個候選x 。最簡單的方式，隨機一個

? ? ? ? ? ? ? ? ? ? （2）Calculate? 計算接受比 alpha 這里定義的接受比等于 新產生的x在某種分布中的概率/x0在某種分布中的概率。這里的某種分布跟前面的高斯或者平均不是一回事，這里的分布更多的要體現問題的方面。這也是問題在這個算法中唯一的一次介入。

? ? ? ? ? ? ? ? ? ? （3）Accepted?

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(3_1) 產生一個隨機數u,其范圍在[0,1]

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(3_2) if u<=alpha? ?那么候選的x就去掉候選狀況，正式替代

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（3_3）if u>alpha 那么候選的x無效

說明:1,為什么u<alpha才有效替換？

我們可以看到，當p(before)=p(after)時，這時，按照我們的常人觀點，此時變不變是一半對一半！因為兩者的概率一樣，所以發生的概率應該一樣的。但是這里并不是，如果發生了這種事，直接替換。

當p(after)<p(before)時，此時變化后的狀態存在的概率要低于不變的概率，此時并不是一定不變，仍然以一定的概率變。

說明2，第（2）步的計算概率是最跟問題有關的。前面的都是隨機分布變化，跟問題一點關系沒有。請問，問題的概率分布一定存在嗎？如果存在，一般怎么構造這種概率表達式？

我也沒接觸太多。目前，只見到過這個伊辛模型。在熱力學統計中，確實有一個物理量對應這種概率分布是多少，稱為配分函數。這里是配分函數的統計意義解釋。

代碼實現

? import random import matplotlib.pyplot as plt import numpy as np import copy import math import time '''這個就是MCMC模擬，用來模擬某個溫度下XY Ising模型分布， 幾點注意： 注意1，二維伊辛模型，我們用矩陣來模擬。 注意2，旋轉的方向，我們用0到2pi表示吧 算法過程： 一，用一個對稱的分布，高斯分布初始化矩陣 二，下面是循環（1）產生一個候選的自旋方向，為了連續，假設新產生的自旋方向變化最多比原來變化pi/2也就是旋轉90度（2）計算兩個概率，這里熱統中的配分函數正比于概率，因此我們用配分函數Z的比值Z(變化后)/Z(變化前)=exp(-(E后-E前)/kT) ,這里k是玻爾茲曼常數，T是開爾文溫度。令這個值是alpha(3)判斷是都接受*********這個規則有進一步討論*************（3_1)產生一個隨機數u在【0,1】之間（3_2）如果u<=alhpa,接受候選的，改變此時的自旋狀態(3_3)如果u>alpha，不接受候選的，不改變此時的自旋狀態 inputdata: S :matrix 隨機分布,假設已經產生 param: T 溫度delta 最大的變化度數，默認是90度，也可以調整為其他 outputdata:S '''def MetropolisHastings(S,T,numsOfItera):deltamax=0.5k = 1 #玻爾茲曼常數 for sdw in range(numsOfItera): # k2 = np.random.randint(0,S.shape[0]**2)i = np.random.randint(0,S.shape[0])j = np.random.randint(0,S.shape[0])# print('產生的隨機位置是：',i,j)#time.sleep(0.1)for m in range(1): delta = (2*np.random.random()-1)*deltamaxnewAngle = S[i][j]+delta# print(delta)energyBefore = getEnergy(i=i,j=j,S=S,angle=None)energyLater = getEnergy(i,j,S=S,angle=newAngle)alpha = math.exp(-(energyLater-energyBefore)/(k*T))#print(alpha)# if alpha>=1:# print('大于1的哦')if alpha >=1:S[i][j]=newAngleelif np.random.uniform(0.0,1.0)<=1.0*alpha:S[i][j]=newAnglereturn S#計算i,j位置的能量 = 與周圍四個的相互能之和 def getEnergy(i,j,S,angle=None):width = S.shape[0]height = S.shape[1]# print('矩陣的寬和高是',width,height)top_i = i-1 if i>0 else width-1bottom_i = i+1 if i<(width-1) else 0left_j = j-1 if j>0 else height-1right_j = j+1 if j<(height-1) else 0enviroment = [[top_i,j],[bottom_i,j],[i,left_j],[i,right_j]]# print(i,j,enviroment)#print(enviroment)energy = 0if angle ==None:# print('angle==None') for num_i in range(0,4,1): energy += -np.cos(S[i][j]-S[enviroment[num_i][0]][enviroment[num_i][1]])else:# print('Angle')for num_i in range(0,4,1):energy += -np.cos(angle-S[enviroment[num_i][0]][enviroment[num_i][1]])return energy#S=2*np.pi*np.random.rand(30,30) #計算整個格子的能力，為了求平均內能 def calculateAllEnergy(S):energy =0for i in range(len(S)):for j in range(len(S[0])):energy +=getEnergy(i,j,S)averageEnergy = energy/(len(S[0])*len(S))return averageEnergy/2#print(S) #for j in range(1000):# print(j)# MetropolisHastings(S,10) #這個是輸入樣本的多少，格子的尺寸，溫度。中間那個循環，是隨機取迭代的過程 def getWeightValue(numsOfSample,sizeOfSample,temperature):for i in range(numsOfSample): #產生個數print('+++++++正在計算第%s個樣本++++++++++'%i)S=2*np.pi*np.random.random((sizeOfSample,sizeOfSample))initialEnergy = calculateAllEnergy(S)print('系統的初始能量是:',initialEnergy)newS = np.array(copy.deepcopy(S))for nseeps in range(100): newS = MetropolisHastings(newS,temperature,sizeOfSample**2)aveEnergy = calculateAllEnergy(newS)plot(newS)print('系統的平均能量是:',aveEnergy)reshaped = np.reshape(newS,(1,sizeOfSample**2))if i ==0:s = copy.deepcopy(reshaped)continueelse:s = np.row_stack((s,reshaped))return s #運行getweightValue函數，中間已經把結果會成圖了 res = getWeightValue(1,40,2) #print(len(res)) #畫成箭頭圖表示出現 def plot(S): X, Y = np.meshgrid(np.arange(0,S.shape[0]),np.arange(0,S.shape[0]))U = np.cos(S)V = np.sin(S)plt.figure()#plt.title('Arrows scale with plot width, not view')Q = plt.quiver(X, Y, U, V, units='inches') #qk = plt.quiverkey(Q, 0.3, 0.3, 1, r'$2 \frac{m}{s}$', labelpos='E',# coordinates='figure')plt.show()

?

?

結果展示：

?

?

系統值是：

+++++++正在計算第0個樣本++++++++++ 系統的初始能量是: 0.006919117595 系統的平均能量是: -0.494289314266 +++++++正在計算第0個樣本++++++++++ 系統的初始能量是: 0.154688600481 系統的平均能量是: -0.315447987184 +++++++正在計算第0個樣本++++++++++ 系統的初始能量是: 0.00142378493144 系統的平均能量是: -0.588001265121

跟點擊打開鏈接對比，

我們看到在溫度等于2的時候，還是很類似的。

其他溫度的結果：

?

?

總結

以上是生活随笔為你收集整理的蒙特卡洛模拟Ising模型的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。