《计算材料学》PPT课件.ppt

计算材料学ComputationalMaterialsScience第二讲 蒙特卡罗（MonteCarlo）方法主讲：张晖电话：Email:本堂课 主要内容Monte Carlo模拟发展简介Monte Carlo模拟基本原理Monte Carlo模拟典型算法Monte Carlo模拟典型应用蒙特卡洛法是什么？蒙特卡洛法是什么？蒙特卡洛蒙特卡洛(Monte Carlo)方法，方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。蒙特蒙特 卡罗方法卡罗方法（Monte Carlo method），也称统计模拟方法统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。MC的基本思想 MCMC基本思想很早以前就被人们所发现和利用。基本思想很早以前就被人们所发现和利用。1717世纪，人们就知道用事件发生的世纪，人们就知道用事件发生的“频率频率”来决定事来决定事件的件的“概率概率”。但要真正实现随机抽样是很困难的，。但要真正实现随机抽样是很困难的，甚至几乎是不可能的。甚至几乎是不可能的。高速计算机的出现，使得用数学方法在计算机高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。上大量、快速地模拟这样的试验成为可能。确定性系统确定性系统随机性系统随机性系统模拟模拟自然界自然界Monte-Carlo模拟，即随机模拟（重复模拟，即随机模拟（重复“试验试验”）重复试验重复试验计算机模拟计算机模拟Monte Carlo方法：方法：亦称统计模拟方法，亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte Carlo名字的由来：名字的由来：是由是由Metropolis在二次世界大战期间提出的：在二次世界大战期间提出的：Manhattan计计划，研究与原子弹有关的中子输运过程；划，研究与原子弹有关的中子输运过程；Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博闻名的首都，该城以赌博闻名Nicholas Metropolis(1915-1999)Monte-Carlo,MonacoMonte CarloMonte Carlo方法简史方法简史简单地介绍一下简单地介绍一下Monte Carlo方法的发展历史方法的发展历史1 1、BuffonBuffon投针实验：投针实验：18世纪，法国数学家世纪，法国数学家ComtedeBuffon利用投针实验估计利用投针实验估计 的值的值dL1777年法国科学家布丰提出的一种计算圆周率的方法年法国科学家布丰提出的一种计算圆周率的方法随机投针法，即著名的随机投针法，即著名的布丰投针布丰投针问题。这一方法的问题。这一方法的步骤是：步骤是：1)取一张白纸，在上面画上许多条间距为取一张白纸，在上面画上许多条间距为d的平行线。的平行线。2)取一根长度为取一根长度为l（ld）的针，随机地向画有平行的针，随机地向画有平行直线的纸上掷直线的纸上掷n次，观察针与直线相交的次数，记为次，观察针与直线相交的次数，记为m 3）计算针与直线相交的概率）计算针与直线相交的概率 布丰本人证明了，这个概率是布丰本人证明了，这个概率是:p=2l/(d)，为圆周率为圆周率：利用这个公式可以用概率的方法得到圆周率的近似利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料值。下面是一些资料实验者年代投掷次数相交次数圆周率估计值沃尔夫1850500025313.1596史密斯1855320412193.1554德摩根18806003833.137福克斯188410304893.1595拉泽里尼1901340818083.1415929赖纳192525208593.1795布丰投针实验是第一个用几何形式表达概率问题的例子，他首次使用随机实验处理确定性数学问题，为概率论和蒙特卡罗方法的发展起到一定的推动作用。Monte Carlo方法之方法之随机数的产生许多计算机系统都有随机数生成函数F90:callrandom_seedcallrandom_number(a)2、ISEED=RTC()X=RAN(ISEED)Y=RAN(ISEED）Matlab:x=rand(N)产生元素在(0,1)间随机分布的N*N矩阵s=rand(state,0)重设该生成函数到初始状态注意：上述随机数序列均具周期性，如上页random子程序的周期约230。实例一、计算值计算过程：计算过程：1 1、构造或描述问题的概率过程、构造或描述问题的概率过程2 2、从概率密度函数出发进行随机抽样，实、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的现从已知概率分布的抽样，得到特征量的一些模拟结果一些模拟结果计算均值计算均值Monte Carlo方法之方法之典型算法与应用考虑平面上的一个边长为考虑平面上的一个边长为1 1的正方形及其内部的一个形的正方形及其内部的一个形状不规则的状不规则的“图形图形”，如何求出这个，如何求出这个“图形图形”的面积的面积呢？呢？Monte CarloMonte Carlo方法是这样一种方法是这样一种“随机化随机化”的方法：的方法：向该正方形向该正方形“随机地随机地”投掷投掷N N个点，若有个点，若有M M个点落于个点落于“图形图形”内，则该内，则该“图形图形”的面积近似为的面积近似为M/NM/N。用该方法计算的基本思路是：1、根据圆面积的公式：s=R2,当R=1时，S=。2、由于圆的方程是：x2+y2=1(x2为x的平方的意思),因此1/4圆面积为x轴、y轴和上述方程所包围的部分。3、如果在1*1的正方形中均匀地落入随机点，则落入1/4圆中的点的概率就是1/4圆的面积。其4倍，就是圆面积。由于半径为1，该面积的值为的值。REALR,R1,R2,PIISEED=RTC()N0=0N=300000DOI=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R1.0)N0=N0+1ENDDOPI=4.0*N0/NWRITE(*,*)PIEND面积的计算面积的计算f(x)x辛普逊方法辛普逊方法I=Sn蒙特蒙特-卡洛方法卡洛方法f(x)x在长方形中均匀投在长方形中均匀投N0组组(x,y)如如 yf(x),则则 N=N+1I=(N/N0)S0SS011MC的优点MCMC与传统数学方法相比，具有与传统数学方法相比，具有直观性强，简便易行的优点，直观性强，简便易行的优点，该方法能处理一些其他方法无法解决的负责问题，并且容该方法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以易在计算机上实现，在很大程度上可以代替许多大型、难代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。实验误差。Monte Carlo方法利用方法利用随机抽样随机抽样的方法来求解物理问题的方法来求解物理问题;数值解法数值解法:从一个物理系统的数学模型出发从一个物理系统的数学模型出发,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态;Monte Carlo方法方法并非只能用来解决包含随机过程并非只能用来解决包含随机过程的问题的问题:例例如如:用用Monte Carlo方法计算定积分方法计算定积分.对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程,然后用然后用Monte Carlo方法进行求解方法进行求解注意以下两点：注意以下两点：MC的应用自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；数值分析：数值分析：数学问题，求积分，求逆矩阵，解线性代数等数学问题，求积分，求逆矩阵，解线性代数等经济学模拟：经济学模拟：库存问题，随机服务系统中排队问题库存问题，随机服务系统中排队问题 人口问题：人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争人口的出生，传染病的蔓延；乃至动物的生态竞争金属学：扩散、组织长大、相金属学：扩散、组织长大、相变过程变过程蒙特-卡洛模拟的意义能研究不同边界、不同材料的影响能研究不同边界、不同材料的影响 理论不可能、实验耗费太大理论不可能、实验耗费太大 用于实验设计用于实验设计无污染无污染 反应堆防护反应堆防护 核弹爆炸核弹爆炸能摆脱实验误差能摆脱实验误差 作理论和实验的桥梁作理论和实验的桥梁Monte Carlo模拟的步骤：模拟的步骤：1.根据欲研究的物理系统的性质，建立能够描述该根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量系统特性的理论模型，导出该模型的某些特征量的概率密度函数；的概率密度函数；2.（即构造或描述问题的概率过程）（即构造或描述问题的概率过程）2.从概率密度函数出发进行随机抽样，实现从已知从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；概率分布的抽样，得到特征量的一些模拟结果；3.有了明确的概率过程后，为了实现过程的有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的数值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。变量的大量试验值。产生已知概率分布的随机变产生已知概率分布的随机变量，是实现量，是实现MC方法的关键步骤方法的关键步骤，其中最基本的是，其中最基本的是（0，1）均匀分布。）均匀分布。3.对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。4.模拟结果的检验模拟结果的检验Monte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数(pdf)必须给出描述一个物理系统的一组概必须给出描述一个物理系统的一组概率密度函数率密度函数;随机数产生器随机数产生器能够产生在区间能够产生在区间0,1上（均匀）分布的随机数上（均匀）分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,1上均匀分布的随机数出发上均匀分布的随机数出发,随随机抽取服从给定的机抽取服从给定的pdf的随机变量的随机变量;模拟结果记录模拟结果记录记录一些感兴趣的量的模拟结果记录一些感兴趣的量的模拟结果误差估计误差估计必须确定统计误差（或方差）随模拟次数以及其必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；它一些量的变化；减少方差的技术减少方差的技术利用该技术可减少模拟过程中计算的次数；利用该技术可减少模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法实例二 定积分计算事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分!计算计算x*2在在(0,1)上积分上积分计算过程：计算过程：1、构造或描述问题的概率过程、构造或描述问题的概率过程:产生服从分产生服从分布布f（x）的随机变量）的随机变量Xi()（i=1,2,N）2、从概率密度函数出发进行随机抽样，实、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的现从已知概率分布的抽样，得到特征量的一些模拟结果一些模拟结果计算均值计算均值()REALYY=0N=300000ISEED=RTC()DOI=1,NX=RAN(ISEED)Y=Y+X*2/NENDDOWRITE(*,*)YENDlimx2dx (dx 0)Monte Carlo方法另一个重要问题：方法另一个重要问题：随机数随机数 随机数：由单位矩阵分布中所产生的简单子样称为随机数序列，其中的每一个个体称为随机数。但真正的随机数的不适合电子计算机上使用，因为它需要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，且其产生的序列无法重复实现，使程序无法复算，结果无法验证，同时需要增添随机数发生器和电路联系等附加设备。伪随机数：是有数学递推公式所产生的随机数。（近似的具备随机数的性质。）An+1=T(A);An+1=An+k+1伪随机的优点和缺点：判断伪随机数好坏的方法：1、它能够有较好的均匀性和独立性；2、它的费用大小，即指所消耗计算机的时间；3、容量要求尽可能大。随机数产生的办法产生均匀分布随机数的几种方法；（1）物理方法；（2）数学方法。伪随机数产生方法：加同余法乘同余法乘加同余法取中方法逆变换法合成法筛选法。关于随机数的几点注意关于随机数的几点注意注注1 1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注注2 2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。FORTRAN语言产生随机数的实例random_number(x)random_number(x)产生一个产生一个0 0到到1 1之间的随机数（之间的随机数（x x可以是向可以是向量），但是每次总是那几个数。量），但是每次总是那几个数。用了用了random_seed()random_seed()后，系统根据日期和时间随机地提供种后，系统根据日期和时间随机地提供种子，使得随机数更随机了。子，使得随机数更随机了。programrandomprogramrandomreal:xreal:xcallrandom_seed()!callrandom_seed()!系统根据日期和时间随机地提供种系统根据日期和时间随机地提供种子子 callrandom_number(x)!callrandom_number(x)!每次的随机数就都不一样了每次的随机数就都不一样了 write(*,*)xwrite(*,*)xstopstopendprogramrandom endprogramrandom Monte Carlo方法之方法之随机数的产生许多计算机系统都有随机数生成函数F90:callrandom_seedcallrandom_number(a)2、ISEED=RTC()X=RAN(ISEED)Y=RAN(ISEED）Matlab:x=rand(N)产生元素在(0,1)间随机分布的N*N矩阵s=rand(state,0)重设该生成函数到初始状态注意：上述随机数序列均具周期性，如上页random子程序的周期约230。FinitedifferenceapproximationofdifferentialequationsAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtaime-CentralspaceFICK第二定律Fick第二定律稳态扩散解REALC0(0:1000+1),C(0:1000+1)!C(DISTANCE)C0=0.1C0(0)=0.8!BOUNDARYCONDITIONC0(1001)=0.1!BOUNDARYCONDITIONC=C0OPEN(1,FILE=F:DIF.dat)DOJT=1,50000！TimeDOIX=1,1000!distanceC(IX)=C0(IX)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX)!C(0)=0.8!C(1001)=0.1IF(JT=50000)WRITE(1,*)IX,C(IX)ENDDOC0=CENDDOEND应用之二生日问题MC模拟模拟假设有假设有n个人在一起，各自的生日为个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直天之一，根据概率理论，与很多人的直觉相反，只需觉相反，只需23个人便有大于个人便有大于50的几的几率人群中至少有率人群中至少有2个人生日相同。个人生日相同。n理论几率模拟几率100.1170.110200.4110.412230.5270.520300.7060.692450.9410.93646500.9860.987INTEGERM(1:10000),NUMBER1(0:364),NUMBER2REALX,YISEED=RTC()DOJ=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1DOI=1,365Y=RAN(ISEED)NUMBER2=INT(365*Y+1)ETR=COUNT(NUMBER1.EQ.NUMBER2)IF(ETR=1)THENEXITELSEJJJ=JJJ+1M(J)=JJJNUMBER1(I)=NUMBER2ENDIFENDDOENDDODOI=1,10000IF(M(I).LE.23)SUM=SUM+1ENDDOPRINT*,SUM/10000ENDMC在材料学领域的应用 随机行走背景背景如，布朗运动最简单、无限制随机行如，布朗运动最简单、无限制随机行走走(Unrestricted randon walk,RW)Unrestricted randon walk,RW)startend平均平方端平均平方端-端位移端位移:,自然科学和社会生活中很多现象都与随机运动有关自然科学和社会生活中很多现象都与随机运动有关可以模拟的内容？扩散；分子运动；。如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了所有分子距原点的方均距离：！Monte Carlo Simulation of One Dimensional Diffusion INTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XX(J,I)：X*X,J:第几天实验,I:第几步跳跃!XXM(I):THEMEANOFXXWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验X=0!DOI=1,IMAX!第几步跳跃RN=RAN(ISEED)IF(RN0.5)THENX=X+1ELSEX=X-1ENDIFXX(J,I)=X*XENDDOENDDOOPEN(1,FILE=“f:DIF1.DAT)DOI=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX!WRITE(1,*)I,XXM(I)ENDDOCLOSE(1)END!Monte Carlo Simulation of Two Dimensional DiffusionINTEGERX,Y,XY(1:1000,1:1000)REALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验X=0!Y=0!DOI=1,IMAX!第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THENx=xy=y-1ENDIFIF(RN.LT.0.5.AND.RN.GE.0.25)THENx=xy=y+1ENDIFIF(RN.LT.0.75.AND.RN.GE.0.5)THENx=x-1y=yENDIFIF(RN.GE.0.75)THENx=x+1y=yENDIFXY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=“f:DIF2.DAT)DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END!MonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DOJ=1,JMAX!第几天实验X=0!Y=0!DOI=1,IMAX!第几步跳跃RN=4*RAN(ISEED)+1X=X+XN(RN)Y=Y+YN(RN)XY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=C:DIF2.DAT)DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)！做三维空间随机行走？！做三维空间随机行走？END实际环境是复杂的实际环境是复杂的:可变，缺陷，风，可变，缺陷，风，季节，电场等。季节，电场等。空间中有随机性缺陷空间中有随机性缺陷不退行走不退行走自回避行走自回避行走例如模拟高分子的位形用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形，行走过的位置代表的是构成分子的原子或官能团，因此，无限制随机行走忽略了体斥效应。不退行走就是禁止在每一步行走后立即倒退，可以解决刚走的一步与上一步重叠的问题。但不退行走没有完全解决高分子的体斥效应。自回避行走就是所有已走过的位置不能再走，这样就完全解决了体斥效应问题。气体分子的随机行走假设在一个空旷封闭的房间中心滴上一滴香水，挥发的香水分子随机的与空气中的粒子发生碰撞，最终会扩散到整个房间。这里我们将通过计算机模拟400个分子在二维平面内的随机运动讨论熵变现在来讨论一下熵在我们这个气体扩散模型中是什么意思。在刚开始的时候，我们的系统处在一个非常有序的状态：所有的分子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随着时间的推移，分子开始不断的向外扩散，这时系统出现了无序状态，熵开始逐渐增加。系统的熵和我们预测的一样在随时间增大，但当所有分子布满整个边界以内区域时，系统的熵开始趋于稳定。从这一结果中我们可以了解：当所有的分子随机的布满整个区域时，虽然当我们跟踪某一个确定分子时，它还是在区域内到处乱窜，但每个小区域内分子的密度却不会再变化了。所以，一旦气体分子扩散到整个区域以后，不管我们再等上多少时间，系统的熵都不会再有太大的起伏。换句话说，让系统自动回到开始的状态，即所有分子都在原点的状态，已经不可能了。画一个圆晶粒USEMSFLIB!INTEGERXR,YR!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(1:XR,1:YR)X0=XR/2!圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态（颜色）DOI=1,XRDOJ=1,YRIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J)+3)IER=SETPIXEL(I,J)ENDDOENDDOEND画晶界!画一个圆USEMSFLIBINTEGERXR，YR!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNSXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2!圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态（颜色）DOI=1,XRDOJ=1,YRIF(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)ENDIFIER=SETPIXEL(I,J)ENDDOENDDOEND如何画有一定宽度的晶界？