《计算物理(研究生用)》篇

《《计算物理(研究生用)》篇》由会员分享，可在线阅读，更多相关《《计算物理(研究生用)》篇（89页珍藏版）》请在金锄头文库上搜索。

1、第 三 讲 蒙 特 卡 罗 方 法 第三篇 蒙特卡罗方法 3.1 蒙特卡罗方法概述 蒙特卡罗(Monte Carlo)方法：利用随机数统计 地去计算和模拟给定的问题。 说明： 1、也称统计模拟法、随机抽样法或统计试验法。 2、Monte Carlo方法的命名：世界上著名的赌城摩 洛哥的Monte Carlo。 3、方法：先构造一个与物理问题等价的随机过程 ，当完成大量的随机试验后，结果由多次事件的 平均值给出。 第 三 讲 蒙 特 卡 罗 方 法 5、求解确定性物理问题：如微分方程、定积分、线 性方程等。同其它数值计算方法相比是速度慢。 6、求解复杂物理问题：如果物理学问题的严格算法 不知道，

2、或非常复杂，则蒙特卡罗方法有意想不 到的成功。特别在分子运动学、输运现象、布朗 运动、放射性衰变等问题中，由于本身有一定的 统计规律性，这种方法很奏效。 7、计算机模拟：蒙特卡罗方法广泛应用于“计算机 模拟”，在计算机上模拟真实过程。 4、随机数的抽样：现在都用计算机程序产生， 一般不用物理方法抽样。计算机产生的是伪随机 数。 第 三 讲 蒙 特 卡 罗 方 法 3.2 蒙特卡罗方法的原理 用蒙特卡罗方法寻找某个未知量x时，利用计算 机产生的随机变量，得到的期望值E() x= E() 3.1 为了估算的平均值，构造随机变量的若干独立 的实验数序列i|i=1，2，3，n，得 由统计理论随机变量是

3、有限的均方差D() 3.3 第 三 讲 蒙 特 卡 罗 方 法 实际上我们不知道变量的均方差，但在计算过程 中可以根据（3.3）式进行估计 3.5 式 ( 3.4) 可清楚地看出，蒙特卡罗方法以 的速 率收敛，这是它的主要缺点。但与其它数值方法比 ，它能解决其它数值方法不能解决的问题，这正是 它有强大生命力的原因所在。 然后引入概率误差 3.4 第 三 讲 蒙 特 卡 罗 方 法 3.3 伪（赝）随机数的产生 随机数产生：有两种，一是物理方法，二是数学 方法。 物理方法：产生的随机数仅遵循统计规律，无任 何其它规律。如骰子、抽奖机等。理论上Monte Carlo方法可以利用物理方法产生随机数来

4、研究 。但实际应用不可行，运算量太大，速度上不去 。 数学方法：产生的随机数是利用计算机程序运算 得到，称它为伪（赝）随机数。它虽然不是真正 的物理随机数，但它能通过一系列的统计检验， 可以放心使用。 第 三 讲 蒙 特 卡 罗 方 法 （1）伪随机数的“伪”性：伪随机数有一定周期性 ，只不过周期很长。若在一次数值模拟中所用的 随机数都在一个周期内，模拟效果与真正的随机 数是一致的。 （2）统计检验：在使用前，应对伪随机数进行统 计检验，以确定其统计性能的优劣。主要是均匀 性和独立性。 （3）程序函数：多数程序设计语言都有此函数。 如C、Matlab等用rand()函数来实现，此函数产 生的是

5、均匀分布的随机数。 说明： 第 三 讲 蒙 特 卡 罗 方 法 其中a，c，N为给定的三个整数， (c* Xn-1 )(mod N )是指c* Xn-1被N整除后取余数。 上述关系可推出一组X1、 X2 、 X3 伪随机数 。 a，c，N选择适当，这一序列伪随机数具有足够长 的周期性。 (4) 伪随机数的产生方法：主要是乘同余法，即下面 的伪随机数递推关系： Xn=(c* Xn-1 )(mod N ) ， X0 =a 第 三 讲 蒙 特 卡 罗 方 法 例如：产生032767之间均匀分布的随机数递推公式 Xn =(889* Xn -1 )(mod32768 ) ， X0 =13 若产生01之间

6、均匀分布的随机数，上式改为 Xn =(889* Xn -1 )(mod32768 ) /32767 ， X0 =13 这个递推公式产生的伪随机数周期很长，是比 较实用的一个产生伪随机数的公式。 第 三 讲 蒙 特 卡 罗 方 法 3.4 伪（赝）随机变量的抽样 实际上，大多数伪随机变量不满足0，1均匀分 布，而是具有符合分布密度函数为f(x)的分布。 抽取符合f(x) 随机数的步骤如下： （2）从上面随机数总体中抽取一个简单子样，使它 满足分布密度函数f(x) ： 1，2，3 （1）在0，1之间抽取均匀分布的伪随机数序列 ： 1，2，3 第 三 讲 蒙 特 卡 罗 方 法 一、离散型分布随机变

7、量的直接抽样法 随机变量的概率表 其中xi为离散型随机变量的跳跃点，Pi为概率。 q变换方法： 第一步 构造一个矢量 第 三 讲 蒙 特 卡 罗 方 法 第二步 产生一个0，1间随机数，找到区间 使恰好落在这个区间内。 第三步 P取对应的值xj，其表达式为 第 三 讲 蒙 特 卡 罗 方 法 例1.光子与物质相互作用的抽样问题。 物理过程：光子与物质作用有康普顿效应、光 电效应和电子对效应三种类型 。其中光电效应和 电子对效应为光子吸收过程。 设三种过程的碰撞截面分别s 、e和p，则总截 面T = s +e+p 。 根据给定的0，1之间均匀分布的随机数 ，求应 产生那种效应？ 第 三 讲 蒙

8、特 卡 罗 方 法 （3）产生0,1间均匀随机数： 若 s /T，则发生康普顿散射； 若s/T(s+e)/T ， 则发生光电过程； 若(s + e ) /T ，则发生电子对产生过程。 （2）构造一个矢量 解：（1）随机变量的概率表 第 三 讲 蒙 特 卡 罗 方 法 概率表 例2 Possion分布： Possion分布是离散型分布，其概率函数为 求符合此分布的随机变量 第 三 讲 蒙 特 卡 罗 方 法 二、连续型分布随机变量的直接抽样法 一个对应一个。如果积分式能求出解析形式的 反函数，则得到变换公式 F=F-1() f(x) x a 1- b 对连续型随机变量的值，可直接解连续分布方程

9、其中f(x)为在（a，b）区间的归一化分布密度函数。 第 三 讲 蒙 特 卡 罗 方 法 连续型分布的困难性：在实际问题中，连续型分 布是很复杂的，有的只能给出分布函数F解析表 达式，但给不出其反函数F-1() ，如分布。有 的则连分布函数都给不出解析式，如正态分布。 由于这一情况，对于连续型分布采用求反函数的 方法就及其困难，这必须寻找其它方法。 第 三 讲 蒙 特 卡 罗 方 法 例1.均匀分布： 在区间（a，b）内找出符合均匀分布的随机变量 。 解法2：线性变换，令 解法1： 第 三 讲 蒙 特 卡 罗 方 法 例2指数分布： 指数分布的一般形式为 f(x)=e-x x 0 其中0 ，求

10、符合此分布的随机变量 求反函数 =-1/ln(1-F() 其中0F 1，由上式得 =-1/ln(1-) 因为1-0,1 和0,1是等价，所以其式等价于 =-1/ln 第 三 讲 蒙 特 卡 罗 方 法 三、连续型分布随机变量的变换抽样法 变换抽样法：变换积分形式得到反函数的抽样法。 变换抽样的基本思想： （1）寻找到一个直接抽样的(y) ，记(y) 。 （2）再寻找一个适当的变换关系x=g(y)。且g(y)反函 数存在，一阶导数连续，记g-1(x)=h(x)。 （3）根据概率论，x满足密度函数h(x)*|h(x)| 。 （4）如果选择合适的g(y) ，使之 f(x)= h(x)*|h(x)|

11、则通过对(y)的抽样 ，再变换=g()得到满足f(x)的 抽样值。 第 三 讲 蒙 特 卡 罗 方 法 例1正态分布(高斯分布)： 正态分布的形式为 求符合此分布的随机变量。 求积分 第 三 讲 蒙 特 卡 罗 方 法 正态分布抽样的方法有许多种。正态分布在统计 物理学中是最重要的分布之一，也是可用最多的方 法来产生随机数的分布之一。 取其中一个即可 另外是在0，2之间的数， 必须=22 可用替代 。 第 三 讲 蒙 特 卡 罗 方 法 乘分布舍选抽样法的基本思想：若随机变量 满足密度函数f(x)，直接抽样不方便或不可求出 时，可设法把f(x)写成 f(x)=H(x) f1(x) 四、连续型分

12、布的乘分布舍选抽样法(第二类) 其中H(x), f1(x)是归一化密度函数 。 A） H(x) 有上界M，即H(x) M B） f1(x)可直接抽样h=F 1() C）对f(x)的乘分布舍选抽样过程 如左图。 其中f是符合f(x)的随机变量 取 第 三 讲 蒙 特 卡 罗 方 法 例1Maxwell分布： Maxwell分布的概率密度为 当 时，即为统计物理学的速率分布。 求符合此分布的随机变量 解：用乘分布舍选抽样法。把f(x)写成 其中 第 三 讲 蒙 特 卡 罗 方 法 取 B）f1(x)可直接抽样h=F1() C）对f(x)的抽样过程是挑选h，使 h符合f(x)的随机分布。 A）H(x

13、)是有上界M，即H(x)M。由H(x)取极值 第 三 讲 蒙 特 卡 罗 方 法 3.5 Monte Carlo方法求圆周率 求的物理模型如图所示，边长 为1的正方形内切一个圆。于是 正方形和圆的面积分别为 S正=1，S圆= /4 T正：计算机掷出的 随机数总数； T圆：随机点(x,y)落 在圆内的总数.即， 当 时的随机点总数 。 y 0.5 0.5 x 若正方形内有T正个随机点(x,y) x=1-0.5 ； y=2-0.5 其中为0，1之间的随机数。 圆内的随机点数为T圆，那么有 S正/S圆=T正/T圆 得： =4T正/T圆 第 三 讲 蒙 特 卡 罗 方 法 K=10 , M=1000

14、, P=0,i=1 j=1 j=j+1 P=P+1 PI=4P/(i*(j-1) , printf(PI) i=i+1 开始 结束 R2=0.25 i=K j=M 第 三 讲 蒙 特 卡 罗 方 法 %用Monte Carlo方法对圆周率进行计算的Matlab程序，文件名 PI_val.m% K=100;M=5000;p=0; for i=1:K for j=1:M X=rand(1,2); R2=(X(1,1)-0.5)2+(X(1,2)-0.5)2; if R20.25|R2=0.25 F1F0 （4）如果连续产生m个随机数还不满足|F1|F0|， 则将b减半再进行第三步运算。 （5）重复

15、上述过程，直到|F0|0，在求根区域内，模拟给 出方程近似解（x1*, x2*, , xn*）使之满足 （x1*, x2*, , xn*）0，在-b,b之间，反复产生均匀分布的 随机数组(1,2,n)，计算 F1=（x1*+1, x2*+2, , xn*+n） (3)直到发现一组(1,2,n)使得F1F0为止，然后令 xi*+ixi* i=1,2,n; F1F0 (4)若连续产生m组随机数(1,2,n)，还不满足F1F0 ，则将b减半再重复2、3过程，直到F0 1表示裂变数增加，发生链式反应； f 1表示反应逐渐停止； f =1表示处于临界状态，用临界质量Mc表示。 第 三 讲 蒙 特 卡 罗 方 法 用计算机模拟具有一定大小及形状的体积内发生 的大量随机裂变，然后计算放出中子被吸收引起 的裂变数Nin，求出相应的f 值。 为简单起见，考虑 铀的几何形状为长 方形a*a*b，如图所 示。 蒙特卡罗法研究这 种形状的核材料临 界问题： 第 三 讲 蒙 特 卡 罗 方 法 （1）核裂变的随机位置