蒙特卡罗(Monte Carlo)模拟

1、17. 蒙特卡罗(Monte Carlo)模拟蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。该方法源于美国在WWI后研制原子弹的“曼哈顿计划”。 S. M. Ulam (1908-1984)和该计划的主持人之一、数学家冯诺伊曼（J.von Neumann ）用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法。其基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定。高速计算机的出现，使得用数学方法在计算机上大量模拟这样的试验成为可能。科技计算及社会中的问题比计算要复杂得多。比如金融

2、衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(Curse of Dimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。Monte Carlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。另一类形式与Monte Carlo方法相似，但理论基础不同的方法“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)近年来也

3、获得迅速发展。我国数学家华罗庚、王元提出的“华王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提出高数百倍，并可计算精确度。http:/ call random_seed call random_number(a)Matlab: x=rand(N) 产生元素在(0, 1)间随机分布的N*N矩阵 s= rand(state,0) 重设该生成函数到初始状态计算程序产生的随机数不是真正的随机数，它们是确定的，但

4、看上去是随机的，且能通过一些随机性的检验，故常称为伪随机数。如对32位字长的计算机real procedure random(xi)integer array (li)nreal array (xi)nl0 any integer that 1 l0 231-1for i=1 to n doli =(231-1)除以16807 li-1的余数xi li/ (231-1)endfor其它算法1. 取初值x0 (0,1)，let xi be the fractional part of (+ xi-1)5 2. Let u0 =1. For i=1,2,3, let ui be the remai

5、nder of (8t-3) ui-1 divided by 28, and xi = ui/2i. Here t can be any large integer注意：上述随机数序列均具周期性，如上页random子程序的周期约230。Maple has a collection of random number generators.如：With(stats)x :=uniform(0.1):seq(x(),i=1.10);Matlab: x=rand(10,1)产生10个随机数。如 (b-a)r+a x (a,b) integer(n+1)r) x 0,1,2,n试产生1000个在椭圆x2

6、+4y2=4内均匀分布的随机点：方法：在 中均匀产生足够 多随机点，当位于椭圆内的点数为1000时停止。可由随机数序列 r (0,1)构造一系列随机数序列1122yx，xy有时随机数产生函数通不过严格的随机性测试。而在一些计算（如MC积分）中随机性非常重要。故使用更长字节的计算机更好。应用之一估计面积和体积Numerical integration选取(0, 1)中随机数序列x1， x2， x3， xn。则niixfnxxf110)(1d)(误差约 ，它并不能和一些高级的数值积分算法比拟，但对多维情况，MC方法却很有吸引力。n1 niiiizyxfnzyxzyxf1101010),(1ddd)

7、,(我们可产生一系列随机数可简单取3个随机数构成一个随机点，即,.,654321),.,(),(654321相应地，niibaxfnxxfab1)(1d)(1 niiidcbayxfnyxyxfabcd1),(1dd),()(1一般地，A) in points random nover of (average) of measure(d)(fAxfA例如： 求其中25. 0)5 . 0()5 . 0( : ),(22yxyxyxyxdd) 1ln(sin在该区域中取5000个点，可得该积分约0.57.计算体积Pseudo-code1sin10 10 102yezxzyxzyxProgram v

8、olume_regionInteger parameter n 5000, iprt 1000Real array integer i, mReal vol,x,y,zCall random(rij)For I=1 to n dox ri1y ri2z ri3If then m m+1If mod(I,iprt)=0 thenVol real(m)/real(I)Output I, volEndifEndforEnd program3)(nijr1 sin2yezxzyx计算体积之作业：冰淇淋的体积2222221) 1(yxzzyxyxz应用之二MC模拟生日问题生日问题假设有假设有n个人在一起

9、，各自的生日为个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直天之一，根据概率理论，与很多人的直觉相反，只需觉相反，只需23个人便有大于个人便有大于50的几的几率人群中至少有率人群中至少有2个人生日相同。个人生日相同。365)1(365.3653633653643653651nn 理论几率 模拟几率100.117 0.110200.411 0.412230.507 0.520300.706 0.692450.941 0.93650 0.986 0.987Real function prob(n,npts)Integer i, nptsLogical birthReal sumS

10、um=0For I=1 to npts doIf birth(n) then sum sum+1EndforProb sum/real(npts)End function蒲丰的投针问题蒲丰的投针问题Buffons Needle problem蒲丰证明在相距蒲丰证明在相距d 的一系列平行线上的一系列平行线上投长为投长为 l 针。针与线的交点数除以投针。针与线的交点数除以投的次数等于的次数等于2 l / (d)。Real function prob(n,npts)Integer i, nptsLogical birthReal sumSum=0For I=1 to npts doIf birth(

11、n) then sum sum+1EndforProb sum/real(npts)End function不妨取不妨取d l。产生随机数。产生随机数u和和 ，其中，其中u为针中心距最近的线的距离，故为针中心距最近的线的距离，故u小小于等于于等于0.5， 为夹角，根据对称性，为夹角，根据对称性，可取可取0到到/2。usin21即为相交usin21算时需要用上是否不自洽（姜冰）作业中子屏蔽问题中子屏蔽问题Neutron Shielding problem假设假设铅墙长为5，中子在铅中的平均自由程为1，中子与铅原子碰撞，中子与铅原子碰撞后各向同性散射。令碰撞后各向同性散射。令碰撞8次后中次后中子能

12、量耗尽，试求穿透铅墙的中子能量耗尽，试求穿透铅墙的中子的比例。子的比例。暂不考虑垂直纸面的运动，则中暂不考虑垂直纸面的运动，则中子的水平位移是子的水平位移是 。1721cos.coscos1入口铅墙（长为5）221点问题点问题Blackjack 21 problem参数拟合问题参数拟合问题Parameter fitting problem辐射转移问题辐射转移问题Radiative transfer problem1 Ahn,S.-H., Lee, H.-W. & Lee, H.-M. 2000, JKAS, 33, p29.2 Zheng Zheng & Jordi Miral

13、da-Escude. 2002, ApJ, 578, p33.3Watson, A.M. & Henney, W.J. 2001, RMxAA, 37, p221.4Lorne W. Avery & Lewis L. H. , 1968, ApJ, 152, p493.5Lawrence J. C., Peter D. N. & Jeffery D.S. 1972, ApJ, 176, 439.6David A. Neufeld, 1990, ApJ, 350, p216.辐射转移问题对三维辐射转移问题,蒙特卡罗模拟几乎是必须的计算方法1. 将空间划分为Nx*Ny*Nz个网格2. 选取空间中的随机点产生光子包(photon package), 并此向随机方向传播3. 考虑光子与粒子的随机碰撞产生具有某种分布的随机数)(xfy dxxfxF)()()(1xfFdxdxFdF即正比于即分布密度试验粒子数值模拟研究粒子在电磁场中的运动研究粒子加速等全粒子试验粒