计算材料学概述-之--蒙特卡洛方法.详解优秀PPT.ppt

《计算材料学概述-之--蒙特卡洛方法.详解优秀PPT.ppt》由会员分享，可在线阅读，更多相关《计算材料学概述-之--蒙特卡洛方法.详解优秀PPT.ppt（112页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、计算材料学概述计算材料学概述第三章第三章蒙特卡罗方法蒙特卡罗方法（Monte Carlo）主要内容主要内容oMonte Carlo模拟发展简介模拟发展简介oMonte Carlo模拟基本原理模拟基本原理oMonte Carlo模拟典型算法模拟典型算法oMonte Carlo模拟典型应用模拟典型应用 蒙特卡洛法是什么？蒙特卡洛法是什么？o蒙特卡洛蒙特卡洛(Monte Carlo)方法，是在简洁的理论准方法，是在简洁的理论准则基础上，接受反复随即抽样的方法，解决困难系则基础上，接受反复随即抽样的方法，解决困难系统的问题。其实质是一种概率和统计的问题。统的问题。其实质是一种概率和统计的问题。o 蒙特

2、蒙特卡罗方法（卡罗方法（Monte Carlo method），），也称统计模拟方法，是二十世纪四十年头中期由于也称统计模拟方法，是二十世纪四十年头中期由于科学技术的发展和电子计算机的独创，而被提出的科学技术的发展和电子计算机的独创，而被提出的一种以概率统计理论为指导的一类特别重要的数值一种以概率统计理论为指导的一类特别重要的数值计算方法。是指运用随机数（或更常见的伪随机数）计算方法。是指运用随机数（或更常见的伪随机数）来解决很多计算问题的方法。来解决很多计算问题的方法。MC的基本思想 MC MC基本思想很早以前就被人们所发觉和利基本思想很早以前就被人们所发觉和利用。用。1717世纪，人们就知

3、道用事务发生的世纪，人们就知道用事务发生的“频率频率”来确定事务的来确定事务的“概率概率”。但要真正实现随机。但要真正实现随机抽样是很困难的，甚至几乎是不行能的。抽样是很困难的，甚至几乎是不行能的。高速计算机的出现，使得用数学方法在高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可计算机上大量、快速地模拟这样的试验成为可能。能。确定性系统确定性系统随机性系统随机性系统模拟模拟自然界自然界Monte-Carlo模拟，即随机模拟（重复模拟，即随机模拟（重复“试验试验”）重复试验重复试验计算机模拟计算机模拟Monte Carlo方法：方法：亦称统计模拟方法，亦称统计模拟方法，

4、statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte Carlo名字的由来：名字的由来：是由是由Metropolis在二次世界大战期间提出的：在二次世界大战期间提出的：Manhattan支支配，探讨与原子弹有关的中子输运过程；配，探讨与原子弹有关的中子输运过程；Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博著名的首都，该城以赌博著名Nicholas Metropolis(1915-1999)Monte-Carlo,MonacoMonte CarloMonte Carlo方法简史方法简史简洁地介

5、绍一下简洁地介绍一下Monte CarloMonte Carlo方法的发展历史方法的发展历史1 1、BuffonBuffon投针试验：投针试验：18世纪，法国数学家世纪，法国数学家Comte de Buffon利用投针试验估计利用投针试验估计 的值的值dL1777年法国科学家布丰提出的一种计算圆周率的方法年法国科学家布丰提出的一种计算圆周率的方法随机投随机投针法，即著名的布丰投针问题。这一方法的步骤是：针法，即著名的布丰投针问题。这一方法的步骤是：1)取一张白纸，在上面画上很多条间距为取一张白纸，在上面画上很多条间距为d的平行线。的平行线。2)取一根长度为取一根长度为l（ld）的针，随机地向画

6、有平行直线的纸上的针，随机地向画有平行直线的纸上掷掷n次，视察针与直线相交的次数，记为次，视察针与直线相交的次数，记为m 3）计算针与直线相交的概率）计算针与直线相交的概率 布丰本人证明白，这个概率是布丰本人证明白，这个概率是:p=2l/(d)，为圆周率为圆周率：利用这个公式可以用概率的方法得到圆周率的近似值。下面是利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料一些资料 试验者年头投掷次数相交次数圆周率估计值沃尔夫1850500025313.1596史密斯1855320412193.1554德摩根18806003833.137福克斯188410304893.1595拉泽里尼19

7、01340818083.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)重设该生

8、成函数到初始状态留意：上述随机数序列均具周期性，如上页random子程序的周期约230。实例一、计算值计算过程：计算过程：1、构造或描述问题的概率过程、构造或描述问题的概率过程2、从概率密度函数动身进行随机抽样，实现、从概率密度函数动身进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些从已知概率分布的抽样，得到特征量的一些模拟结果模拟结果计算均值计算均值Monte Carlo方法之方法之典型算法与应用o考虑平面上的一个边长为考虑平面上的一个边长为1 1的正方形及其内部的一个形态不规的正方形及其内部的一个形态不规则的则的“图形图形”，如何求出这个，如何求出这个“图形图形”的面积呢？的面积呢

9、？Monte Monte CarloCarlo方法是这样一种方法是这样一种“随机化随机化”的方法：向该正方形的方法：向该正方形“随机随机地地”投掷投掷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，该面积的值为

10、的值。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(*,*)PIENDo面积的计算面积的计算f(x)x辛普逊方法辛普逊方法I=Sn蒙特蒙特-卡洛方法卡洛方法f(x)x在长方形中均匀投在长方形中均匀投N0组组(x,y)如如 yf(x),则则 N=N+1I=(N/N0)S0SS011MC的优点oMCMC与传统数学方法相比，具有直观性强，简便易行的优点，该方与传统数学方法相比，具有直观性强，简便易

11、行的优点，该方法能处理一些其他方法无法解决的负责问题，并且简洁在计算机法能处理一些其他方法无法解决的负责问题，并且简洁在计算机上实现，在很大程度上可以代替很多大型、难以实现的困难试验上实现，在很大程度上可以代替很多大型、难以实现的困难试验和社会行为。无污染、无危急、能摆脱试验误差。和社会行为。无污染、无危急、能摆脱试验误差。Monte Carlo方法利用方法利用随机抽样随机抽样的方法来求解物理问题的方法来求解物理问题;数值解法数值解法:从一个物理系统的数学模型动身从一个物理系统的数学模型动身,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态;Mont

12、e Carlo方法方法并非只能用来解决包含随机过程并非只能用来解决包含随机过程的问题的问题:例例如如:用用Monte Carlo方法计算定积分方法计算定积分.对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程,然后用然后用Monte Carlo方法进行求解方法进行求解留意以下两点：留意以下两点：MC的应用自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；数值分析：数值分析：数学问题，求积分，求逆矩阵，解线性代数等数学问题，求积分，求逆矩阵，解线性代数等经济学模

13、拟：经济学模拟：库存问题，随机服务系统中排队问题库存问题，随机服务系统中排队问题 人口问题：人口的诞生，传染病的扩散；乃至动物的生态竞争人口问题：人口的诞生，传染病的扩散；乃至动物的生态竞争金属学：扩散、组织长大、相金属学：扩散、组织长大、相变过程变过程蒙特-卡洛模拟的意义o能探讨不同边界、不同材料的影响能探讨不同边界、不同材料的影响o 理论不行能、试验耗费太大理论不行能、试验耗费太大o 用于试验设计用于试验设计o无污染无污染o 反应堆防护反应堆防护o 核弹爆炸核弹爆炸o能摆脱试验误差能摆脱试验误差o 作理论和试验的桥梁作理论和试验的桥梁Monte Carlo模拟的步骤：模拟的步骤：1.依据欲

14、探讨的物理系统的性质，建立能够描述该系依据欲探讨的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概统特性的理论模型，导出该模型的某些特征量的概率密度函数；率密度函数；2.（即构造或描述问题的概率过程）（即构造或描述问题的概率过程）2.从概率密度函数动身进行随机抽样，实现从已知概从概率密度函数动身进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；率分布的抽样，得到特征量的一些模拟结果；3.有了明确的概率过程后，为了实现过程的有了明确的概率过程后，为了实现过程的数值模拟，必需实现从已知概率分布的随机数的抽数值模拟，必需实现从已知概率分布的随机数的抽样，进

15、行大量的随机模拟试验，从中获得随机变量样，进行大量的随机模拟试验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是的大量试验值。产生已知概率分布的随机变量，是实现实现MC方法的关键步骤，其中最基本的是（方法的关键步骤，其中最基本的是（0，1）匀整分布。）匀整分布。3.对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。4.模拟结果的检验模拟结果的检验Monte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数(pdf)必需给出描述一个物理系统的一组概必需给出描述一个物理系统的一组概率密度函数率密度函数;随机数产生器随

16、机数产生器能够产生在区间能够产生在区间0,1上（匀整）分布的随机数上（匀整）分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,1上匀整分布的随机数动身上匀整分布的随机数动身,随随机抽取听从给定的机抽取听从给定的pdf的随机变量的随机变量;模拟结果记录模拟结果记录记录一些感爱好的量的模拟结果记录一些感爱好的量的模拟结果误差估计误差估计必需确定统计误差（或方差）随模拟次数以及其必需确定统计误差（或方差）随模拟次数以及其它一些量的变更；它一些量的变更；削减方差的技术削减方差的技术利用该技术可削减模拟过程中计算的次数；利用该技术可削减模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行

17、计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法实例二 定积分计算o事实上，不少的统计问题，如计算概率、各阶距等，最终都归结为定积分的近似计算问题。o下面考虑一个简洁的定积分o!计算计算x*2在在(0,1)上积分上积分o计算过程：计算过程：o1、构造或描述问题的概率过程、构造或描述问题的概率过程:产生听从分布产生听从分布f（x）的随机变量）的随机变量Xi()（i=1,2,N）o2、从概率密度函数动身进行随机抽样，实现、从概率密度函数动身进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些从已知概率分布的抽样，得到特征量的一些模拟结果模拟结果计算均值计算均值()REALYY=0N

18、=300000ISEED=RTC()DOI=1,NX=RAN(ISEED)Y=Y+X*2/NENDDOWRITE(*,*)YENDlimx2dx (dx 0)Monte Carlo方法另一个重要问题：方法另一个重要问题：随机数随机数 随机数：由单位矩阵分布中所产生的简洁子样称为随机数序列，其中的每一个个体称为随机数。但真正的随机数的不适合电子计算机上运用，因为它须要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，且其产生的序列无法重复实现，使程序无法复算，结果无法验证，同时须要增加随机数发生器和电路联系等附加设备。伪随机数：o是有数学递推公式所产生的随机数。（近似的具备随机数的性质

19、。）oAn+1=T(A);An+1=An+k+1o伪随机的优点和缺点：o推断伪随机数好坏的方法：o1、它能够有较好的匀整性和独立性；o2、它的费用大小，即指所消耗计算机的时间；o3、容量要求尽可能大。随机数产生的方法o产生匀整分布随机数的几种方法；o（1）物理方法；o（2）数学方法。o伪随机数产生方法：o加同余法o乘同余法o乘加同余法o取中方法o逆变换法o合成法o筛选法 o。关于随机数的几点留意关于随机数的几点留意o注注1 1 由由于于匀匀整整分分布布的的随随机机数数的的产产生生总总是是接接受受某某个个确确定定的的模模型型进进行行的的，从从理理论论上上讲讲，总总会会有有周周期期现现象象出出现现

20、的的。初初值值确确定定后后，全全部部随随机机数数也也随随之之确确定定，并并不不满满足足真真正正随随机机数数的的要要求求。因因此此通通常常把把由由数数学学方方法法产产生生的的随随机机数数成为伪随机数。成为伪随机数。o 但但其其周周期期又又相相当当长长，在在实实际际应应用用中中几几乎乎不不行行能能出出现现。因因此此，这这种种由由计计算算机机产产生生的的伪伪随随机机数数可可以以当当作作真正的随机数来处理。真正的随机数来处理。o注注2 2 应应对对所所产产生生的的伪伪随随机机数数作作各各种种统统计计检检验验，如如独独立性检验，分布检验，功率谱检验等等。立性检验，分布检验，功率谱检验等等。FORTRAN

21、语言产生随机数的实例orandom_number(x)random_number(x)产生一个产生一个0 0到到1 1之间的随机数（之间的随机数（x x可以是向量），可以是向量），但是每次总是那几个数。但是每次总是那几个数。用了用了random_seed()random_seed()后，系统依据日期和时间随机地供应种子，后，系统依据日期和时间随机地供应种子，使得随机数更随机了。使得随机数更随机了。programrandomprogramrandomreal:xreal:xcallrandom_seed()!callrandom_seed()!系统依据日期和时间随机地供应种子系统依据日期和时间随

22、机地供应种子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)

23、间随机分布的N*N矩阵s=rand(state,0)重设该生成函数到初始状态留意：上述随机数序列均具周期性，如上页random子程序的周期约230。FinitedifferenceapproximationofdifferentialequationsoAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtaime-CentralspaceFICK其次定律Fick其次定律稳态扩散解REALC0(0:1000+1),

24、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个人在一起，各自的生日为个人在

25、一起，各自的生日为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

26、+1)JJJ=1DOI=1,365 Y=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)(Unrestri

27、cted randon walk,RW)startend平均平方端平均平方端-端位移端位移:,自然科学和社会生活中很多现象都与随机运动有关自然科学和社会生活中很多现象都与随机运动有关可以模拟的内容？扩散；分子运动；。o如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了全部分子距原点的方均距离：！MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XX(J,I)：X

28、*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)

29、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

30、=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)ENDo!MonteCarloSimulationofTwoDimensiona

31、lDiffusionooINTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNoREALXYM(1:1000)o!X:INSTANTANEOUSPOSITIONOFATOMo!XY(J,I)：X*Y,J:第几天试验,I:第几步跳动o!XYM(I):THEMEANOFXYoWRITE(*,*)试验天数JMAX，试验次数IMAXoREAD(*,*)JMAX,IMAXoXN=(/0,0,-1,1/)oYN=(/-1,1,0,0/)oISEED=RTC()oDOJ=1,JMAX!第几天试验oX=0!oY=0!oDOI=1,IMAX!第几步跳动oRN=4*RAN(

32、ISEED)+1oX=X+XN(RN)oY=Y+YN(RN)oXY(J,I)=X*X+Y*YoENDDOoENDDOoOPEN(1,FILE=C:DIF2.DAT)oDOI=1,IMAXoXYM=0.0oXYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX!oWRITE(1,*)I,XYM(I)oENDDOoCLOSE(1)！做三维空间随机行走？oEND实际环境是困难的实际环境是困难的:可变，缺陷，风，可变，缺陷，风，季节，电场等。季节，电场等。空间中有随机性缺陷空间中有随机性缺陷不退行走不退行走自回避行走自回避行走o例如模拟高分子的位形o用随机行走方法模拟高分子位形是用随机行走的

33、轨迹代表高分子的位形，行走过的位置代表的是构成分子的原子或官能团，因此，无限制随机行走忽视了体斥效应。o不退行走就是禁止在每一步行走后马上倒退，可以解决刚走的一步与上一步重叠的问题。但不退行走没有完全解决高分子的体斥效应。o自回避行走就是全部已走过的位置不能再走，这样就完全解决了体斥效应问题。o气体分子的随机行走假设在一个空旷封闭的房间中心滴上一滴香水，挥发的香水分子随机的与空气中的粒子发生碰撞，最终会扩散到整个房间。这里我们将通过计算机模拟400个分子在二维平面内的随机运动探讨熵变o现在来探讨一下熵在我们这个气体扩散模型中是什么意思。在刚起先的时候，我们的系统处在一个特别有序的状态：全部的分

34、子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随着时间的推移，分子起先不断的向外扩散，这时系统出现了无序状态，熵起先渐渐增加。o系统的熵和我们预料的一样在随时间增大，但当全部分子布满整个边界以内区域时，系统的熵起先趋于稳定。从这一结果中我们可以了解：当全部的分子随机的布满整个区域时，虽然当我们跟踪某一个确定分子时，它还是在区域内到处乱窜，但每个小区域内分子的密度却不会再变更了。所以，一旦气体分子扩散到整个区域以后，不管我们再等上多少时间，系统的熵都不会再有太大的起伏。换句话说，让系统自动回到起先的状态，即全部分子都在原点的状态，已经不行能了。画一个圆晶粒ooUSEMSFLIBo!IN

35、TEGERXR,YR!在的区域中画一个圆oPARAMETERXR=400,YR=400oINTEGERR,S(1:XR,1:YR)oX0=XR/2!圆心位置X0，YOoY0=YR/2oR=MIN(X0-10,Y0-10)!圆半径oS=0!像素的初始状态（颜色）oDOI=1,XRoDOJ=1,YRoIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10oIER=SETCOLOR(S(I,J)+3)oIER=SETPIXEL(I,J)oENDDOoENDDOoEND画晶界o!画一个圆oUSEMSFLIBoINTEGERXR，YR!在的区域中画一个圆oPARAMETERXR=400,YR

36、=400oINTEGERR,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNSoXN=(/0,0,-1,1/)oYN=(/-1,1,0,0/)oX0=XR/2!圆心位置X0，YOoY0=YR/2oR=MIN(X0-10,Y0-10)!圆半径oS=0!像素的初始状态（颜色）oDOI=1,XRoDOJ=1,YRoIF(I-X0)*2+(J-Y0)*20)THENoIER=SETCOLOR(9)oELSEoIER=SETCOLOR(8)oENDIFoIER=SETPIXEL(I,J)oENDDOoENDDOoENDo如何画有确定宽度的晶界？例题例求前N个自然数倒数和。即求 程序

37、须要几个变量：累加的项数I(整型)，存放通项值的T(实型)，存放总和的S(实型)，待输入的总项数N(整型)。PROGRAM MAINPROGRAM MAIN IMPLICIT NONE IMPLICIT NONE INTEGER N,I INTEGER N,I REAL T,S REAL T,S READ*,N READ*,N S=0.0;I=1 S=0.0;I=1 DO DO T=1.0/I T=1.0/I S=S+T S=S+T I=I+1 I=I+1 IF(IN)EXIT IF(IN)EXIT END DO END DO PRINT*,SUM=,S PRINT*,SUM=,SEND PR

38、OGRAM MAINEND PROGRAM MAIN 程序的运行结果如下15 SUM=3.182290用DO结构处理循环问题，最关键的是要处理好初值、循环终止条件和循环体语句的关系。FinitedifferenceapproximationofdifferentialequationsoAdifferentialequationcanbeapproximatedbyafinitedifferencescheme.ForexampleForwardtimeBackwardspaceCentralspaceForwardtime-CentralspaceFICK其次定律Fick其次定律稳态扩散解R

39、EALC0(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.8C(1001)=0.1IF(JT=50000)WRITE(1,*)IX,C(IX)ENDDOC0=CENDDOPRINT*,PROGRAMMEISFINISHE

40、DCLOSE(1)END蒙特卡罗蒙特卡罗(Monte Carlo)方法方法1、简介2、相关几个例子Monte Carlo方法：方法：亦称统计模拟方法，亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte Carlo名字的由来：名字的由来：是由是由Metropolis在二次世界大战期间提出的：在二次世界大战期间提出的：Manhattan支支配，探讨与原子弹有关的中子输运过程；配，探讨与原子弹有关的中子输运过程；Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博著名的首都，该城以赌博

41、著名Nicholas Metropolis(1915-1999)Monte-Carlo,MonacooMonte Carlo方方法法的的基基本本思思想想很很早早以以前前就就被被人人们们所所发发觉觉和和利利用用。早早在在17世世纪纪，人人们们就就知知道道用用事事务务产产生生的的“频频率率”来来近近似似事事务务的的“概概 率率”。18世世纪纪人人们们用用投投针针试试验验的的方方法法来来确确定定圆圆周周率率。本本世世纪纪40年年头头电电子子计计算算机机的的出出现现，特特殊殊是是近近年年来来高高速速电电子子计计算算机机的的出出现现，使使得得用用数数学学方方法法在在计计算算机机上上大大量量、快快速速地地

42、模模拟这样的试验成为可能。拟这样的试验成为可能。1930年，年，Enrico Fermi利用利用Monte Carlo方法探讨中子的扩散，方法探讨中子的扩散，并设计了一个并设计了一个Monte Carlo机械装置，机械装置，Fermiac,用于计算核用于计算核反应堆的临界状态反应堆的临界状态Von Neumann是是Monte Carlo方法的正式奠基者方法的正式奠基者,他与他与Stanislaw Ulam合作建立了概率密度函数、反累积分布函数的数学基础，合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中，以及伪随机数产生器。在这些工作中，Stanislaw U

43、lam意识意识到了数字计算机的重要性到了数字计算机的重要性合作起源于合作起源于Manhattan工程：利用工程：利用ENIAC(Electronic Numerical Integrator and Computer)计算产额计算产额Monte Carlo模拟的应用：模拟的应用：自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；试验探测器的模拟试验探测器的模拟数值分析：数值分析：利用利用Monte Carlo方法求积分方法求积分材料学中的应用：材料学中的应用：扩散、晶粒的长大、再结晶等

44、扩散、晶粒的长大、再结晶等Monte Carlo模拟在物理探讨中的作用模拟在物理探讨中的作用Monte Carlo模拟的步骤：模拟的步骤：1.依据欲探讨的物理系统的性质，建立能够描述该系统特性依据欲探讨的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；的理论模型，导出该模型的某些特征量的概率密度函数；2.从概率密度函数动身进行随机抽样，得到特征量的一些模从概率密度函数动身进行随机抽样，得到特征量的一些模拟结果；拟结果；3.对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。留意以下两点：留意以下两点：Monte

45、 Carlo方法与数值解法的不同方法与数值解法的不同:Monte Carlo方法利用随机抽样的方法来求解物理问题方法利用随机抽样的方法来求解物理问题;数值解法数值解法:从一个物理系统的数学模型动身从一个物理系统的数学模型动身,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态;Monte Carlo方法并非只能用来解决包含随机的过程的问题方法并非只能用来解决包含随机的过程的问题:很多利用很多利用Monte Carlo方法进行求解的问题中并不包含随方法进行求解的问题中并不包含随机过程机过程 例如例如:用用Monte Carlo方法计算定积分方法计算定积分

46、.对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程,然后用然后用Monte Carlo方法进行求解方法进行求解Monte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数 必需给出描述一个物理系统的一组概率密必需给出描述一个物理系统的一组概率密度函数度函数;随机数产生器随机数产生器能够产生在区间能够产生在区间0,1上匀整分布的随机数上匀整分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,1上匀整分布的随机数动身上匀整分布的随机数动身,随随机抽取听从给定的机抽取听从给定的pdf的随机变量的随机变量;模拟结果记录模拟结果记录记录一些感爱好的量

47、的模拟结果记录一些感爱好的量的模拟结果误差估计误差估计必需确定统计误差（或方差）随模拟次数以及其必需确定统计误差（或方差）随模拟次数以及其它一些量的变更；它一些量的变更；削减方差的技术削减方差的技术利用该技术可削减模拟过程中计算的次数；利用该技术可削减模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法确定性系统确定性系统随机性系统随机性系统模拟模拟自然界自然界Monte-Carlo模拟，即随机模拟（重复模拟，即随机模拟（重复“试验试验”）重复试验重复试验计算机模拟计算机模拟蒙特卡罗方法的基本原理及思想o当所要求解的问题是某

48、种事务出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事务出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。o把蒙特卡罗解题归结为三个主要步骤：o构造或描述概率过程；o实现从已知概率分布抽样；o建立各种估计量。o对于本身就具有随机性质的问题，主要是正确描述和模拟这个概率过程；对于不是随机性质的确定性问题，比如计算定积分，就必需事先构造一个人为的概率过程。最简洁、最基本、最重要的一个概率分布是(0,1)上的匀整分布（或称矩形分布）。o实现从已知概率分布抽样(模拟试验)：构造了概率模型以后，产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法

49、模拟试验的基本手段，这也是蒙特卡罗方法被称为随机抽样的缘由。o建立各种估计量：构造了概率模型并能从中抽样后，即实现模拟试验后，就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。o收敛性收敛性:由大数定律,Monte-Carlo模拟的收敛是以概率而言的o误差误差:用频率估计概率时误差的估计，可由中心极限定理，给定置信水平 的条件下，有：oo模拟次数模拟次数:由误差公式得随机数的产生原理随机数的产生原理-计算机抽样试验计算机抽样试验o 随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于匀整分布U(0,1的随机数。o首先，介绍听从匀整分布U(0,1的随机数的产

50、生方法。o其次，介绍听从其他各种分布的随机数的产生方法。以及听从正态分布的随机数的产生方法。o最终，关于随机数的几点注。(0,1匀整分布随机数的产生匀整分布随机数的产生o最常用的是在(0,1区间内匀整分布的随机数，其他分布的随机数可利用匀整分布随机数来产生。o匀整分布随机数的产生：乘同余法o产生区间(0,1内匀整分布的随机数的递推公式是：oo其中，是乘因子，M是模数。第一式称做以M为模数（modulus）的同余式，即以M除后得到的余数记为。当给定了一个初值（称为种子），计算出的即在(0,1上匀整分布的随机数。例1：oR=RAND()oISEED=RTC()R=RAN(ISEED)其他各种分布的