1、合肥工业大学电子科学与应用物理学院第十一章第十一章 蒙特卡罗方法蒙特卡罗方法11.1 蒙特卡罗方法概述蒙特卡罗方法概述 蒙特卡罗蒙特卡罗(Monte Carlo)方法方法：利用随机数统计地：利用随机数统计地去计算和模拟给定的问题。去计算和模拟给定的问题。q说明：说明：1、也称也称统计模拟法、随机抽样法或统计试验法统计模拟法、随机抽样法或统计试验法。2、Monte Carlo方法的命名方法的命名：世界上著名的赌城摩世界上著名的赌城摩洛哥的洛哥的Monte Carlo。3、方法：方法：先构造一个与物理问题等价的随机过程，先构造一个与物理问题等价的随机过程，当完成大量的随机试验后，结果由多次事件的平
2、当完成大量的随机试验后，结果由多次事件的平均值给出。均值给出。合肥工业大学电子科学与应用物理学院5、求解确定性物理问题求解确定性物理问题：如微分方程、定积分、线：如微分方程、定积分、线性方程等。同其它数值计算方法相比是速度慢。性方程等。同其它数值计算方法相比是速度慢。6、求解复杂物理问题：求解复杂物理问题：如果物理学问题的严格算法如果物理学问题的严格算法不知道，或非常复杂，则蒙特卡罗方法有意想不不知道，或非常复杂，则蒙特卡罗方法有意想不到的成功。特别在分子运动学、输运现象、布朗到的成功。特别在分子运动学、输运现象、布朗运动、放射性衰变等问题中，由于本身有一定的运动、放射性衰变等问题中，由于本身
3、有一定的统计规律性，这种方法很奏效。统计规律性，这种方法很奏效。7、计算机模拟：计算机模拟：蒙特卡罗方法广泛应用于蒙特卡罗方法广泛应用于“计算机计算机模拟模拟”，在计算机上模拟真实过程。，在计算机上模拟真实过程。4、随机数的抽样：随机数的抽样：现都用计算机程序来产生，一般现都用计算机程序来产生，一般不用物理方法抽样。计算机产生的是伪随机数。不用物理方法抽样。计算机产生的是伪随机数。合肥工业大学电子科学与应用物理学院11.2 蒙特卡罗方法的原理蒙特卡罗方法的原理用蒙特卡罗方法寻找某个未知量用蒙特卡罗方法寻找某个未知量x时，利用计算机时，利用计算机产生的随机变量产生的随机变量，得到的期望值，得到的
4、期望值E()x=E()为了估算为了估算的平均值，构造随机变量的平均值，构造随机变量的若干独立的若干独立的实验数序列的实验数序列i|i=1，2，3，n，得，得 niinx11合肥工业大学电子科学与应用物理学院 11.3 伪（赝）随机变量的抽样伪（赝）随机变量的抽样实际上，大多数伪随机变量不满足实际上，大多数伪随机变量不满足0，1均匀分布，均匀分布，而是具有符合分布密度函数为而是具有符合分布密度函数为f(x)的分布。的分布。抽取符合抽取符合f(x)随机数的步骤如下：随机数的步骤如下：（2）从上面随机数总体中抽取一个简单子样，使它）从上面随机数总体中抽取一个简单子样，使它满足分布密度函数满足分布密度
5、函数f(x)：1，2，3（1）在）在0，1之间抽取均匀分布的伪随机数序列之间抽取均匀分布的伪随机数序列 ：1，2，3 合肥工业大学电子科学与应用物理学院一、离散型分布随机变量的直接抽样法一、离散型分布随机变量的直接抽样法nnPxPxPxPx,.,.,332211随机变量随机变量的概率表的概率表 其中其中xi为离散型随机变量为离散型随机变量的跳跃点的跳跃点,Pi为相应概率。为相应概率。q变换方法变换方法：第一步第一步 构造一个矢量构造一个矢量 ).,.,(321321211nPPPPPPPPPPP合肥工业大学电子科学与应用物理学院.,.3211321jjPPPPPPPP第二步第二步 产生一个产生
6、一个0，1间随机数间随机数，找到区间，找到区间使使恰好落在这个区间内。恰好落在这个区间内。j1i1-j1i iijPPPx当第三步第三步 P取对应的值取对应的值xj，其表达式为，其表达式为 合肥工业大学电子科学与应用物理学院例例1.光子与物质相互作用的抽样问题。光子与物质相互作用的抽样问题。物理过程物理过程：光子与物质作用有康普顿效应、光电光子与物质作用有康普顿效应、光电效应和电子对效应三种类型效应和电子对效应三种类型。其中光电效应和电。其中光电效应和电子对效应为光子吸收过程。子对效应为光子吸收过程。设三种过程的碰撞截面分别设三种过程的碰撞截面分别s、e和和p，则总截面，则总截面T=s+e+p
7、。根据给定的根据给定的0，1之间均匀分布的随机数之间均匀分布的随机数，求应，求应产生那种效应？产生那种效应？合肥工业大学电子科学与应用物理学院 （3）产生）产生0,1间均匀随机数间均匀随机数：若若 s/T，则发生康普顿散射；，则发生康普顿散射；若若s/T(s+e)/T，则发生光电过程；则发生光电过程；若若(s+e)/T，则发生电子对产生过程。，则发生电子对产生过程。（2）构造一个矢量构造一个矢量),(),(321211TpesTesTsPPPPPPP解：（解：（1）随机变量的概率表随机变量的概率表TpTeTspes合肥工业大学电子科学与应用物理学院二、连续型分布随机变量的直接抽样法二、连续型分
8、布随机变量的直接抽样法一个一个对应一个对应一个。如果积分式能求出解析形式的。如果积分式能求出解析形式的反函数，则得到变换公式反函数，则得到变换公式 F F=F F-1-1()11.12)11.12f(x)xa1-b对连续型随机变量对连续型随机变量的值，可直接解连续分布方程的值，可直接解连续分布方程11.11 )()(FaFdxxf其中其中f(x)为在（为在（a，b）区间的归一化分布密度函数。）区间的归一化分布密度函数。合肥工业大学电子科学与应用物理学院例例1.均匀分布：均匀分布：在区间（在区间（a，b）内找出符合均匀分布的随机变量）内找出符合均匀分布的随机变量。其它 0bxa 1)(abxf解
9、法解法2：线性变换，令线性变换，令 aabaAabAbaAA)(;|101010得由aabFabadxabdxxfFaa)()(1)(1)(解法解法1：合肥工业大学电子科学与应用物理学院例例2指数分布：指数分布：指数分布的一般形式为指数分布的一般形式为 f(x)=e-x x 0其中其中0，求符合此分布的随机变量，求符合此分布的随机变量求反函数求反函数 =-1/ln(1-F()其中其中0F 1，由上式得，由上式得 =-1/ln(1-)1|)()()(：00eedxxfdxxfFx解因为因为1-0,1 和和0,1是等价，所以其式等价于是等价，所以其式等价于 =-1/ln合肥工业大学电子科学与应用物
10、理学院 11.4 Monte Carlo方法求圆周率方法求圆周率求求的物理模型如图所示，边长的物理模型如图所示，边长为为1的正方形内切一个圆。于是的正方形内切一个圆。于是正方形和圆的面积分别为正方形和圆的面积分别为 S正正=1，S圆圆=/45.0;22RyxRT正正：计算机掷出的：计算机掷出的随机数总数；随机数总数；T圆圆：随机点：随机点(x,y)落落在圆内的总数。即，在圆内的总数。即，当当 时的随机点总数。时的随机点总数。y0.50.5x若正方形内有若正方形内有T正正个随机点个随机点(x,y)x=1-0.5 ；y=2-0.5其中其中为为0，1之间的随机数。之间的随机数。圆内的随机点数为圆内的
11、随机点数为T圆圆，那么有，那么有 S正正/S圆圆=T正正/T圆圆 得得:=4T正正/T圆圆合肥工业大学电子科学与应用物理学院K=10,M=1000,P=0,i=1j=1j=j+1P=P+1PI=4P/(i*(j-1),printf(PI)i=i+122221,5.0,5.0yxRyx开始开始结束结束R R2 2=0.25=0.25i=Kj=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)
12、;R2=(X(1,1)-0.5)2+(X(1,2)-0.5)2;if R2 1表示裂变数增加，发生链式反应；表示裂变数增加，发生链式反应；f 1表示反应逐渐停止；表示反应逐渐停止；f=1表示处于临界状态，用临界质量表示处于临界状态，用临界质量Mc表示。表示。合肥工业大学电子科学与应用物理学院用计算机模拟具有一定大小及形状的体积内发生的用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然后计算放出中子被吸收引起的裂大量随机裂变，然后计算放出中子被吸收引起的裂变数变数Nin，求出相应的，求出相应的f 值。值。为简单起见，考虑铀为简单起见，考虑铀的的几何几何形状为长方形形状为长方形a*a*b
13、，如图所示。，如图所示。蒙特卡罗法研究这种蒙特卡罗法研究这种形状的核材料临界问形状的核材料临界问题：题：ZYX(X0Y0Z0)(X1Y1Z1)Oaab合肥工业大学电子科学与应用物理学院（1）核裂变的随机位置核裂变的随机位置（X0，Y0，Z0）随机点的坐标范围为随机点的坐标范围为-0.5a X0 0.5a；-0.5a Y0 0.5a；-0.5b Z0 0.5b（2）裂变的两个中子出射方向裂变的两个中子出射方向(,)裂变产生的中子出射方向用裂变产生的中子出射方向用和和来描述。从（来描述。从（X0，Y0，Z0）点放出的中子与以此为圆心的单位球上某一）点放出的中子与以此为圆心的单位球上某一面积发生碰撞
14、，它的碰撞几率只与球面上被碰撞的面面积发生碰撞，它的碰撞几率只与球面上被碰撞的面积大小成正比。或者说积大小成正比。或者说碰撞按立体角均匀分布碰撞按立体角均匀分布。一次裂变需要一次裂变需要3个随机数来确定个随机数来确定裂变点的位置裂变点的位置。合肥工业大学电子科学与应用物理学院立体角：立体角：d=sindd=-ddcos按立体角内均匀分布是指按立体角内均匀分布是指:角在角在02之间均匀分布；之间均匀分布；在在之间，以之间，以cos的值的值在在-11之间均匀分布。之间均匀分布。dd（X0,Y0,Z0）ZXY注意：注意：cos均匀分布，不是均匀分布，不是角均匀分布。如果按照角均匀分布。如果按照为为均
15、匀分布抽样，则在均匀分布抽样，则在=/2附近有较大的权重。附近有较大的权重。一次裂变需要一次裂变需要4个随机数来确定两个中子的方向。个随机数来确定两个中子的方向。合肥工业大学电子科学与应用物理学院(3)(3)中子平均自由程中子平均自由程 平均自由程平均自由程：中子与铀核碰撞所走的平均路程。：中子与铀核碰撞所走的平均路程。在平均自由程内与核碰撞几率相同：在平均自由程内与核碰撞几率相同：例如平均例如平均自自由程为由程为1cm，中子在核块内飞越的距离为，中子在核块内飞越的距离为0.3cm，则中子有则中子有30%几率的与铀核发生碰撞被吸收。几率的与铀核发生碰撞被吸收。中子飞越距离中子飞越距离d可以用可
16、以用0,1之间的随机数表示。之间的随机数表示。一次裂变需要一次裂变需要2个个随机数来确定两个随机数来确定两个中子的平均自由程。中子的平均自由程。cossinsincossin010101dZZdYYdXX中子的位置为中子的位置为 合肥工业大学电子科学与应用物理学院判断中子与铀核发生碰撞判断中子与铀核发生碰撞：若（：若（X1，Y1，Z1）点）点在体积内就认为可以发生碰撞被吸收，反之中子在体积内就认为可以发生碰撞被吸收，反之中子就不会引起下一次裂变。就不会引起下一次裂变。NNfin若有若有N个随机点裂变，个随机点裂变，Nin次总碰撞数（有效中子次总碰撞数（有效中子数），则残存因子数），则残存因子合
17、肥工业大学电子科学与应用物理学院M,N,S,Ep;Nin=0,b=(M/S2)1/3,a=(MS)1/3,J=1(k)=rand,k=1,2,9;X0=a*(1)-0.5),Y0=a*(2)-0.5),Z0=b*(3)-0.5)k =1=2（2k+2）cos=2(2k+3)-0.5,d=(k+7)X1=X0+dsincosY1=Y0+dsinsinZ1=Z0+dcos开始开始J=N N Y 第一个中子第一个中子 k=2F=Nin/N,输出输出F,M,SNin=Nin+1结束结束N Y 第二个中子第二个中子(|X1|-a/20).and.(|Y1|-a/2)0.and.(|Z1|-b/20)k=
18、2N Y J=J+1N Y|F-1|Ep 计算残存因子的流程图计算残存因子的流程图合肥工业大学电子科学与应用物理学院以下给出计算残存因子的流程图的说明：以下给出计算残存因子的流程图的说明：M=V：铀核块的质量（假设密度均匀且数值为：铀核块的质量（假设密度均匀且数值为1）s=a/b：核块边长比：核块边长比N：总随机裂变次数总随机裂变次数Ep=10-4：计算：计算f=1时要求的误差（临界状态）时要求的误差（临界状态）（2 2）总共有）总共有N个随机裂变点，每个点要求九个个随机裂变点，每个点要求九个0，1之间的随机数，之间的随机数，i(i=1,2,=1,2,，9 9）3/123/12332)/()(
19、SMbMSasbsabaVM（1 1）由）由M和和s求长方体的边长求长方体的边长 合肥工业大学电子科学与应用物理学院（3）随机裂变点的坐标三个随机裂变点的坐标三个（4 4）裂变后两个中子的出射方向四个）裂变后两个中子的出射方向四个 第一个中子两个方向第一个中子两个方向 第二个中子两个方向第二个中子两个方向)5.0()5.0()5.0(302010aZaYaX)5.0(2cos25141)5.0(2cos27262（5 5）裂变后两个中子的飞越距离两个）裂变后两个中子的飞越距离两个 d1=8 d2=9合肥工业大学电子科学与应用物理学院取九个随机数，按式计算第一个中子的（取九个随机数，按式计算第一
20、个中子的（X1，Y1，Z1），判断此坐标点是否在铀块内，若是在内部，判断此坐标点是否在铀块内，若是在内部则则Nin增加增加1 1，接着同样方法计算第二个中子；，接着同样方法计算第二个中子；重复上步重复上步N次；次；计算计算 f=Nin/N；调节常数调节常数M和和s使使 f=1，得到临界时的，得到临界时的Mc和和s的的值。值。（6）模拟计算：模拟计算：合肥工业大学电子科学与应用物理学院11.7 11.7 蒙特卡罗方法对趋向平衡态的模拟蒙特卡罗方法对趋向平衡态的模拟 势力学第二定律势力学第二定律：孤立体系，在平衡态时熵不变。：孤立体系，在平衡态时熵不变。在非平衡态时熵趋向增加，最后熵达到最大值而在
21、非平衡态时熵趋向增加，最后熵达到最大值而趋向平衡态。趋向平衡态。典型实例典型实例：正中间用隔板隔开的一个密闭盒子，：正中间用隔板隔开的一个密闭盒子，开始时左边充满某种气体，当隔板中间开一小孔开始时左边充满某种气体，当隔板中间开一小孔后，由于分子无规则的热运动，一些分子通过小后，由于分子无规则的热运动，一些分子通过小孔进入隔板的右边，以后隔板两边分子来回运动，孔进入隔板的右边，以后隔板两边分子来回运动，最后达到两边的压强相等。从统计观点来看，时最后达到两边的压强相等。从统计观点来看，时间很长后分子全部位于左边的情况几率很小。间很长后分子全部位于左边的情况几率很小。合肥工业大学电子科学与应用物理学
22、院一、一、趋向平衡态趋向平衡态数学模型数学模型利用利用N个硬币可以描述上个硬币可以描述上述问题的模型述问题的模型（Ehrenfest模型）。模型）。国徽面国徽面分子在左边；分值面分子在左边；分值面分子在右边。分子在右边。初始时初始时N个硬币都是国徽面。个硬币都是国徽面。随机取硬币改变状态，国徽面随机取硬币改变状态，国徽面分值面或分值面分值面或分值面国徽面，研究硬币状态。国徽面，研究硬币状态。经多次取硬币改变状态，结果国徽面和分值面的经多次取硬币改变状态，结果国徽面和分值面的硬币数基本相等（相当于硬币数基本相等（相当于平衡态）。平衡态）。AB合肥工业大学电子科学与应用物理学院二、理论分析二、理论
23、分析经无规地改变经无规地改变X次后，有次后，有n个分值面向上的硬币。个分值面向上的硬币。此时再取下一个硬币，随机取到分值面向上的几此时再取下一个硬币，随机取到分值面向上的几率为率为n/N，相当于有，相当于有n/N的几率分值面的几率分值面国徽面，国徽面，设设PB=n/N。对国徽面向上的硬币而言，随机取到的几率为对国徽面向上的硬币而言，随机取到的几率为1-PB，表示有，表示有1-1-PB的几率使国徽面的几率使国徽面分值面，设分值面，设PA=1-PB。平均来说，下一次试验分值面向上几率净增加为平均来说，下一次试验分值面向上几率净增加为 n=PA-PB=1-2n/N合肥工业大学电子科学与应用物理学院对
24、于每一次状态变化用对于每一次状态变化用X表示，即表示，即 X=X(ti+1)-X(ti)=1 于是前式变为于是前式变为 n=(1-2n/N)X 利用初始条件在利用初始条件在X=0时时n=0，对上式积分得：，对上式积分得：)1(21/2NXBePPB1/2XdXNndn21把把n和和X当成连续处理，当成连续处理，上上式变为式变为 合肥工业大学电子科学与应用物理学院（1）经过经过X次变化后，平均来说硬币分值面向上的次变化后，平均来说硬币分值面向上的几率是几率是PB ；（2 2）当）当X时，时，PB0.5，这就是平衡态。，这就是平衡态。（3 3）按具体过程画出）按具体过程画出PB=n/N随随X变化的
25、曲线，可以变化的曲线，可以看出看出PB是平均地以指数形式趋向平衡值是平均地以指数形式趋向平衡值0.5，PB曲曲线在指数形式的曲线附近进行涨落。线在指数形式的曲线附近进行涨落。说明：说明：合肥工业大学电子科学与应用物理学院三、蒙特卡罗方法模拟三、蒙特卡罗方法模拟 模拟模拟方法方法：初始时：初始时N个硬币的都是国徽面，经过个硬币的都是国徽面，经过多次随机取硬币，每取一次硬币改变一次状态，多次随机取硬币，每取一次硬币改变一次状态，即国徽面即国徽面分值面或分值面分值面或分值面国徽面。最后统计国徽面。最后统计一下硬币国徽面或分值面所占比例。结果应该是一下硬币国徽面或分值面所占比例。结果应该是各占各占1/
26、2，此时也就是，此时也就是平衡态。平衡态。模拟模拟曲线曲线：如果按模拟过程画出：如果按模拟过程画出PB=n/N随随X变化变化的曲线，可以看出的曲线，可以看出PB是平均地以指数形式趋向平是平均地以指数形式趋向平衡值衡值0.5，并在曲线附近涨落变化。，并在曲线附近涨落变化。合肥工业大学电子科学与应用物理学院NU(I)=1结束结束I=1，2，N输入：输入：N，MX，n=0=rand（）；（）；K=INTN*+1；NU(K)=-NU(K)n=n-NU(K)输出：输出：J，n，PB开始开始Y J=1，2，,MXPB=n/NN：硬币的总数：硬币的总数 MX：无规地改变硬币状：无规地改变硬币状态的总次数态的
27、总次数NU(I)：一维数组，用来：一维数组，用来描述每一个硬币的状态。描述每一个硬币的状态。定义定义分值面国徽面 1 1)(INUn：分值面的硬币数：分值面的硬币数K：1,N间的随机数；间的随机数；利用利用N*rand()+1取整取整（MATLAB语言：语言：fix()；C语言：语言：(int)()。合肥工业大学电子科学与应用物理学院实验十实验十 蒙特卡罗方法的应用蒙特卡罗方法的应用一、实验目的一、实验目的1了解蒙特卡罗原理和方法；了解蒙特卡罗原理和方法；2计算计算光子与物质相互作用的抽样问题；光子与物质相互作用的抽样问题；3蒙特卡罗方法对趋向平衡态的模拟。蒙特卡罗方法对趋向平衡态的模拟。二、
28、实验内容二、实验内容1实现伪（赝）随机变量的各种抽样。实现伪（赝）随机变量的各种抽样。2用用MATLAB计算计算光子与物质相互作用。光子与物质相互作用。3蒙特卡罗方法模拟热力学趋向平衡态的物理过程。蒙特卡罗方法模拟热力学趋向平衡态的物理过程。合肥工业大学电子科学与应用物理学院三、练习题三、练习题1.简述蒙特卡罗方法的基本原理。简述蒙特卡罗方法的基本原理。2.铀块的核临界安全尺寸的模拟计算：铀块的核临界安全尺寸的模拟计算：用用MonteCarlo法计算长方体法计算长方体aab铀块的核临界铀块的核临界安全尺寸。安全尺寸。要求（要求（1）说明计算方法。（）说明计算方法。（2）画出流程图。（）画出流程图。（3）编程上机运行。编程上机运行。已知：中子飞射平均自由程为已知：中子飞射平均自由程为=1。
