《《计算物理(研究生用)》[第3篇].ppt》由会员分享，可在线阅读，更多相关《《计算物理(研究生用)》[第3篇].ppt（89页珍藏版）》请在三一办公上搜索。
1、合肥工业大学理学院,第三篇 蒙特卡罗方法,3.1 蒙特卡罗方法概述 蒙特卡罗(Monte Carlo)方法：利用随机数统计地去计算和模拟给定的问题。说明：1、也称统计模拟法、随机抽样法或统计试验法。2、Monte Carlo方法的命名：世界上著名的赌城摩洛哥的Monte Carlo。3、方法：先构造一个与物理问题等价的随机过程，当完成大量的随机试验后，结果由多次事件的平均值给出。,合肥工业大学理学院,5、求解确定性物理问题：如微分方程、定积分、线性方程等。同其它数值计算方法相比是速度慢。6、求解复杂物理问题：如果物理学问题的严格算法不知道，或非常复杂，则蒙特卡罗方法有意想不到的成功。特别在分子
2、运动学、输运现象、布朗运动、放射性衰变等问题中，由于本身有一定的统计规律性，这种方法很奏效。7、计算机模拟：蒙特卡罗方法广泛应用于“计算机模拟”，在计算机上模拟真实过程。,4、随机数的抽样：现在都用计算机程序产生，一般不用物理方法抽样。计算机产生的是伪随机数。,合肥工业大学理学院,3.2 蒙特卡罗方法的原理,用蒙特卡罗方法寻找某个未知量x时，利用计算机产生的随机变量，得到的期望值E()x=E()3.1为了估算的平均值，构造随机变量的若干独立的实验数序列i|i=1，2，3，n，得,由统计理论随机变量是有限的均方差D()3.3,合肥工业大学理学院,实际上我们不知道变量的均方差，但在计算过程中可以根
3、据（3.3）式进行估计 3.5,式(3.4)可清楚地看出，蒙特卡罗方法以 的速率收敛，这是它的主要缺点。但与其它数值方法比，它能解决其它数值方法不能解决的问题，这正是它有强大生命力的原因所在。,然后引入概率误差 3.4,合肥工业大学理学院,3.3 伪（赝）随机数的产生,随机数产生：有两种，一是物理方法，二是数学方法。物理方法：产生的随机数仅遵循统计规律，无任何其它规律。如骰子、抽奖机等。理论上Monte Carlo方法可以利用物理方法产生随机数来研究。但实际应用不可行，运算量太大，速度上不去。数学方法：产生的随机数是利用计算机程序运算得到，称它为伪（赝）随机数。它虽然不是真正的物理随机数，但它
4、能通过一系列的统计检验，可以放心使用。,合肥工业大学理学院,（1）伪随机数的“伪”性：伪随机数有一定周期性，只不过周期很长。若在一次数值模拟中所用的随机数都在一个周期内，模拟效果与真正的随机数是一致的。（2）统计检验：在使用前，应对伪随机数进行统计检验，以确定其统计性能的优劣。主要是均匀性和独立性。（3）程序函数：多数程序设计语言都有此函数。如C、Matlab等用rand()函数来实现，此函数产生的是均匀分布的随机数。,说明：,合肥工业大学理学院,其中a，c，N为给定的三个整数，(c*Xn-1)(mod N)是指c*Xn-1被N整除后取余数。上述关系可推出一组X1、X2、X3 伪随机数。a，c
5、，N选择适当，这一序列伪随机数具有足够长的周期性。,(4)伪随机数的产生方法：主要是乘同余法，即下面的伪随机数递推关系：Xn=(c*Xn-1)(mod N)，X0=a,合肥工业大学理学院,例如：产生032767之间均匀分布的随机数递推公式 Xn=(889*Xn-1)(mod32768)，X0=13,若产生01之间均匀分布的随机数，上式改为 Xn=(889*Xn-1)(mod32768)/32767，X0=13,这个递推公式产生的伪随机数周期很长，是比较实用的一个产生伪随机数的公式。,合肥工业大学理学院,3.4 伪（赝）随机变量的抽样,实际上，大多数伪随机变量不满足0，1均匀分布，而是具有符合分
6、布密度函数为f(x)的分布。抽取符合f(x)随机数的步骤如下：,（2）从上面随机数总体中抽取一个简单子样，使它满足分布密度函数f(x)：1，2，3,（1）在0，1之间抽取均匀分布的伪随机数序列：1，2，3,合肥工业大学理学院,一、离散型分布随机变量的直接抽样法,随机变量的概率表其中xi为离散型随机变量的跳跃点，Pi为概率。,变换方法：第一步 构造一个矢量,合肥工业大学理学院,第二步 产生一个0，1间随机数，找到区间使恰好落在这个区间内。,第三步 P取对应的值xj，其表达式为,合肥工业大学理学院,例1.光子与物质相互作用的抽样问题。,物理过程：光子与物质作用有康普顿效应、光电效应和电子对效应三种
7、类型。其中光电效应和电子对效应为光子吸收过程。设三种过程的碰撞截面分别s、e和p，则总截面T=s+e+p。根据给定的0，1之间均匀分布的随机数，求应产生那种效应？,合肥工业大学理学院,（3）产生0,1间均匀随机数：若 s/T，则发生康普顿散射；若s/T(s+e)/T，则发生光电过程；若(s+e)/T，则发生电子对产生过程。,（2）构造一个矢量,解：（1）随机变量的概率表,合肥工业大学理学院,概率表,例2 Possion分布：Possion分布是离散型分布，其概率函数为求符合此分布的随机变量,合肥工业大学理学院,二、连续型分布随机变量的直接抽样法,一个对应一个。如果积分式能求出解析形式的反函数，
8、则得到变换公式 F=F-1(),合肥工业大学理学院,连续型分布的困难性：在实际问题中，连续型分布是很复杂的，有的只能给出分布函数F解析表达式，但给不出其反函数F-1()，如分布。有的则连分布函数都给不出解析式，如正态分布。由于这一情况，对于连续型分布采用求反函数的方法就及其困难，这必须寻找其它方法。,合肥工业大学理学院,例1.均匀分布：在区间（a，b）内找出符合均匀分布的随机变量。,解法2：线性变换，令,解法1：,合肥工业大学理学院,例2指数分布：指数分布的一般形式为 f(x)=e-x x 0其中0，求符合此分布的随机变量,求反函数=-1/ln(1-F()其中0F 1，由上式得=-1/ln(1
9、-),因为1-0,1 和0,1是等价，所以其式等价于=-1/ln,合肥工业大学理学院,三、连续型分布随机变量的变换抽样法,变换抽样法：变换积分形式得到反函数的抽样法。变换抽样的基本思想：（1）寻找到一个直接抽样的(y)，记(y)。（2）再寻找一个适当的变换关系x=g(y)。且g(y)反函数存在，一阶导数连续，记g-1(x)=h(x)。（3）根据概率论，x满足密度函数h(x)*|h(x)|。,（4）如果选择合适的g(y)，使之 f(x)=h(x)*|h(x)|则通过对(y)的抽样，再变换=g()得到满足f(x)的抽样值。,合肥工业大学理学院,例1正态分布(高斯分布)：正态分布的形式为 求符合此分
10、布的随机变量。,求积分,合肥工业大学理学院,正态分布抽样的方法有许多种。正态分布在统计物理学中是最重要的分布之一，也是可用最多的方法来产生随机数的分布之一。,取其中一个即可,另外是在0，2之间的数，必须=22,可用替代。,合肥工业大学理学院,乘分布舍选抽样法的基本思想：若随机变量满足密度函数f(x)，直接抽样不方便或不可求出时，可设法把f(x)写成 f(x)=H(x)f1(x),四、连续型分布的乘分布舍选抽样法(第二类),其中H(x),f1(x)是归一化密度函数。A）H(x)有上界M，即H(x)MB）f1(x)可直接抽样h=F 1()C）对f(x)的乘分布舍选抽样过程如左图。,合肥工业大学理学
11、院,例1Maxwell分布：Maxwell分布的概率密度为 当 时，即为统计物理学的速率分布。求符合此分布的随机变量,解：用乘分布舍选抽样法。把f(x)写成其中,合肥工业大学理学院,B）f1(x)可直接抽样h=F1(),C）对f(x)的抽样过程是挑选h，使h符合f(x)的随机分布。,A）H(x)是有上界M，即H(x)M。由H(x)取极值,合肥工业大学理学院,3.5 Monte Carlo方法求圆周率,求的物理模型如图所示，边长为1的正方形内切一个圆。于是正方形和圆的面积分别为 S正=1，S圆=/4,T正：计算机掷出的随机数总数；T圆：随机点(x,y)落在圆内的总数.即，当 时的随机点总数。,若
12、正方形内有T正个随机点(x,y)x=1-0.5；y=2-0.5其中为0，1之间的随机数。,圆内的随机点数为T圆，那么有 S正/S圆=T正/T圆 得：=4T正/T圆,合肥工业大学理学院,合肥工业大学理学院,%用Monte Carlo方法对圆周率进行计算的Matlab程序，文件名PI_val.m%K=100;M=5000;p=0;for i=1:K for j=1:M X=rand(1,2);R2=(X(1,1)-0.5)2+(X(1,2)-0.5)2;if R20.25|R2=0.25 end end i PI01=4*p/(i*M)end,合肥工业大学理学院,3.6蒙特卡罗方法求解确定性问题(
13、定积分),物理问题是以数学物理的经典方程表示出来的，如玻尔兹曼方程、拉普拉斯方程、麦克斯韦方程等。首先从解析解分析，若不行就用数值分析法求数值解。若还解决不了，就设法用蒙特卡罗方法求解了。蒙特卡罗方法求解确定性问题的基本思想是：对于给定的确定性问题，设计一个概率模型，然后采用一定的抽样方法按照所设计的概率模型抽样，最后把由这个模型产生的一个数学特征作为该确定问题的近似解。,合肥工业大学理学院,一一维定积分计算,1平均值法设定积分为,思想：随机分割n个条，在什么位置分割由随机数xi确定取0,1间均匀分布的随机数序列i(i=1,2,3，n)，并令,只要n足够大，则有,合肥工业大学理学院,2 随机投
14、点法,设积分,若，且g(x)不满足0g(x)1。可以设法将积分转换成其中，M和L为g(x)的最大值和最小值，如图所示。,合肥工业大学理学院,（1）在单位正方形产生两个随机1,20,1,（2）检验随机点(1,2)是否落入G区，如果条件2g(1)满足，则记录点(1,2)落入G内一次(记录器m=m+1)。（3）设在N次试验中，有m个随机点满足上式，那么 Im/N,合肥工业大学理学院,二多重定积分计算,设多重定积分为,取0，1间均匀分布随机点列为,并令,只要m足够大，则有,合肥工业大学理学院,3.7 蒙特卡罗方法求解确定性问题(非线性方程和非线性方程组求根),给定非线性方程和非线性方程组 f(x)=0
15、 fi(x1,x2,xn)=0 i=1,2,n其中f 和 fi 分别是 x 和 x1,x2,xn 的非线性函数。在指定的求根范围内求解x*和x1*,x2*,xn*。解决上述问题可采用数值解法。最常用的有牛顿下山法、埃特金迭代法等；还有线性化方法和求函数极小值方法。,合肥工业大学理学院,在使用上述数值算法时，有时会遇到一些困难，甚至方法失效，近似解都不能获得。在这种情况下，可以用蒙特卡罗方法求解。虽然该方法耗时较多，但方法简单，比较实用。因此对要求精度较高的求解问题，可先用蒙特卡罗方法求得一个初解，再进行其它方法的运算求精度更高的解。,合肥工业大学理学院,一蒙特卡罗方法求非线性方程的一个实根,设
16、函数方程 f(x)=0（1）给定常数b0、产生随机数的个数m 和精度。（2）选取初值x*，计算出F0=f(x*)。（3）在-b,b上反复产生均匀分布的随机数，对于每个计算出F1=f(x*+)，直到发现|F1|F0|为止。然后令 x*+x*;F1F0（4）如果连续产生m个随机数还不满足|F1|F0|，则将b减半再进行第三步运算。（5）重复上述过程，直到|F0|为止。若遇到不收敛时，则可以适当调整b和m的值。,合肥工业大学理学院,二蒙特卡罗方法求非线性方程组的一组实根,设非线性方程组 fi(x1,x2,xn)=0 i=1,2,n 定义目标函数,目的是以下降法求出函数(x)的极小值点，此极小值点可以
17、认为是非线性方程的一个近似解。具体 对预先给定精度0，在求根区域内，模拟给出方程近似解（x1*,x2*,xn*）使之满足（x1*,x2*,xn*）求解过程采用随机搜索算法：,合肥工业大学理学院,(2)给定常数b0，在-b,b之间，反复产生均匀分布的随机数组(1,2,n)，计算 F1=（x1*+1,x2*+2,xn*+n）(3)直到发现一组(1,2,n)使得F1F0为止，然后令 xi*+ixi*i=1,2,n;F1F0(4)若连续产生m组随机数(1,2,n)，还不满足F1F0，则将b减半再重复2、3过程，直到F0为止。求出非线性方程组的近似根(x1*,x2*,xn*)。若遇到不收敛时，则适当调整
18、b与m，再进行运算。,(1)选取初值(x1*,x2*,xn*)，计算 F0=(x1*,x2*,xn*),合肥工业大学理学院,3.8 蒙特卡罗方法求解确定性问题(椭圆型差分方程),实际应用中求解偏微分方程可用差分法。例如：Laplace方程第一边值问题。在迭代过程中要占用大量的存贮单元。但是在许多情况下，只需计算某些特殊位置的函数值就行了。,例如：燃气轮机上叶片的截面如图形状，主要研究叶片尖端P附近的物理现象，而其它位置并不重要。因为当叶片表面温度很高时，可能使叶片尖端P点附近的内部物理结构发生较大变化。,合肥工业大学理学院,一、蒙特卡罗方法解 Poisson方程第一边值问题,问题的提出：二维P
19、oisson方程第一边值问题是,D：二维域；：D的边界；P：D内点集；Q：线上点集。,差分方程：在D内，以步长h为方形网边长。网格有两类，第一类是四个相邻网格位于D内或边界上，称内部网格点D，图中用“”表示。第二类是相邻网格点小于四个，这类点称为边界网格点，图中用“”表示。,合肥工业大学理学院,二维Poisson方程相应的差分方程为,其中P是D 内的节点，P1j(j=1,2,3,4)是与P相邻的四个节点，如图所示。Q是边界上的节点。,构造Monte Carlo方法：为了求函数u(p)在域D内节点P处的值u(p)，构造随机游动模型：,u(Q)=f(Q),Q,合肥工业大学理学院,设想有一个质点自P
20、=P0点出发，以等概率1/4向与P相邻的四个点P11,P12,P13,P14某一点随机游动一步；然后再按同样的方式自新的位置向与之邻接的四点中某一点随机游动一步。如此下去，直到质点第一次到达边界上的一个子节点处时，游动即告终止。,(2)四个方向是等概率的原因是差方程中存在，即Laplace算子 的X，Y偏微分前的系数相同。若不相同，就不是等概率。,说明：(1)实际上是构造马尔可夫链；,合肥工业大学理学院,设一条游动路线p是p:P=P0P1P2P3Pk-1Q,计算随机变量：游动刚开始=0，质点每游一步，应把游动之前的值增加一项，这项值为该点Pi的函数值，即h2q(Pi)/4；质点到达边界时的节点
21、Q处时，游动停止，此时值再最后增加一项f(Q)。因此与路线p相应的随机变量为,合肥工业大学理学院,当质点是从上的Q出发时，质点等于不游动，那么,说明：（1）可以证明的数学期望E()=u(p)(证明略)（2）由P点开始而在边界上终止的游动称为一次随机游动，并有一个随机变量对应。,=(p)=f(Q),经过N次随机游动，即获得N个随机变量1,2N。其平均值,合肥工业大学理学院,二、蒙特卡罗方法解一般椭圆型方程第一边值问题,对于一般椭圆型方程的第一边值问题，除能构造出等概率模型外，还可以构造变概率模型。,其中,问题提出：一般三维椭圆型方程第一边值问题的可写成,合肥工业大学理学院,差分形式：上式可化为差
22、分方程形式：,D,网格点域,其中,合肥工业大学理学院,h的条件 网格步长h应满足条件,游动模型 下面构造随机游动概率模型：设质点自p=p0D处出发，按照如下概率表,其中min,max取最小值和最大值。,合肥工业大学理学院,分别向与P相邻的四个节点P11,P12,P13,P14处随机游动一步。若质点到达P1j(j=1,2,3,4)，则再按概率,其中X1j=X(P1j),X:a,b,c,d 分别向P1j(j=1,2,3,4)相邻的四个节点随机游动一步。如此重复下去，直到第一次到达边界点Q处时为止。,游动路线：p:P=P0P1P2P3Pk-1Q,合肥工业大学理学院,质点出发前=0，每游动一次，在中需
23、增加一项-(h2/2)(j)Zj。最后离开Pk-1D到达边界节点Q时，在中再增加一项(k-1)f(Q)，那么对于p的路线定义随机变量值为：,其中Wi=W(Pi),Zi=Z(Pi),i=0,1,2,k-1.。,如果质点游动的出发点是边界节点Q，则它就不动，此时为=(Q)=f(Q),合肥工业大学理学院,这样把P点开始游动到上Q点终止称为一次随机游动，经过N次随机游动，即得1,2N,合肥工业大学理学院,实际链式反应过程是一个随机过程，可以利用蒙特卡罗方法来研究链式反应。,3.9 蒙特卡罗方法对链式反应的模拟（核临界安全模拟计算）,链式反应模型：某些元素(如235U、239U)是不稳定的，它能自发地裂
24、变，变成碎片，同时放出中子。中子被其他核材料吸收，使其更不稳，导致迅速裂变，放出更多的中子。,合肥工业大学理学院,核爆炸：一般235U核的半衰期很长，约7.07*109 年。因此在某一时刻内，只有非常少量的235U发生裂变，释放的能量很小。但是当裂变一旦引起链式反应后，将释放出大量能量，导致核爆炸。,链式反应要求：裂变产生的二个中子中，平均来说被其他235U核吸收的中子要大于1，以保证增殖过程。若小于1，则是熄灭的过程。能否链式反应，与235U核团的体积(质量)有关，裂变中子在小质量的235U核团中，吸收几率很小，大质量235U核团中才能持续裂变。临界质量：能维持链式反应的核材料最小质量。,合
25、肥工业大学理学院,临界条件定量分析,残存因子：N次核裂变，Nin个中子被铀核吸收,核爆炸：当铀块的质量大于MC时，发生链式反应。原子弹原理：取两块小于临界质量的235U，将两者合在一起，质量超过临界质量。安全运输：理论计算临界质量十分重要。,f 1表示裂变数增加，发生链式反应；f 1表示反应逐渐停止；f=1表示处于临界状态，用临界质量Mc表示。,合肥工业大学理学院,用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然后计算放出中子被吸收引起的裂变数Nin，求出相应的f 值。,为简单起见，考虑铀的几何形状为长方形a*a*b，如图所示。蒙特卡罗法研究这种形状的核材料临界问题：,合肥工业大学
26、理学院,（1）核裂变的随机位置（X0，Y0，Z0）,随机点的坐标范围为-0.5a X0 0.5a；-0.5a Y0 0.5a；-0.5b Z0 0.5b,（2）裂变的两个中子出射方向(,)裂变产生的中子出射方向用和来描述。从（X0，Y0，Z0）点放出的中子与以此为圆心的单位球上某一面积发生碰撞，它的碰撞几率只与球面上被碰撞的面积大小成正比。或者说碰撞按立体角均匀分布。,一次裂变需要3个随机数来确定裂变点的位置。,合肥工业大学理学院,立体角：d=sindd=-ddcos,按立体角内均匀分布是指:角在02之间均匀分布；在之间，以cos的值在-11之间均匀分布。,注意：cos均匀分布，不是角均匀分布
27、。如果按照为均匀分布抽样，则在=/2附近有较大的权重。,一次裂变需要4个随机数来确定两个中子的方向。,合肥工业大学理学院,(3)中子平均自由程,平均自由程：中子与铀核碰撞所走的平均路程。在平均自由程内与核碰撞几率相同：例如平均自由程为1cm，中子在核块内飞越的距离为0.3cm，则中子有30%几率的与铀核发生碰撞被吸收。中子飞越距离d可以用0,1之间的随机数表示。,一次裂变需要2个随机数来确定两个中子的平均自由程。,中子的位置为,合肥工业大学理学院,判断中子与铀核发生碰撞：若（X1，Y1，Z1）点在体积内就认为可以发生碰撞被吸收，反之中子就不会引起下一次裂变。,若有N个随机点裂变，Nin次总碰撞
28、数（有效中子数），则残存因子,合肥工业大学理学院,计算残存因子的流程图,合肥工业大学理学院,以下给出计算残存因子的流程图的说明：,M=V：铀核块的质量（假设密度均匀且数值为1）s=a/b：核块边长比N：总随机裂变次数Ep=10-4：计算f=1时要求的误差（临界状态）,（2）总共有N个随机裂变点，每个点要求九个0，1之间的随机数，i(i=1,2,，9）,（1）由M和s求长方体的边长,合肥工业大学理学院,（3）随机裂变点的坐标三个,（4）裂变后两个中子的出射方向四个 第一个中子两个方向 第二个中子两个方向,（5）裂变后两个中子的飞越距离两个 d1=8 d2=9,合肥工业大学理学院,取九个随机数，按
29、式计算第一个中子的（X1，Y1，Z1），判断此坐标点是否在铀块内，若是在内部则Nin增加1，接着同样方法计算第二个中子；重复上步N次；计算 f=Nin/N；调节常数M和s使 f=1，得到临界时的Mc和s的值。,（6）模拟计算：,合肥工业大学理学院,3.10 蒙特卡罗方法对趋向平衡态的模拟,势力学第二定律：孤立体系，在平衡态时熵不变。在非平衡态时熵趋向增加，最后熵达到最大值而趋向平衡态。典型实例：正中间用隔板隔开的一个密闭盒子，开始时左边充满某种气体，当隔板中间开一小孔后，由于分子无规则的热运动，一些分子通过小孔进入隔板的右边，以后隔板两边分子来回运动，最后达到两边的压强相等。从统计观点来看，时
30、间很长后分子全部位于左边的情况几率很小。,合肥工业大学理学院,一、趋向平衡态数学模型,利用N个硬币可以描述上述问题的模型（Ehrenfest模型）。,国徽面分子在左边；分值面分子在右边。初始时N个硬币都是国徽面。随机取硬币改变状态，国徽面分值面或分值面国徽面，研究硬币状态。经多次取硬币改变状态，结果国徽面和分值面的硬币数基本相等（相当于平衡态）。,合肥工业大学理学院,二、理论分析,经无规地改变X次后，有n个分值面向上的硬币。此时再取下一个硬币，随机取到分值面向上的几率为n/N，相当于有n/N的几率分值面国徽面，设PB=n/N。对国徽面向上的硬币而言，随机取到的几率为1-PB，表示有1-PB的几
31、率使国徽面分值面，设PA=1-PB。平均说，下一次试验分值面向上的几率净增加为 n=PA-PB=1-2n/N,合肥工业大学理学院,对于每一次状态变化用X表示，即 X=X(ti+1)-X(ti)=1 于是前式变为 n=(1-2n/N)X,把n和X当成连续处理，上式变为,合肥工业大学理学院,（1）经过X次变化后，平均来说硬币分值面向上的几率是PB；（2）当X时，PB0.5，这就是平衡态。（3）按具体过程画出PB=n/N随X变化的曲线，可以看出PB是平均地以指数形式趋向平衡值0.5，PB曲线在指数形式的曲线附近进行涨落。,说明：,合肥工业大学理学院,三、蒙特卡罗方法模拟,模拟方法：初始时N个硬币的都
32、是国徽面，经过多次随机取硬币，每取一次硬币改变一次状态，即国徽面分值面或分值面国徽面。最后统计一下硬币国徽面或分值面所占比例。结果应该是各占1/2，此时也就是平衡态。模拟曲线：如果按模拟过程画出PB=n/N随X变化的曲线，可以看出PB是平均地以指数形式趋向平衡值0.5，并在曲线附近涨落变化。,合肥工业大学理学院,N：硬币的总数 MX：无规地改变硬币状态的总次数NU(I)：一维数组，用来描述每一个硬币的状态。定义,n：分值面的硬币数K：1,N间的随机数；利用N*rand()+1取整（MATLAB语言：fix()；C语言：(int)()。,合肥工业大学理学院,3.11 麦曲罗保利斯（Metropo
33、lis）方法,一、统计物理中的理论问题,积分是对整个体系的相空间进行，若是N个单原子分子构成的体系则有6N重积分(双原子，有12N重积分)。一般上式不能进行解析积分计算。,合肥工业大学理学院,如果能对体系的微观量A进行随机抽样，则上式变为,其中Xi是第i个微观状态(N个坐标和动量组成)，A(Xi)为微观状态Xi所对应的微观量，M为抽样的次数。,说明：（1）原来有权重的统计平均值，现变为算术平均。（2）计算时，由于因子e-H存在，若均匀抽样，则对平均值贡献大的微观状态被抽到次数少，得到的平均值偏低。,合肥工业大学理学院,(X)是高维数的分布，又有一个未知的归一化常数(即分母)，按(X)抽样实际上
34、仍然不可能。,解决的办法：选择按分布来进行抽样,合肥工业大学理学院,二、Metropolis方法,Metropolis思想：设想微观量的抽样通过Markov过程来进行。具体就是不要彼此独立无关地抽取微观状态Xj，而由前一个状态Xi通过适当的跃迁概率W(XiXj)得到，当跃迁次达到足够大后再抽取。开始以(Xi)分布的Xi态是随机取的，当步数M时，随机分布(Xi)趋向平衡态分布。对Xi到Xj的Markov过程的跃迁几率用W(XiXj)表示，在步数M足够大后，Markov过程上的两个微观状态Xi和Xj 之间满足细致平衡条件：(Xi)W(XiXj)=(Xj)W(Xj Xi),合肥工业大学理学院,于是有
35、,上述Markov过程常称为Markov链。,合肥工业大学理学院,确定Markov链进行的方向：当H0时，则进一步用下式判断,能量增加较少允许跃迁，增加较多不允许跃迁，判断依据是。物理原因是体系存在涨落。经N步(N很大)后，体系从随机的初态达到了平衡态附近，抽样计算平均值,其中是0，1的随机数。,合肥工业大学理学院,相应的偏差计算公式为,上面的Markov过程，画出图解如图：,合肥工业大学理学院,三、实例：XY模型的二维自旋系统,1.二维自旋系统的XY模型该体系的哈密顿量为,其中J为耦合常数；表示只对近邻的元素自旋求和；每对只作用一次；i是第i个自旋方向与XY平面法向的夹角，范围在0,2之间。
36、,合肥工业大学理学院,2.周期边界条件,体系有N*N个自旋点，自旋点处在正方形的格点上。边界取周期性边界条件，相当于多个N*N的正方形连接。,3.模拟计算能量平均值 体系的初始状态为随机态,i取0,2间随机值。用Metropolis方法计算某个温度（如J=1）时二维自旋体系的能量平均值和它的偏差。,程序流程图：,合肥工业大学理学院,Nx：数组能接受的Nx*NxN：由人输入的 N*NN*N自旋粒子初始随机角度，值的范围在0,2周期边界条件：N+i=i,计算时可能用到N+2=2M：跃迁步数，每调用一次函数change()就是Markov过程跃迁一次。1000：非平衡态向平衡态预跃迁次数。,主函数（
37、一）,合肥工业大学理学院,求能量的平均值和它的偏差K：计算能量时抽样次数Ex：计算最邻近的哈密顿量Hi，主要有四个最邻近的,主函数,每对只作用一次，即仅求I+1，J+1的和,不用求I-1，J-1的和,合肥工业大学理学院,跃迁函数,如果跃迁，则重整周期性，否则返回,（1）I，J从2开始到N+1结束，仍为N*N个点计算，主要避免下标J-1或I-1有0的情况（2）对N*N个格点的每个(I,J)粒子都要计算，判断是否跃迁到新的状态A，即(ij),Y,N,合肥工业大学理学院,3.12 经典体系自由能的蒙特卡罗计算,一、自由能F的引入,体系的哈密顿量H：具有相互作用的多体系统，其哈密顿量为,其中i、Pi：
38、第i个粒子的广义坐标、广义动量；N：体系的总粒子数；V(i-j)：粒子间相互作用势能；ij下标：所有粒子间求和。,合肥工业大学理学院,自由能：由统计物理知，体系自由能：,其中=1/kT、Z：体系的配分函数 式中积分是对广义坐标、广义动量进行。,热力学量：由热力学第一定律：dH=TdS-PdV-Mdh和自由能的定义：F=H-TS 得 dF=-SdT-PdV-Mdh由此可知,其中P：压强;M：磁距；h：磁场；S：熵；CV：比热注意符号表示不要混淆：内能是H(非焓)，不是U，磁场是h，不是H。,合肥工业大学理学院,由上式可知，只要计算出自由能F，就可以得到体系的全部热力学平衡性质了；直接计算自由能F
39、非常困难；计算自由能F主要是计算配分函数Z，由于体系内存在相互作用，实际上配分函数的解析准确计算也非常困难。,分析：,合肥工业大学理学院,二、计算方法,计算自由能差值：为了避免直接计算自由能，可以利用计算自由能的差值来计算热力学其它物理量。现在定义,合肥工业大学理学院,积分公式离散化：微观状态一般是分裂的情况，将上式离散化，得,合肥工业大学理学院,变化过大的解决：一般来说，变化范围都很大。通常用Markov链到平衡以后的能量HE作为参考点进行解决。,现在只要求出Q(T)，就可以求出体系中其它热力学量了。,合肥工业大学理学院,热力学物理量的计算：,(2)比热,合肥工业大学理学院,(3)某一温度下磁化强度M和磁化率的计算：,若体系存在外磁场h时，则同样计算自由能的差值。,得,令,合肥工业大学理学院,分析：(1)综上所述，可以利用固定某一温度和磁场，作为热力学状态参考点，采用蒙特卡罗法计算不同温度和磁场下的自由能差值，来得到宏观热力学物理量。,(2)考虑二自旋格子系统，存在外磁场。若取XY模型，则哈密顿量为,利用上节的程序，只需将哈密顿量改变为现在的形式表示，其中计算能量平均值和偏差改为不同温度和不同磁场时的式子计算即可。,
