2019年12月25日1数学建模专题之蒙特卡罗方法重庆理工大学主讲:肖汉光Email:simenxiao1211@163.com数学建模专题之MonteCarlo方法2019年12月25日2内容提纲引言MonteCarlo模拟基本思想随机数生成函数应用小实例排队论模拟2009年B题病床合理安排蒙特卡罗求解积分2010年A题储油罐的变位识别与罐容表标定MonteCarlo方法预测搜救区域MonteCarlo方法计算零件可靠度MonteCarlo方法求解规划问题数学建模专题之MonteCarlo方法2019年12月25日3引言(Introduction)MonteCarlo方法:蒙特卡罗方法,又称随机模拟方法,属于计算数学的一个分支,它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。亦称统计模拟方法,statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来:•MonteCarlo是摩纳哥(monaco)最大的城市,该城以赌博闻名。Monte-Carlo,Monaco数学建模专题之MonteCarlo方法2019年12月25日4引言(Introduction)MonteCarlo方法的起源:二十世纪四十年代中期,由于科学技术的发展和电子计算机的发明,蒙特卡罗方法作为一种独立的方法被提出来,并首先在核武器的试验与研制中得到了应用。NicholasMetropolis(1915-1999)JohnVonNeumann(1903-1957)数学建模专题之MonteCarlo方法2019年12月25日5引言(Introduction)MonteCarlo方法的应用:物理:核物理,热力学与统计物理,粒子输运问题等数学:多重积分、解微分方程、非线性方程组求解等工程领域:真空技术,水力学,激光技术等经济学领域:期权定价、项目管理、投资风险决策等其他领域:化学、医学,生物,生产管理、系统科学、公用事业等方面。随着科学技术的发展,其应用范围将更加广泛。数学建模专题之MonteCarlo方法2019年12月25日6MonteCarlo方法的基本思想Buffon投针实验1768年,法国数学家ComtedeBuffon利用投针实验估计的值dLp2dL数学建模专题之MonteCarlo方法2019年12月25日7SolutionThepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:),[),[dATherandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.TheprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegraldldAdpld20sin01数学建模专题之MonteCarlo方法2019年12月25日8例1.蒲丰投针问题利用关系式:求出π值其中N为投计次数,n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。alP2)(22nNalaPl数学建模专题之MonteCarlo方法2019年12月25日9一些人进行了实验,其结果列于下表:实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929数学建模专题之MonteCarlo方法2019年12月25日10基本思想由上面的例子可以看出,当所求问题的解是某个事件的概率,或者是某个随机变量的数学期望,或者是与之有关的量时,通过某种试验的方法,得出该事件发生的频率,再通过它得到问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法的基本思想虽然早已被人们提出,却很少被使用。直到电子计算机出现后,使得人们可以通过电子计算机来模拟巨大数目的随机试验过程,使得蒙特卡罗方法得以广泛地应用,在现代化的科学技术中发挥应有的作用。数学建模专题之MonteCarlo方法2019年12月25日11例1在我方某前沿防守地域,敌人以一个炮排(含两门火炮)为单位对我方进行干扰和破坏.为躲避我方打击,敌方对其阵地进行了伪装并经常变换射击地点.经过长期观察发现,我方指挥所对敌方目标的指示有50%是准确的,而我方火力单位,在指示正确时,有1/3的射击效果能毁伤敌人一门火炮,有1/6的射击效果能全部毁伤敌人火炮.现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来,确定有效射击的比率及毁伤敌方火炮的平均值。分析:这是一个概率问题,可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果,而显示不出作战行动的动态过程.为了能显示我方20次射击的过程,现采用模拟的方式。举例数学建模专题之MonteCarlo方法2019年12月25日12需要模拟出以下两件事:[2]当指示正确时,我方火力单位的射击结果情况[1]观察所对目标的指示正确与否模拟试验有两种结果,每一种结果出现的概率都是1/2.因此,可用投掷一枚硬币的方式予以确定,当硬币出现正面时为指示正确,反之为不正确.模拟试验有三种结果:毁伤一门火炮的可能性为1/3(即2/6),毁伤两门的可能性为1/6,没能毁伤敌火炮的可能性为1/2(即3/6).这时可用投掷骰子的方法来确定:如果出现的是1、2、3三个点:则认为没能击中敌人;如果出现的是4、5点:则认为毁伤敌人一门火炮;若出现的是6点:则认为毁伤敌人两门火炮.问题分析数学建模专题之MonteCarlo方法2019年12月25日13i:要模拟的打击次数;k1:没击中敌人火炮的射击总数;k2:击中敌人一门火炮的射击总数;k3:击中敌人两门火炮的射击总数;E:有效射击比率;E1:20次射击平均每次毁伤敌人的火炮数.符号说明数学建模专题之MonteCarlo方法2019年12月25日14模拟框图初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子点数?k1=k1+1k2=k2+1k3=k3+1k1=k1+1i<20?E=(k2+k3)/20E1=(0*k1+1*k2+2*k3)/20停止硬币正面?YNNY1,2,34,56数学建模专题之MonteCarlo方法2019年12月25日15模拟结果消灭敌人火炮数试验序号投硬币结 果指示正确指 示不正确掷骰子结 果0121正∨4∨2正∨4∨3反∨∨4正∨1∨5正∨2∨6反∨∨7正∨3∨8正∨6∨9反∨∨10反∨∨数学建模专题之MonteCarlo方法2019年12月25日16消灭敌人火炮数试验序号投硬币结 果指示正确指 示不正确掷骰子结 果01211正∨2∨12反∨∨13正∨3∨14反∨∨15正∨6∨16正∨4∨17正∨2∨18正∨4∨19反∨∨20正∨6∨从以上模拟结果可计算出:E=7/20=0.35/()E=0.5数学建模专题之MonteCarlo方法2019年12月25日17理论计算设:观察所对目标指示正确确观察所对目标指示不正10jA0:射中敌方火炮的事件;A1:射中敌方一门火炮的事件;A2:射中敌方两门火炮的事件.则由全概率公式:E=P(A0)=P(j=0)P(A0∣j=0)+P(j=1)P(A0∣j=1)=25.02121021P(A1)=P(j=0)P(A1∣j=0)+P(j=1)P(A1∣j=1)=613121021P(A2)=P(j=0)P(A2∣j=0)+P(j=1)P(A2∣j=1)=1216121021E1=33.01212611数学建模专题之MonteCarlo方法2019年12月25日18结果比较理论计算和模拟结果的比较分类项目无效射击有效射击平均值模拟0.650.350.5理论0.750.250.33虽然模拟结果与理论计算不完全一致,但它却能更加真实地表达实际战斗动态过程.要使结果接近理论计算值,必须加大模拟次数,这就要求使用计算机模拟了。用蒙特卡洛方法进行计算机模拟的步骤:[1]设计一个数学模型.[2]根据流程图编写程序,模拟随机现象.可通过具有各种概率分布的模拟随机数来模拟随机现象.[3]分析模拟结果,计算所需要结果.数学建模专题之MonteCarlo方法2019年12月25日19Matlab中的随机数生成函数rand(m,n)生成一个满足均匀分布的mn随机矩阵,矩阵的每个元素都在(0,1)之间。rand(4,5)ans=0.44690.11600.22340.87430.16250.15280.25090.67330.39370.30980.88620.75970.81880.93700.68110.03140.89830.94890.43690.9341数学建模专题之MonteCarlo方法2019年12月25日20Matlab中的随机数生成函数unifrnd(a,b,m,n)生成一个满足均匀分布的mn随机矩阵,矩阵的每个元素都在(a,b)之间。unifrnd(-1,1,5,4)ans=0.8948-0.4123-0.7244-0.31060.1982-0.9361-0.5161-0.23040.89770.7290-0.55400.1897-0.1920-0.13490.73550.0702-0.9179-0.81440.5284-0.3329数学建模专题之MonteCarlo方法2019年12月25日21Matlab中的随机数生成函数randn(m,n)生成一个满足正态mn的随机矩阵hist(randn(1,1000),50)数学建模专题之MonteCarlo方法2019年12月25日22Matlab中的随机数生成函数exprnd(mu,sigma,m,n)生成一个满足N(mu,sigma)正态mn的随机矩阵hist(normrnd(-1,3,1,1000),50)数学建模专题之MonteCarlo方法2019年12月25日23Matlab中的随机数生成函数normrnd(lamda,m,n)生成一个指数分布为lamda的mn随机矩阵hist(exprnd(10,1,1000),50)数学建模专题之MonteCarlo方法2019年12月25日24Matlab中的随机数生成函数poissrnd(lamda,m,n)生成一个泊松分布为lamda的mn随机矩阵hist(poissrnd(3,1,1000),50)数学建模专题之MonteCarlo方法2019年12月25日25Matlab中的随机数生成函数randperm(10)ans=21013875694randperm(m)生成一个由1:m组成的随机排列perms(1:n)生成由1:n组成的全排列,共n!个数学建模专题之MonteCarlo方法2019年12月25日26name的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......random('name',A1,A2,A3,M,N)Matlab中的随机数生成函数数学建模专题之MonteCarlo方法2019年12月25日27随机投掷均匀硬币,验证国徽朝上与朝下的概率是否都是1/2n=10000;%给定试验次数m=0;fori=1:nx=randperm(2)