馬爾可夫鏈理論和Monte Carlo 取樣的實現(xiàn)

1、馬爾可夫鏈理論和Monte Carlo 取樣的實現(xiàn) Monte Carlo 取樣不直接使用 P(X) ,而是以某種方式取樣, 大量的樣本最終符合所需的分布. Monte Carlo利用轉(zhuǎn)移矩陣 , 從當(dāng)前 x 生成一個狀態(tài) y . 馬爾可夫鏈和Monte Carlo 馬爾可夫鏈?zhǔn)且粋€簡單的隨機過程, 給定轉(zhuǎn)移矩陣 W,可以得到平衡分布 P. Monte Carlo 是馬爾可夫鏈的計算機實現(xiàn). Monte Carlo中, P 給定, 我們需要尋找 W 使得 P = P W. 各態(tài)歷經(jīng) 對所有 n nmax, 所有 x 和 x 細(xì)致平衡馬爾可夫鏈的收斂條件(充要)證明考慮很多個平行的Markov鏈

2、, 在一個給定的某一步, 有Nr個鏈處于第r個態(tài), Ns個鏈處于第s個態(tài). 于是在下一步從r態(tài)到s態(tài)的數(shù)目為從s態(tài)到r態(tài)的數(shù)目為從r態(tài)到s態(tài)的凈轉(zhuǎn)移的數(shù)目為若w( xr! xs)滿足細(xì)致平衡條件, 則上式成為這是一個十分重要的結(jié)果, 上式表明, 如果二個狀態(tài)之間不滿足分布P, 則這一Markov 過程的演化結(jié)果將總是使其趨于滿足. 這樣, 就證明了我們的論斷.Metropolis 算法 (1953) Metropolis 算法取 T 為一個對稱的轉(zhuǎn)移矩陣 The Paper (7500 citations from 1988 to 2003)THE JOURNAL OF CHEMICAL PH

3、YSICS VOLUME 21, NUMBER 6 JUNE, 1953Equation of State Calculations by Fast Computing MachinesNICHOLAS METROPOLIS, ARIANNA W. ROSENBLUTH, MARSHALL N. ROSENBLUTH, AND AUGUSTA H. TELLER,Los Alamos Scientific Laboratory, Los Alamos, New MexicoANDEDWARD TELLER, * Department of Physics, University of Chic

4、ago, Chicago, Illinois(Received March 6, 1953) A general method, suitable for fast computing machines, for investigating such properties as equations of state for substances consisting of interacting individual molecules is described. The method consists of a modified Monte Carlo integration over co

5、nfiguration space. Results for the two-dimensional rigid-sphere system have been obtained on the Los Alamos MANIAC and are presented here. These results are compared to the free volume equation of state and to a four-term virial coefficient expansion.1087The Calculation Number of particles N = 224 M

6、onte Carlo sweep 60 Each sweep took 3 minutes on MANIAC Each data point took 5 hoursMANIAC the Computer and the ManSeated is Nick Metropolis, the background is the MANIAC vacuum tube computer正則分布的抽樣方法正則分布的抽樣方法:選擇一個滿足細(xì)致平衡條件的轉(zhuǎn)移幾率;產(chǎn)生一個Markov 鏈, 丟掉鏈的前而面M個狀態(tài);用其余狀態(tài)進(jìn)行物理量的計算. 考慮從r態(tài)到s態(tài)的轉(zhuǎn)移, 若二狀態(tài)的能量差為則:當(dāng)年Metro

7、polis 選擇 :目前常用的另一種選擇是:應(yīng)當(dāng)注意的是, w的選擇并不唯一, 只要滿足細(xì)致平衡條件的要求即可, 但不同的w收斂速度往往差別很大, 如何選擇合適的w以達(dá)到盡可能快的收斂速度和盡可能高的計算精度仍然是當(dāng)前Monte Carlo算法研究的前沿課題之一.The Ising Model-+- - -The energy of configuration s isE(s) = - J si sjwhere i and j run over a lattice, denotes nearest neighbors, s = 1s = s1, s2, , si, 例題例題, Ising模型的

8、模擬模型的模擬Ising 模型:式中J稱為交換積分, h為外場, si 可取值(1, -1), 稱為自旋變量. Ising 模型是最簡單的非平庸統(tǒng)計物理模型, 它是由德國物理學(xué)家 Lenz 在二十年代提出的, 這一模型可用來描述單軸各向異性磁性系統(tǒng), 合金等物理體系, 同時也是一個十分有興趣的理論模型.例題例題, Ising模型的模擬模型的模擬Ising 最早給出了這一模型在一維情況下的嚴(yán)格解, 證明了在一維下這一模型不存在相變. Onsager 于1944 年做出了零場下這一模型在二維空間的嚴(yán)格解并計算了它的相變溫度, 比熱在相變點的行為等熱力學(xué)量. 楊振寧在1952 年解出了外場很小時二維

9、空間的 Ising 模型, 求出了序參量的臨界行為. 由于對這一模型的很多形為目前了解的比較透徹, 因此它經(jīng)常被用來做為檢驗各種數(shù)值方法或解析近似方法的標(biāo)準(zhǔn). 1969年，A. E. Ferdinand和M. E. Fisher求得了有限尺寸2D Ising 模型在周期性邊界條件下的嚴(yán)格解，成為檢驗?zāi)M結(jié)果的一個有效標(biāo)準(zhǔn)。感興趣的物理量平均能量 比熱，利用下面的公式:3. 磁化 = 磁化率 5.Binder 4 階累積量 自旋相關(guān)函數(shù)時間相關(guān)函數(shù)2D Ising 模型的比熱From D P Landau, Phys Rev B 13 (1976) 2997.一個算法一個算法 選擇一個格點 i,

10、 其自旋將考慮作翻轉(zhuǎn) si! -si. 計算與此翻轉(zhuǎn)相聯(lián)系的能量變化 H. 計算這一翻轉(zhuǎn)的轉(zhuǎn)移幾率 w. 產(chǎn)生一在 0,1 之間均勻分布的隨機數(shù) . 如果w, 則翻轉(zhuǎn)該自旋, 否則, 保持不變. 不論何種情況, 其結(jié)果都作為一新的狀態(tài). 分析該狀態(tài), 為計算平均值收集數(shù)據(jù).有限尺寸標(biāo)度與相變漲落:h A2i 通常與某種響應(yīng)函數(shù)相聯(lián)系: 漲落耗散定理!若A(x)為系統(tǒng)的哈密頓量, 則平均每個粒子的能量 hEi為:比熱是一個強度量, 由此可以推斷:自平均效應(yīng). 當(dāng)系統(tǒng)的尺度趨于無限時, 其漲落趨于0!有限尺寸標(biāo)度與相變這實際上是中心極限定理的一個結(jié)果. 考慮把系統(tǒng)分成一些小的系統(tǒng), 小系統(tǒng)之間的相

11、互作用可以忽略時(這總可以做到?!), 就是一系列獨立, 同分布的樣本, 從而中心極限定理成立.由此, 能量的分布可以寫為:P(E), N=1, 5, 10, 20 , 30h E i =1KB T2 c=1分布寬度 / 1/N1/2Ising 模型:平均磁化:自發(fā)磁化:注意求極限的順序!有限尺寸標(biāo)度與相變由于Ising模型的哈密頓量在h=0時具有 Si $ Si 的對稱性, 所以,如果先取h=0, 則得不到自發(fā)磁化. 自發(fā)磁化是熱力學(xué)極限下的產(chǎn)物, 由各態(tài)歷經(jīng)破壞而得到.有限尺寸標(biāo)度與相變由于:磁化率:有限尺寸標(biāo)度與相變由此:基于在討論能量時的同樣的分析, 我們得到序參量(磁化)的分布.當(dāng)h

12、=0時, 由于對稱性, M和-M是對稱的. 在臨界溫度以上, M=0, 所以當(dāng)TTc 時, 為了反映這種對稱性, 分布函數(shù)寫為:有限尺寸標(biāo)度與相變P(S), N=1, 5, 10, 20 , 30h Si =0.8KBT=1有限尺寸標(biāo)度與相變?nèi)绾斡嬎阈騾⒘?注意到前面的分析, 對于有限系統(tǒng), 如果直接計算 h S i, 則在任何溫度都得到的是0!注意到自平均效應(yīng), 有下述結(jié)果:可以使用上述任何一種進(jìn)行計算, 然后外推到N ! 1.有限尺寸標(biāo)度與相變磁化隨MCS的變化，有限系統(tǒng)的各態(tài)歷經(jīng)演示。2D Ising 模型，L64有限尺寸標(biāo)度與相變相變點會發(fā)生什么?在二級相變點(現(xiàn)一般稱為臨界點, 二級

13、相變也稱為連續(xù)相變), 系統(tǒng)的比熱, 磁化率等物理量發(fā)散, 發(fā)散的根源在于相關(guān)長度發(fā)散. 這樣, 前面的分析將不再成立. 在臨界點, 熱力學(xué)量具有很有趣的標(biāo)度行為: 以單軸各向異性磁體為例(理論上近似以Ising模型描述). 在T! Tc 時, 引進(jìn)約化溫度:有限尺寸標(biāo)度與相變比熱: 磁化率: 自發(fā)磁化: 關(guān)聯(lián)長度: 當(dāng)t=0時, 磁化與外場之間有:有限尺寸標(biāo)度與相變 , , , , 稱為臨界指數(shù), 是反映臨界點本質(zhì)的量. 這些指數(shù)的來源與關(guān)聯(lián)長度的發(fā)散有關(guān). 還有一個指數(shù), 這里不做介紹.當(dāng)系統(tǒng)有限時, 不會有發(fā)散! Monte Carlo只能計算有限體系, 如何從Monte Carlo計算

14、獲得相變和臨界現(xiàn)象的信息?有限尺寸標(biāo)度與相變?nèi)绾握遗R界點?如何算臨界指數(shù)指數(shù)?為此, 我們考慮有限體系的標(biāo)度理論.基本出發(fā)點基本出發(fā)點: 設(shè)有限系統(tǒng)的線度為L, 則關(guān)聯(lián)長度最多為L, TC 時的序參量分布序參量的有限尺寸效應(yīng)有限尺寸標(biāo)度與相變磁化率的有限尺寸效應(yīng)T TC 時的序參量分布TTC 時的序參量分布考慮磁化率: 在臨界點附近有一個峰, 最高處的溫度定義為有限系統(tǒng)的臨界溫度Tc(L), 峰寬為 T. 當(dāng) L ! 1 時, Tc(L)! Tc(1), T! 0假定: 磁化率的極大值: 磁化在Tc(L)不為0, 而是:Fisher 假定假定: 在臨界點只有一個重要的特征長度, 即關(guān)聯(lián)長度.有

15、限尺寸標(biāo)度與相變由于 , 磁化率峰的寬度對應(yīng)于 , 于是 有限尺寸標(biāo)度與相變由此得到:Finite-Size ScalingShift of TcFrom A M Ferrenberg and D P Landau, Phys Rev B 44 (1991) 5081Accurate Exponent RatioFrom J S Wang, R H Swendsen, and R Koteck, Phys Rev. B, 42 (1990) 2465.有限尺寸標(biāo)度與相變Binder 累積量:性質(zhì)1: T , 遠(yuǎn)離臨界區(qū)域有限尺寸標(biāo)度與相變性質(zhì)2: T Tc, L, 遠(yuǎn)離臨界區(qū)域有限尺寸標(biāo)度與相變性質(zhì)3: T Tc, L, 在臨界區(qū)域U* 是一個與L關(guān)系很小的普適值.應(yīng)用: 確定Tc, 對不同的L計算UL, 對各種對(L, L), 計算比值 UL/UL, 并對T作圖, 這些比值的曲線將這些比值的