计算材料学-2

计算材料学ComputationalMaterialsScience第二讲蒙特卡罗（MonteCarlo）方法主讲：张晖电话：[phone]Email:[e-mail]本堂课主要内容MonteCarlo模拟发展简介MonteCarlo模拟基本原理MonteCarlo模拟典型算法MonteCarlo模拟典型应用蒙特卡洛法是什么？蒙特卡洛(MonteCarlo)方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。蒙特·卡罗方法（MonteCarlomethod），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。MC的基本思想MC基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。但要真正实现随机抽样是很困难的，甚至几乎是不可能的。高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。确定性系统随机性系统模拟自然界Monte-Carlo模拟，即随机模拟（重复“试验”）重复试验计算机模拟MonteCarlo方法：亦称统计模拟方法，statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来：•是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；•MonteCarlo是摩纳哥（monaco)的首都，该城以赌博闻名NicholasMetropolis(1915-1999)Monte-Carlo,MonacoMonteCarlo方法简史简单地介绍一下MonteCarlo方法的发展历史1、Buffon投针实验：18世纪，法国数学家ComtedeBuffon利用投针实验估计的值dLp2dL1777年法国科学家布丰提出的一种计算圆周率的方法——随机投针法，即著名的布丰投针问题。这一方法的步骤是：1)取一张白纸，在上面画上许多条间距为d的平行线。2)取一根长度为l（ld）的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m3）计算针与直线相交的概率．布丰本人证明了，这个概率是:p=2l/(πd)，π为圆周率：利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料实验者年代投掷次数相交次数圆周率估计值沃尔夫1850500025313.1596史密斯1855320412193.1554德摩根18806003833.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=πR^2,当R=1时，S=π。2、由于圆的方程是：x^2+y^2=1(x^2为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的优点MC与传统数学方法相比，具有直观性强，简便易行的优点，该方法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。MonteCarlo方法利用随机抽样的方法来求解物理问题;数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;MonteCarlo方法并非只能用来解决包含随机过程的问题:例如:用MonteCarlo方法计算定积分.对这样的问题可将其转换成相关的随机过程,然后用MonteCarlo方法进行求解注意以下两点：MC的应用自然现象的模拟：宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；数值分析：数学问题，求积分，求逆矩阵，解线性代数等经济学模拟：库存问题，随机服务系统中排队问题人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争金属学：扩散、组织长大、相变过程蒙特-卡洛模拟的意义能研究不同边界、不同材料的影响理论不可能、实验耗费太大用于实验设计无污染反应堆防护核弹爆炸能摆脱实验误差作理论和实验的桥梁MonteCarlo模拟的步骤：1.根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；（即构造或描述问题的概率过程）2.从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是实现MC方法的关键步骤，其中最基本的是（0，1）均匀分布。3.对模拟结果进行分析总结，预言物理系统的某些特性。4.模拟结果的检验MonteCarlo算法的主要组成部分概率密度函数(pdf)—必须给出描述一个物理系统的一组概率密度函数;随机数产生器—能够产生在区间[0,1]上（均匀）分布的随机数抽样规则—如何从在区间[0,1]上均匀分布的随机数出发,随机抽取服从给定的pdf的随机变量;模拟结果记录—记录一些感兴趣的量的模拟结果误差估计—必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；减少方差的技术—利用该技术可减少模拟过程中计算的次数；并行和矢量化—可以在先进的并行计算机上运行的有效算法实例二定积分计算事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分dxxfba!计算x**2在(0,1)上积分计算过程：1、构造或描述问题的概率过程:产生服从分布f（x）的随机变量Xi()（i=1,2,···,N）2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果——计算均值()niix12nixiy12REALYY=0N=300000ISEED=RTC()DOI=1,NX=RAN(ISEED)Y=Y+X**2/NENDDOWRITE(*,*)YENDlim∑x2dx(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(*,*)xstopendprogramrandomMonteCarlo方法之随机数的产生许多计算机系统都有随机数生成函数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.ForexampletuudxdutttForwardtimexuudxdutxxtx2222xuuudxudtxxtxtxxBackwardspaceCentralspace21222xuuutuuxututxxtxtxxtxtxForwardtaime-CentralspaceFICK第二定律21222xuuutuuxututxxtxtxxtxtxvxCCCDtCCxCDttxxtxtxxtxtx21222C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)=