1、1蒙 特 卡 罗 随 机 模 拟 投 点 法 在 数 字 积 分 中 的应 用数 学 与 应 用 数 学 0901 班 ： 张 瑞 宸指 导 老 师 ： 任 明 慧摘 要 ：本 文 首 先 介 绍 了 蒙 特 卡 罗 方 法 的 产 生 和 发 展 ，然 后 分 析 了 蒙 特 卡 罗 方法 计 算 数 值 积 分 的 理 论 原 理 ，最 后 给 出 了 蒙 特 卡 罗 方 法 计 算 数 值 积 分 的MATLAB 编 程 实 现 ，全 文 主 要 是 讨 论 了 蒙 特 卡 罗 方 法 在 定 积 分 计 算 的 应 用 。而 蒙 特 卡 罗 的 优 点 ：可以计算被积函数非常复杂的定积
2、分、重积分，并且维数没有限制，这是别的数值积分方法还未达到的。 蒙 特 卡 罗 的 缺 点 ：收 敛 速 度 慢 ，误差 一 般 较 大 ，且 是 概 率 的 误 差 ，不 是 真 正 的 误 差 。关 键 词 ：蒙 特 卡 罗 方 法 ，均 值 估 计 法 ，数 值 积 分 ，Matlab 编 程 Abstract：This paper first introduces the emergence and development of the Monte Carlo method, and then analyze the theoretical principles of Monte Ca
3、rlo numerical integration method, Full-text mainly discussed the application of the Monte Carlo method in the definite integral. The advantages of Monte Carlo: can be calculated the integrable functions very complex definite integral, Multiple integrals, and dimension no limit, other numerical integ
4、ration methods have not yet reached. Monte Carlo Disadvantages: slow convergence speed, error generally higher, and the probability of error, not a real error.Keywords: Monte Carlo method，Mean estimation method，numerical integral，Matlab programming 0 引言历史上有记载的蒙特卡罗试验始于十八世纪末期（约 1777 年） ，当时布丰（Buffon）为了
5、计算圆周率，设计了一个“投针试验” ，后文会给出。虽然方法已经存在了 200 多年，此方法命名为蒙特卡罗则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用蒙特卡罗方法模拟中子对某种特殊材料的穿透作用。出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一 von Neumann 灵犀一点选择摩洛哥著名赌城 蒙特卡罗（Monte 蒙特卡罗随机模拟投点法在数值积分中的应用2Carlo） 作为该项目名称，自此这种方法也就被命名为 Monte Carlo 方法广为流传。蒙特卡罗方法，又名随机模拟法或统计实验法它是以概率统计理论为基础，依据大数定律（样本均值替代总体均值）利用电子计算机
6、数字模拟技术，解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的求出数值解的方法。一般蒙特卡罗法在数学中最常见的应用就是蒙特卡罗随机投点和蒙特卡罗数值积分。1 蒙特卡罗方法的产生与发展蒙特卡罗方法是在二战期间产生和发展起来的他的奠基者是美籍匈牙利人数学家冯诺伊曼（J.Von Neuman
7、n 1903-1957）由于通常计算量相当大而电子计算机在当时还没有出现，所有运算只能用手工进行，故而相当长的时间里蒙特卡罗方法难以推广。1.1 蒙特卡罗方法的产生作为蒙特卡罗方法的最初应用，是解决蒲丰氏问题1777 年，法国数学家Buffon 提出利用投针实验求解的问题：设平面上有无数多条距离为 1 的等距平行线，现向该平面随机的投掷一根长度为 的针。随机投针是指针的中心点于最近的平行线间的距离 均匀l x分布在0,1/2上，针与平行线的夹角 （不管相交与否）均匀分布在 上，0,如图一所示，故而，故其 ， 概率密度函数分别为1(0,)2xU(0,)。1(),()2px故我们得到针与线相交的充
8、要条件是 2sinlx蒙特卡罗随机模拟投点法在数值积分中的应用3图1.1 针与线相交的几种情况则针与线相交的概率是 (sin)2lpx= si sin200(l lpdx= 2.sinlld所以得到圆周率 。2(si)llpxp假如我们能做大量的投针实验并记录下针与线的相交次数，则可以根据大数定律估计出针线相交的概率 P。投针实验 N 次可能有 n 次使针与任意平行线相交, 那么 ，显然，试验次数 N 越多，P 的近似程度越好。历史上曾有npN几位学者做过这样的投针试验，并用手工计算出 值，结果参见表1。表 1实验者 时间（年） 针长 l投针次数 相交次数 估计值Wolf 1850 0.80
9、5000 2532 3.15956Smith 1855 0.60 3024 1218 3.15665Fox 1884 0.75 1030 489 3.15951Lazzarini 1925 0.83 3408 1808 3.14159292应该 指出，上述试验的精度 一般不会很高。譬如，假 设21,pl则。则由中心极限定理， 如果试验次数为 N，则 p 的估计值 渐进服从 ，近似为 。因(1),pN0.231(.6,)N蒙特卡罗随机模拟投点法在数值积分中的应用4此，若要以95%的概率保证 p 的精确到三位有效数字，即 与 p 的差距小于0.001，则 N 必须满足 N 。2251.960.3/
10、.18.90重复进行上千次的投针实验和手工计算，要消耗大量的人力、财力，蒙特卡罗方法虽然能解决此类问题，但得不到推广应用。1.2 蒙特卡罗方法的发展20世纪40年代以后，随着电子计算机的出现和发展，人们有可能用计算机来模拟这类实验和计算。计算机具有计算速度高和存储容量大的特点，采用数字模拟技术可以代替许多实际上非常庞大而复杂的实验，并迅速将实验结果进行运算处理，于是 Monte Carlo 方法重新被提起，引起世人重视，应用日渐广泛。实际上，采用 Monte Carlo 方法在计算机上建立模型来解决 Buffon 问题是非常简单的。我们在计算机上进行模拟试验，给定 ，我们可以在计算机上随机产生
11、 和l x，然后判断 是否成立。若成立，则针线相交，否则不交。假如我们2sinlx在计算机上独立的产生 N 对这样的 和 ，并记录下 成立的次数,记为x2sinlxn。则 估计值可取为 ，这就是随机模拟计算的结果。2l借助 Matlab 在计算机上进行投针实验， （相应程序见附录1）所取得的计算结果参见表2。表 2N =0.4l=0.6l=0.8l=0.9l=1.0l50000 3.1252442 3.1583934 3.1392246 3.14642309 3.1464351100000 3.1241457 3.1211978 3.1269543 3.1324504 3.1286664200
12、000 3.1414435 3.1396353 3.1373472 3.1366262 3.1371811500000 3.1519392 3.1470787 3.1455443 3.1422826 3.1425492800000 3.1483358 3.1465609 3.1440825 3.1423004 3.14302941000000 3.1435296 3.1459072 3.1412646 3.1419370 3.1423368蒙特卡罗随机模拟投点法在数值积分中的应用52000000 3.1497197 3.1462907 3.1418722 3.1410570 3.1407306
13、5000000 3.1393330 3.1398735 3.1410784 3.1416836 3.1412106其中 = 3.1416。不难看出，利用计算机运算不仅结果精确，而且迅速。随着电子计算机的普及，蒙特卡罗方法作为一种独到的方法得到开发，并首先应用到核武器的试验与研制中。尤其是各种可是编程方法的不断涌现，更显示出蒙特卡罗方法的最独到的优点，即形象直观地用数学方法在电子计算机上实现数字模拟实验。2 蒙特卡罗方法计算数值积分的理论原理随机模拟是一种随机试验的方法，也称为蒙特卡罗方法（Monte Carlo）。这种方法利用随机试验，根据频率与概率、平均值与期望值等之间的关系，1推断出预期的
14、结果。2.1 蒙特卡罗算法的理论原理随机模拟法的基本原理非常简单，下面给以直观的说明。首先构造一个定积分 ，其中积分 （即 1/41024dx102dxS单位圆的面积）用随机模拟方法求得。图 2.1 中粗线是 1/4 单位圆，如果向图2.1 中边长为 1 的正方形里随机投 块小石头，n当 很大时小石头会大致均匀地分布在正方形n中，数一下落在 1/4 单位圆内的小石头，假定有 个，那么 / 就可以看作 1/4 单位圆面积kkn/4 的近似值。小石头的位置坐标可以用产生均匀分布随机数的程序得到，记为是0,1区间均匀分布随机数（ 相互独立） ，记录满足),.21(,niyxi iyx,的数量 ，即得
15、 =4 / 。2iikkn我们用概率论中的大数定律来说明这个直观认识的原理。图2.1 xy0 11蒙特卡罗随机模拟投点法在数值积分中的应用6大数定律（伯努利（Bernoull）定理） 设 是 次独立重复试验中事2kn件 A 发生的次数， 是事件 A 每次试验中发生的概率，则对任意的正数 ，有p 1limpnkPn（1）若规定“向图中正方形随机投一块小石头落在四分之一单位圆里”为事件A 发生，则 A 发生的概率 应该等于四分之一单位圆面积，随机投 块石头就p n是独立重复做 次试验，事件 A 发生 次，由（1）式， 无限变大时 / 与nkk之差小于任意一个数 的概率趋于 1。p用计算机算一下就会
16、发现，即使 很大结果也不好，并且很不稳定，远不n如常规数值积分的几种方法。实际上，随机模拟法很少用来做定积分，它的特点是能够方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。2.2 蒙特卡罗方法计算定积分的理论原理通常有两类办法计算定积分。（1）随机投点法在“投石算面积”的例子中，事件 A 在每次试验中发生的概率 是四分p之一单位圆面积，即（2）dxdyxp102102次试验由计算机完成，采用0,1区间上的均匀分布产生相互独立的随机n数。记这样产生的 个点的坐标为 。事件 A 发生等价于n),.21),(niyxi满足 ， A 发生的次数是满足 点的个),(iyx21i
17、ix 2iix),.1(n数 。由伯努利定理， 可以用 / 近似替代。kpkn这种方法可以推广如下。对任意区间 内的连续函数 ，满足,ba)(f，为计算定积分 ，先由计算机产生 个点的坐标dxfc)(,c)0adxf)(n，其中 分别为 和 区间上的均匀分布随机数，,iy.21niyx,c蒙特卡罗随机模拟投点法在数值积分中的应用7然后记录 个点中满足 的数目 ，则n)(iixfykcabncdabba )()(（3）（2）均值估计法这种方法依据概率论的一下两个定理：大数定理（辛钦定理） 若随机变量 相互独立，服从同一2 nY,.21个分布，且具有数学期望 ,则对任意的正数 ，有iEY).,1(
18、n（4）lm1iinP随机变量函数的期望 若随机变量 的概率分布密度是2 X，则随机变量的函数 的期望为)(bxap)(XfYdxpfEba)(（5）于是，当 为 区间均匀分布的随机变量时，X,ba )/(1)ab， （ 5）式给出)(bxa)()(xfEabdxfba（6）只要产生 区间相互独立、均匀分布的随机数 ， 就,ba i),.21(n)(iixfy相互独立。根据辛钦定理，当 很大时期望 就可以用 的平nXfEYi均值近似，所以由（6）式得到niiba xfadxf1)()()(（7）与随机投点法相比，均值估计法没有 的限制，只需dxfc)(,c)0蒙特卡罗随机模拟投点法在数值积分中
19、的应用8计算 ，不需要产生随机数 ，也不需要作 的比较，显然大为方)(ixf iyiy)(ixf便。2.3 蒙特卡罗方法计算重积分的理论原理用随机模拟作数值积分的优点不仅在于计算简单，尤其是它可以方便地推广到计算多重积分，而不少多重积分是其他方法很难或者根本无法计算的。而通过上文对蒙特卡罗计算定积分的两种方法的比较，我们可知均值估计法比投点法更好操作，经过计算还可知均值估计法精确度（若选择的密度函数适合）会更高。重积分本文讨论均值估计法。考虑 维积分：k kkdxxfI .),.(. 2121其中 为 维积分区域， ( ， )。kViiiba:,.选取 上的一个概率密度函数 ， ( ，),.2
20、1kxpViiibxa:)，且 满足 ，ki,.21),.(21kxp 1,(. 21 kkd，故：,.(21kx kkk dxxpxpfI .),.(),.(. 212121其中, ),.(21kxfE ),.(21k其中 , mi kxpf12),.( ),.(21kx（8）这里 是落在 区域中的随机投点数。m对于二重积分： dxyfI),(其中 ， 。记 ，设 ，bxa: gyxc)()(21 ,dcbaVS分别是 ， 的测度，若选取 上的均匀分布，则VS是,蒙特卡罗随机模拟投点法在数值积分中的应用9。由几何概率得 ，即 ，于是：y)(x, 1),(Syxg nmSVnSVdxygyxf
21、I),(,mif1),(iyxfS1,miVyxfnS1),(综上可得， （9）miVfnI1),(3 蒙特卡罗方法 MATLAB 编程实例用随机模拟作数值积分主要利用均匀分布产生的随机数及相应的判断、计算和简单运算。MATLAB 提供的命令unifrnd(a,b,m,n)产生 行 列 区间上均匀分布随机数。当 ， 时，可用 rand(m,n)。mn,ba0a1b例 1 用蒙特卡罗方法计算 。解(1) 随机投点法。 MATLAB 程序如下：n=10000;x=rand(2,n);k=0;for i=1:nif x(1,i)2+x(2,i)2=feval(fai1,x(k),y(k)=feval
22、(fai1,x),ytic;s1=mtc(inline(exp(-x.2/2).*sin(x.2+y),inline(-sqrt(1-x.2/2),inline(sqrt(1-x.2/2),-0.5,1,-1,1,100000),tocs1=0.409 010 343 471 83Elapsed time is 40.094 000 seconds.tic;s2=mtcx(inline(exp(-x.2/2).*sin(x.2+y),inline(-sqrt(1-x.2/2),inline(sqrt(1-x.2/2),-0.5,1,-1,1,100000),tocs2=0.411 839 62
23、5 630 18Elapsed time is 0.047 000 seconds.在运用了向量化编程后运行时间从40.094s 缩短到 0.047s，是原来的1/853，效果好；误差为0.00 008 992 054 612 000 366 ，也比前者好。若将上面程序稍作如下修改，即输入参数减少2个，那么当 =1 000 000时，n运行时间大约增加0.45s，精度不变。function s=mtcxg(f,fai1,fai2,a,b,c,d,n)if nargin=feval(fai1,x),ytic;s3=mtcxg(inline(exp(-x.2/2).*sin(x.2+y),inli
24、ne(-sqrt(1-x.2/2),inline(sqrt(1-x.2/2),-0.5,1,-1,1,100000),tocs3=0.412 005 502 840 08Elapsed time is 0.922 000 seconds.对比常规数值积分和蒙特卡罗方法数值积分，同样数量的 值也就n意味这几乎相同的计算量常规数值积分结果的精确度要高于蒙特卡罗数值积分的结果。那么，我们为何还需要用蒙特卡罗来算数值积分呢？答案的关键在于，常规数值积分的精度直接取决于每个维度上取点数量，维度增加了，但是每个维度上要取的点却不能减少。在多重积分中，随着被积函数维度增加，需要计算的函数值数量以指数速度递增
25、。例如在一重积分 中，只要沿着 x 轴取 个点；要达到相同abdxf)(n大小的精确度，在 重积分 中，仍然需要k kk dxf .).() 212121在每个维度上取 个点， 个纬度的坐标相组合，共需要计算 个坐标对应n的 函数值。取点越多，会占用计算机大量内存，也需要更长运算时间，最)(xf终导致这种计算方法不可行！ 蒙特卡罗方法却不同，不管是积分有多少重，取 个点计算的结果精确度都n差不多。因此，即使在一重积分的情形下，蒙特卡罗方法的效率比不过常规数值积分，但随着积分维度增加，常规数值积分的速度呈指数下降，蒙特卡罗方法的效率却基本不变。经验表明，当积分重数达到 4 重积分甚至更高时，蒙特
26、卡罗方法将远远优于常规数值积分方法。4 总结归纳起来本论文主要论述了：蒙特卡罗算法在数值积分中的运用。本篇论文给出了蒙特卡罗方法应用在数值积分中的原理，讨论蒙特卡罗算法运用在数蒙特卡罗随机模拟投点法在数值积分中的应用13值积分中的两种计算方法，随机投点法和均值估计法，具体蒙特卡罗算法的实现给出了 MATLAB 编程的实例。最后对蒙特卡罗算法的优缺点进行了讨论。蒙特卡罗理论的应用非常广泛,两种蒙特卡罗数值计算方法的应用也很全面,由于作者的知识水平,研究能力有限,蒙特卡罗理论在各行各业中都有广泛使用,思想、理论等方面仍存在欠缺之处,但这些都是以后的努力研究的方向,有侍继续发掘。 参考文献：1曲双石
27、，王会娟.MonteCarlo 方法及其应用J.统计教育，2009.1.2茆诗松，程依明，濮 晓龙.概率论与数理统计教程M. 高等教育出版社，2004 年 7 月.3高雷阜. Monte-Carlo 理论与优化方法的研究J.辽宁工程技术大学学报( 自然科学版).2002 年 03 期.4李世刚,刘 辉,陈标华. 素数的一个特殊性质及其用于伪随机数生成的方法J.北京化工大学学报( 自然科学版 ).2003 年 03 期.5张韵华，奚梅成，陈效群.数值计算方法与计算J-2 版.-北京科学出版社，2006.6孙维君，秦华.Monte Carlo 方法在三重积分中的应用J. 山东理工大学学报（自然科学
28、报），2008 第 1 期.7柴中林， 银俊成.蒙特卡罗方法计算定积分的进一步讨论J.应用数学与计算数学学报.2008.第 1 期.8张韵华，奚梅成，陈效群.数值计算方法与算法M. 科学出版社 2006.第 2版.9 JONATHANJ.FORSTER ,JOHNW.M CDONALD and 蒙特卡罗随机模拟投点法在数值积分中的应用14PETERW.F.SMITH.Monte Carlo methods to calculateJ.2005.42,37-46.spring-verlag.致谢:首先要感谢以梅院长等领导对我们的精心统筹安排,科学培养,使我们能够完成所学的知识积累,完成论文。我更
29、要感谢在论文完成过程中一直给与悉心指导的任老师，任老师在指导论文的过程中,耐心的讲解和教诲，使我能掌握写作本文的基础知识并理解掌握参考文献的写作的框架思路,尤其感谢姚芬芳等多位同学，是你们的帮助让我顺利完成论文的完善和校对。再一次对给予我帮助的老师、同学们表示最诚挚的感谢！蒙特卡罗随机模拟投点法在数值积分中的应用15附录附录 1：a=1; % 设置两条平行线之间的距离l=0.6; % 投针的长度counter=0; % 针与平行线相交的次数n=10000000; % 投掷次数x=unifrnd(0,a/2,1,n); %产生 n 个(0,a/2)之间均匀分布的随机数，这里 a/2 是投针的中点到最近的平行线的距离phi=unifrnd(0,pi,1,n); % 产生 n 个(0,pi)之间均匀分布的随机数，这里 pi 是投针到最近的平行线的角度for i=1:nif x(i)l*sin(phi(i)/2 % 只要 x 小于 l*sin(phi(i)/2,则相交counter=counter+1;endendfrequency=counter/n; % 计算相交的频率，即相交次数比总次数Pi=2*l/(a*frequency) % 从相交的频率总求的 pi%运行结果
