二维伊辛模型的蒙特卡罗数值模拟

1、. 毕 业 论 文题目：二维伊辛模型的蒙特卡罗数值模拟二一一年五月目 录摘要.1abstract.2 第一章 引言.3第二章 伊辛模型.42.1 伊辛模型的意义.42.2 伊辛模型的历史与发展.42.3 二维伊辛模型的基本结构.5第三章 蒙特卡罗方法.63.1 蒙特卡罗方法的产生与发展.63.2 蒙特卡罗方法的基本思想.6 3.3 文中采用蒙特卡罗方法时用到的的基本理论.7第四章 模拟的理论和过程.8 4.1 计算理论.8 4.2 模拟过程.9第五章 数据分析和讨论.10参考文献.14附录.15致谢.19.二维伊辛模型的蒙特卡罗数值模拟摘要：本文主要介绍应用蒙特卡罗方法对二维伊辛模型进行数值模

2、拟的基本思路和方法，然后具体应用fortran语言进行数值模拟的过程。文中首先介绍了伊辛模型对解决相变问题的意义及其历史与发展，然后介绍了蒙特卡罗方法的产生与发展及其基本内容和思想，最后具体应用计算机fortran语言用蒙特卡罗方法对二维伊辛模型的相关物理量进行数值模拟的方法和过程，并对所得数据进行分析和讨论。关键字：二维伊辛模型；蒙特卡罗方法；哈密顿量；磁化强度；fortran语言 the simulation of monte carlo for the two dimensional ising modelabstract：this paper describes the monte c

3、arlo method to simulate two dimensional ising model of the basic ideas and methods, and then specific applications fortran language simulation process. the paper first describes the ising model to solve the problem of phase change the meaning and history and development, and then introduced the mont

4、e carlo method and its basic contents and development of production and ideas, and finally specific applications the computer language fortran and monte carlo method to simulate the related physical of the two dimensional ising model of method and process, and the data were analyzed and discussed.ke

5、ywords：the two dimensional ising model; monte carlo method; hamiltonian; fortran language第一章 引言统计物理学中，物质相变方面的研究是一个重要的领域，物质经过相变，要出现新的结构和物性，而得到物质的相变过程是非常重要且意义深远的，伊辛模型就是试图去解决这些问题的一个重要模型。然而，应用蒙特卡罗方法去作用于伊辛模型，又是一个重要的思想，并且确实可以模拟一些过程，得到一些重要结论。本文主要是以蒙特卡罗方法为基础，应用计算机程序来模拟二维伊辛的一些基本物理量，并结合理论，对所得数据进行分析，完成最基础的学习研究

6、。伊辛模型和蒙特卡罗方法都是非常重要而且被广泛使用的，这里，通过对二维伊辛模型的蒙特卡罗数值模拟，可以了解甚至解决热力学统计物理中，相变领域的一些问题，它可以简化这些问题的处理方法，并求出精确地解，从而可以对现在得到广泛关注的三维伊辛模型的精确求解的理论打下基础进而可以对其进行思考，如果有可能，可以更深入理解二维伊辛模型的精确求解，并去了解三维伊辛模型的相关理论。第二章 伊辛模型伊辛模型是描述物质相变的一种模型。物质经过相变，要出现新的结构和物性。发生相变的系统一般是在分子之间有较强相互作用的系统，又称合作系统。它是一个简单且被广泛使用的物理模型，下面就对从伊辛模型的意义及其历史与发展两方面对

7、其做简要的介绍，并且简单介绍二维伊辛模型的基本结构。2.1 伊辛模型的意义“伊辛 (ising)模型是一个最简单且可以提供非常丰富的物理内容的模1,它可以被用来帮助我们发现物理世界的原则。它不仅可以用来描述晶体的磁性,还可以用来描述非常广泛的现象,如合金中的有序无序转变、液氦到超流态的转变、液体的冻结和蒸发、晶格气体、玻璃物质的性质、森林火灾、城市交通、蛋白质分子进入它们的活性形式的折叠等。科学家对伊辛模型的广泛兴趣还源于它是描述相互作用的粒子 (或自旋 )最简单的模型。它可以用来测试研究相互作用的粒子在多体系统 (特别是理解在临界点及其附近的合作现象和临界行为)任何近似方法的理想工具。”2以

8、上这段文字正是对伊辛模型的意义的最好诠释，而且三维伊辛模型的意义更为重大，它能解决更为广泛的问题，最为重要的应用是它能很好地显示连续相变过程,特别是在相变的临界温度附近的临界现象。对伊辛模型的研究，使其成为成熟的模型，对于许多领域都有重大意义。2.2 伊辛模型的历史与发展统计物理学通过对配分函数求导得系统的热力学函数，从而确定系统的全部平衡性质。但是在相变点某些热力学量会发生突变或者出现无穷尖峰，那么从单一的配分函数表达式能否同时描述各种相和相的转变呢？回答这个疑难的一个方法是建立包含系统最本质特征的简化模型，严格导出其在相变点的宏观特性。从20世纪30年代中叶开始，经过半个世纪以来对统计模型

9、的大量研究，已形成统计物理学的一个专门领域，建立了很多模型，伊辛模型就是其中的一个模型。经过数位科学家的研究和探讨，伊辛模型经过几十年的发展，有了很大的突破，现在，一维，二维伊辛模型理论已经比较成熟，科学家已经得他们在一些条件下的精确解，但是二维伊辛模型在其他一些条件下还有得到精确解，比如具有次紧邻相互作用的二维伊辛模型和磁场下的伊辛模型至今没有精确解，而且三维伊辛模型也没有得到精确解，还没有达到成熟的阶段，还处于探索和研究中。以上只是一个简单的介绍，对伊辛模型的详细介绍，可以参考文献2,这是一篇具有代表性的论文，它对伊辛模型做了完整而详细的介绍，对于了解伊辛模型，很有参考价值。2.3 二维伊

10、辛模型的基本结构“伊辛模型是一个非常简单的模型1,在一维、二维或三维的每个格点上占据一个自旋。自旋是电子的一个内禀的性质,每个自旋在空间有两个量子化的方向,即其指向可以向上或向下。”2如上所说，一个n维的伊辛模型，在每个格点上占据一个自旋，指向向上或向下，并与周围的自旋相互作用，伊辛模型只考虑最近邻自旋的作用，因而每个自旋与周围自旋相互作用的强弱与周围最近邻子自旋数有关，即随着维度的增加,每个自旋的最近邻自旋数增加,与周围自旋的相互作用也在增强2。下面简要介绍正方二维伊辛模型的基本结构。如图1，一个正方格子，每个格点代表一个自旋（向上(+1)或向下(-1)），仅与周围上下左右最近邻的四个自旋有

11、相互作用，这就是，一个简单的正方二维伊辛模型。在求解一些物理量或者解决一些问题时，就以这个模型为基础，根据需要确定格子大小，初始自旋指向，然后以仅考虑最近邻相互作用为基础，根据需要，设定条件，处理问题。这就是最基本的思想。 图1其中空心圆点代表自旋向上，实心圆点代表自旋向下。 第三章 蒙特卡罗方法monte-carlo(蒙特-卡罗,又译蒙特-卡洛)模拟法,是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统

12、计特征，故借用赌城蒙特卡罗命名。 下面就蒙特卡罗方法的产生与发展及其基本思想以及文中所涉及到的相关理论。3.1 蒙特卡罗方法的产生与发展蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员s.m.乌拉姆和j.冯诺伊曼首先提出。数学家冯诺伊曼用驰名世界的赌城摩纳哥monte carlo来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777年，法国数学家布丰（georges louis leclere de buffon，17071788）提出用投针实验的方法求圆周率。这被认为是蒙特卡罗方法的起源和最初应用。由于蒙特卡罗方法的计算量很

13、大，使得它起初得不到广泛应用，但随着20世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟庞大且复杂的试验成为可能，从而蒙特卡罗方法被重新提起并得到重视，随着时间的推移，它已成为一种被广泛使用的方法。3.2 蒙特卡罗方法的基本思想蒙特卡罗方法的基本思想是当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件

14、的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗方法正是基于此思路进行分析的。蒙特卡罗方法回避了结构可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多，就可得到一个比较精确的失效概率和可靠度指标，于是蒙特卡罗方法为很多以前难以解决和分析的问题提供了一种思想，使他们得到解决和进展，是蒙特卡罗方法称为一种基础而重要的思想。3.3 文中采用蒙特卡罗方法时用到的的基本理论在伊辛模型的蒙特卡罗数值模拟方法中，会用到metropolis准则和马

15、尔可夫链的一些基本内容，在这里简单介绍下。metropolis3准则，是粒子在温度t时趋于平衡的概率为 ，其中e为温度t时的内能，e为其改变量，k为boltzmann常数。这是在进行蒙特卡罗循环时，应用metropolis抽样技术的条件，根据概率 的大小，决定格点上的自旋状态是否发生变化。马尔可夫链，因安德烈马尔可夫（a. a. markov，18561922）得名，是数学中具有马尔可夫性质的离散时间随机过程。该过程中，在给定当前知识或信息的情况下，过去（即当期以前的历史状态）对于预测将来（即当期以后的未来状态）是无关的。在数值模拟过程中，随机选取一个自旋状态以后，构造出马尔可夫链，由metr

16、opolis转移概率决定马尔可夫链过程进行的方向3,4,5，进行蒙特卡罗循环,从而完成模拟过程。mcmc方法，即马尔可夫链蒙特卡罗(markov chain monte carlo)方法，和metropolis蒙特卡罗原理，是本文数值模拟用到的俩个基本原理方法。对它们的详细论述与分析请参考见本文末尾的参考文献。第四章 模拟的理论和过程下面就从一个二维伊辛模型开始，先从理论方面简要介绍介绍数值模拟中的理论思想，及涉及到的物理量，再在此基础上，介绍模拟过程，即用蒙特卡罗方法模拟伊辛模型的步骤，这里就会用到前面提到的metropolis准则和mcmc方法。4.1 计算理论 如图25,是某一线度为n的

17、方形二维伊辛模型的某一格点的自旋（s0）及其周围最近邻的四个自旋（s1,s2,s3,s4）,格子上伊辛自旋模型自旋取向分布的一个状态， 或叫自旋组态，si取值为+1或-1，即自旋向上或向下，称为伊辛变量4。由于只考虑每个自旋周围最近邻的四个自旋，那么si与其周围最近邻自旋sj相互作用能为5，那么伊辛模型系统的哈密顿量为： （1）其中j为耦合常数，与物质有关。同时系统的磁化强度为： （2） 图2 二维方形伊辛模型的某一格点的自旋及其最近邻自旋 在给定的温度t下，伊辛模型系统从初始状态开始，在metropolis准则和mcmc方法二者的条件下，经过一定时间，最终达到相应温度下的平衡态，进而可以根据

18、进行相应物理量的数值模拟。这就是计算理论的基本思想4.2 模拟过程及步骤模拟过程主要分三个部分，第一部分为设定初始状态,如温度t，耦合常数j，各格点的自旋状态（可以随机，随机更接近实际情况），以及格子线度n等数据；第二部分为模拟过程，随机选定的某个自旋状态，在metropolis准则和mcmc方法的判断下，确定是否改变现有状态；第三部分为循环过程，即重复第二部分的内容，得到足够多的数据，求解磁化强度的平均值，分析哈密顿量的变化。重复以上三部分内容，即可得到不同初始值下的物理量的值。 具体模拟步骤4,5为：(1) 给定初始值，温度t，耦合常数j，格子线度n，这里只列出主要的数据变量，对于其他的物

19、理量，在数据分析中在具体设定。(2) 随机选取二维伊辛模型的各格点的初始状态，算出系统初始能量。(3) 随机选取一个自旋状态，构造马尔科夫链，求解翻转后的系统能量，然后求解能量差,计算metropolis转移概率函数的值，转移概率函数为： （3） 设,其中r称为接受率3。于是转移概率变为： (4)上式确定马尔可夫链的行进方向，如果出于边界上，则采用周期性边界条件9。（4）当时，直接发生偏转；当时，需要系统产生随机数x来决定 是否发生偏转，x在0,1之间，当wx ,发生偏转，当wx , 不发生偏转，但这时它已是新的状态。用接受率r来表示是否翻转的条件3要更为简单，当或者时，发生偏转，否则不发生偏

20、转。编程时用这种思路，会简化问题。（5）重复（2）到（4）步，当进行足够多次后，是系统达到平衡状态，系统愈大，所需时间愈长。（6）根据数据求得所需物理量，改变初始值，重复上述过程，结合每种初始值计算所得的物理量的数据，根据数据作图，用图分析讨论，得出结论。第五章 数据分析和讨论 在这一章中将对数值模拟过程中的物理量进行具体设定和说明，然后对所得的物理量的数据进行作图，再分析数据。本文进行数值模拟的计算机语言为fortran语言，数值模拟的程序见附录。其中需要设定的物理量值有（可以参见参考文献8）： 1.耦合常数j,它与物质的性质有关，在模拟时，设定j=0.5。2.二维伊辛模型的线度，即n,n与

21、模拟效果有关，是一个重要的变量，模拟时取n=10,20,50,100。3.温度t,也是影响模拟效果的一个重要变量,模拟时温度范围在0到6.0之间。4.玻耳兹曼常量,即boltzmann constant，取为1.0。5.自旋状态取值为+1或-1。6.模拟时，假设系统经过50000次metropolis蒙特卡罗循环达到平衡状态（n=100时，是经过500000次，因为它达到平衡时需要的时间更长），此时求解各物理量的初始值，然后开始求单个自旋物理量绝对值的平均值，在此，进行200000次蒙特卡罗循环，即求解次循环过程的单个自旋物理量绝对值的平均值，其他的物理量的模拟可以参见参考文献7,8。以上即为

22、模拟过程中的主要参量。 下面将结合数据作出的图表就行分析，主要从俩个方面，一个方面是不同的线度随温度t的变化趋势，另一个方面是在一定温度下，不同的线度n随时间t的变化趋势这两方面又分别从哈密顿量和磁化强度来进行分析。 首先分析n随t的变化趋势，如下图3，图中横坐标是t,纵坐标是系统中单个自旋在平衡态（假设系统经过50000次（n=100时是500000）metropolis蒙特卡罗循环达到平衡状态）下的磁化强度的绝对值的平均值|m|（由于单个自旋的状态值是正负对称的，所以磁化强度直接求和基本是趋于零的，从未无法分析数据，所以取绝对值）,由图可看出，对于不同的n来说|m|在t=1.0到2.0之间

23、都发生很大变化，由1.0很快趋于0，实际上这就是二维伊辛模型的相变，由此可知相变就发生在这个温度范围内，仔细观察后，发现，随着n的增大,|m|的突变愈加明显和尖锐，从而由图中对于临界温度的确定就愈加精确，更接近理论值，这种情况和理论上相符，说明模拟思想和过程是对的。因此，在模拟时，n，尽量使n足够大，模拟时间尽可能长，所得结果就会更接近理论值，在j=0.5，其他参量给定时，由图中|m|开始变化的温度大致看出，此次模拟的临界温度大概在t=1.2左右的范围内。 图 3 图 4 分析图4，横坐标是温度，纵坐标是系统平衡时单个自旋的能量e的平均值，由图可知，随着温度的变化，不同的n，其曲线基本重合，临

24、界温度的值基本和图3情况相同。图 5分析图5，横坐标是时间，纵坐标是系统平衡后系统磁化强度，选择t=2.5，j=0.5进行模拟，可知，系统达到平衡后，系统磁化强度随时间以m=0时为中心轴，在它上下进行波动，可以看出，在一段时间内系统磁化强度的平均值趋于零，系统大致处于平衡态，说明蒙特卡罗方法模拟的准确性，正确的反映了二维伊辛模型平衡时的磁化强度。分析图6，横坐标是时间，纵坐标是系统平衡时的系统哈密顿量，选择t=2.5，j=0.5进行模拟，与图5相吻合，随着时间的推移，不同n值的系统哈密顿量以各自平衡时的平均哈密顿量为中心轴，进行上下波动，也可看出，经过一段时间，系统哈密顿量的平均值基等于平均哈

25、密顿量，与理论相符合，正确反映了二维伊辛模型平衡时的系统哈密顿量。 图6最后将由图所得的结果与二维伊辛模型在理论上的精确解进行比较，由图3与图4可以得到，无论线度n的值为多少，磁化强度或者哈密顿量发生变化最为尖锐时的温度大概在t=1.2左右，其中磁化强度的变化更为显著，而理论上所得到的精确解，即有温度t=1.1364，二者相比较可以得到模拟的误差大约在5.6%的左右范围内，误差最大为29.6%，再次说明编写的程序的正确性，及该方法处理二维伊辛模型相变问题的可行性。 由以上分析可知，蒙特卡罗方法对二维伊辛模型的数值模拟所得结果与理论结果相符合，说明了该方法的可行性，蒙特卡罗方法可以较好的模拟二维

26、伊辛模型。同时从数据所得到的图表可以看出，在进行数值模拟时，要想较好的模拟出相变过程和得到临界温度，就需要使系统线度n在计算能力允许的范围内尽可能的大，而且模拟时间尽可能的长，这样得到的结果会更加接近理论值，使模拟更加可靠。 参考文献：1. ising e. eeitrag zur theorie des ferromagnetismus j. z phys,1925, 31: 253-258；2. 张志东，伊辛模型的研究进展简介，chinese journal of nature，vol.30 no.2：98-101；3nikolay prokofev， classical monte ca

27、rlo & metropolis algorithm，lecture notes：theory of monte carlo methods by nikolay prokofev；4刘策军，郑有因，二维伊辛模型自旋状态图样的蒙特卡罗模拟，华南农业大学学报，j6 (4)，1995：101-105；5. 李晓寒，王宗笠，宁旭，二维ising模型临界相变的monte-carlo数值模拟，安徽大学学报(自然科学版)，10002162(2008)03005604：56-59； 6. w. selke, saarbrficken, m. e. fisher, two-dimensional ising

28、models with competing interactions-a monte carlo study，z. physik b-condensed matter 40, 71-77 , may 16, 1980；7. robert knegjens, simulation of the 2d ising mode, may 13,2008；8. indrek mandre, the ising model course project in simulation of physical processes, tallinn 2008；9. 北京大学物理系，量子统计物理学m，北京大学出版社

29、，1987：380-460.附 录以下是本文进行模拟所编写的fortran语言程序。real, parameter : n=50 !定义系统线度real, parameter : j=0.5 !定义耦合常数和温度 integer, parameter : l=10000 !数据显示频率，即每l步显示一次数据 integer : s(n,n),z,c,d,e,a,b !二维伊辛模型自旋参数的定义real : m,x1,x2,x3,h,de,t real : r !r为接受率 do t=0,6.0,0.05 m=n*n call random_seed () !定义随机数 s=1 !自旋初始值的给定 z=0 !俩物理量在假设达到平衡时初始值的求解do while(z0)then de=de+s(i,k2) else de=de+s(n,k2) endif i=k1+1; if(i0)then de=de+s(k1,i) elsede=de+s(k1,n) endif i=k2+1; if(i1.or.x30)then e=e+s(c,d)*s(i,d) else e=e+s(c,d)*s(n,d) endif i=c+1; if(i0)then e=e+s(c