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

计算材料学概述计算材料学概述第三章第三章蒙特卡罗方法蒙特卡罗方法（Monte Carlo）主要内容主要内容o Monte Carlo模拟发展简介模拟发展简介o Monte Carlo模拟基本原理模拟基本原理o Monte Carlo模拟典型算法模拟典型算法o Monte Carlo模拟典型应用模拟典型应用 蒙特卡洛法是什么？蒙特卡洛法是什么？o 蒙特卡洛蒙特卡洛(Monte Carlo)方法，方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。o 蒙特蒙特卡罗方法卡罗方法（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世纪，法国数学家世纪，法国数学家Comte de Buffon利用投针实验估计利用投针实验估计 的值的值dLp2dL1777年法国科学家布丰提出的一种计算圆周率的方法年法国科学家布丰提出的一种计算圆周率的方法随机投随机投针法，即著名的针法，即著名的布丰投针布丰投针问题。这一方法的步骤是：问题。这一方法的步骤是： 1) 取一张白纸，在上面画上许多条间距为取一张白纸，在上面画上许多条间距为d的平行线。的平行线。 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.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子程序的周期约230。实例一、计算值o 计算过程：计算过程：o 1 1、构造或描述问题的概率过程、构造或描述问题的概率过程o 2 2、从概率密度函数出发进行随机抽样，实、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一现从已知概率分布的抽样，得到特征量的一些模拟结果些模拟结果计算均值计算均值Monte Carlo方法之方法之典型算法与应用o考虑平面上的一个边长为考虑平面上的一个边长为1 1的正方形及其内部的一个形状不规的正方形及其内部的一个形状不规则的则的“图形图形”，如何求出这个，如何求出这个“图形图形”的面积呢？的面积呢？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，该面积的值为的值。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(*,*)PIENDo面积的计算面积的计算f (x)x辛普逊方法辛普逊方法I = Sn蒙特蒙特-卡洛方法卡洛方法f (x)x在长方形中均匀投在长方形中均匀投N0组组(x,y)如如 yf(x), 则则 N=N+1I =(N/N0)S0SS011MC的优点oMCMC与传统数学方法相比，具有与传统数学方法相比，具有直观性强，简便易行的优点，直观性强，简便易行的优点，该方该方法能处理一些其他方法无法解决的负责问题，并且容易在计算机法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以上实现，在很大程度上可以代替许多大型、难以实现的复杂实验代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。和社会行为。无污染、无危险、能摆脱实验误差。Monte Carlo方法利用方法利用随机抽样随机抽样的方法来求解物理问题的方法来求解物理问题;数值解法数值解法:从一个物理系统的数学模型出发从一个物理系统的数学模型出发,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态;Monte Carlo方法方法并非只能用来解决包含随机过程并非只能用来解决包含随机过程的问题的问题: 例例如如:用用Monte Carlo方法计算定积分方法计算定积分. 对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程, 然后用然后用Monte Carlo方法进行求解方法进行求解注意以下两点：注意以下两点：MC的应用自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；数值分析：数值分析：数学问题，求积分，求逆矩阵，解线性代数等数学问题，求积分，求逆矩阵，解线性代数等经济学模拟：经济学模拟：库存问题，随机服务系统中排队问题库存问题，随机服务系统中排队问题 人口问题：人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争人口的出生，传染病的蔓延；乃至动物的生态竞争金属学：扩散、组织长大、相金属学：扩散、组织长大、相变过程变过程蒙特-卡洛模拟的意义o能研究不同边界、不同材料的影响能研究不同边界、不同材料的影响 理论不可能、实验耗费太大理论不可能、实验耗费太大 用于实验设计用于实验设计o无污染无污染 反应堆防护反应堆防护 核弹爆炸核弹爆炸o能摆脱实验误差能摆脱实验误差 作理论和实验的桥梁作理论和实验的桥梁Monte Carlo模拟的步骤：模拟的步骤：1. 根据欲研究的物理系统的性质，建立能够描述该根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量系统特性的理论模型，导出该模型的某些特征量的概率密度函数；的概率密度函数； （即构造或描述问题的概率过程）（即构造或描述问题的概率过程）2. 从概率密度函数出发进行随机抽样，实现从已知从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；概率分布的抽样，得到特征量的一些模拟结果； 有了明确的概率过程后，为了实现过程的数有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的抽值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。量的大量试验值。产生已知概率分布的随机变量，产生已知概率分布的随机变量，是实现是实现MC方法的关键步骤方法的关键步骤，其中最基本的是（，其中最基本的是（0，1）均匀分布。）均匀分布。3. 对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。4. 模拟结果的检验模拟结果的检验Monte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数(pdf) 必须给出描述一个物理系统的一组概必须给出描述一个物理系统的一组概率密度函数率密度函数;随机数产生器随机数产生器能够产生在区间能够产生在区间0,1上（均匀）分布的随机数上（均匀）分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,1上均匀分布的随机数出发上均匀分布的随机数出发,随随机抽取服从给定的机抽取服从给定的pdf的随机变量的随机变量;模拟结果记录模拟结果记录记录一些感兴趣的量的模拟结果记录一些感兴趣的量的模拟结果误差估计误差估计必须确定统计误差（或方差）随模拟次数以及其必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；它一些量的变化；减少方差的技术减少方差的技术利用该技术可减少模拟过程中计算的次数；利用该技术可减少模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法实例二 定积分计算o 事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。o 下面考虑一个简单的定积分 dxxfbao ! 计算计算x*2在在(0,1)上积分上积分o 计算过程：计算过程：o 1、构造或描述问题的概率过程、构造或描述问题的概率过程:产生服从分产生服从分布布f（x）的随机变量）的随机变量Xi( )（i=1,2, ,N）o 2、从概率密度函数出发进行随机抽样，实现、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些从已知概率分布的抽样，得到特征量的一些模拟结果模拟结果计算均值计算均值( )niix12nixiy12REAL YY=0N=300000ISEED=RTC()DO I=1,N X=RAN(ISEED) Y=Y+X*2/N END DOWRITE(*,*)YEND limx2dx (dx 0)Monte Carlo方法另一个重要问题：方法另一个重要问题：随机数随机数 随机数：由单位矩阵分布中所产生的简单子样称为随机数序列，其中的每一个个体称为随机数。 但真正的随机数的不适合电子计算机上使用，因为它需要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，且其产生的序列无法重复实现，使程序无法复算，结果无法验证，同时需要增添随机数发生器和电路联系等附加设备。伪随机数：o 是有数学递推公式所产生的随机数。（近似的具备随机数的性质。） An+1=T(A);An+1= An+k+1o 伪随机的优点和缺点：o 判断伪随机数好坏的方法： 1、它能够有较好的均匀性和独立性； 2、它的费用大小，即指所消耗计算机的时间； 3、容量要求尽可能大。随机数产生的办法o产生均匀分布随机数的几种方法； （1）物理方法；（2）数学方法。o伪随机数产生方法：加同余法乘同余法乘加同余法取中方法逆变换法合成法筛选法 。关于随机数的几点注意关于随机数的几点注意o 注注1 1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。o 注注2 2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。 FORTRAN 语言产生随机数的实例orandom_number(xrandom_number(x) )产生一个产生一个0 0到到1 1之间的随机数（之间的随机数（x x可以是向量），可以是向量），但是每次总是那几个数。但是每次总是那几个数。用了用了random_seedrandom_seed()()后，系统根据日期和时间随机地提供种子，后，系统根据日期和时间随机地提供种子，使得随机数更随机了。使得随机数更随机了。programrandomprogramrandomreal:xreal:xcallcallrandom_seedrandom_seed()!()!系统根据日期和时间随机地提供种子系统根据日期和时间随机地提供种子 callcallrandom_numberrandom_number(x)!(x)!每次的随机数就都不一样了每次的随机数就都不一样了 write(write(* *, ,* *)x)xstopstopendprogramrandom endprogramrandom 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子程序的周期约230。Finite difference approximation of differential equationso A differential equation can be approximated by a finite difference scheme. For exampletuudxdutttForward timexuudxdutxxtx2222xuuudxudtxxtxtxxBackward spaceCentral space21222 xuuutuuxututxxtxtxxtxtxForward t aime-Central spaceFICK 第二定律21222 xuuutuuxututxxtxtxxtxtxvxCCCDtCCxCDttxxtxtxxtxtx21222 C 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.8!C(1001)=0.1 IF(JT=50000) WRITE(1,*) IX,C(IX) END DO C0=C END DOEND21222 C xCCCDtCCxCDttxxtxtxxtxtx应用之二 生日问题MC模拟模拟假设有假设有n个人在一起，各自的生日为个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直天之一，根据概率理论，与很多人的直觉相反，只需觉相反，只需23个人便有大于个人便有大于50的几的几率人群中至少有率人群中至少有2个人生日相同。个人生日相同。365)1(365.3653633653643653651nn 理论几率 模拟几率100.117 0.110200.411 0.412230.527 0.520300.706 0.692450.941 0.93650 0.986 0.987INTEGER 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,10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END MC在材料学领域的应用 随机行走背景背景如，布朗运动最简单、无限制随机如，布朗运动最简单、无限制随机行走行走( (Unrestricted Unrestricted randonrandon walk,RWwalk,RW ) )startendtR平均平方端平均平方端-端位移端位移: vttR2221v,自然科学和社会生活中很多现象都与随机运动有关自然科学和社会生活中很多现象都与随机运动有关可以模拟的内容？扩散；分子运动；。o 如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了所有分子距原点的方均距离：！Monte Carlo Simulation of One Dimensional Diffusion INTEGER 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 DOOPEN(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 WRITE(*,*) 实验天数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 DOEND 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)ENDo!Monte Carlo Simulation of Two Dimensional Diffusion ooINTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNoREAL XYM(1:1000)o!X:INSTANTANEOUS POSITION OF ATOMo!XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃o!XYM(I): THE MEAN OF XY oWRITE(*,*) 实验天数JMAX，实验次数 IMAXoREAD(*,*) JMAX,IMAXoXN=(/0,0,-1,1/)oYN=(/-1,1,0,0/)oISEED=RTC()oDO J=1,JMAX !第几天实验oX=0 !oY=0 !oDO I=1,IMAX !第几步跳跃oRN=4*RAN(ISEED)+1oX=X+XN(RN)o Y=Y+YN(RN)oXY(J,I)=X*X+Y*YoEND DOoEND DOoOPEN(1,FILE=C:DIF2.DAT)oDO I=1,IMAXoXYM=0.0oXYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !oWRITE(1,*) I, XYM(I)oEND DOoCLOSE(1) ！做三维空间随机行走？！做三维空间随机行走？oEND实际环境是复杂的实际环境是复杂的:可变，缺陷，风，可变，缺陷，风， 季节，电场等。季节，电场等。 空间中有随机性缺陷空间中有随机性缺陷不退行走不退行走自回避行走自回避行走o 例如模拟高分子的位形 用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形，行走过的位置代表的是构成分子的原子或官能团，因此，无限制随机行走忽略了体斥效应。 不退行走就是禁止在每一步行走后立即倒退，可以解决刚走的一步与上一步重叠的问题。但不退行走没有完全解决高分子的体斥效应。 自回避行走就是所有已走过的位置不能再走，这样就完全解决了体斥效应问题。o 气体分子的随机行走 假设在一个空旷封闭的房间中心滴上一滴香水，挥发的香水分子随机的与空气中的粒子发生碰撞，最终会扩散到整个房间。这里我们将通过计算机模拟400个分子在二维平面内的随机运动 讨论熵变o 现在来讨论一下熵在我们这个气体扩散模型中是什么意思。在刚开始的时候，我们的系统处在一个非常有序的状态：所有的分子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随着时间的推移，分子开始不断的向外扩散，这时系统出现了无序状态，熵开始逐渐增加。 o 系统的熵和我们预测的一样在随时间增大，但当所有分子布满整个边界以内区域时，系统的熵开始趋于稳定。从这一结果中我们可以了解：当所有的分子随机的布满整个区域时，虽然当我们跟踪某一个确定分子时，它还是在区域内到处乱窜，但每个小区域内分子的密度却不会再变化了。所以，一旦气体分子扩散到整个区域以后，不管我们再等上多少时间，系统的熵都不会再有太大的起伏。换句话说，让系统自动回到开始的状态，即所有分子都在原点的状态，已经不可能了。 画一个圆晶粒ooUSE MSFLIB o!INTEGER XR,YR !在的区域中画一个圆oPARAMETER XR=400,YR=400 oINTEGER R,S(1:XR,1:YR) oX0=XR/2 ! 圆心位置X0，YOoY0=YR/2oR=MIN(X0-10,Y0-10) !圆半径oS=0 !像素的初始状态（颜色）oDO I=1,XRoDO J=1,YRoIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10oIER=SETCOLOR(S(I,J)+3)o IER=SETPIXEL(I,J)oEND DOoEND DOoEND画晶界o!画一个圆oUSE MSFLIB oINTEGER XR，YR !在的区域中画一个圆oPARAMETER XR=400,YR=400 oINTEGER R,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNS oXN=(/0,0,-1,1/)oYN=(/-1,1,0,0/)oX0=XR/2 ! 圆心位置X0，YOoY0=YR/2oR=MIN(X0-10,Y0-10) !圆半径oS=0 !像素的初始状态（颜色）oDO I=1,XRoDO J=1,YRoIF(I-X0)*2+(J-Y0)*20)THENoIER=SETCOLOR(9)oELSEoIER=SETCOLOR(8)oEND IFo IER=SETPIXEL(I,J)oEND DOoEND DOoENDo 如何画有一定宽度的晶界？例题例 求前N个自然数倒数和。即求ns1.31211 程序需要几个变量：累加的项数I(整型)，存放通项值的T(实型)，存放总和的S(实型)，待输入的总项数N(整型)。PROGRAM MAINPROGRAM MAIN IMPLICIT NONE IMPLICIT NONE INTEGER N,I INTEGER N,I REAL T,S REAL T,S READ READ * *,N,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 PRINT * *,SUM=,S,SUM=,SEND PROGRAM MAINEND PROGRAM MAIN 程序的运行结果如下15 SUM= 3.182290 用DO结构处理循环问题，最关键的是要处理好初值、循环终止条件和循环体语句的关系。Finite difference approximation of differential equationso A differential equation can be approximated by a finite difference scheme. For exampletuudtdutttForward timexuudxdutxxtx2222xuuudxudtxxtxtxxBackward spaceCentral space21222 xuuutuuxututxxtxtxxtxtxForward time-Central spaceFICK 第二定律21222 xuuutuuxututxxtxtxxtxtx21222 C xCCCDtCCxCDttxxtxtxxtxtxFick 第二定律稳态扩散解 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,*) IX,C(IX) END DO C0=C END DOPRINT *, PROGRAMME IS FINISHED CLOSE(1) END21222 C xCCCDtCCxCDttxxtxtxxtxtx蒙特卡罗蒙特卡罗(Monte Carlo)方法方法1、简介2、相关几个例子Monte Carlo方法：方法：亦称统计模拟方法，亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte Carlo名字的由来：名字的由来： 是由是由Metropolis在二次世界大战期间提出的：在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；计划，研究与原子弹有关的中子输运过程； Monte Carlo是摩纳哥（是摩纳哥（monaco)的首都，该城以赌博闻名的首都，该城以赌博闻名Nicholas Metropolis (1915-1999)Monte-Carlo, Monacoo Monte Carlo方法的基本思想很早以前就被人们方法的基本思想很早以前就被人们所发现和利用。早在所发现和利用。早在17世纪，人们就知道用事件产世纪，人们就知道用事件产生的生的“频率频率”来近似事件的来近似事件的“概概 率率”。 18世纪人世纪人们用投针试验的方法来决定圆周率们用投针试验的方法来决定圆周率。本世纪本世纪40年年代电子计算机的出现，特别是近年来高速电子计算代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。地模拟这样的试验成为可能。1930年，年，Enrico Fermi利用利用Monte Carlo方法研究中子的扩散，方法研究中子的扩散，并设计了一个并设计了一个Monte Carlo机械装置，机械装置，Fermiac,用于计算核用于计算核反应堆的临界状态反应堆的临界状态Von Neumann是是Monte Carlo方法的正式奠基者方法的正式奠基者,他与他与Stanislaw Ulam合作建立了概率密度函数、反累积分布函数的数学基础，合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中，以及伪随机数产生器。在这些工作中， Stanislaw Ulam意识意识到了数字计算机的重要性到了数字计算机的重要性合作起源于合作起源于Manhattan工程：利用工程：利用ENIAC(Electronic Numerical Integrator and Computer)计算产额计算产额Monte Carlo模拟的应用：模拟的应用：自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；实验探测器的模拟实验探测器的模拟数值分析：数值分析：利用利用Monte Carlo方法求积分方法求积分材料学中的应用：材料学中的应用：扩散、晶粒的长大、再结晶等扩散、晶粒的长大、再结晶等Monte Carlo模拟在物理研究中的作用模拟在物理研究中的作用Monte Carlo模拟的步骤：模拟的步骤：1. 根据欲研究的物理系统的性质，建立能够描述该系统特性根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；的理论模型，导出该模型的某些特征量的概率密度函数；2. 从概率密度函数出发进行随机抽样，得到特征量的一些模从概率密度函数出发进行随机抽样，得到特征量的一些模拟结果；拟结果；3. 对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。注意以下两点：注意以下两点： Monte Carlo方法与数值解法的不同方法与数值解法的不同:Monte Carlo方法利用随机抽样的方法来求解物理问题方法利用随机抽样的方法来求解物理问题;数值解法数值解法:从一个物理系统的数学模型出发从一个物理系统的数学模型出发,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态; Monte Carlo方法并非只能用来解决包含随机的过程的问题方法并非只能用来解决包含随机的过程的问题:许多利用许多利用Monte Carlo方法进行求解的问题中并不包含随方法进行求解的问题中并不包含随机过程机过程 例如例如:用用Monte Carlo方法计算定积分方法计算定积分. 对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程, 然后用然后用Monte Carlo方法进行求解方法进行求解Monte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数 必须给出描述一个物理系统的一组概率密必须给出描述一个物理系统的一组概率密度函数度函数;随机数产生器随机数产生器能够产生在区间能够产生在区间0,1上均匀分布的随机数上均匀分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,1上均匀分布的随机数出发上均匀分布的随机数出发,随随机抽取服从给定的机抽取服从给定的pdf的随机变量的随机变量;模拟结果记录模拟结果记录记录一些感兴趣的量的模拟结果记录一些感兴趣的量的模拟结果误差估计误差估计必须确定统计误差（或方差）随模拟次数以及其必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；它一些量的变化；减少方差的技术减少方差的技术利用该技术可减少模拟过程中计算的次数；利用该技术可减少模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法确定性系统确定性系统随机性系统随机性系统模拟模拟自然界自然界Monte-Carlo模拟，即随机模拟（重复模拟，即随机模拟（重复“试验试验”）重复试验重复试验计算机模拟计算机模拟蒙特卡罗方法的基本原理及思想o 当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。o 把蒙特卡罗解题归结为三个主要步骤： 构造或描述概率过程； 实现从已知概率分布抽样； 建立各种估计量。 o对于本身就具有随机性质的问题，主要是正确描述和模拟这个概率过程；对于不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。 o实现从已知概率分布抽样(模拟试验)：构造了概率模型以后，产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。o建立各种估计量： 构造了概率模型并能从中抽样后，即实现模拟实验后，就要确定一个随 机变量，作为所要求的问题的解，我们称它为无偏