《计算物理(本科)》第11章课件

第十一章蒙特卡罗方法§11.1蒙特卡罗方法概述蒙特卡罗(MonteCarlo)方法：利用随机数统计地去计算和模拟给定的问题。说明：1、也称统计模拟法、随机抽样法或统计试验法。2、MonteCarlo方法的命名：世界上著名的赌城摩洛哥的MonteCarlo。3、方法：先构造一个与物理问题等价的随机过程，当完成大量的随机试验后，结果由多次事件的平均值给出。合肥工业大学电子科学与应用物理学院第十一章蒙特卡罗方法§11.1蒙特卡罗方法概述合15、求解确定性物理问题：如微分方程、定积分、线性方程等。同其它数值计算方法相比是速度慢。6、求解复杂物理问题：如果物理学问题的严格算法不知道，或非常复杂，则蒙特卡罗方法有意想不到的成功。特别在分子运动学、输运现象、布朗运动、放射性衰变等问题中，由于本身有一定的统计规律性，这种方法很奏效。7、计算机模拟：蒙特卡罗方法广泛应用于“计算机模拟”，在计算机上模拟真实过程。4、随机数的抽样：现都用计算机程序来产生，一般不用物理方法抽样。计算机产生的是伪随机数。合肥工业大学电子科学与应用物理学院5、求解确定性物理问题：如微分方程、定积分、线性方程等。同其2§11.2蒙特卡罗方法的原理用蒙特卡罗方法寻找某个未知量x时，利用计算机产生的随机变量ξ，得到的期望值E(ξ)x=E(ξ)为了估算ξ的平均值，构造随机变量ξ的若干独立的实验数序列{ξi|i=1，2，3，……n}，得合肥工业大学电子科学与应用物理学院§11.2蒙特卡罗方法的原理用蒙特卡罗方法寻找某个未知量3§11.3伪（赝）随机变量的抽样实际上，大多数伪随机变量不满足[0，1]均匀分布，而是具有符合分布密度函数为f(x)的分布。抽取符合f(x)随机数的步骤如下：（2）从上面随机数总体中抽取一个简单子样，使它满足分布密度函数f(x)

ξ：ξ1，ξ2，ξ3……（1）在[0，1]之间抽取均匀分布的伪随机数序列

γ：γ1，γ2，γ3……合肥工业大学电子科学与应用物理学院§11.3伪（赝）随机变量的抽样实际上，大多数伪随机变量4一、离散型分布随机变量的直接抽样法随机变量ξ的概率表

其中xi为离散型随机变量ξ的跳跃点,Pi为相应概率。变换方法：第一步构造一个矢量

合肥工业大学电子科学与应用物理学院一、离散型分布随机变量的直接抽样法随机变量ξ的概率表变换方法5第二步产生一个[0，1]间随机数γ，找到区间使γ恰好落在这个区间内。第三步ξP取对应的值xj，其表达式为合肥工业大学电子科学与应用物理学院第二步产生一个[0，1]间随机数γ，找到区间第三步6例1.γ光子与物质相互作用的抽样问题。物理过程：γ光子与物质作用有康普顿效应、光电效应和电子对效应三种类型。其中光电效应和电子对效应为光子吸收过程。设三种过程的碰撞截面分别σs、σe和σp，则总截面σT=σs+σe+σp。根据给定的[0，1]之间均匀分布的随机数γ，求应产生那种效应？合肥工业大学电子科学与应用物理学院例1.γ光子与物质相互作用的抽样问题。物理过程：γ光子与物质7（3）产生[0,1]间均匀随机数γ：若γ＜σs/σT，则发生康普顿散射；若σs/σT≤γ＜(σs+σe)/σT，则发生光电过程；若γ≥(σs+σe)/σT，则发生电子对产生过程。（2）构造一个矢量解：（1）随机变量的概率表合肥工业大学电子科学与应用物理学院（3）产生[0,1]间均匀随机数γ：8二、连续型分布随机变量的直接抽样法一个γ对应一个ξ。如果积分式能求出解析形式的反函数，则得到变换公式

ξF=F-1(γ)11.12f(x)xaξ1-γγb对连续型随机变量ξ的值，可直接解连续分布方程其中f(x)为在（a，b）区间的归一化分布密度函数。合肥工业大学电子科学与应用物理学院二、连续型分布随机变量的直接抽样法一个γ对应一个ξ。如果积分9例1.均匀分布：

在区间（a，b）内找出符合均匀分布的随机变量。解法2：线性变换，令

解法1：合肥工业大学电子科学与应用物理学院例1.均匀分布：

在区间（a，b）内找出符合均匀分布的随机10例2．指数分布：

指数分布的一般形式为f(x)=λe-λxx≥0

其中λ>0，求符合此分布的随机变量ξ求反函数ξ=-1/λln(1-F(ξ))其中0≤F≤1，由上式得ξ=-1/λln(1-γ)因为1-γ∈[0,1]和γ∈[0,1]是等价，所以其式等价于

ξ=-1/λlnγ合肥工业大学电子科学与应用物理学院例2．指数分布：

指数分布的一般形式为f(x)=λe11§11.4MonteCarlo方法求圆周率π求π的物理模型如图所示，边长为1的正方形内切一个圆。于是正方形和圆的面积分别为S正=1，S圆=π/4T正：计算机掷出的随机数总数；T圆：随机点(x,y)落在圆内的总数。即，当时的随机点总数。y0.50.5x若正方形内有T正个随机点(x,y)x=γ1-0.5；y=γ2-0.5其中γ为[0，1]之间的随机数。圆内的随机点数为T圆，那么有S正/S圆=T正/T圆得:π=4T正/T圆合肥工业大学电子科学与应用物理学院§11.4MonteCarlo方法求圆周率π求π的12K=10,M=1000,P=0,i=1j=1j=j+1P=P+1PI=4P/(i*(j-1)),printf(PI)i=i+1开始结束R2<=0.25i<=Kj<=M合肥工业大学电子科学与应用物理学院K=10,M=1000,P=0,i=1j=1j=j+13%用MonteCarlo方法对圆周率进行计算的Matlab程序，文件名PI_val.m%K=100;M=5000;p=0;fori=1:Kforj=1:MX=rand(1,2);R2=(X(1,1)-0.5)^2+(X(1,2)-0.5)^2;ifR2<0.25p=p+1;endendiPI01=4*p/(i*M)end合肥工业大学电子科学与应用物理学院%用MonteCarlo方法对圆周率进行计算的Matlab14§11.5蒙特卡罗方法求解确定性问题

(定积分)

蒙特卡罗方法求解确定性问题的基本思想是：对于给定的确定性问题，设计一个概率模型，然后采用一定的抽样方法按照所设计的概率模型抽样，最后把由这个模型产生的一个数学特征作为该确定问题的近似解。合肥工业大学电子科学与应用物理学院§11.5蒙特卡罗方法求解确定性问题

15一．一维定积分计算f(x)xoabf(x)1．平均值法设定积分为思想：随机分割n个条，在什么位置分割由随机数xi确定取[0,1]间均匀分布的随机数序列γi(i=1,2,3……，n)，并令只要n足够大，则有合肥工业大学电子科学与应用物理学院一．一维定积分计算f(x)xoabf(x)1．平均值法思想：162．随机投点法g(x)xoabg(x)ML设积分若，且g(x)不满足0≤g(x)≤1。可以设法将积分转换成其中，M和L为g(x)的最大值和最小值，如图所示。

合肥工业大学电子科学与应用物理学院2．随机投点法g(x)xoabg(x)ML设积分若17（1）在单位正方形产生两个随机ξ1,ξ2∈[0,1]（2）检验随机点(ξ1,ξ2)是否落入G区，如果条件ξ2≤g(ξ1)满足，则记录点(ξ1,ξ2)落入G内一次(记录器m=m+1)。（3）设在N次试验中，有m个随机点满足上式，那么I≈m/N合肥工业大学电子科学与应用物理学院（1）在单位正方形产生两个随机ξ1,ξ2∈[0,1]（2）检18二．多重定积分计算设多重定积分为取[0，1]间均匀分布随机点列为并令只要m足够大，则有合肥工业大学电子科学与应用物理学院二．多重定积分计算设多重定积分为取[0，1]间均匀分布随机点19实际链式反应过程是一个随机过程，可以利用蒙特卡罗方法来研究链式反应。§11.6蒙特卡罗方法对链式反应的模拟

（核临界安全模拟计算）

链式反应模型：某些元素(如235U、239U)是不稳定的，它能自发地裂变，变成碎片，同时放出中子。中子被其他核材料吸收，使其更不稳，导致迅速裂变，放出更多的中子。合肥工业大学电子科学与应用物理学院实际链式反应过程是一个随机过程，可以利用蒙特卡罗方法来研究链20核爆炸：一般235U核的半衰期很长，约7.07*109年。因此在某一时刻内，只有非常少量的235U发生裂变，释放的能量很小。但是当裂变一旦引起链式反应后，将释放出大量能量，导致核爆炸。链式反应要求：裂变产生的二个中子中，平均来说被其他235U核吸收的中子要大于1，以保证增殖过程。若小于1，则是熄灭的过程。能否链式反应，与235U核团的体积(质量)有关，裂变中子在小质量的235U核团中，吸收几率很小，大质量235U核团中才能持续裂变。临界质量：能维持链式反应的核材料最小质量。合肥工业大学电子科学与应用物理学院核爆炸：一般235U核的半衰期很长，约7.07*109年。21临界条件定量分析

残存因子：N次核裂变，Nin个中子被铀核吸收核爆炸：当铀块的质量大于MC时，发生链式反应。原子弹原理：取两块小于临界质量的235U，将两者合在一起，质量超过临界质量。安全运输：理论计算临界质量十分重要。

f>1表示裂变数增加，发生链式反应；f<1表示反应逐渐停止；f=1表示处于临界状态，用临界质量Mc表示。合肥工业大学电子科学与应用物理学院临界条件定量分析残存因子：N次核裂变，Nin个中子被铀核吸22用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然后计算放出中子被吸收引起的裂变数Nin，求出相应的f值。

为简单起见，考虑铀的几何形状为长方形a*a*b，如图所示。蒙特卡罗法研究这种形状的核材料临界问题：ZYX(X0Y0Z0)(X1Y1Z1)Oaab合肥工业大学电子科学与应用物理学院用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然23（1）核裂变的随机位置

（X0，Y0，Z0）随机点的坐标范围为

-0.5a<X0<0.5a；-0.5a<Y0<0.5a；-0.5b<Z0<0.5b

（2）裂变的两个中子出射方向(θ,φ)裂变产生的中子出射方向用θ和φ来描述。从（X0，Y0，Z0）点放出的中子与以此为圆心的单位球上某一面积发生碰撞，它的碰撞几率只与球面上被碰撞的面积大小成正比。或者说碰撞按立体角均匀分布。一次裂变需要3个随机数来确定裂变点的位置。合肥工业大学电子科学与应用物理学院（1）核裂变的随机位置（X0，Y0，Z0）随机点的坐标24立体角：dΩ=sinθdθdφ=-dφdcosθ按立体角内均匀分布是指:

φ角在0～2π之间均匀分布；θ在０～π之间，以cosθ的值在-1～1之间均匀分布。dφdθ（X0,Y0,Z0）ZXYθφ注意：

cosθ均匀分布，不是θ角均匀分布。如果按照θ为均匀分布抽样，则在θ=π/2附近有较大的权重。

一次裂变需要4个随机数来确定两个中子的方向。合肥工业大学电子科学与应用物理学院立体角：dΩ=sinθdθdφ=-dφdcosθ按立体角内均25(3)中子平均自由程

平均自由程：中子与铀核碰撞所走的平均路程。在平均自由程内与核碰撞几率相同：例如平均自由程为1cm，中子在核块内飞越的距离为0.3cm，则中子有30%几率的与铀核发生碰撞被吸收。中子飞越距离d可以用[0,1]之间的随机数表示。一次裂变需要2个随机数来确定两个中子的平均自由程。中子的位置为

合肥工业大学电子科学与应用物理学院(3)中子平均自由程平均自由程：中子与铀核碰撞所走的平均26判断中子与铀核发生碰撞：若（X1，Y1，Z1）点在体积内就认为可以发生碰撞被吸收，反之中子就不会引起下一次裂变。若有N个随机点裂变，Nin次总碰撞数（有效中子数），则残存因子合肥工业大学电子科学与应用物理学院判断中子与铀核发生碰撞：若（X1，Y1，Z1）点在体积内就认27M,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+dsinθcosφY1=Y0+dsinθsinφZ1=Z0+dcosθ开始J<=N

N

Y

第一个中子①⑤②③④k’=2F=Nin/N,输出F,M,SNin=Nin+1结束N

Y

第二个中子(|X1|-a/2<0).and.(|Y1|-a/2)<0.and.(|Z1|-b/2<0)k’=2N

Y

J=J+1①②③④⑤N

Y

|F-1|<Ep

计算残存因子的流程图合肥工业大学电子科学与应用物理学院M,N,S,Ep;Nin=0,b=(M/S2)1/3,a=(28以下给出计算残存因子的流程图的说明：

M=V：铀核块的质量（假设密度均匀且数值为1）s=a/b：核块边长比N：总随机裂变次数Ep=10-4：计算f=1时要求的误差（临界状态）（2）总共有N个随机裂变点，每个点要求九个[0，1]之间的随机数，γi(i=1,2,……，9）（1）由M和s求长方体的边长

合肥工业大学电子科学与应用物理学院以下给出计算残存因子的流程图的说明：M=V：铀核块的质量29（3）随机裂变点的坐标三个

（4）裂变后两个中子的出射方向四个第一个中子两个方向第二个中子两个方向

（5）裂变后两个中子的飞越距离两个

d1=γ8d2=γ9合肥工业大学电子科学与应用物理学院（3）随机裂变点的坐标三个（4）裂变后两个中子的出射方向四30取九个随机数，按式计算第一个中子的（X1，Y1，Z1），判断此坐标点是否在铀块内，若是在内部则Nin增加1，接着同样方法计算第二个中子；重复上步N次；计算f=Nin/N；调节常数M和s使f=1，得到临界时的Mc和s的值。

（6）模拟计算：合肥工业大学电子科学与应用物理学院取九个随机数，按式计算第一个中子的（X1，Y1，Z1），31§11.7蒙特卡罗方法对趋向平衡态的模拟

势力学第二定律：孤立体系，在平衡态时熵不变。在非平衡态时熵趋向增加，最后熵达到最大值而趋向平衡态。典型实例：正中间用隔板隔开的一个密闭盒子，开始时左边充满某种气体，当隔板中间开一小孔后，由于分子无规则的热运动，一些分子通过小孔进入隔板的右边，以后隔板两边分子来回运动，最后达到两边的压强相等。从统计观点来看，时间很长后分子全部位于左边的情况几率很小。合肥工业大学电子科学与应用物理学院§11.7蒙特卡罗方法对趋向平衡态的模拟势力学第二定律：32一、趋向平衡态数学模型利用N个硬币可以描述上述问题的模型（Ehrenfest模型）。国徽面→分子在左边；分值面→分子在右边。初始时N个硬币都是国徽面。随机取硬币改变状态，国徽面→分值面或分值面→国徽面，研究硬币状态。经多次取硬币改变状态，结果国徽面和分值面的硬币数基本相等（相当于平衡态）。

AB合肥工业大学电子科学与应用物理学院一、趋向平衡态数学模型利用N个硬币可以描述上述问题的模型国徽33二、理论分析经无规地改变X次后，有n个分值面向上的硬币。此时再取下一个硬币，随机取到分值面向上的几率为n/N，相当于有n/N的几率分值面→国徽面，设PB=n/N。对国徽面向上的硬币而言，随机取到的几率为1-PB，表示有1-PB的几率使国徽面→分值面，设PA=1-PB。平均来说，下一次试验分值面向上几率净增加为Δn=PA-PB=1-2n/N合肥工业大学电子科学与应用物理学院二、理论分析经无规地改变X次后，有n个分值面向上的硬币。合肥34对于每一次状态变化用ΔX表示，即ΔX=X(ti+1)-X(ti)=1于是前式变为Δn=(1-2n/N)ΔX

利用初始条件在X=0时n=0，对上式积分得：PB1/2X把n和X当成连续处理，上式变为

合肥工业大学电子科学与应用物理学院对于每一次状态变化用ΔX表示，即利用初始条件在X=0时n=035（1）经过X次变化后，平均来说硬币分值面向上的几率是PB；（2）当X→∞时，PB→0.5，这就是平衡态。（3）按具体过程画出PB=n/N随X变化的曲线，可以看出PB是平均地以指数形式趋向平衡值0.5，PB曲线在指数形式的曲线附近进行涨落。说明：合肥工业大学电子科学与应用物理学院（1）经过X次变化后，平均来说硬币分值面向上的几率是PB；36三、蒙特卡罗方法模拟

模拟方法：初始时N个硬币的都是国徽面，经过多次随机取硬币，每取一次硬币改变一次状态，即国徽面→分值面或分值面→国徽面。最后统计一下硬币国徽面或分值面所占比例。结果应该是各占1/2，此时也就是平衡态。模拟曲线：如果按模拟过程画出PB=n/N随X变化的曲线，可以看出PB是平均地以指数形式趋向平衡值0.5，并在曲线附近涨落变化。合肥工业大学电子科学与应用物理学院三、蒙特卡罗方法模拟模拟方法：初始时N个硬币的都是国徽面，37NU(I)=1结束I=1，2，……N输入：N，MX，n=0γ=rand（）；K=INT[N*γ+1]；

NU(K)=-NU(K)n=n-NU(K)输出：J，n，PB开始Y

J=1，2，…，,MXPB=n/NN：硬币的总数

MX：无规地改变硬币状态的总次数NU(I)：一维数组，用来描述每一个硬币的状态。定义n：分值面的硬币数K：[1,N]间的随机数；利用N*rand()+1取整（MATLAB语言：fix()；C语言：(int)(……))。合肥工业大学电子科学与应用物理学院NU(I)=1结束I=1，2，……N输入：N，MX，n=0γ38实验十蒙特卡罗方法的应用一、实验目的1．了解蒙特卡罗原理和方法；2．计算γ光子与物质相互作用的抽样问题；3．蒙特卡罗方法对趋向平衡态的模拟。二、实验内容1．实现伪（赝）随机变量的各种抽样。2．用MATLAB计算γ光子与物质相互作用。3．蒙特卡罗方法模拟热力学趋向平衡态的物理过程。合肥工业大学电子科学与应用物理学院实验十蒙特卡罗方法的应用合肥工业大学电子科学与应用物理学39三、练习题1.简述蒙特卡罗方法的基本原理。2.铀块的核临界安全尺寸的模拟计算：用Monte—Carlo法计算长方体a×a×b铀块的核临界安全尺寸。要求（1）说明计算方法。（2）画出流程图。（3）编程上机运行。已知：中子飞射平均自由程为λ=1。合肥工业大学电子科学与应用物理学院三、练习题合肥工业大学电子科学与应用物理学院40第十一章蒙特卡罗方法§11.1蒙特卡罗方法概述蒙特卡罗(MonteCarlo)方法：利用随机数统计地去计算和模拟给定的问题。说明：1、也称统计模拟法、随机抽样法或统计试验法。2、MonteCarlo方法的命名：世界上著名的赌城摩洛哥的MonteCarlo。3、方法：先构造一个与物理问题等价的随机过程，当完成大量的随机试验后，结果由多次事件的平均值给出。合肥工业大学电子科学与应用物理学院第十一章蒙特卡罗方法§11.1蒙特卡罗方法概述合415、求解确定性物理问题：如微分方程、定积分、线性方程等。同其它数值计算方法相比是速度慢。6、求解复杂物理问题：如果物理学问题的严格算法不知道，或非常复杂，则蒙特卡罗方法有意想不到的成功。特别在分子运动学、输运现象、布朗运动、放射性衰变等问题中，由于本身有一定的统计规律性，这种方法很奏效。7、计算机模拟：蒙特卡罗方法广泛应用于“计算机模拟”，在计算机上模拟真实过程。4、随机数的抽样：现都用计算机程序来产生，一般不用物理方法抽样。计算机产生的是伪随机数。合肥工业大学电子科学与应用物理学院5、求解确定性物理问题：如微分方程、定积分、线性方程等。同其42§11.2蒙特卡罗方法的原理用蒙特卡罗方法寻找某个未知量x时，利用计算机产生的随机变量ξ，得到的期望值E(ξ)x=E(ξ)为了估算ξ的平均值，构造随机变量ξ的若干独立的实验数序列{ξi|i=1，2，3，……n}，得合肥工业大学电子科学与应用物理学院§11.2蒙特卡罗方法的原理用蒙特卡罗方法寻找某个未知量43§11.3伪（赝）随机变量的抽样实际上，大多数伪随机变量不满足[0，1]均匀分布，而是具有符合分布密度函数为f(x)的分布。抽取符合f(x)随机数的步骤如下：（2）从上面随机数总体中抽取一个简单子样，使它满足分布密度函数f(x)

ξ：ξ1，ξ2，ξ3……（1）在[0，1]之间抽取均匀分布的伪随机数序列

γ：γ1，γ2，γ3……合肥工业大学电子科学与应用物理学院§11.3伪（赝）随机变量的抽样实际上，大多数伪随机变量44一、离散型分布随机变量的直接抽样法随机变量ξ的概率表

其中xi为离散型随机变量ξ的跳跃点,Pi为相应概率。变换方法：第一步构造一个矢量

合肥工业大学电子科学与应用物理学院一、离散型分布随机变量的直接抽样法随机变量ξ的概率表变换方法45第二步产生一个[0，1]间随机数γ，找到区间使γ恰好落在这个区间内。第三步ξP取对应的值xj，其表达式为合肥工业大学电子科学与应用物理学院第二步产生一个[0，1]间随机数γ，找到区间第三步46例1.γ光子与物质相互作用的抽样问题。物理过程：γ光子与物质作用有康普顿效应、光电效应和电子对效应三种类型。其中光电效应和电子对效应为光子吸收过程。设三种过程的碰撞截面分别σs、σe和σp，则总截面σT=σs+σe+σp。根据给定的[0，1]之间均匀分布的随机数γ，求应产生那种效应？合肥工业大学电子科学与应用物理学院例1.γ光子与物质相互作用的抽样问题。物理过程：γ光子与物质47（3）产生[0,1]间均匀随机数γ：若γ＜σs/σT，则发生康普顿散射；若σs/σT≤γ＜(σs+σe)/σT，则发生光电过程；若γ≥(σs+σe)/σT，则发生电子对产生过程。（2）构造一个矢量解：（1）随机变量的概率表合肥工业大学电子科学与应用物理学院（3）产生[0,1]间均匀随机数γ：48二、连续型分布随机变量的直接抽样法一个γ对应一个ξ。如果积分式能求出解析形式的反函数，则得到变换公式

ξF=F-1(γ)11.12f(x)xaξ1-γγb对连续型随机变量ξ的值，可直接解连续分布方程其中f(x)为在（a，b）区间的归一化分布密度函数。合肥工业大学电子科学与应用物理学院二、连续型分布随机变量的直接抽样法一个γ对应一个ξ。如果积分49例1.均匀分布：

在区间（a，b）内找出符合均匀分布的随机变量。解法2：线性变换，令

解法1：合肥工业大学电子科学与应用物理学院例1.均匀分布：

在区间（a，b）内找出符合均匀分布的随机50例2．指数分布：

指数分布的一般形式为f(x)=λe-λxx≥0

其中λ>0，求符合此分布的随机变量ξ求反函数ξ=-1/λln(1-F(ξ))其中0≤F≤1，由上式得ξ=-1/λln(1-γ)因为1-γ∈[0,1]和γ∈[0,1]是等价，所以其式等价于

ξ=-1/λlnγ合肥工业大学电子科学与应用物理学院例2．指数分布：

指数分布的一般形式为f(x)=λe51§11.4MonteCarlo方法求圆周率π求π的物理模型如图所示，边长为1的正方形内切一个圆。于是正方形和圆的面积分别为S正=1，S圆=π/4T正：计算机掷出的随机数总数；T圆：随机点(x,y)落在圆内的总数。即，当时的随机点总数。y0.50.5x若正方形内有T正个随机点(x,y)x=γ1-0.5；y=γ2-0.5其中γ为[0，1]之间的随机数。圆内的随机点数为T圆，那么有S正/S圆=T正/T圆得:π=4T正/T圆合肥工业大学电子科学与应用物理学院§11.4MonteCarlo方法求圆周率π求π的52K=10,M=1000,P=0,i=1j=1j=j+1P=P+1PI=4P/(i*(j-1)),printf(PI)i=i+1开始结束R2<=0.25i<=Kj<=M合肥工业大学电子科学与应用物理学院K=10,M=1000,P=0,i=1j=1j=j+53%用MonteCarlo方法对圆周率进行计算的Matlab程序，文件名PI_val.m%K=100;M=5000;p=0;fori=1:Kforj=1:MX=rand(1,2);R2=(X(1,1)-0.5)^2+(X(1,2)-0.5)^2;ifR2<0.25p=p+1;endendiPI01=4*p/(i*M)end合肥工业大学电子科学与应用物理学院%用MonteCarlo方法对圆周率进行计算的Matlab54§11.5蒙特卡罗方法求解确定性问题

(定积分)

蒙特卡罗方法求解确定性问题的基本思想是：对于给定的确定性问题，设计一个概率模型，然后采用一定的抽样方法按照所设计的概率模型抽样，最后把由这个模型产生的一个数学特征作为该确定问题的近似解。合肥工业大学电子科学与应用物理学院§11.5蒙特卡罗方法求解确定性问题

55一．一维定积分计算f(x)xoabf(x)1．平均值法设定积分为思想：随机分割n个条，在什么位置分割由随机数xi确定取[0,1]间均匀分布的随机数序列γi(i=1,2,3……，n)，并令只要n足够大，则有合肥工业大学电子科学与应用物理学院一．一维定积分计算f(x)xoabf(x)1．平均值法思想：562．随机投点法g(x)xoabg(x)ML设积分若，且g(x)不满足0≤g(x)≤1。可以设法将积分转换成其中，M和L为g(x)的最大值和最小值，如图所示。

合肥工业大学电子科学与应用物理学院2．随机投点法g(x)xoabg(x)ML设积分若57（1）在单位正方形产生两个随机ξ1,ξ2∈[0,1]（2）检验随机点(ξ1,ξ2)是否落入G区，如果条件ξ2≤g(ξ1)满足，则记录点(ξ1,ξ2)落入G内一次(记录器m=m+1)。（3）设在N次试验中，有m个随机点满足上式，那么I≈m/N合肥工业大学电子科学与应用物理学院（1）在单位正方形产生两个随机ξ1,ξ2∈[0,1]（2）检58二．多重定积分计算设多重定积分为取[0，1]间均匀分布随机点列为并令只要m足够大，则有合肥工业大学电子科学与应用物理学院二．多重定积分计算设多重定积分为取[0，1]间均匀分布随机点59实际链式反应过程是一个随机过程，可以利用蒙特卡罗方法来研究链式反应。§11.6蒙特卡罗方法对链式反应的模拟

（核临界安全模拟计算）

链式反应模型：某些元素(如235U、239U)是不稳定的，它能自发地裂变，变成碎片，同时放出中子。中子被其他核材料吸收，使其更不稳，导致迅速裂变，放出更多的中子。合肥工业大学电子科学与应用物理学院实际链式反应过程是一个随机过程，可以利用蒙特卡罗方法来研究链60核爆炸：一般235U核的半衰期很长，约7.07*109年。因此在某一时刻内，只有非常少量的235U发生裂变，释放的能量很小。但是当裂变一旦引起链式反应后，将释放出大量能量，导致核爆炸。链式反应要求：裂变产生的二个中子中，平均来说被其他235U核吸收的中子要大于1，以保证增殖过程。若小于1，则是熄灭的过程。能否链式反应，与235U核团的体积(质量)有关，裂变中子在小质量的235U核团中，吸收几率很小，大质量235U核团中才能持续裂变。临界质量：能维持链式反应的核材料最小质量。合肥工业大学电子科学与应用物理学院核爆炸：一般235U核的半衰期很长，约7.07*109年。61临界条件定量分析

残存因子：N次核裂变，Nin个中子被铀核吸收核爆炸：当铀块的质量大于MC时，发生链式反应。原子弹原理：取两块小于临界质量的235U，将两者合在一起，质量超过临界质量。安全运输：理论计算临界质量十分重要。

f>1表示裂变数增加，发生链式反应；f<1表示反应逐渐停止；f=1表示处于临界状态，用临界质量Mc表示。合肥工业大学电子科学与应用物理学院临界条件定量分析残存因子：N次核裂变，Nin个中子被铀核吸62用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然后计算放出中子被吸收引起的裂变数Nin，求出相应的f值。

为简单起见，考虑铀的几何形状为长方形a*a*b，如图所示。蒙特卡罗法研究这种形状的核材料临界问题：ZYX(X0Y0Z0)(X1Y1Z1)Oaab合肥工业大学电子科学与应用物理学院用计算机模拟具有一定大小及形状的体积内发生的大量随机裂变，然63（1）核裂变的随机位置

（X0，Y0，Z0）随机点的坐标范围为

-0.5a<X0<0.5a；-0.5a<Y0<0.5a；-0.5b<Z0<0.5b

（2）裂变的两个中子出射方向(θ,φ)裂变产生的中子出射方向用θ和φ来描述。从（X0，Y0，Z0）点放出的中子与以此为圆心的单位球上某一面积发生碰撞，它的碰撞几率只与球面上被碰撞的面积大小成正比。或者说碰撞按立体角均匀分布。一次裂变需要3个随机数来确定裂变点的位置。合肥工业大学电子科学与应用物理学院（1）核裂变的随机位置（X0，Y0，Z0）随机点的坐标64立体角：dΩ=sinθdθdφ=-dφdcosθ按立体角内均匀分布是指:

φ角在0～2π之间均匀分布；θ在０～π之间，以cosθ的值在-1～1之间均匀分布。dφdθ（X0,Y0,Z0）ZXYθφ注意：

cosθ均匀分布，不是θ角均匀分布。如果按照θ为均匀分布抽样，则在θ=π/2附近有较大的权重。

一次裂变需要4个随机数来确定两个中子的方向。合肥工业大学电子科学与应用物理学院立体角：dΩ=sinθdθdφ=-dφdcosθ按立体角内均65(3)中子平均自由程

平均自由程：中子与铀核碰撞所走的平均路程。在平均自由程内与核碰撞几率相同：例如平均自由程为1cm，中子在核块内飞越的距离为0.3cm，则中子有30%几率的与铀核发生碰撞被吸收。中子飞越距离d可以用[0,1]之间的随机数表示。一次裂变需要2个随机数来确定两个中子的平均自由程。中子的位置为

合肥工业大学电子科学与应用物理学院(3)中子平均自由程平均自由程：中子与铀核碰撞所走的平均66判断中子与铀核发生碰撞：若（X1，Y1，Z1）点在体积内就认为可以发生碰撞被吸收，反之中子就不会引起下一次裂变。若有N个随机点裂变，Nin次总碰撞数（有效中子数），则残存因子合肥工业大学电子科学与应用物理学院判断中子与铀核发生碰撞：若（X1，Y1，Z1）点在体积内就认67M,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+dsinθcosφY1=Y0+dsinθsinφZ1=Z0+dcosθ开始J<=N

N

Y

第一个中子①⑤②③④k’=2F=Nin/N,输出F,M,SNin=Nin+1结束N

Y

第二个中子(|X1|-a/2<0).and.(|Y1|-a/2)<0.and.(|Z1|-b/2<0)k’=2N

Y

J=J+1①②③④⑤N

Y

|F-1|<Ep

计算残存因子的流程图合肥工业大学电子科学与应用物理学院M,N,S,Ep;Nin=0,b=(M/S2)1/3,a=(68以下给出计算残存因子的流程图的说明：

M=V：铀核块的质量（假设密度均匀且数值为1）s=a/b：核块边长比N：总随机裂变次数Ep=10-4：计算f=1时要求的误差（临界状态）（2）总共有N个随机裂变点，每个点要求九个[0，1]之间的随机数，γi(i=1,2,……，9）（1）由M和s求长方体的边长

合肥工业大学电子科学与应用物理学院以下给出计算残存因子的流程图的说明：M=V：铀核块的质量69（3）随机裂变点的坐标三个

（4）裂变后两个中子的出射方向四个第一个中子两个方向第二个中子两个方向

（5）裂变后两个中子的飞越距离两个

d1=γ8d2=γ9合肥工业大学电子科学与应用物理学院（3）随机裂变点的坐标三个（4）裂变后两个中子的出射方向四70取九个随机数，按式计算第一个中子的（X1，Y1，Z1），判断此坐标点是否在铀块内，若是在内部则Nin增加1，接着同样方法计算第二个中子；重复上步N次；计算f=Nin/N；调节常数M和s使f=1，得到临界时的Mc和s的值。

（6）模拟计算：合肥工业大学电子科学与应用物理学院取九个随机数，按式计算第一个中子的（X1，Y1，Z1），71§11.7蒙特卡罗方法对趋向平衡态的模拟

势力学第二定律：孤立体系，在平衡态时熵不变。在非平衡态时熵趋向增加，最后熵达到最大值而趋向平衡态。典型实例：正中间用隔板隔开的一个密闭盒子，开始时左边充满某种气体，当隔板中间开一小孔后，由于分子无规则的热运动，一些分子通过小孔进入隔板的右边，以后隔板两边分子来回运动，最后达到两边的压强相等。从统计观点来看，时间很长后分子全部位于左边的情况几率