《计算材料学》PPT课件

计算材料学ComputationalMaterialsScience,第二讲蒙特卡罗（MonteCarlo）方法,主讲：张晖电话mail:huizhang,本堂课主要内容,MonteCarlo模拟发展简介MonteCarlo模拟基本原理MonteCarlo模拟典型算法MonteCarlo模拟典型应用,蒙特卡洛法是什么？,蒙特卡洛(MonteCarlo)方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。蒙特卡罗方法（MonteCarlomethod），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。,MC的基本思想,MC基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。但要真正实现随机抽样是很困难的，甚至几乎是不可能的。高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，即随机模拟（重复“试验”）,重复试验,计算机模拟,MonteCarlo方法：,亦称统计模拟方法，statisticalsimulationmethod利用随机数进行数值模拟的方法,MonteCarlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,MonteCarlo是摩纳哥（monaco)的首都，该城以赌博闻名,NicholasMetropolis(1915-1999),Monte-Carlo,Monaco,MonteCarlo方法简史,简单地介绍一下MonteCarlo方法的发展历史,1、Buffon投针实验：,18世纪，法国数学家ComtedeBuffon利用投针实验估计的值,1777年法国科学家布丰提出的一种计算圆周率的方法随机投针法，即著名的布丰投针问题。这一方法的步骤是：1)取一张白纸，在上面画上许多条间距为d的平行线。2)取一根长度为l（ld）的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m3）计算针与直线相交的概率布丰本人证明了，这个概率是:p=2l/(d)，为圆周率：利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料,实验者年代投掷次数相交次数圆周率估计值沃尔夫1850500025313.1596史密斯1855320412193.1554德摩137福克斯188410304893.1595拉泽里尼1901340818083.1415929赖纳192525208593.1795布丰投针实验是第一个用几何形式表达概率问题的例子，他首次使用随机实验处理确定性数学问题，为概率论和蒙特卡罗方法的发展起到一定的推动作用。,MonteCarlo方法之随机数的产生许多计算机系统都有随机数生成函数,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、构造或描述问题的概率过程2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果计算均值,MonteCarlo方法之典型算法与应用,考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？MonteCarlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点，若有M个点落于“图形”内，则该“图形”的面积近似为M/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,蒙特-卡洛方法,11,MC的优点,MC与传统数学方法相比，具有直观性强，简便易行的优点，该方法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。,MonteCarlo方法利用随机抽样的方法来求解物理问题;,数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;,MonteCarlo方法并非只能用来解决包含随机过程的问题:例如:用MonteCarlo方法计算定积分.对这样的问题可将其转换成相关的随机过程,然后用MonteCarlo方法进行求解,注意以下两点：,MC的应用,自然现象的模拟：,数值分析：,数学问题，求积分，求逆矩阵，解线性代数等,经济学模拟：库存问题，随机服务系统中排队问题,人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争,金属学：扩散、组织长大、相变过程,蒙特-卡洛模拟的意义,能研究不同边界、不同材料的影响理论不可能、实验耗费太大用于实验设计无污染反应堆防护核弹爆炸能摆脱实验误差作理论和实验的桥梁,MonteCarlo模拟的步骤：,根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；（即构造或描述问题的概率过程）,从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是实现MC方法的关键步骤，其中最基本的是（0，1）均匀分布。,对模拟结果进行分析总结，预言物理系统的某些特性。,4.模拟结果的检验,MonteCarlo算法的主要组成部分,概率密度函数(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(*,*)YEND,limx2dx(dx0),MonteCarlo方法另一个重要问题：随机数,随机数：由单位矩阵分布中所产生的简单子样称为随机数序列，其中的每一个个体称为随机数。但真正的随机数的不适合电子计算机上使用，因为它需要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，且其产生的序列无法重复实现，使程序无法复算，结果无法验证，同时需要增添随机数发生器和电路联系等附加设备。,伪随机数：,是有数学递推公式所产生的随机数。（近似的具备随机数的性质。）An+1=T(A);An+1=An+k+1伪随机的优点和缺点：判断伪随机数好坏的方法：1、它能够有较好的均匀性和独立性；2、它的费用大小，即指所消耗计算机的时间；3、容量要求尽可能大。,随机数产生的办法,产生均匀分布随机数的几种方法；（1）物理方法；（2）数学方法。伪随机数产生方法：加同余法乘同余法乘加同余法取中方法逆变换法合成法筛选法。,关于随机数的几点注意,注1由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,FORTRAN语言产生随机数的实例,random_number(x)产生一个0到1之间的随机数（x可以是向量），但是每次总是那几个数。用了random_seed()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。programrandomreal:xcallrandom_seed()!系统根据日期和时间随机地提供种子callrandom_number(x)!每次的随机数就都不一样了write(*,*)xstopendprogramrandom,MonteCarlo方法之随机数的产生许多计算机系统都有随机数生成函数,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。,Finitedifferenceapproximationofdifferentialequations,Adifferentialequationcanbeapproximatedbyafinitedifferencescheme.Forexample,Forwardtime,Backwardspace,Centralspace,Forwardtaime-,Centralspace,FICK第二定律,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理论几率模拟几率0.1170.1100.4110.4120.5270.5200.7060.6920.9410.936500.9860.987,INTEGERM(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/10000END,MC在材料学领域的应用随机行走,背景,如，布朗运动最简单、无限制随机行走(Unrestrictedrandonwalk,RW),平均平方端-端位移:,自然科学和社会生活中很多现象都与随机运动有关,可以模拟的内容？扩散；分子运动；。,如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了所有分子距原点的方均距离：,！MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,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,!MonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,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个分子在二维平面内的随机运动,讨论熵变,现在来讨论一下熵在我们这个气体扩散模型中是什么意思。在刚开始的时候，我们的系统处在一个非常有序的状态：所有的分子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随