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 b
5、e the remainder 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,
6、2,n,试产生1000个在椭圆x2+4y2=4内均匀分布的随机点：方法：在 中均匀产生足够 多随机点，当位于椭圆内的点数为1000时停止。,可由随机数序列 r (0,1)构造一系列随机数序列,有时随机数产生函数通不过严格的随机性测试。而在一些计算（如MC积分）中随机性非常重要。故使用更长字节的计算机更好。,应用之一,估计面积和体积,Numerical integration,选取(0, 1)中随机数序列x1， x2， x3， xn。则,误差约 ，它并不能和一些高级的数值积分算法比拟，但对多维情况，MC方法却很有吸引力。,我们可产生一系列随机数可简单取3个随机数构成一个随机点，即,相应地，,一般
7、地，,例如： 求,其中,在该区域中取5000个点，可得该积分约0.57.,计算体积,Pseudo-code,Program volume_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 program,计算体积之作业：冰淇淋的体积,应用之
8、二,MC模拟,生日问题,假设有n个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50的几率人群中至少有2个人生日相同。,n 理论几率 模拟几率0.117 0.1100.411 0.4120.507 0.5200.706 0.6920.941 0.93650 0.986 0.987,Real function prob(n,npts)Integer i, nptsLogical birthReal sumSum=0For I=1 to npts doIf birth(n) then sum sum+1EndforProb sum/real(npts)E
9、nd 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(n) then sum sum+1EndforProb sum/real(npts)End function,不妨取d l。产生随机数u和，其中u为针中心距最近的线的距离，故u小于等于0.5， 为夹角，根据对称性，可取0到/2
10、。,算时需要用上是否不自洽（姜冰）,作业,中子屏蔽问题Neutron Shielding problem,假设铅墙长为5，中子在铅中的平均自由程为1，中子与铅原子碰撞后各向同性散射。令碰撞8次后中子能量耗尽，试求穿透铅墙的中子的比例。,暂不考虑垂直纸面的运动，则中子的水平位移是 。,入口,铅墙（长为5）,21点问题Blackjack 21 problem,参数拟合问题Parameter fitting problem,辐射转移问题Radiative transfer problem,1 Ahn,S.-H., Lee, H.-W. & Lee, H.-M. 2000, JKAS, 33, p29
11、.2 Zheng Zheng & Jordi Miralda-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.,辐射转移问题,对三维辐射转移问题,蒙特卡罗模拟几乎是必须的,计算
12、方法,将空间划分为Nx*Ny*Nz个网格选取空间中的随机点产生光子包(photon package), 并此向随机方向传播考虑光子与粒子的随机碰撞,产生具有某种分布的随机数,试验粒子数值模拟,研究粒子在电磁场中的运动,研究粒子加速等,全粒子,试验粒子,研究微观不稳定性、反常电阻等等,忽略粒子间相互作用及对电磁场的反馈,也可考虑库仑碰撞。,算法 （常微分）,蛙 跳,龙格库塔,保证了E=0时满足能量守恒。,CUEBIT,如Mori et al. 1998, ApJ, 494, 430,如Liu, W. J. et al. 2007, Adv. Space Res., in press,Hamilton et al. 2003, SP, 214, 339,模拟结果,
