基于蒙特卡洛思想的定积分数值解法.doc

精品海 南 大 学数理统计课程设计题目： 基于蒙特卡洛思想的定积分数值解法 班级： 信息与计算科学 姓名： 体贴的瑾色 学号： 指导教师： 日期： 2017.06 目录基于蒙特卡洛方法的定积分数值解法3摘要3Abstract3一、前言4二、蒙特卡洛求解定积分法42.1 随机投点法42.2平均值法62.3两种方法的比较7三、蒙特卡洛计算二重积分8四、结语10参考文献:10附录：MATLAB程序10基于蒙特卡洛方法的定积分数值解法 摘要微积分是现代数学和现代物理的一个重要基础，而定积分在生活生产上也具有十分广泛的应用，然而，定积分的计算特别是涉及到较复杂的定积分的话便变的十分困难.本文介绍了蒙特卡洛方法解决定积分求解的两种思路：随机取点法和平均值法，并给出了用matlab实现两种方法的算法程序，最后用两种方法分别计算了几个实例.关键词：蒙特卡洛 定积分 matlab AbstractCalculus is an important foundation of modern mathematics and modern physics, and definite integral are also widely used in life production. However, it is very difficult to calculate the definite integral, especially when it comes to more complex. This paper introduces two methods of Monte Carlo method to solve definite integral solution: random point method and mean method, and gives the algorithm to realize two methods by matlab. Finally, we use two methods to calculate several examples respectively.Keys: Monte Carlo definite integral matlab一、前言蒙特卡洛方法（英语：Monte Carlo method），也称统计模拟方法，是1940年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法.是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法.蒙特卡洛方法是属于统计实验的一种方法，主要通过统计抽样实验为各种各样的数学问题提供近似解，所以也被称为随机抽样技术.而正是因为蒙特卡洛方法的原理使得这种方法得到的结果基本上都会存在误差，不过可喜的是这种误差随着取的随机数的数量的增大而减小.所以，随着计算机的发展，蒙特卡洛算法的精确度一直在提高.定积分的求解是蒙特卡洛方法解决的最好的问题之一，特别是在涉及到重积分的情况下，这种方法往往是科研工作者比较喜欢的方法之一.二、蒙特卡洛求解定积分法2.1 随机投点法 伯努利大数定理说明：随着试验次数的增大，事件发生的频率与其频率的偏差大于预定给定的精度的可能性愈来愈小，要多小就多小，这就是说当试验次数足够多的话，频率稳定于概率.所以，随机投点法的方法的原理就是重复进行大量实验，然后用观测到的频率代替概率作为所求结果.比如求定积分，其中，可设二维随机变量服从正方形上的均匀分布，则可知服从上的均匀分布，也服从上的均匀分布且和相互独立，又记事件 则的概率为 我们只需要找到事件出现的频率即可，即将看成是向正方形内的随机投的点，用随机点在区域中的频率作为定积分的值.算法流程是：(1).先用计算机产生上均匀分布的个随机数，组成对随机数 其中应该足够大.(2).对这对随机数记录满足的次数，由此可得到事件发生的频率，则.注意对于一般区间上的定积分，其中为可积函数，文献1给出一种解法是做线性变换将转换成的形式.本文给出另外一种算法来计算，我们知道一重定积分的几何意义是求曲线和坐标轴所包围的面积，但是这里面临一个问题就是包围的面积可负可正，另外此时的总区域面积也不再是1，所以如果按照上述流程(2)显然会出现问题，所以，本文对(2)进行改良后的新的第二步是：(2) 记，记录满足出现的次数，和满足的次数，算得新的频率为.若，则，或者 则或者，则.最后得到.譬如计算,其精确值由matlab函数integral给出（下同），取，运行100次后得到的结果为 模拟后得到的值的分析为（误差计算公式为,下同）：方法平均值方差误差精确解随机投点法0.65015 3.733E-04 3.830E-040.64984当时的投点图为2.2平均值法计算定积分 ，其中可积.设随机变量服从上的均匀分布，则的数学期望为所以估计期望就是估计所求值，由辛钦大数定理，可以用观察值的平均去估计的期望值.具体做法如下：先用计算机产生个在上均匀分布的随机数，然后对每个计算，最后得到的估计值为 .至于一般区间上定积分的计算，按照上述流程因为,所以.譬如计算取，运行100次后得到的结果为模拟后得到的值的分析为方法平均值方差误差精确解平均值法1.912942.378E-032.381E-031.913122.3两种方法的比较本为将通过计算几个例子来比较两种方法的优劣性.比如计算综合的的比较为： 方法平均值方差误差精确解随机投点法18.250364.732E-01 4.732E-01 18.25077平均值法18.253138.831E-02 8.886E-02 18.25077又比如计算结果为：综合比较结果为： 方法平均值方差误差精确解随机投点法 9.164791.898E-01 1.901E-01 9.16667平均值法 9.165898.122E-02 8.128E-02 9.16667显然在一重定积分的计算上，虽然两种方法平均值都接近于精确解，但是不管是平均值法方差还是误差都较之随机投点法更小，所以平均值法性能更好一些.三、蒙特卡洛计算二重积分计算二重积分其实和计算一重积分差距并不是很大。在随机投点法中与只需要将生成的随机数变成三维并把投点区域换成一个正方体即可，比如求区间在上的二重积分，计算公式应该为其中与上文含义一致.而平均值法除了生成的随机数也要变成三维之外，最大的不同在于计算一般积分区域上的积分，计算公式应该变成.另外值得一提的是，因为matlab在最优化上效果并不是十分的出色 ，所以本文在随机投点法中为了求投点区域，本文采用lingo求函数最大值和最小值.譬如求 取，使用两种方法得到的结果为方法平均值方差误差精确解随机投点法5.20314 4.70790E-044.71132E-04 5.20320平均值法5.20321 2.14826E-052.15035E-05 5.20320又比如求积分 模拟后得到的结果为：方法平均值方差误差精确解随机投点法 5.20836 3.64836E-013.67495E-01 5.20320平均值法5.20268 2.12750E-022.13024E-02 5.20320四、结语通过上文给出的结果，可以看到蒙特卡洛方法用来求定积分的误差十分的小，每次得到的结果也十分的稳定，所以是求解定积分的一种很好的方法.但是，显然不管在计算一重还是二重定积分上，平均值法还是比随机投点法性能要更好一些.参考文献:1茆诗松,程依明.概率论与数理统计M. 高等教育出版社.20112朱陆陆. 蒙特卡洛方法及应用D.华中师范大学,2014.附录：MATLAB程序感谢下载载function gailv(time,num,dlim,ulim)%time 次数%num 生成的点数%min 最大值%max 最小值%随机投点法fun=(x)exp(x.2)-x;gun=(x)-(exp(x.2)-x);val=integral(fun,dlim,ulim)I=zeros(time,1);,min=fminbnd(fun,dlim,ulim);,max=fminbnd(gun,dlim,ulim);max=-max;for i=1:time i1=0;i2=0;xi=unifrnd(dlim,ulim,num,1);if max0 yi=max*rand(num,1); i1=1; else yi=min+(max-min)*rand(num,1); endendy=fun(xi);k=sum(yi0&y0);p=sum(yi=y&yi0&y0);if i=100 a=(yi0&y0);b=(yi=y&yi0&y0);c=(a+b); hold on plot(xi(a),yi(a),ro) plot(xi(b),yi(b),bo) plot(xi(c),yi(c),co)% fplot(exp(-x.3)-x),dlim,ulim) hold offendI(i,1)=(k-p)/num*(max-min)+i1*min-i2*max)*(ulim-dlim)*(ulim-dlim); ave_r=sum(I)/time;wucha_r=sum(I-val).2);var_r=sum(I-ave_r).2);endfigure plot(1:time,I,bo);%平均值法for i=1:time xi=dlim+(ulim-dlim)*rand(num,1);y=fun(xi);I(i,1)=sum(y)/num;endhold onplot(1:time,zeros(1,time)+val,c)title(随机投点法计算定积分结果);plot(1:time,I,ro);hold offave_a=sum(I)/time;var_a=sum(I-ave_a).2);wucha_a=sum(I-val).2);fprintf(方法ttt平均值ttt方差tttt误差n)fprintf(随机投点法t%1.5ftt%1.5Ett%1.5En,ave_r,var_r wucha_r);fprintf(平均值法tt%1.5ftt%1.5Ett%1.5En,ave_a,var_a wucha_a);二重定积分程序:function gailv2(time,num,dlim,ulim,ydlim,yulim)%time 次数%num 生成的点数%min 最大值%max 最小值%随机投点法fun=(x,y)(-x.2+y+5);val=dblquad(fun,dlim,ulim,ydlim,yulim);I=zeros(time,1);max=6.3;min=4.06;for i=1:time i1=0;i2=0;xi=dlim+(ulim-dlim)*rand(num,1);yi=ydlim+(yulim-ydlim)*rand(num,1);if max0 zi=max*rand(num,1); i1=1; else zi=min+(max-min)*rand(num,1); endendy=fun(xi,yi);k=sum(zi0&y=0);p=sum(zi=y&zi0&y0); plot3(xi,yi,y,b); hold off endI(i,1)=(k-p)/num*(max-min)+i1*min-i2*max)*(ulim-dlim)*(yulim-ydlim); endave_r=sum(I)/time;wucha_r=sum(I-val).2);var_r=sum(I-ave_r).2);figureplot(1:time,I,bo);%平均值法for i=1:time xi=dlim+(ulim-dlim)*rand(num,1);yi=unifrnd(ydlim,yulim,num,1);% zi=min+(max-min)*rand(num,1);y=fun(xi,yi);I(i,1)=sum(y)/num*(ulim-dlim)*(yulim-ydlim);endhold onplot(1:time,zeros(1,time)+val,b)plot(1:time,I,r+);title(两种方法的比较结果);legend(随机投点法,精确解,平均值法);xlabel(次数);