数学建模-蒙特卡罗方法

2017美赛培训：蒙特卡罗方法主讲人：老教练报告人：陈蓉蓉陈雨路杨园园什么叫蒙特卡罗方法？•蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。基本思想•当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率，数学期望有关的量时，通过某种试验的方法，得出该事件发生的概率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。•当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。蒙特卡罗方法的特点•优点：•1、能够比较逼真地描述具有随机性质的事物的特点及物理实验过程•2、受几何条件限制小•3、收敛速度与问题的维数无关•4、具有同时计算多个方案与多个未知量的能力•5、误差容易确定•6、程序结构简单，易于实现•缺点：•1收敛速度慢•2误差具有概率性•3在粒子输运问题中，计算结果与系统大小有关•所以在使用蒙特卡罗方法时，要“扬长避短”，只对问题中难以用解析（或数值）方法处理的部分，使用蒙特卡罗方法计算，对那些能用解析（或数值）方法处理的部分，应当尽量使用解析方法主要应用范围•粒子输运问题（实验物理，反应堆物理，高能物理）•统计物理•典型数学问题•真空技术•激光技术以及医学•生物•探矿什么是随机数？•在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称为随机数序列，其中每一个体称为随机数•符号：•两个特点：独立性，均匀性21，产生随机数•随机数表方法•物理方法随机数表•随机数表是由0,1,2,3,4,5,6,7,8,9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立，这些数字序列叫作随机数字序列。•（如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为763425891...，要想得到三位有效数字的随机数一次为0.763，0.425，0.891...）物理方法•利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。•作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。•一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：•其中或者为0，或者为1。因此，利用物理方法在计算机产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5),,1(imim222m-22-11缺点•随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。•物理方法产生的随机数序列无法重复实现，不能进行程序复算。给验证结果带来很大困难。而且增加随机数发生器和电路联接等附加设备，费用昂贵。因此该方法也不适合在计算机上使用。伪随机数•用递推公式产生随机数序列。,2,1,,,11nnTknnnkk21，，，伪随机数存在的两个问题•递推公式和初始值确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。•由于随机数序列是由递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现重复。随机数序列出现周期性的循环现象。解决方案•第一个问题：不能从本质上加以改变，但只要递推公式选的比较好，随机数间的相互独立性是可以近似满足的。•第二个问题：因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。应用：蒙特卡罗方法计算积分可以通俗地说，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量ｇ(r)的数学期望通过某种试验，得到Ｎ观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f(r)中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值•作为积分的估计值（近似值）。0)()(drrfrggNiiNrgNg1)(1•求积分（2.1）ba)(fdxx蒙特卡罗方法步骤如下：•1、在区间【a,b】上利用计算机均匀产生n个随机数x1，x2·····xn，这个可以在MATLAB软件中用unifrnd命令实现。•2、计算每一个随机数相应的被积函数值f（x1），f（x2）····f（xn）。•3、计算被积函数值的平均值•4、所以2.1式的值≈n)(1niixfn)(1niixf简单定积分例子：用蒙特卡罗方法求•首先我们进行严格的数学计算，便于后面与蒙特卡洛计算方法所得结果形成对比：已知的原函数是，那么定积分值就是：•我们可以在Matlab中输入以下代码进行精确计算：exp（2）-exp（0），这个值是此定积分的真实值.•下面进行蒙特卡洛计算上述定积分，其MATLAB代码如下：•N=500;•x=unifrnd(0,2,N,1);•mean(2*exp(3*x.^2))dx203x2e389056099.6e-e02xexe解释语句：•上述三条语句完整实现了蒙特卡洛计算上述定积分步骤.第一条语句是设定了停止条件，共做N次MonteCarlo模拟.第二条语句实现了在积分区间上均匀产生N个随机数.第三条语句实现蒙特卡洛计算方法的面积逼近.对N设置不同的值，观察所得蒙特卡洛计算方法定积分值，如表1所示，我们可以发现：当不断增大N值时，所得结果越来越接近真实值.蒙特卡罗计算方法性质：•蒙特卡洛计算方法依据概率统计理论，具有统计特性，主要表现在以下三个方面.•1.蒙特卡洛计算方法具有随机性、不确定性.即每次运行结果都会不一样，因为计算机产生随机数并不是可以重现的.•2.蒙特卡洛计算方法具有统计稳定性。虽然每次运行产生随机数是不一样的，但是随机数的概率分布是一样的，所以蒙特卡洛计算方法可以满足概率统计的稳定性.•3.随着随机数数量增加，蒙特卡洛计算方法所得结果会更加逼近真实值，这就是我们所讲的依概率收敛到真实值意义.蒙特卡罗方法实验面积、体积计算问题冰淇淋锥的体积计算思考题与练习题蒙特卡罗方法——随机投点试验求近似解引例.给定曲线y=2–x2和曲线y3=x2，曲线的交点为：P1(–1，1)、P2(1，1)。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积。P=rand(10000,2);x=2*P(:,1)-1;y=2*P(:,2);II=find(y=2-x.^2&y.^3=x.^2);M=length(II);S=4*M/10000plot(x(II),y(II),'g.')S=2.11361532)2(11322dxxxS例5.14计算其中D为y=x–2与y2=x所围D的边界曲线交点为：(–1，1)，(4，2)，被积函数在求积区域内的最大值为16。积分值是三维体积，该三维图形位于立方体区域0≤x≤4，–1≤y≤2，0≤z≤16内，立方体区域的体积为192。data=rand(10000,3);x=4*data(:,1);y=-1+3*data(:,2);z=16*data(:,3);II=find(x=y.^2&x=y+2&z=x.*(y.^2));M=length(II);V=192*M/10000Ddxdyxy2例5.15用蒙特卡罗方法计算dxdydzzyx)(222其中,积分区域是由和z=1所围成。22yxz被积函数在求积区域上的最大值为2。所以有四维超立方体–1≤x≤1，–1≤y≤1，0≤z≤1，0≤u≤2P=rand(10000,4);x=-1+2*P(:,1);y=-1+2*P(:,2);z=P(:,3);u=2*P(:,4);II=find(zsqrt(x.^2+y.^2)&z=1&u=x.^2+y.^2+z.^2);M=length(II);V=8*M/10000实验:蒙特卡罗方法计算体积22yxz2211yxz&}20,11,11|),,{(zyxzyxx=2*rand-1产生–1到1之间的随机数y=2*rand-1产生–1到1之间的随机数z=2*rand;产生0到2之间的随机数冰淇淋锥含于体积=8的六面体NmV822由于rand产生0到1之间的随机数,所以N个点均匀分布于六面体中,锥体中占有m个,则锥体与六面体体积之比近似为m:Nfunction[q,error]=MonteC(L)ifnargin==0,L=7;endN=10000;fork=1:LP=rand(N,3);x=2*P(:,1)-1;y=2*P(:,2)-1;z=2*P(:,3);R2=x.^2+y.^2;R=sqrt(R2);II=find(z=R&z=1+sqrt(1-R2));m=length(II);q(k)=8*m/N;enderror=q-pi;实验参考程序蒙特卡罗方法计算体积半球体积332R圆锥体积hR231实验任务一：记录L次实验的实验数据及误差实验任务二：修改实验程序MonteC计算L次实验数据均值及均值误差(mean计算平均值)序号1234567数据误差L8163264128256均值误差functionicecream(m,n)ifnargin==0,m=20;n=100;endt=linspace(0,2*pi,n);r=linspace(0,1,m);x=r'*cos(t);y=r'*sin(t);z1=sqrt(x.^2+y.^2);z2=1+sqrt(1+eps-x.^2-y.^2);X=[x;x];Y=[y;y];Z=[z1;z2];mesh(X,Y,Z)view(0,-18)colormap([001]),axisoff冰淇淋锥体积3R冰淇淋锥图形绘制程序5.下面程序绘出二维图形填充图(右图)。分析每条语句功能,给程序中语句写注记x1=-1:0.1:1;y1=x1.^2.^(1/3);x2=1:-0.1:-1;y2=2-x2.^2;fill([x1,x2],[y1,y2],'g')axisoff