UA MATH575B 數值分析下 統計物理的隨機模擬方法5

• Ising Model
• Gibbs Sampling
• • Glauber Dynamics

這一講介紹Ising Model，它是MCMC與Gibbs Sampler在統計物理中的一個經典應用之一。

Ising Model

Ising Model是用來描述物質的鐵磁性的。考慮一個$d$維的晶格，假設$\Lambda$是所有晶格點的集合，每一個晶格點有自旋態${\sigma }_i\in \{?1,+1\},?i\in \Lambda$，所有晶格點的自旋態的集合$\sigma =\{{\sigma }_i,?i\in \Lambda \}$是這個晶格的自旋組態。定義第$i$個晶格點與第$j$個晶格點的相互作用參數為$J_{ij}$（$J_{ij}=0$如果兩個晶格不相鄰），外加的磁場對第$j$個晶格的作用是$h_j$，這個晶格的Hamilton量是
$$H(\sigma )=?\sum _{i,j\in \Lambda }{J}_{ij}{\sigma }_i{\sigma }_j?\mu \sum _{j\in \Lambda }{h}_j{\sigma }_j$$
沒有外場時，晶格的Hamilton量是
$$H(\sigma )=?\sum _{i,j\in \Lambda }{J}_{ij}{\sigma }_i{\sigma }_j$$

$J_{ij}>0$，這個晶格具有強磁性（鐵磁性），相鄰電子自旋同向；

$J_{ij}<0$，這個晶格具有反強磁性（反鐵磁性），相鄰電子自旋反向；

$\mu$是晶格點磁矩的值。定義自旋組態的分布為組態幾率，記為$P(\sigma )$，熱平衡下自旋組態服從Boltzmann分布：
$$P_{\beta }(\sigma )=\frac{e^{?\beta H(\sigma )}}{Z_{\beta }}$$
其中$\beta =\frac{1}{kT}$，$k$是Boltzmann常數， $Z_{\beta }$是配分函數：
$$Z_{\beta }=\sum _{?\sigma }{e}^{?\beta H(\sigma )}$$
假設物理量$f(\sigma )$是自旋組態的函數，則它的期望是
$$?f{?}_{\beta }=\sum _{?\sigma }f(\sigma )P_{\beta }(\sigma )$$

理論求解Ising Model比較困難，$d=1、2$的被研究透徹了，但是$d>2$的目前還沒有解析解。數值計算的難度也很大，因為自選組態的可能性是隨晶格的個數指數增加的。所以一般數值上用Monte Carlo方法模擬Ising Model。

Gibbs Sampling

要用Monte Carlo方法模擬Ising Model先要生成$P(\sigma )$的隨機數。生成多元分布的隨機數比生成一元分布的隨機數更復雜得多，但如果每一元的條件分布比較好求的話就可以用Gibbs采樣來做，Gibbs采樣就是通過條件分布sequentially來做采樣，相當于把生成多元分布的隨機數問題又變成了一元的。用Gibbs Sampling模擬Ising Model的物理方法叫做Glauber Dynamics。

Glauber Dynamics

考慮$|\Lambda |=N$，$d=2$的晶格，把晶格點排列在一個網格上，隨機選擇一個晶格點，假設坐標為$(x,y)$：

計算與$(x,y)$相鄰的四個晶格點的自旋之和：
$$S={\sigma }_{x+1,y}+{\sigma }_{x?1,y}+{\sigma }_{x,y+1}+{\sigma }_{x,y?1}$$

如果$(x,y)$的自旋方向反轉，計算能量變化量：
$$\Delta E=2{\sigma }_{x,y}S$$

如果能量變化量小于0，那就反轉一下；如果能量變化量不小于0，那就按如下的概率接受自旋方向反轉：
$$\alpha =e^{?\Delta E/T}$$

重復這個過程超過$N$次。這個算法的作用是模擬晶格自旋組態，接受率的另一種寫法是用Boltzmann分布：
$$\alpha =\frac{e^{?{\sigma }_{x,y}S/T}}{e^{ST}+e^{?ST}}=\frac{e^{?({\sigma }_{x,y}+1)S/T}}{1+e^{?2ST}}$$

下面是我老師課件例子的一個代碼，我們來簡單讀一下。兩個define定義了我們模擬的是$256\times 256$的晶格的2D網格，主函數前四行都是一些變量定義，有兩個比較關鍵的信息：這段代碼模擬的晶格溫度是$T=3$，自旋組態存在二維數組$S$中。$n$和$p$表示的是某個晶格點的相鄰點，網格同一行最左邊與最右邊的晶格點、同一列最上邊與最下邊的晶格點被視為相鄰點，老師說這是某種神秘的量子糾纏。主函數第五行在隨機決定晶格的初始自旋組態。

從主函數第六行開始到倒數第二行，這些代碼就是在做輸出以及做Glauber Dynamics了。這個for循環前三行寫輸出，后面那個for循環第一行是隨機選擇一個晶格點；第二行是計算反轉的能量變化，只是相鄰晶格點的自旋和減2再乘2做了一個標準化，這樣做的好處是讓第三行做接受與否的判斷的時候可以用1和$e^{?\Delta E/T}$來歸一化概率，因為這段代碼用的0表示反向自旋，所以需要這樣做，如果用-1的話就不需要這樣做，可以直接按上面描述的步驟。我們看一下模擬結果的輸出：

這個輸出把晶格點用像素點表示，自旋用1和0表示，對應像素的亮度。從這三輸出我們可以觀察到，隨著溫度升高，相鄰晶格點的自旋狀態就越不那么統一。另外定義
$$M=\frac{1}{N}\sum _{x,y}{\sigma }_{x,y}$$
為磁化強度。我們可以計算磁化強度的自相關函數：

總結

以上是生活随笔為你收集整理的UA MATH575B 数值分析下 统计物理的随机模拟方法5的全部內容，希望文章能夠幫你解決所遇到的問題。

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