MonteCarlo模拟

1、 毕 业 论 文（设计） 题目：Monte Carlo方法研究及应用学 院：数理与信息学院专 业：物理学班 级：B12物理学 号：120603130学生姓名：张夏燕指导教师：钦小平二一六年五月Monte Carlo方法研究及应用摘 要：Monte Carlo方法又称统计实验方法，将事物的特点和物理过程较逼真的描述出来。随着计算机的快速发展，Monte Carlo方法广泛得到应用，在统计物理、反应堆临界安全分析、辐射输运计算、石油测井、燃料循环、金融和航天等领域应用广泛。本文的主要内容安排如下：第一章阐述Monte Carlo方法的发展背景及应用。第二章介绍Monte Carlo方法在具体物理问

2、题中的两个重要算法：Metropolis和Heat bath算法。第三章介绍Monte Carlo方法在统计物理学中的应用：相变、临界现象、临界慢化，介绍平衡态和短时动力学。第四章对Monte Carlo方法进行总论述和展望。关键词：Monte Carlo方法、Metropolis算法、Heat bath算法、铁磁材料相变- - method and application of Monte CarloAbstract:Monte Carlo method, also known as statistical method, can be realistic to describe the c

3、haracteristics of physical experiment process. In recent years, with the rapid development of the computer, the Monte Carlo method has been widely used. In statistical physics、weapons physics、 reactor criticality safety analysis、fuel cycle、radiation transport calculations、the oil well logging、biomed

4、icine、 finance、information and space , Monte Carlo method is widely applied.The main contents of this paper are as follows:：The first chapter describes the present of Monte Carlo methods application in various fields in the results.The second chapter introduces the two important algorithms of Monte

5、Carlo method in compele physics problems: Metropolis and bath-Heat algorithm.The third chapter, the application of Monte Carlo method in statistic physics is introduced: phase transition, critical phenomena, critical slowing down. The equilibrium and short-time dynamics are introduced.In the fourth

6、chapter, the general discussion and Prospect of Monte Carlo method are discussed.Keywords: Monte Carlo method、Metropolis algorithm、Heat - bath algorithm、phase transition 目录目录3引 言3一、Monte Carlo方法的背景与现状51.1发展背景61.2应用现状6二、Monte Carlo方法的重要算法72.1简单抽样和重要抽样22.2Metropolis算法72.3Heat-bath算法8三、Monte Carlo方法的应用

7、123.1磁性材料相变123.2临界慢化123.3短时动力学和平衡态13四、Monte Carlo方法与展望14参考文献15 引 言在摩纳哥公国有这样一座城市，该城市以赌博出名Monte Carlo。因为在赌博游戏中随机性是特别突显的，所以以该城来命名Monte Carlo方法，其实预示了这是一种随机方法。18世纪法国出现一种随机投针实验，设计者是蒲丰，通过该实验可以计算圆周率的值1。普遍认为Monte Carlo方法起源于那时，但正式命名是在“曼哈顿计划”中。追溯到第二次世界大战，美国秘密组织由奥本海默主持的原子弹研制计划-“曼哈顿计划”，原子弹的研制成功至此刷新世界高端武器排行。在研制原子

8、弹的过程中，有一个问题，即中子在不均匀的介质中的散射和输运的行为难以精确求解；根据已知理论可以写出一系列关于中子行为的微分-积分方程，却难以精确求解。研究员放弃精确求解方程组，转而求助于一种基于计算机的模拟实验。利用这种模拟实验解决了当时所遇到的问题，但由于曼哈顿计划的保密性，研究员暂时称这之为Monte Carlo方法，后来这个名称延用至今。Monte Carlo方法是一种数值模拟方法，在概率统计理论为基奠。由于Monte Carlo方法能将物理实验过程和事物的特征描述的较为准确，所以人们越来越广泛地使用该法。Monte Carlo方法的工作原理是对大量的随机样本进行各种统计，从而了解整个系

9、统，得到我们感兴趣的值。基本的思想是，当随机样本足够多时，事件发生的频率数值上等于事件发生的概率。从下面的两个例子中，简单应用Monte Carlo方法。关于圆周率的计算，不必采用蒲丰投针实验，也可以考虑一个更简单的方法。 图0.1（用于计算）如图所示，边长为2r的正方形，半径为r的圆，令前者与后者相互内切。在圆心处为原点水平竖直方向延伸、建立xoy平面，在x.y属于-r,r的区域内，用matlab随机产生M个点。根据计算随机点到圆心的距离判断该点是否属于圆内（上）还是圆外，继而得到点在圆内（上）的频率；另一方面圆面积/正方形面积=/4.令前者与后者在数值上相等，即得到的估算值，随机产生的点数

10、M越多，估计值越接近真实值。同样的思想可以用在计算积分上，考虑如图所示的积分 图0.2（用于计算积分值）在如图所示的区域内随机的产生M个点，根据”f(x)小于等于y“来判断该点是否在阴影内（上），得到点落在阴影部分的频率；当M足够大时，该值应等于阴影面积与正方形面积之比，又因为正方形面积为1，继可得阴影面积即为积分值。上述两个例子都体现了，Monte Carlo方法是先给出大量随机样本，再根据这些样本从而了解系统，继计算我们感兴趣的值。例如上述积分，可以很复杂，复杂到无法笔算得解，但此方法仍有效。因此可以看到，Monte Carlo方法实际上解决了数学上无法解析求解的困境。的确，Monte C

11、arlo方法不是直接解决方程而是绕过去，寻找新的解决方案，大大降低解决难度，在其余领域中同样也可以解决无法突破的瓶颈。本文主要从以下几个方面阐述Monte Carlo方法的研究及应用。第一章主要介绍Monte Carlo方法起源和应用领域所获得成果。第二章主要介绍Monte Carlo方法两种重要的基本算法。第三章主要介绍Monte Carlo方法在统计物理学方面的应用。在讨论相变、临界慢化等处。第四章对Monte Carlo方法的综述和展望。 一、 Monte Carlo方法的起源和现状物理学自18世纪至20世纪以来一直作为自然的主导科学，作出了巨大的贡献，在人类社会亦或是其它科学的发展中。

12、经典的物理研究主要依赖实验和理论。物理学是门很古老的学科，从牛顿定律开始就已经有四百年，但是研究问题的手段不古老。伴随着计算机技术的高速发展，研究物理问题很自然的依靠起计算机，所以计算模拟成为新物理研究手段。计算物理是建立物理问题的模型，框定物理定律后，利用计算机进行演算；如果辅助使用绘图软件还能模拟出直观图像。1.1 起源背景18世纪法国出现一种随机投针实验，设计者是蒲丰，通过该实验可以计算圆周率的值1。常认为Monte Carlo方法即起源于此。 图1.1.1（用于投针实验）投针步骤：1） 在白纸上面画出许多条平行线，令其间距均为d2） 取一根针，长度为L，随机地在纸上掷n次。针与白纸上的

13、平行线相交，观察相交的次数，记为m，得到针与直线相交的事件发生频率3）从概率论上，推得针与直线相交的概率蒲丰本人证明了这个概率是： (1.1.1)试验者时间投掷次数相交次数圆周率估计值Wolf1850500025323.1596Smith1855320412183.1554C.De Morga137Fox188410304893.1595Lazzerini1901340818083.1415929Reina192525208593.1795图1.1.2（投针实验表）类似于投针实验，从概率论的角度理论上推得一个事件发生的概率，再做实验得到从实际操作上得到的该事件发生的频

14、率，当实验次数接近无限时两者数值相等，通过这样来估算一个未知的量，称为Monte Carlo方法。相当一部分问题，用传统解决方法很难解决，但可以用概率模型进行描述。由于概率模型具备不确定的随机要素，那么如果用确定性手段解决难度会很大。特别是有的模型，第一是难以作定量分析，第二是即便有了定量分析但计算代价太大，也无法计算。考虑到这种情况，就需要应用Monte Carlo 方法。二战时期，Monte Carlo 方法又对“曼哈顿计划”做出重大贡献2。Monte Carlo 方法是一种随机的数值实验，而数值实验是基于计算机的一种模拟实验。可以看到，Monte Carlo 方法跳过了精确求解的尴尬，也

15、能得到有重大意义的数据支持原子弹的研制。1.2 国内外的现状从20世纪40年代开始，在快速发展的计算机制造和科技的背景下，Monte Carlo方法被抽出来作为独立的分支，进一步继续得到发展，后来在原子弹研制中首次得到应用。直到80年代，二维算法被美国研究组开发，对晶粒生长的过程重新进行模拟，例如，再结晶、多晶材料的晶粒生长，有序-无序磁畴转变过程等物理学过程。1983年，一个关于Monte Carlo 方法的新算法被Anderson提出，被应用于模拟晶粒生长过程。1992年，Anderson结合晶粒间相互作用，通过Monte Carlo 方法模拟晶粒边界能量和点缺陷对微观结构的影响，最后前者

16、与实验真实值符合的相当理想。1994年，修正后的算法被Zacharia和Radhakrishnan所合作提出，此算法即在模拟晶粒长大的动力学曲线上考虑了模拟时间和实际时间的线性关系。1995年，晶粒尺寸与真实参数之间的关系被他们所得到，是通过使用修正的模型研究晶粒边界的钉扎作用在焊接热影响区的行为。1999年，Sjahanian等人考虑到晶粒边界迁移，进一步研究焊接热影响区晶粒尺寸生长模拟基础3。近段时间以来，在高速发展的电子计算机背景下，随机数值实验被人们开始大范围的用来解决各种问题。另一方面，在不计其数的科学工作者的不懈努力下，Monte Carlo 方法不断完善。因为临界慢化等现象出现，

17、传统Monte Carlo 方法不适用，大量新算法得到提出，例如Cluster算法、Multi-canonical算法、reweighting算法、Hybrid算法、Monte Carlo 重整化群方法等。而后，随着拟随机序列的出现，Monte Carlo 方法已经发展到了Quasi-Monte Carlo 方法，收敛速度和精确度都有提升。目前，三个关于Monte Carlo 方法的研究问题是：（1）产生随机数器，特别是在（0，1）区间上均匀伪随机数的产生。（2）误差问题，下文后提到临界慢化现象，这使得传统Monte Carlo 方法不适用，需要寻找收敛速度更快，精确度也更高的修正方法。（3）

18、应用问题二、 Monte Carlo方法的主要算法2.1 简单抽样和重要抽样考虑到Monte Carlo方法本身就是一种随机抽样的方法，所以对随机数的选取是均匀的还是非均匀的，就有两种分类，即简单抽样7和重要抽样。考虑此积分 图2.1.1（用于计算积分）令 （2.1.1）用Monte Carlo方法简单抽样是这样计算积分：在中均匀地随机地选取，设，所以各点间距为，显然有 （2.1.2）通过分析，可知假如取得点足够多，那么面积元近似为长方形的的准确度足够高，最终结果就足够准确。另一方面，从式（2.1.2）得知，只要M越大，误差越小。误差会即又可以由涨落理论得知，关于S的误差为 均匀取M个点，将所

19、有点所对应的量求和。计算机会取很多个点，想要越准确的答案就需要越长时间。考虑下面这E(x)的图像，再要求其积分。 图2.1.2（用于说明“平均主义”）可以发现，在0,15和118,130区间内，函数值近似于0。简单抽样是在x所属的全区间内均匀随机的取M个点，那么其实取在0,15和118,130区间内的点对应的量对S贡献很小很小，但是又在计算过程中占用了一部分的计算时间，这个问题称为“平均主义”。考虑一下，如果可以这样做：在对S贡献很小很小的区间取少一些的点，在对S贡献大的区间多取一些点，那么即可以减少计算时间又可以保证计算准确度。上述办法就是区别于简单抽样的重要抽样。设积分 （2.1.3）其中

20、就是用来调制的函数，那么我们来按照的分布来得到不均匀分布的，这分布有个直观的特点，那就是：在0,15和118,130区间内的点较少，在其余区间内的点会较多。则 （2.1.4）对照简单抽样和重要抽样的S的形式，是完全一样。按照重要抽样的计算原理可知，它克服了“平均主义”的缺点。但并不是说重要抽样就比简单抽样高级，只是在处理不均匀函数积分时，重要抽样的效率更高；如果是均匀函数积分，可以看见两者效率并没有太大的差别。2.2 Metropolis算法先来介绍Markov过程，再介绍Metropolis算法。2.2.1 Markov过程Markov过程7是一个随机过程。在已知目前的状态下，其之后的演变独

21、立于它以前的演变，这是Markov过程最大特点；像树林里兔子的数量、候车站等车的人数等等都属于Markov过程。虽然对于决定论者这点很难接受，像经典的牛顿力学都是给定初态和受力情况，其将来如何演变就与以往状态关系密切；但是到了微观领域，概率论成为主导，此时Markov过程就容易接受的多些。其实，很多的物理研究过程就是属于Markov过程。下面先来构造一个Markov过程：产生的方法：为了构造某个Markov过程，先得给出一个初态和演变过程所遵循的动力学规则，由随机地产生。那么，给一个，可产生如果随机过程满足某些条件，则当足够大时，即达到“平衡态”时，按分布。这里的“一定条件”，我们在后面的章节

22、会继续介绍。2.2.2 Metropolis算法Metropolis算法8是将markov过程实现于计算机里的重要算法之一，是随机抽样的方法对现实进行数值模拟。需要满足以下两个条件，即（1）各态历经从概率上说，可取任何允许存在的态。（2）细致平衡 Markov过程由从到的转移概率定义设为过程从转移到的概率，为从到的概率。则细致平衡条件为 (2.2.1)Metropolis算法是这样处理的：令 (2.2.2)设 (2.2.3)设 (2.2.4)即满足细致平衡条件。但是，注意到是转移概率，所以应有归一化条件 (2.2.5)上面的还不满足归一化条件。所以，完整的Metropolis算法的转移概率是

23、(2.2.6)其中是从选中的概率，同时。该转移概率仍满足细致平衡条件。同时，归一化条件也能保证。一般系统遵循波尔兹曼能量分布律即 (2.2.7)根据Metropolis算法，则有转移概率为 (2.2.8)此时我们再取一个随机数，将其与转移概率做比较，就可以得到状态是否改变的结论，具体应用在下节展开。2.2.3 Metropolis算法应用Metropolis算法在物理中很常用，下面考虑具体的应用8图2.2.1（二维伊辛模型）图为二维平面铁磁材料，考虑二维伊辛模型：每个格点代表一个粒子，粒子自旋只有两个方向，用表示，方向向上时用“1”表示，方向向下时用“-1”表示。初态如图所示，即高温完全随机态

24、，考虑近邻粒子相互作用，则系统的哈密顿量为 (2.2.9)其中是处的自旋，表示最近邻之和,是近邻作用能参数（1） 选择正中间格点，考虑是否翻转（2） 计算初、末态构型的能量、，继而计算的能量差（3） 如果那么就直接从态变成态，如果，则在之间均匀抽取随机数，若则自旋构型从态变成态，否则仍处于态（4） 完成操作后，换一个格点，重复步骤（2）、（3）说明：每完成一次操作，就形成一个新的自旋构型。用Metropolis算法模拟出系统自旋构型的随机变换过程。 2.3 Heat-bath算法Heat-bath 算法12虽然没有Metropolis算法那么普遍，但也是多体系统研究中的典型方法之一。Heat-

25、bath 算法直接取转移概率为 (2.3.1)这算法自然是各态历经并且满足细致平衡条件的。显然，Heat-bath算法的特点是与无关，也就是与初态无关。取 (2.3.2) (2.3.3)我们认为系统遵循波尔兹曼能量分布律即 (2.3.4)根据Heat-bath算法，则有转移概率为自旋取正的概率为： (2.3.5)自旋取负的概率为 (2.3.6)此时我们再取一个随机数，将其与两个转移概率做比较，就可以得到状态是否改变的结论。2.3.1 Heat-bath算法应用Heat-bath算法在物理中也较常用10图2.2.1为平面铁磁材料，现同样考虑二维伊辛模型初态如图2.2.1所示，图中X表示所需要考虑

26、格点，只考虑近邻粒子相互作用，则哈密顿量为式子2.2.9所示其中是处的自旋，表示最近邻之和，是近邻作用能参数（1） 考虑某一格点：若它自旋向上则计算得系统的能量，若自旋向下则计算得系统的能量（2） 则自旋取正的几率为，式2.3.5；则自旋取负的几率为，式2.3.6（3）在之间均匀抽取随机数，若，则自旋向上，反之自旋向下（4）完成操作后，换一个格点，重复（2）、(3)说明：每完成一次操作，就形成一个新的自旋构型。用Heat-bath算法模拟出系统自旋构型的随机变换过程。2.4 Metropolis和Heat-bath算法与Markov过程的关系很多的物理过程是属于Markov过程的，那么要用计算

27、机解析这些过程首先要能将这个过程模拟出来，Metropolis和Heat-bath算法就是在解决这件事。根据Metropolis算法和Heat-bath算法对转移概率的定义，可以看出，两者都满足上述的两个条件，并且这种随机取样的方式也满足Markov过程的特点。所以，这两套算法是在计算机上实现Monte Carlo方法的第一步大量随机样本。例如上述的二维伊辛模型中所考虑的粒子自旋问题，粒子自旋翻不翻转是个随机的过程，可以利用这两种基本算法，仿真出系统粒子自旋翻转情况。三、 Monte Carlo方法的应用3.1 磁性材料相变磁性材料的相变是一个典型的Markov过程，上一章讲过可以用两种基本算

28、法实现Markov过程，所以磁性材料的相变可以由Monte Carlo方法在计算机上实现。通常，磁性材料处于无序态（无磁性）时是高温，处于又序态（有磁性）时是低温。为了描述磁性材料的有序与无序，我们用一个重要的物理量序参数Q；当系统处于有序态时，即认为Q非零值；当系统处于无序态时，即认为Q为零。同一种方法定义的序参数可能在不同物理系统中作用不同，所以这种定义序参数的方式并不是普适的。例如，在铁磁材料中，常把序参数定义为磁化强度的绝对值M。体系的宏观状态称为相。随着控制参数的改变，比如温度，体系宏观状态可能会发生突变，即称为相变。9接下来介绍用Monte Carlo方法对相变过程进行模拟。考虑一

29、铁磁材料，采用最简单的二维伊辛模型，初态如图所示 图3.1（二维伊辛模型）考虑近邻粒子相互作用，则系统的哈密顿量为 (3.1.1)其中是处的自旋，表示最近邻之和,是近邻作用能参数利用Heat-bath算法，则有(3.1.2) (3.1.3)我们考虑单自旋偏转法，与无关的自旋的能量不必计算。设 (3.1.4) (3.1.5) (3.1.6) (3.1.7)（1） 选定某个格点，计算（2） 产生均匀分布的随机数， 如果，否则（3） 换一个格点，重复步骤（1）、（2）（4） 令T（t）=T(t)+1 ，重复步骤（1）-（3）将各个点都进行扫描，形成一个个新的自旋构型，这是一个由Heatbath算法所

30、模拟出来的关于磁性材料的markov过程。每一个新的自旋构型，可以由平均磁化强度绝对值M定义写出对应的序参量 (3.1.8)一系列的M值来描述磁性材料的有无序，当M为非零值时，认为系统处于有序态；当M为零时，认为系统处于无序态。通过观察这一系列的M值，可以看到相变的过程。3.2 临界慢化 在上文中，我们用单自旋翻转法，一次只考虑一个自旋是否翻转。因此，在某一段演化时间内，各个自旋构型的大部分状态是没有发生改变的。这些大部分相同的自旋构型，几乎可以看成同一自旋构型。假设在 时间内自旋构型没有改变，则称为关联时间。显然，当存在时间关联时，Monte Carlo模拟的误差应当重新估计，修正为这里是抽

31、样总数。当不大时，这对数值模拟的影响并不大。但是，在一些重要而特殊的情形，例如，相变点附近，发散。这会极大地降低Monte Carlo方法的效率，称之为动力学慢化。如果是在二级相变点附近，便是所谓的临界慢化。后来改进Monte Carlo算法，克服临界慢化困难，一直是物理学家和计算科学家的一个重要课题。3.3 短时动力学和平衡态一般，动力学的问题被分成两大类，即（非）平衡态动力学。非平衡态动力学是未达到平衡态的动力学行为11；平衡态动力学，简单说来，是平衡态时的动力学行为或者说是 的动力学行为。研究远未达到平衡态时的行为称为短时动力学，短时动力学研究的是系统动力学的早期行为。虽然临界指数、Z等

32、物理参量是在平衡态中被定义的，但是这些参量也很大程度上主导了动力学的行为。所以，通过研究动力学的早期行为，我们便可以测量这些平衡态的物理参量，预言系统达到平衡态时的行为。四、 Monte Carlo方法的综述与展望这几十年来，Monte Carlo方法不仅应用在数学、物理中，还应用在工程、金融、电网设计等领域。很多物理学的研究对象，常常是一个Markov过程。随着计算机的快速普及与发展，利用计算机解决物理问题是大势所趋。首先需要把所研究对象在计算机里实现，然后进行各种物理量的计算。例如，单自旋翻转的二维Ising模型。模型中每一个新的自旋构型的产生是随机的，是一个典型的Markov过程。已介绍的两种重要算法，Metropolis和Heat-bath算法，能够在计算机上实现这一随机过程。对每一个新的自旋构型计算其序参量，这就在计算机上实现相变过程。但是，例如单自旋翻转这种迭代方式，会造成这样一个问题：接下去的许多新的自旋构型中有一大部分自旋是没有发生改变的。这导致Monte Carlo方法误差值中出现一个量-称为关联时间。不在相变附近时，值趋于0，那么对误差影响就小，可以忽略；在相变附近时，值发散，那么此时误差很大，Monte Carlo方法效