蒙特卡罗(Monte Carlo)模拟讲义课件

17.

蒙特卡罗(MonteCarlo)模拟蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。该方法源于美国在WWI后研制原子弹的“曼哈顿计划”。S.M.Ulam(1908-1984)和该计划的主持人之一、数学家冯·诺伊曼（J.vonNeumann）用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法。蒙特卡罗模拟其基本思想很早以前就被人们所发现和利用：17世纪，用事件发生的“频率”来决定事件的“概率”；19世纪人们用投针试验的方法来决定π。高速计算机：用数学方法在计算机上大量模拟这样的试验必要性科技计算及社会中的问题比计算π要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(CurseofDimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。MonteCarlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。另一类形式与MonteCarlo方法相似，但理论基础不同的方法—“拟蒙特卡罗方法”(Quasi-MonteCarlo方法)—近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华—王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为LowDiscrepancySequences)代替MonteCarlo方法中的随机数序列。对某些问题该方法的实际速度一般可比MonteCarlo方法提出高数百倍，并可计算精确度。/gkjqy/rkx/rd2.htm拟蒙特卡罗模拟许多计算机系统都有随机数生成函数F90:callrandom_seedcallrandom_number(a)计算程序产生的随机数不是真正的随机数，它们是确定的，但看上去是随机的，且能通过一些随机性的检验，故常称为伪随机数若采用当前时刻作随机数种子，则也是不可重复的随机数如对32位字长的计算机realprocedurerandom((xi))integerarray(li)nrealarray

(xi)nl0anyintegerthat1<l0<231-1fori=1tondoli=(231-1)除以16807li-1的余数xi

＝li/(231-1)endfor其它算法1.取初值x0(0,1)，letxi

bethefractionalpartof(π+xi-1)52.Letu0=1.Fori=1,2,3,…,letui

betheremainderof(8t-3)ui-1dividedby28,andxi

=ui/2i.Heretcanbeanylargeinteger注意：上述随机数序列均具周期性，如上页random子程序的周期约230。Maplehasacollectionofrandomnumbergenerators.如：With(stats)x

:=uniform(0..1):seq(x(),i=1..10);Matlab:x=rand(10,1)产生10个随机数

x=rand(N)产生元素在(0,1)间随机分布的N*N矩阵

s=rand(‘state’,0)重设该生成函数到初始状态如(b-a)r+a

x

(a,b)integer((n+1)r)x

{0,1,2,……n}试产生1000个在椭圆x2+4y2=4内均匀分布的随机点：方法：在中均匀产生足够多随机点，当位于椭圆内的点数为1000时停止。可由随机数序列r

(0,1)构造一系列随机数序列有时随机数产生函数通不过严格的随机性测试。而在一些计算（如MC积分）中随机性非常重要。故使用更长字节的计算机更好。应用之一估计面积和体积Numericalintegration选取(0,1)中随机数序列x1，x2，x3，……

xn。则误差约，它并不能和一些高级的数值积分算法比拟，但对多维情况，MC方法却很有吸引力。我们可产生一系列随机数可简单取3个随机数构成一个随机点，即相应地，一般地，例如：求其中在该区域中取5000个点，可得该积分约0.57.计算体积Pseudo-codeProgramvolume_regionIntegerparametern5000,iprt1000Realarrayintegeri,mRealvol,x,y,zCallrandom((rij))ForI=1tondoxri1yri2zri3Ifthenmm+1Ifmod(I,iprt)=0thenVolreal(m)/real(I)OutputI,volEndifEndforEndprogram计算体积之作业：冰淇淋的体积应用之二MC模拟生日问题假设有n个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50％的几率人群中至少有2个人生日相同。n理论几率模拟几率0.1170.1100.4110.4120.5070.5200.7060.6920.9410.936500.9860.987Realfunctionprob(n,npts)Integeri,nptsLogicalbirthRealsumSum=0ForI=1tonptsdoIfbirth(n)thensumsum+1EndforProbsum/real(npts)EndfunctionLogicalfunctionbirth(n)……蒲丰的投针问题Buffon’sNeedleproblem蒲丰证明在相距d的一系列平行线上投长为l针。针与线的交点数除以投的次数等于2l/(πd)。Realfunctionprob(n,npts)Integeri,nptsLogicalbirthRealsumSum=0ForI=1tonptsdoIfbirth(n)thensumsum+1EndforProbsum/real(npts)Endfunction不妨取d＝l。产生随机数u和，其中u为针中心距最近的线的距离，故u小于等于0.5，为夹角，根据对称性，可取0到π/2。算时需要用上是否不自洽（姜冰）作业中子屏蔽问题NeutronShieldingproblem假设铅墙长为5，中子在铅中的平均自由程为1，中子与铅原子碰撞后各向同性散射。令碰撞8次后中子能量耗尽，试求穿透铅墙的中子的比例。暂不考虑垂直纸面的运动，则中子的水平位移是。入口铅墙（长为5）21点问题Blackjack21problem牛刀小试参数拟合问题Parameterfittingproblem辐射转移问题Radiativetransferproblem[1]Ahn,S.-H.,Lee,H.-W.&Lee,H.-M.2000,JKAS,33,p29.[2]

ZhengZheng&JordiMiralda-Escude.2002,ApJ,578,p33.[3]Watson,A.M.&Henney,W.J.2001,RMxAA,37,p221.[4]LorneW.Avery&LewisL.H.,1968,ApJ,152,p493.[5]LawrenceJ.C.,PeterD.N.&JefferyD.S.1972,ApJ,176,439.[6]DavidA.Neufeld,1990,ApJ,350,p216.辐射转移问题对三维辐射转移问题,蒙特卡罗模拟几乎是必须的计算方法将空间划分为Nx*Ny*Nz个网格选取空间中的随机点产生光子包(photonpackage),并此向随机方向传播考虑光子与粒子的随机碰撞产生具有某种分布的随机数试验粒子数值模拟研究粒子在电磁场中的运动研究粒子加速等全粒子试验粒子研究微观