蒙特卡罗模拟

1、 蒙特卡罗蒙特卡罗( (MonteCarlo)模拟，又称蒙特卡模拟，又称蒙特卡罗方法、统计试验法等罗方法、统计试验法等. . MC模拟是静态模拟，描述特定时间点上模拟是静态模拟，描述特定时间点上的系统行为的系统行为.模拟过程中模拟过程中不出现时间不出现时间参数。参数。基本思想基本思想: :把随机事件把随机事件（变量）的概率特征与（变量）的概率特征与数学分析的解联系起来数学分析的解联系起来. . 概率特征：概率特征：随机事件的概率和随机变量的随机事件的概率和随机变量的数学期望等数学期望等. . 用试验方法确定用试验方法确定 例例7.3.1 用用MC 模拟求圆周率模拟求圆周率的估计值的估计值. 1

2、10 设二维随机变量设二维随机变量(X, Y)在正方形内在正方形内服从均匀分布服从均匀分布. . (X, Y)落在圆内的概率为落在圆内的概率为: : 4122 YXP计算机上做计算机上做n 次掷点试验：次掷点试验： 产生产生n 对二维随机点对二维随机点(xi，yi) ，i1 ,2, , n .xi 和和yi 是是RND 随机数对随机数对. . 检查每对随机数是否满足检查每对随机数是否满足:122 iiyx相当于第相当于第i个随机点落个随机点落在在1/4圆内圆内.若有若有k 个点落在个点落在l4圆内圆内 随机事件随机事件“点落入点落入1/4圆内圆内”的的频率为频率为 k/n 根据概率论中的大数定

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

4、值法平均值法的算法如下：的算法如下：(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 足够大时，得近似公式：足够大时，得近似公式： niinbafabdxxfI11)()

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

6、是事件A出现的频率作为概率出现的频率作为概率p的估计的估计 nkpn 答案：答案：22)1( zppn其中其中, z是正态分布的临界值是正态分布的临界值.n次独立试验中次独立试验中A出现的次数出现的次数knB(n, p).由中由中心极限定理知心极限定理知 )()( pnkpnPppPn )1()1()1(pnpnpnpnpkpnpnPn)1()1(pnpnpnpn 1)1(2 pnpn 11)1(2pnpn令令平均值法平均值法 在给定在给定和和下所需的试验次数下所需的试验次数的估计式为的估计式为 01202222)(11/niixxnSSzn查得正态分布的临界值查得正态分布的临界值z, ,可解

7、得可解得 22)1( zppn试验次数估计式的分析试验次数估计式的分析 01202222)(11/niixxnSSzn为估计概率为估计概率p做做模拟模拟, ,却又需要却又需要用用p去估计模拟去估计模拟次数次数n. 如何计算如何计算S2 ？ 解决方法解决方法：先做：先做n0 次模拟（称为学习样本次模拟（称为学习样本) )，根，根据学习样本据学习样本. .（1）先求出先求出p的估计，再估计模拟次数的估计，再估计模拟次数n : : 22)1(.1 zppn22)1( zppn（2）计算出的样本方差计算出的样本方差S2 , ,用来估计用来估计n. 2. M C模拟的估计精度模拟的估计精度与试验次数与试

8、验次数n的平的平方根成反比方根成反比, 若精度若精度提高提高10倍倍, ,则试验次数则试验次数n要增大要增大100倍倍. . P197表表8.2中列出了置信度为中列出了置信度为0.95 时时, , 在不同在不同精度精度及概率及概率p条件下频率法所需试验次数。条件下频率法所需试验次数。 对该表进行分析，能得到什么结论？对该表进行分析，能得到什么结论？1. 精度提高，试验次数大幅提高；精度提高，试验次数大幅提高；2. 事件发生概率越接近事件发生概率越接近0.5，试验次数越高；，试验次数越高； 核反应堆屏蔽层是用一定厚度的铅包围反应核反应堆屏蔽层是用一定厚度的铅包围反应堆，用以阻挡或减弱反应堆发出的

9、各种射线堆，用以阻挡或减弱反应堆发出的各种射线.在各种射线中在各种射线中, 中子对人体伤害极大，因此，中子对人体伤害极大，因此，在屏蔽层的设计中在屏蔽层的设计中, 了解中子穿透屏蔽层的概了解中子穿透屏蔽层的概率对反应堆的安全运行至关重要率对反应堆的安全运行至关重要.例例7.3.2 核反应堆屏蔽层设计问题核反应堆屏蔽层设计问题1.问题背景问题背景假定屏蔽层是理想的均匀平板假定屏蔽层是理想的均匀平板 一个中子进入屏蔽层后运动的物理过程一个中子进入屏蔽层后运动的物理过程: :中中子以初速度子以初速度v0和方向角和方向角射入屏蔽层射入屏蔽层, ,运动一运动一段距离后与铅核发生碰撞，中子获得新的速度段距

10、离后与铅核发生碰撞，中子获得新的速度及方向及方向(v1,1). 再游动一段距离后再游动一段距离后, ,与铅核发生与铅核发生第二次碰撞第二次碰撞, ,并获得新的状态并获得新的状态(v2,2),如此等等如此等等, ,经过若干次碰撞后经过若干次碰撞后, ,出现下述情况之一时中子出现下述情况之一时中子终止运动过程终止运动过程 1）中子被弹回反应堆；中子被弹回反应堆；2）中子穿透屏蔽层；）中子穿透屏蔽层；3）第第n次碰撞后次碰撞后, ,中子被屏蔽层吸收中子被屏蔽层吸收. .D返回返回穿透穿透吸收吸收三种三种状态状态为使屏蔽层的厚度达到安全设计要求为使屏蔽层的厚度达到安全设计要求, ,在计在计算机上对中子

11、在屏蔽层的运动过程进行模拟算机上对中子在屏蔽层的运动过程进行模拟 阐述中子的运动，阐述中子的运动，为模拟做理论准备为模拟做理论准备2. 简化假设：简化假设： *1 假定屏蔽层平行板厚度为假定屏蔽层平行板厚度为D=3d，其中其中d 为两次碰撞之间中子的平均游动距离；为两次碰撞之间中子的平均游动距离； *2 假设在第假设在第10 次碰撞以后，中子速度下降次碰撞以后，中子速度下降到为某一很小数值而终止运动（被引收）到为某一很小数值而终止运动（被引收）. . 因每次碰撞后因每次碰撞后, 中子因损失一部分能量而速中子因损失一部分能量而速度下降度下降. *3 假定中子在屏蔽层内相继两次碰撞之间假定中子在屏

12、蔽层内相继两次碰撞之间游动的距离服从指数分布游动的距离服从指数分布； *4 中子经碰撞后的弹射角中子经碰撞后的弹射角 U(0, 2).思考思考：请仔细分析以上假设的合理性：请仔细分析以上假设的合理性.3. 中子运动的数学描述中子运动的数学描述引进变量引进变量： 弹射角弹射角i 第第i 次碰撞后中子的运动方向与次碰撞后中子的运动方向与x 轴正向的夹角轴正向的夹角. .xi 第第i 次碰撞后中子所处位置与屏蔽层内次碰撞后中子所处位置与屏蔽层内 壁的距离壁的距离.R i 中子在第中子在第 i 次碰撞前后的游动距离次碰撞前后的游动距离.DD0 xxi i三个变量三个变量均为随机均为随机变量变量 中子在

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

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

15、数( (个个) ) 穿透穿透(%) 吸收吸收(%) 返回返回( (%) 100100030005000 30.026.026.526.3 28.023.421.8 22.0 42.050.651.7 51.7 中子穿透屏蔽层的百分比超过了中子穿透屏蔽层的百分比超过了1/4, 模拟模拟结果表明屏蔽层厚度结果表明屏蔽层厚度D=3d不合适不合适. 多厚的屏蔽层才能使穿透的概率多厚的屏蔽层才能使穿透的概率W106？问问 题题： 如何解决这个问题？如何解决这个问题？思路？思路？1. 计算机收索法计算机收索法 增大屏蔽层的厚度，如增大屏蔽层的厚度，如D=6d、12d、24d、36d，交由计算机进行模拟交由

16、计算机进行模拟, 并搜索到所并搜索到所求解求解.2. 分析法分析法 12mDDD 设计屏蔽层的厚度：设计屏蔽层的厚度： x=mD 将屏蔽层视为将屏蔽层视为m层厚度均为层厚度均为D的平行板的平行板. . 由于碰撞的能量损失由于碰撞的能量损失, ,中子穿过屏蔽层的平中子穿过屏蔽层的平均速度会逐层下降均速度会逐层下降. . 设设WD 是中子穿过厚度为是中子穿过厚度为D 屏蔽层的概率屏蔽层的概率,则则穿过整个屏蔽层的概率穿过整个屏蔽层的概率W 满足满足 mDWW 利用模拟结果利用模拟结果：当当D=3d，WD0.25，令令(WD)m106， 或或 (m)410 6967. 9602. 064log10log6 m取屏蔽层的厚度取屏蔽层的厚度 x=10D=30d，可使穿透屏蔽层的概率可使穿透屏蔽层的概率w106 注注：模拟：模拟5000个中子个中子的运动的运动, ,用穿透屏蔽的用穿透屏蔽的频率估计穿透概率频率估计穿透概率, ,由表由表8.2可知精度大约只有可知精度大约只有1 1, ,模拟精度太低模拟精度太低, ,应适当增大模拟次数应适当增大模拟次数. .总结：总结：模拟的意义？模拟的意义？ 1.模拟方法本质上是模拟方法本质上是试验性试验性的的, ,模拟系统是模拟系统是 现实系