计算材料学——3.ppt

1、例题,例 求前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 *,SUM=,S END PROGRAM MAIN,程序的运行结果如下 15 SUM= 3.182290,用DO结构处理循环问题，最关键的是要处理好初值、循环终止条件和循环体语句的关系。,Finite differen

2、ce 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-,Central space,FICK 第二定律,Fick 第二定律稳态扩散解,REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE) C0=0.1 C0(0)=0.8 !BOUNDARY CONDI

3、TION C0(1001)=0.1 !BOUNDARY CONDITION C=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distance C(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 DO PRINT *, PROGRAMME IS FINISHED CLOSE(1) END,蒙特卡罗(Monte Carlo)方法,1、简介 2

4、、相关几个例子,Monte Carlo方法：,亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件产生的“频率”来近似事件的“概

5、率”。 18世纪人们用投针试验的方法来决定圆周率。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,1930年，Enrico Fermi利用Monte Carlo方法研究中子的扩散，并设计了一个Monte Carlo机械装置，Fermiac,用于计算核反应堆的临界状态,Von Neumann是Monte Carlo方法的正式奠基者,他与Stanislaw Ulam合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中， Stanislaw Ulam意识到了数字计算机的重要性,合作起源于Ma

6、nhattan工程：利用ENIAC(Electronic Numerical Integrator and Computer)计算产额,Monte Carlo模拟的应用：,自然现象的模拟：,实验探测器的模拟,数值分析：,利用Monte Carlo方法求积分,材料学中的应用：,扩散、晶粒的长大、再结晶等,Monte Carlo模拟在物理研究中的作用,Monte Carlo模拟的步骤：,根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；,从概率密度函数出发进行随机抽样，得到特征量的一些模拟结果；,对模拟结果进行分析总结，预言物理系统的某些特性。,注

7、意以下两点：,Monte Carlo方法与数值解法的不同:,Monte Carlo方法利用随机抽样的方法来求解物理问题;,数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;,Monte Carlo方法并非只能用来解决包含随机的过程的问题:,许多利用Monte Carlo方法进行求解的问题中并不包含随机过程 例如:用Monte Carlo方法计算定积分. 对这样的问题可将其转换成相关的随机过程, 然后用Monte Carlo方法进行求解,Monte Carlo算法的主要组成部分,概率密度函数 必须给出描述一个物理系统的一组概率密度函数;,随机数产生器能够产生

8、在区间0,1上均匀分布的随机数,抽样规则如何从在区间0,1上均匀分布的随机数出发,随机抽取服从给定的pdf的随机变量;,模拟结果记录记录一些感兴趣的量的模拟结果,误差估计必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；,减少方差的技术利用该技术可减少模拟过程中计算的次数；,并行和矢量化可以在先进的并行计算机上运行的有效算法,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，即随机模拟（重复“试验”）,重复试验,计算机模拟,蒙特卡罗方法的基本原理及思想,当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的

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

10、从中抽样后，即实现模拟实验后，就要确定一个随 机变量，作为所要求的问题的解，我们称它为无偏估计。,收敛性: 由大数定律, Monte-Carlo模拟的收敛是以概率而言的 误差: 用频率估计概率时误差的估计，可由中心极限定理，给定置信水平 的条件下，有： 模拟次数:由误差公式得,随机数的产生原理-计算机抽样实验,随机数的产生是实现MC计算的先决条件。而大多数概率分布的随机数的产生都是基于均匀分布U(0,1的随机数。 首先，介绍服从均匀分布U(0,1的随机数的产生方法。 其次，介绍服从其他各种分布的随机数的产生方法。以及服从正态分布的随机数的产生方法。 最后，关于随机数的几点注。,(0,1均匀分布

11、随机数的产生,最常用的是在(0,1区间内均匀分布的随机数，其他分布的随机数可利用均匀分布随机数来产生。 均匀分布随机数的产生：乘同余法 产生区间(0,1内均匀分布的随机数的递推公式是： 其中，是乘因子，M是模数。第一式称做以M为模数（modulus）的同余式，即以M 除后得到的余数记为。当给定了一个初值（称为种子）， 计算出的 即在(0,1上均匀分布的随机数。,例1：,R=RAND() ISEED=RTC() R=RAN(ISEED),其他各种分布的随机数的产生,基本方法有如下三种： 逆变换法 合成法 筛选法,逆变换法,设随机变量 的分布函数为 ，定义 定理 设随机变量 服从 上的均匀分布，则

12、 的分布函数为 。 因此，要产生来自 的随机数，只要先产生来自 的随机数，然后计算 即可。 其步骤为,合成法,合成法的应用最早见于Butlter 的书中。 构思如下： 如果 的密度函数 难于抽样，而 关于 的条件密度函数 以及 的密度函数 均易于抽样，则 的随机数可如下产生： 可以证明由此得到 的服从 。,筛选抽样,假设我们要从 抽样，如果可以将 表示成 ，其中 是一个密度函数且易于抽样，而 ， 是常数，则 的抽样可如下进行： 定理 设 的密度函数 ，且 ，其中 ， ， 是一个密度函数。令 和 分别服从 和 ，则在 的条件 下， 的条件密度为,标准正态分布的随机数,的随机数产生方法很多。简要介

13、绍三种。 法1、 变换法（Box 和Muller 1958） 设 ， 是独立同分布的 变量，令 则 与 独立，均服从标准正态分布。 法2、 结合合成法与筛选法。（略） 法3、 近似方法（利用中心极限定理） 即用 个 变量产生一个 变量。 其中 是抽自 的随机数， 可近似为一个 变量。,由于正态随机数的独立性，当 很小时，按(a)、(b)所构成的过程的相关时间很短，从而具有很高的截止频率。 当然， 过小，样本的计算量过大。因此，选取 适当小即可。,关于随机数的几点注,注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，

14、并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。 注2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,Monte Carlo方法相关引例:,1、Buffon投针实验：,18世纪，法国数学家Comte de Buffon利用投针实验估计的值,Problem of Buffons needle:,If a needle of length l is dropped at random on the middle of a horizon

15、tal 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 random vector which has components:,The random vector is uniformly distributed on t

16、he 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 69399 3751058209 74944 59230 78164 06286 20899 86280 34825 34211 7067982148 08651 32

17、823 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 4127372458 70066 06315 58817 48815 20920 96282 92540 91715 3643678925 90360 01133 0

18、5305 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 2179860943 70277 05392 17176 29317 67523 84674 81846 76694 0513200056 81271 45263 56082

19、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 实验探测器,蒙特-卡洛模拟 用计算机取数据,一个例子 道尔顿板,理论 分布曲线 实验 丢铁球 模拟 计算机丢球,6,模拟的出发点 模型,根据球和钉子碰撞的规律 追踪每个球的运动 根据落入每个格子的几率 给出这一几率

20、的每次实现,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, m s=s+p (i) r(i+1)=s enddo,do 3 k=1, m np(k)=0 do 3 j=1,n if(ra( j ).ge.r(k).and.ra( j ).lt.r(k+1) r(k) ra (j ) r(k+1) np(k)=np(k)+1 Enddo np(k) 是落入第 k 个盒子中的粒子数。,r(m) p(i)=1 r(3)

21、 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方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点，若有M个点落于“图形”内，则该“图形”的面积近似为M/N。,用该方法计算的基本思路是： 1根据圆面积的

22、公式：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,PI ISEED=RTC() N=3000000 DO 10 I=1,N X=RAN(ISEED) Y=RAN(ISEED) R=SQRT(X*X+Y*Y) IF(R.LT.1.0)N0=N0+1 10CONTINUE PI=4.0*N0/N WRITE(*,*)YUAN ZHO

23、U LV IS,PI END,例、 定积分（面积）的MC计算 随机投点法 样本平均值法,定积分的MC计算,事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。 下面考虑一个简单的定积分 为了说明问题，我们首先介绍两种求 的简单的MC方法，然后给出几种较为复杂而更有效的MC方法。,在计算积分上，MC的实用场合是计算重积分 其中 是 维空间的点，当 较大时，用MC方法比一般的数值方法有优点，主要是它的误差与维数 无关。,面积正比于点数。总点数n。点（x,y）,那么我们可以得到 的一个估计 具体试验步骤为,注1 随机投点法的思想简单明了，且每次投点结果服从二项分布，故 ，其中 注2 可证 是 的无偏估计。若用估计的标准差来衡量其精度，则估计 的精度的阶为 。 注3 这里，定积分的解，就对应我们选定的随机变量的概率值。,样本平均值法,基本原理：对积分 ，设 是 上的一个密度函数，改写 可见，任一积分均可以表示为某个随机变量 （函数）的期望。由矩法，若有n个来自a,b的观测值，则可给出 的一个矩估计。 最简单的，若n有限，可取 。,设 是来自 的随机数，则 的一个估计为 具体步骤为 注 可证 是 的无偏估计