蒙特卡罗模拟及例课件

1、 例例1 用用MC 模拟求圆周率模拟求圆周率的估计值的估计值. 设二维随机变量设二维随机变量(X, Y)在正方形内服从均匀分布在正方形内服从均匀分布. . (X, Y)落在圆内的落在圆内的概率为概率为: : 110PX2+Y21=4计算机上做计算机上做n次掷点试验：次掷点试验： 产生产生n 对二维随机点对二维随机点(xi，yi) ，i1 ,2, , n .若有若有k 个点落在个点落在l4圆内圆内 其中其中，xi 和和yi 是是RND 随机数对随机数对. . 相当于第相当于第i个随机点个随机点落在落在1/4圆内圆内.检查每对数是否满足检查每对数是否满足: 随机事件随机事件“点落入点落入1/4圆内

2、圆内”的的频率为频率为 k/n 根据概率论中的大数定律，事件发生的频率依根据概率论中的大数定律，事件发生的频率依概率收敛于事件发生的概率概率收敛于事件发生的概率p，即有即有1lim pPnkn122 iiyx得圆周率得圆周率的估计值为的估计值为=4k/n且当且当试验次数足够大试验次数足够大时，其精度也随之提高时，其精度也随之提高。 分析：分析：实际上概率值为实际上概率值为41102 dxx恰为恰为1/4圆的面积圆的面积 频率法：频率法： 利用随机变量落进指定区域内的频率来计算定积分。利用随机变量落进指定区域内的频率来计算定积分。 平均值法：平均值法： 利用随机变量的平均值（数学期望利用随机变量

3、的平均值（数学期望 )来计算定积分。来计算定积分。 badxxfI)(平均值法平均值法的算法如下：的算法如下：(1)(1) 产生产生RND 随机数：随机数：r1，r2，rn；（2）令令 ui=a(ba)ri，i=1，2，n； （3）计算计算 作为作为I的估计值的估计值。 niiufnab1)(原理分析：原理分析： 设随机变量设随机变量1，2，n相互独立相互独立，且且iU(0,1)f(i)，i=1，2，n 相互独立同分布相互独立同分布abIdxxfabfEbai )(1)( 由由(强强)大数定律知大数定律知.)(1lim1eaabIfnniin 即以概率为即以概率为1 成立。成立。 当当n 足够

4、大时，得近似公式：足够大时，得近似公式： niinbafabdxxfI11)()()(注：注： 平均值法本质上是用样本平均值作为平均值法本质上是用样本平均值作为总体教学期望的估计。总体教学期望的估计。 MC 模拟是一种试验近似方法模拟是一种试验近似方法 , 试验次数如何确定？试验次数如何确定？希望：模拟次数较少、希望：模拟次数较少、 模拟精度较高模拟精度较高频率法的讨论频率法的讨论 频率法是用事件频率法是用事件A出现的频率作为概率出现的频率作为概率p 的估计的估计: nkpn 问题：问题：试验次数试验次数 n 多大时，对给定的置信度多大时，对给定的置信度1（01），），估计精度达到估计精度达到

5、. 即问：取多大的即问：取多大的n 使使 1pnkPppPn成立？成立？答案：答案：22)1( zppn其中，其中， z是正态分布的临界值。是正态分布的临界值。证明证明平均值法平均值法 在给定在给定和和下所需的试验次数下所需的试验次数的估计式为的估计式为 01202222)(11/niixxnSSzn试验次数估计式的分析试验次数估计式的分析 01202222)(11/niixxnSSzn为估计概率为估计概率p做模拟，却又需做模拟，却又需要用要用p去估计模拟次数去估计模拟次数n. 如何计算如何计算S2 ？ 解决方法解决方法：先做：先做n0 次模拟（称为学习样本次模拟（称为学习样本) )，根据学习

6、样本，根据学习样本（1）先求出先求出p的估计，再估计模拟次数的估计，再估计模拟次数n : : 22)1( zppn1.22)1( zppn（2）计算出的样本方差计算出的样本方差S2 , ,用来估计用来估计n. 2. M C模拟的估计精度模拟的估计精度与试验次数与试验次数n的平的平方根成反比方根成反比. 核反应堆屏蔽层是用一定厚度的铅包围反应核反应堆屏蔽层是用一定厚度的铅包围反应堆，用以阻挡或减弱反应堆发出的各种射线。在堆，用以阻挡或减弱反应堆发出的各种射线。在各种射线中，中子对人体伤害极大，因此，在屏各种射线中，中子对人体伤害极大，因此，在屏蔽层的设计中，了解中子穿透屏蔽层的概率，对蔽层的设计

7、中，了解中子穿透屏蔽层的概率，对反应堆的安全运行至关重要反应堆的安全运行至关重要。例例2 核反应堆屏蔽层设计问题核反应堆屏蔽层设计问题1）问题背景）问题背景D返回返回穿透穿透吸收吸收三种状态三种状态阐述中子的运动，阐述中子的运动，为模拟做理论准备为模拟做理论准备2）简化假设：）简化假设： *1 假定屏蔽层平行板厚度为假定屏蔽层平行板厚度为D=3d，其中其中d 为两次碰撞之间中子的平均游动距离；为两次碰撞之间中子的平均游动距离； *2 假设在第假设在第10 次碰撞以后，中子速度下降到次碰撞以后，中子速度下降到为某一很小数值而终止运动（被引收）为某一很小数值而终止运动（被引收）. . *3 假定中

8、子在屏蔽层内相继两次碰撞之间游假定中子在屏蔽层内相继两次碰撞之间游动的距离服从指数分布动的距离服从指数分布； *4 中子经碰撞后的弹射角中子经碰撞后的弹射角 U(0, 2).3）中子运动的数学描述）中子运动的数学描述引进变量引进变量： 弹射角弹射角i 第第i 次碰撞后中子的运动方向与次碰撞后中子的运动方向与x 轴正向的夹角轴正向的夹角. . xi 第第i 次碰撞后中子所处位置与屏蔽层内壁的距离次碰撞后中子所处位置与屏蔽层内壁的距离.DD0 xxi iR i 中子在第中子在第 i 次碰撞前后的游动距离次碰撞前后的游动距离.以上三个变量均为以上三个变量均为随机变量随机变量（1） xi0，中子返回反

9、应堆中子返回反应堆；（2） xiD，中子穿透屏蔽层；中子穿透屏蔽层； 经过第经过第i 次碰撞，中子在屏蔽层内的位置是次碰撞，中子在屏蔽层内的位置是 xi=xi1+Ricosi ，i=1，2，10 ， 中子在屏蔽层里随机游动，第中子在屏蔽层里随机游动，第 i 次碰撞以后，按照它的位置坐标次碰撞以后，按照它的位置坐标 xi，可能有以下三可能有以下三种情况发生：种情况发生：4）模拟过程）模拟过程(1) 产生产生RND随机数对随机数对(ri, ui )；（3） 0 xiD，若若i10，中子中子在屏蔽层内继续运动，在屏蔽层内继续运动，若若 i=10，中子被屏蔽层吸收。中子被屏蔽层吸收。中子三状态中子三状

10、态判别准则判别准则(2) 将将(ri, ui )代入公式计算代入公式计算 ,2,lniiiiurdR (i=1,2,3，10) (3) 将将(Ri, i) 代入公式代入公式 xi = xi1+Ricosi ，i=1，2，10 (4) 判断中子是否穿透屏蔽层判断中子是否穿透屏蔽层.5） 模拟结果分析模拟结果分析 要求穿透屏蔽层的概率数量级为要求穿透屏蔽层的概率数量级为1061010，按假设条件得到一次按假设条件得到一次模拟结果模拟结果如下如下： 第第i 次中子的移动距离和弹射角。次中子的移动距离和弹射角。计算出第计算出第i 次碰撞中子与内壁的距离次碰撞中子与内壁的距离xi . 中子数中子数( (

11、个个) ) 穿透穿透(%) 吸收吸收(%) 返回返回( (%) 100100030005000 30.026.026.526.3 28.023.421.8 22.0 42.050.651.7 51.7 中子穿透屏蔽层的百分比超过了中子穿透屏蔽层的百分比超过了1/4，模拟结模拟结果表明屏蔽层厚度果表明屏蔽层厚度D=3d不合适不合适. 多厚的屏蔽层才能使穿透的概率多厚的屏蔽层才能使穿透的概率W106？问问 题题： 如何解决这个问题？如何解决这个问题？思路？思路？1.计算机收索法计算机收索法 2. 分析法分析法 12mDDD 设计屏蔽层的厚度：设计屏蔽层的厚度： x=mD 将屏蔽层视为将屏蔽层视为m层厚度均为层厚度均为D的平行板的平行板. . 增大屏蔽层的厚度，如增大屏蔽层的厚度，如D=6d、12d、24d、36d，交由计算机进行模拟，并搜索到所交由计算机进行模拟，并搜索到所求解求解. 由于碰撞的能量损失，中子穿过屏蔽层的平由于碰撞的能量损失，中子穿过屏蔽层的平均速度会逐层下降均速度会逐层下降. . 设设WD 是中子穿过厚度为是中子穿过厚度为D 屏蔽层的概率，则屏蔽层的概率，则穿过整个屏蔽层的概率穿过整个屏蔽层的概率W满足满足 mDWW 利用模拟结果利用模拟结果：当当D=3d，WD0.25，令令(WD)m106， 或或 (m