计算材料学概述 之蒙特卡洛方法ppt课件.ppt

《计算材料学概述 之蒙特卡洛方法ppt课件.ppt》由会员分享，可在线阅读，更多相关《计算材料学概述 之蒙特卡洛方法ppt课件.ppt（112页珍藏版）》请在三一办公上搜索。

1、计算材料学概述,第三章蒙特卡罗方法（Monte Carlo）,主要内容,Monte Carlo模拟发展简介Monte Carlo模拟基本原理Monte Carlo模拟典型算法Monte Carlo模拟典型应用,蒙特卡洛法是什么？,蒙特卡洛(Monte Carlo)方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。 蒙特卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机

2、数）来解决很多计算问题的方法。,MC的基本思想,MC基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。但要真正实现随机抽样是很困难的，甚至几乎是不可能的。 高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，即随机模拟（重复“试验”）,重复试验,计算机模拟,Monte Carlo方法：,亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis

3、在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法简史,简单地介绍一下Monte Carlo方法的发展历史,1、Buffon投针实验：,18世纪，法国数学家Comte de Buffon利用投针实验估计的值,1777年法国科学家布丰提出的一种计算圆周率的方法随机投针法，即著名的布丰投针问题。这一方法的步骤是： 1) 取一张白纸，在上面画上许多条间距为d的平行线。

4、2) 取一根长度为l（ld） 的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m 3）计算针与直线相交的概率 布丰本人证明了，这个概率是: p=2l/(d) ，为圆周率 ：利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料,实验者 年代 投掷次数 相交次数 圆周率估计值 沃尔夫 1850 5000 2531 3.1596 史密斯 1855 3204 1219 3.1554 德摩根 1880 600 383 3.137 福克斯 1884 1030 489 3.1595 拉泽里尼 1901 3408 1808 3.1415929 赖纳 1925 2520 859 3

5、.1795 布丰投针实验是第一个用几何形式表达概率问题的例子，他首次使用随机实验处理确定性数学问题，为概率论和蒙特卡罗方法的发展起到一定的推动作用。,Monte Carlo方法之随机数的产生 许多计算机系统都有随机数生成函数,F90: call random_seed call random_number(a)2、ISEED=RTC() X=RAN(ISEED)Y=RAN(ISEED）,Matlab: x=rand(N) 产生元素在(0, 1)间随机分布的N*N矩阵 s= rand(state,0) 重设该生成函数到初始状态,注意：上述随机数序列均具周期性，如上页random子程序的周期约23

6、0。,实例一、计算值,计算过程：1、构造或描述问题的概率过程2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果计算均值,Monte Carlo方法之典型算法与应用,考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点，若有M个点落于“图形”内，则该“图形”的面积近似为M/N。,用该方法计算的基本思路是：1、根据圆面积的公式：s=R2,当R=1时，S=。2、由于圆的方程是：x2+y2=1(x2为x的平方的意思),因此1/4圆面积为x

7、轴、y轴和上述方程所包围的部分。3、如果在1*1的正方形中均匀地落入随机点，则落入1/4圆中的点的概率就是1/4圆的面积。其4倍，就是圆面积。由于半径为1，该面积的值为的值。,REAL R,R1,R2,PIISEED=RTC()N0=0N=300000DO I=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R1.0)N0=N0+1END DOPI=4.0*N0/NWRITE(*,*)PIEND,面积的计算,f (x),x,辛普逊方法,I = Sn,蒙特-卡洛方法,11,MC的优点,MC与传统数学方法相比，具有直观性强，简便易行的优点，该方

8、法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。,Monte Carlo方法利用随机抽样的方法来求解物理问题;,数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;,Monte Carlo方法并非只能用来解决包含随机过程的问题: 例如:用Monte Carlo方法计算定积分. 对这样的问题可将其转换成相关的随机过程, 然后用Monte Carlo方法进行求解,注意以下两点：,MC的应用,自然现象的模拟：,数值分析：,数学问题，求积分，求逆矩阵，解线性代

9、数等,经济学模拟：库存问题，随机服务系统中排队问题,人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争,金属学：扩散、组织长大、相变过程,蒙特-卡洛模拟的意义,能研究不同边界、不同材料的影响 理论不可能、实验耗费太大 用于实验设计无污染 反应堆防护 核弹爆炸能摆脱实验误差 作理论和实验的桥梁,Monte Carlo模拟的步骤：,根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数； （即构造或描述问题的概率过程）,从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果； 有了明确的概率过程后，为了实现过程的数值模拟，

10、必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是实现MC方法的关键步骤，其中最基本的是（0，1）均匀分布。,对模拟结果进行分析总结，预言物理系统的某些特性。,4. 模拟结果的检验,Monte Carlo算法的主要组成部分,概率密度函数(pdf) 必须给出描述一个物理系统的一组概率密度函数;,随机数产生器能够产生在区间0,1上（均匀）分布的随机数,抽样规则如何从在区间0,1上均匀分布的随机数出发,随机抽取服从给定的pdf的随机变量;,模拟结果记录记录一些感兴趣的量的模拟结果,误差估计必须确定统计误差（或方差）随模拟次数以及

11、其它一些量的变化；,减少方差的技术利用该技术可减少模拟过程中计算的次数；,并行和矢量化可以在先进的并行计算机上运行的有效算法,实例二 定积分计算,事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分,! 计算x*2在(0,1)上积分计算过程：1、构造或描述问题的概率过程:产生服从分布f（x）的随机变量Xi( )（i=1,2, ,N）2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果计算均值( ),REAL YY=0N=300000ISEED=RTC()DO I=1,N X=RAN(ISEED) Y=Y+X*

12、2/N END DOWRITE(*,*)YEND,limx2dx (dx 0),Monte Carlo方法另一个重要问题：随机数,随机数：由单位矩阵分布中所产生的简单子样称为随机数序列，其中的每一个个体称为随机数。 但真正的随机数的不适合电子计算机上使用，因为它需要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，且其产生的序列无法重复实现，使程序无法复算，结果无法验证，同时需要增添随机数发生器和电路联系等附加设备。,伪随机数：,是有数学递推公式所产生的随机数。（近似的具备随机数的性质。） An+1=T(A);An+1= An+k+1伪随机的优点和缺点：判断伪随机数好坏的方法： 1、

13、它能够有较好的均匀性和独立性； 2、它的费用大小，即指所消耗计算机的时间； 3、容量要求尽可能大。,随机数产生的办法,产生均匀分布随机数的几种方法； （1）物理方法；（2）数学方法。伪随机数产生方法：加同余法乘同余法乘加同余法取中方法逆变换法合成法筛选法 。,关于随机数的几点注意,注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2 应

14、对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,FORTRAN 语言产生随机数的实例,random_number(x)产生一个0到1之间的随机数（x可以是向量），但是每次总是那几个数。用了random_seed()后，系统根据日期和时间随机地提供种子，使得随机数更随机了。programrandomreal:xcallrandom_seed()!系统根据日期和时间随机地提供种子callrandom_number(x)!每次的随机数就都不一样了write(*,*)xstopendprogramrandom,Monte Carlo方法之随机数的产生 许多计算机系统都有随机

15、数生成函数,F90: call random_seed call random_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。,Finite difference approximation of differential equations,A differential equation can be approximated by

16、 a finite difference scheme. For example,Forward time,Backward space,Central space,Forward t aime-,Central space,FICK 第二定律,Fick 第二定律稳态扩散解,REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE)C0=0.1C0(0)=0.8 !BOUNDARY CONDITIONC0(1001)=0.1 !BOUNDARY CONDITIONC=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX

17、=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.1 IF(JT=50000) WRITE(1,*) IX,C(IX) END DO C0=C END DOEND,应用之二 生日问题,MC模拟,假设有n个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50的几率人群中至少有2个人生日相同。,n 理论几率 模拟几率0.117 0.1100.411 0.4120.527 0.5200.706 0.6920.941 0.93650 0.98

18、6 0.987,INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC()DO J=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1 DO I=1,365 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR=1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DOEND DO DO I=1,

19、10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END,MC在材料学领域的应用 随机行走,背景,如，布朗运动最简单、无限制随机行走(Unrestricted randon walk,RW ),平均平方端-端位移:,自然科学和社会生活中很多现象都与随机运动有关,可以模拟的内容？扩散；分子运动；。,如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了所有分子距原点的方均距离：,！Monte Carlo Simulation of One Dimensional Diffusion INT

20、EGER X,XX(1:1000,1:1000)REAL XXM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃!XXM(I): THE MEAN OF XX WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN0.5)THEN X=X+1ELSE X=X-1 END IF XX(J,I)=X*XEND DOEND DOO

21、PEN(1,FILE=“f:DIF1.DAT)DO I=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX !WRITE(1,*) I, XXM(I)END DOCLOSE(1)END,!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,Y,XY(1:1000,1:1000)REAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃!XYM(I): THE MEAN OF XY WRI

22、TE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THEN x=xy=y-1 END IFIF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=xy=y+1 END IFIF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1y=y END IFIF(RN.GE.0.75)THEN x=x+1y=y END IFXY(J,I)=X*X+Y*YEND DO

23、END DOOPEN(1,FILE=“f:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END,!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃!

24、XYM(I): THE MEAN OF XY WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DO J=1,JMAX !第几天实验X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=4*RAN(ISEED)+1X=X+XN(RN) Y=Y+YN(RN)XY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WR

25、ITE(1,*) I, XYM(I)END DOCLOSE(1) ！做三维空间随机行走？END,实际环境是复杂的:可变，缺陷，风， 季节，电场等。,空间中有随机性缺陷,不退行走,自回避行走,例如模拟高分子的位形 用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形，行走过的位置代表的是构成分子的原子或官能团，因此，无限制随机行走忽略了体斥效应。 不退行走就是禁止在每一步行走后立即倒退，可以解决刚走的一步与上一步重叠的问题。但不退行走没有完全解决高分子的体斥效应。 自回避行走就是所有已走过的位置不能再走，这样就完全解决了体斥效应问题。,气体分子的随机行走 假设在一个空旷封闭的房间中心滴

26、上一滴香水，挥发的香水分子随机的与空气中的粒子发生碰撞，最终会扩散到整个房间。这里我们将通过计算机模拟400个分子在二维平面内的随机运动,讨论熵变,现在来讨论一下熵在我们这个气体扩散模型中是什么意思。在刚开始的时候，我们的系统处在一个非常有序的状态：所有的分子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随着时间的推移，分子开始不断的向外扩散，这时系统出现了无序状态，熵开始逐渐增加。,系统的熵和我们预测的一样在随时间增大，但当所有分子布满整个边界以内区域时，系统的熵开始趋于稳定。从这一结果中我们可以了解：当所有的分子随机的布满整个区域时，虽然当我们跟踪某一个确定分子时，它还是在区域内到

27、处乱窜，但每个小区域内分子的密度却不会再变化了。所以，一旦气体分子扩散到整个区域以后，不管我们再等上多少时间，系统的熵都不会再有太大的起伏。换句话说，让系统自动回到开始的状态，即所有分子都在原点的状态，已经不可能了。,画一个圆晶粒,USE MSFLIB !INTEGER XR,YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R,S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*

28、2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J)+3) IER=SETPIXEL(I,J)END DOEND DOEND,画晶界,!画一个圆USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNS XN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2 ! 圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRI

29、F(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IF IER=SETPIXEL(I,J)END DOEND DOEND如何画有一定宽度的晶界？,例题,例 求前N个自然数倒数和。即求,程序需要几个变量：累加的项数I(整型)，存放通项值的T(实型)，存放总和的S(实型)，待输入的总项数N(整型)。,PROGRAM MAIN IMPLICIT NONE INTEGER N,I REAL T,S READ *,N S=0.0; I=1 DO T=1.0/I S=S+T I=I+1 IF(IN)EXIT END DO PRINT

30、 *,SUM=,SEND PROGRAM MAIN,程序的运行结果如下15 SUM= 3.182290,用DO结构处理循环问题，最关键的是要处理好初值、循环终止条件和循环体语句的关系。,Finite difference approximation of differential equations,A differential equation can be approximated by a finite difference scheme. For example,Forward time,Backward space,Central space,Forward time-,Centra

31、l space,FICK 第二定律,Fick 第二定律稳态扩散解,REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE)C0=0.1C0(0)=0.8 !BOUNDARY CONDITIONC0(1001)=0.1 !BOUNDARY CONDITIONC=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distanceC(IX)=C0(IX)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX)C(0)=0.8C(1001)=0.1 IF(JT=50000)WRITE(1,

32、*) IX,C(IX) END DO C0=C END DOPRINT *, PROGRAMME IS FINISHED CLOSE(1) END,蒙特卡罗(Monte Carlo)方法,1、简介2、相关几个例子,Monte Carlo方法：,亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropol

33、is (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件产生的“频率”来近似事件的“概 率”。 18世纪人们用投针试验的方法来决定圆周率。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,1930年，Enrico Fermi利用Monte Carlo方法研究中子的扩散，并设计了一个Monte Carlo机械装置，Fermiac,用于计算核反应堆的临界状态,Von Neumann是Monte Carlo方法的正

34、式奠基者,他与Stanislaw Ulam合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中， Stanislaw Ulam意识到了数字计算机的重要性,合作起源于Manhattan工程：利用ENIAC(Electronic Numerical Integrator and Computer)计算产额,Monte Carlo模拟的应用：,自然现象的模拟：,实验探测器的模拟,数值分析：,利用Monte Carlo方法求积分,材料学中的应用：,扩散、晶粒的长大、再结晶等,Monte Carlo模拟在物理研究中的作用,Monte Carlo模拟的步骤：,根据欲研究的物理

35、系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；,从概率密度函数出发进行随机抽样，得到特征量的一些模拟结果；,对模拟结果进行分析总结，预言物理系统的某些特性。,注意以下两点：,Monte Carlo方法与数值解法的不同:,Monte Carlo方法利用随机抽样的方法来求解物理问题;,数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;,Monte Carlo方法并非只能用来解决包含随机的过程的问题:,许多利用Monte Carlo方法进行求解的问题中并不包含随机过程 例如:用Monte Carlo方法计算定积分. 对这样的

36、问题可将其转换成相关的随机过程, 然后用Monte Carlo方法进行求解,Monte Carlo算法的主要组成部分,概率密度函数 必须给出描述一个物理系统的一组概率密度函数;,随机数产生器能够产生在区间0,1上均匀分布的随机数,抽样规则如何从在区间0,1上均匀分布的随机数出发,随机抽取服从给定的pdf的随机变量;,模拟结果记录记录一些感兴趣的量的模拟结果,误差估计必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；,减少方差的技术利用该技术可减少模拟过程中计算的次数；,并行和矢量化可以在先进的并行计算机上运行的有效算法,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，

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

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

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

40、R=RAN(ISEED),其他各种分布的随机数的产生,基本方法有如下三种： 逆变换法 合成法 筛选法,逆变换法,设随机变量 的分布函数为 ，定义 定理 设随机变量 服从 上的均匀分布，则 的分布函数为 。因此，要产生来自 的随机数，只要先产生来自 的随机数，然后计算 即可。其步骤为,合成法,合成法的应用最早见于Butlter 的书中。构思如下： 如果 的密度函数 难于抽样，而 关于 的条件密度函数 以及 的密度函数 均易于抽样，则 的随机数可如下产生：可以证明由此得到 的服从 。,筛选抽样,假设我们要从 抽样，如果可以将 表示成 ，其中 是一个密度函数且易于抽样，而 ， 是常数，则 的抽样可如

41、下进行： 定理 设 的密度函数 ，且 ，其中 ， ， 是一个密度函数。令 和 分别服从 和 ，则在 的条件下， 的条件密度为,标准正态分布的随机数,的随机数产生方法很多。简要介绍三种。 法1、 变换法（Box 和Muller 1958）设 ， 是独立同分布的 变量，令 则 与 独立，均服从标准正态分布。法2、 结合合成法与筛选法。（略）法3、 近似方法（利用中心极限定理）即用 个 变量产生一个 变量。其中 是抽自 的随机数， 可近似为一个 变量。,由于正态随机数的独立性，当 很小时，按(a)、(b)所构成的过程的相关时间很短，从而具有很高的截止频率。当然， 过小，样本的计算量过大。因此，选取

42、适当小即可。,关于随机数的几点注,注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,Monte Carlo方法相关引例:,1、Buffon投针实验：,18世纪，法国数学家Comte de Buffon利用投针实验估计的值,Problem of Buffon

43、s needle:,If a needle of length l is dropped at random on the middle of a horizontal surface ruled with parallel lines a distance dl apart, what is the probability that the needle will cross one of the lines?,Solution:,The positioning of the needle relative to nearby lines can be described with a ra

44、ndom vector which has components:,The random vector is uniformly distributed on the region 0,d)0,). Accordingly, it has probability density function 1/d.,The probability that the needle will cross one of the lines is given by the integral,圆周率的值 = 3.14159 26535 89793 23846 26433 83279 50288 41971 693

45、99 3751058209 74944 59230 78164 06286 20899 86280 34825 34211 7067982148 08651 32823 06647 09384 46095 50582 23172 53594 0812848111 74502 84102 70193 85211 05559 64462 29489 54930 3819644288 10975 66593 34461 28475 64823 37867 83165 27120 1909145648 56692 34603 48610 45432 66482 13393 60726 02491 41

46、27372458 70066 06315 58817 48815 20920 96282 92540 91715 3643678925 90360 01133 05305 48820 46652 13841 46951 94151 1609433057 27036 57595 91953 09218 61173 81932 61179 31051 1854807446 23799 62749 56735 18857 52724 89122 79381 83011 9491298336 73362 44065 66430 86021 39494 63952 24737 19070 2179860

47、943 70277 05392 17176 29317 67523 84674 81846 76694 0513200056 81271 45263 56082 77857 71342 75778 96091 73637 1787214684 40901 22495 34301 46549 58537 10507 92279 68925 8923542019 95611 21290 21960 86403 44181 59813 62977 47713 .,什么是蒙特-卡洛模拟？,计算机实验,物理实验 用探测器取数据,5,RHIC-STAR实验探测器,蒙特-卡洛模拟 用计算机取数据,一个例子

48、道尔顿板,理论 分布曲线实验 丢铁球模拟 计算机丢球,6,模拟的出发点 模型,根据球和钉子碰撞的规律 追踪每个球的运动根据落入每个格子的几率 给出这一几率的每次实现,8,9,令 n 为铁球数，m 为格子数。p(i) 是铁球落入第 i 个格子中的概率,再产生 n 个均匀分布的随机数：ra( j ) , j=1, 2 ,., n,将 p(i) 累加s=0.r(1)=0.do 2 i=1, ms=s+p (i)r(i+1)=senddo,do 3 k=1, mnp(k)=0do 3 j=1,nif(ra( j ).ge.r(k).and.ra( j ).lt.r(k+1) r(k) ra (j )

49、r(k+1) np(k)=np(k)+1Enddonp(k) 是落入第 k 个盒子中的粒子数。,r(m) p(i)=1r(3) p(1)+p(2)r(2) p(1)r(1) 0,i=1,m,通过比较第 j 个铁球的 ra(j) 和 r(i) 来决定这个球落入那一格子中：,例一：道尔顿板的蒙特-卡洛模拟,p3,p4,p5,p2,p1,p6,p7,例二 面积的计算,f (x),x,辛普逊方法,I = Sn,蒙特-卡洛方法,11,考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”

50、投掷N个点，若有M个点落于“图形”内，则该“图形”的面积近似为M/N。,用该方法计算的基本思路是：1根据圆面积的公式：s=R2,当R=1时，S=。由于圆的方程是：x2+y2=1(x2为x的平方的意思),因此1/4圆面积为x轴、y轴和上述方程所包围的部分。如果在1*1的正方形中均匀地落入随机点，则落入1/4圆中的点的概率就是1/4圆的面积。其4倍，就是圆面积。由于半径为1，该面积的值为的值。,例1 算圆周率,REAL X,Y,R,PIISEED=RTC()N=3000000DO 10 I=1,N X=RAN(ISEED) Y=RAN(ISEED) R=SQRT(X*X+Y*Y) IF(R.LT.