蒙特卡罗与计算机模拟有代码.ppt

第六讲 蒙特卡罗与计算机模拟,内容：计算机模拟(或称仿真)是一种广义数值计算 方法，适合解决一些规模大、难以解析化 以及受随机因素影响的不确定数学模型目的：了解蒙特卡罗方法的基本思想，掌握利用 Matlab对离散/连续系统进行模拟的方法要求：掌握Matlab随机数函数，处理应用问题了解蒙特卡罗方法的起源和基本思想掌握离散系统和连续系统计算机模拟实例掌握随机函数 rand unifrnd normrnd exprnd了解Matlab仿真模块 Simulink,蒙特卡罗方法的起源和基本思想,蒙特卡罗方法(Monte Carlo method)，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。源于美国在第二次世界大战研制原子弹的“曼哈顿计划”，该计划的主持人之一数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。蒙特卡罗方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用蒲丰投针的方法来计算圆周率，上世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,蒲丰投针实验近似计算圆周率,蒲丰投针实验：法国科学家蒲丰(Buffon)在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子。蒲丰投针实验的重要性并非是为了求得比其它方法更精确的值，而是它开创了使用随机数处理确定性数学问题的先河，是用偶然性方法去解决确定性计算的前导，由此可以领略到从“概率土壤”上开出的一朵瑰丽的鲜花蒙特卡罗方法(MC)蒲丰投针实验可归结为下面的数学问题：平面上画有距离为a的一些平行线，向平面上任意投一根长为l(la)的针，假设针落在任意位置的可能性相同，试求针与平行线相交的概率P(从而求),蒲丰投针实验近似计算圆周率,蒲丰投针实验：如右图所示，以M表示针落下后的中点，以x表示M到最近一条平行线的距离，以表示针与此线的交角：针落地的所有可能结果满足：其样本空间视作矩形区域,面积是:针与平行线相交的条件：它是样本空间子集A，面积是：syms l phi;int(l/2*sin(phi),phi,0,pi);%ans=l因此，针与平行线相交的概率为：从而有：特别当 时,蒲丰投针实验近似计算圆周率,蒲丰投针实验的计算机模拟：format long;%设置15位显示精度a=1;l=0.6;%两平行线间的宽度和针长figure;axis(0,pi,0,a/2);%初始化绘图板set(gca,nextplot,add);%初始化绘图方式为叠加counter=0;n=2010;%初始化计数器和设定投针次数x=unifrnd(0,a/2,1,n);phi=unifrnd(0,pi,1,n);%样本空间for i=1:n if x(i)l*sin(phi(i)/2%满足此条件表示针与线的相交 plot(phi(i),x(i),r.);frame(i)=getframe;%描点并取帧 counter=counter+1;%统计针与线相交的次数 endendfren=counter/n;pihat=2*l/(a*fren)%用频率近似计算%movie(frame,1)%播放帧动画1次,蒲丰投针实验近似计算圆周率,蒲丰投针实验的计算机模拟：,意大利数学家拉泽里尼得到了准确到6位小数的值，不过他的实验因为太准确而受到了质疑,蒲丰投针实验计算圆周率,蒙特卡罗投点法是蒲丰投针实验的推广：在一个边长为a的正方形内随机投点，该点落在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即n=10000;a=2;m=0;for i=1:n x=rand(1)*a;y=rand(1)*a;if(x-a/2)2+(y-a/2)2=(a/2)2)m=m+1;endenddisp(投点法近似计算的为:,num2str(4*m/n);,常见分布的随机数产生语句,蒙特卡罗方法的关键步骤在于随机数的产生，计算机产生的随机数都不是真正的随机数(由算法确定的缘故)，如果伪随机数能够通过一系列统计检验，我们也可以将其当作真正的随机数使用：,常见分布的随机数产生语句,MATLAB可以直接产生满足各种分布的随机数具体命令如下:产生mn阶0,1上均匀分布的随机数矩阵 rand(m,n)产生一个0,1上均匀分布的随机数 rand 产生mn阶a,b上均匀分布的随机数矩阵 unifrnd(a,b,m,n)产生一个a,b上均匀分布的随机数 unifrnd(a,b)产生一个1:n的随机排列(元素均出现且不重复)p=randperm(n)注意:randperm(6)与unifrnd(1,6,1,6)的区别,常见分布的随机数产生语句,产生mn阶均值为mu方差为sigma的正态分布的随机数矩阵 normrnd(mu,sigma,m,n)产生一个均值为mu方差为sigma的正态分布的随机数 normrnd(mu,sigma)产生mn阶期望值为mu(mu=1/)的指数分布的随机数矩阵 exprnd(mu,m,n)产生一个期望值为mu的指数分布的随机数 exprnd(mu)注意:产生一个参数为的指数分布的随机数应输入 exprnd(1/),常见分布的随机数产生语句,产生mn阶参数为A1,A2,A3的指定分布name的随机数矩阵 random(name,A1,A2,A3,m,n)产生一个参数为为A1,A2,A3的指定分布name的随机数 random(name,A1,A2,A3)举例:产生24阶的均值为0方差为1的正态分布的随机数矩阵 random(Normal,0,1,2,4)name的取值可以是(详情参见help random)：norm or Normal/unif or Uniformpoiss or Poisson/beta or Betaexp or Exponential/gam or Gammageo or Geometric/unid or Discrete Uniform,MATLAB随机数的“重置”问题,Matlab的随机数是伪随机数，但在一定的信度之下可以看作真正的随机数。问题是rand函数产生的随机数从一个随机数序列中取出来，而每次启动Matlab时，rand的状态都会被重置(相当于把序列的指针移到了随机数序列的开始)，换言之第一次启动Matlab调用的第n次rand函数与下一次启动调用的第n个rand函数产生相同的数值。如果想打乱这种状态，可以为rand指定一个与当前时间相关的初始状态，而不用默认状态：rand(state,sum(100*clock);或者rand(state,sum(100*clock)*rand);,非常见分布的随机数的产生,对于常见分布随机数，可由相应Matlab函数直接产生，对于非常见分布随机数可如下处理：1 连续型随机变量(以p116指数分布为例)：,syms t x lambda;Fx=int(lambda*exp(-lambda*t),t,0,x)%分布函数syms r;Fxinv=finverse(Fx,x);%求反函数Fxinv=subs(Fxinv,x,r)%替换反函数变量x为rFxinv=inline(Fxinv)x=Fxinv(3,rand)%产生参数 lambda=3 指数分布的随机数%指数分布随机数产生函数已经提供 exprnd(1/3,1,1),非常见分布的随机数的产生,2 离散型随机变量(以p117离散分布为例)：,x=2,4,6,8;px=0.1,0.4,0.3,0.2;%以下为程序片段Fx=0;for n=1:length(px),Fx=Fx,sum(px(1:n);endr=rand;index=find(rFx);x(index(1)-1)%已编写通用离散分布随机数产生程序 scatrnd(x,px,n),离散系统的计算机模拟实例,范例 海港系统的卸载货物问题(p110-111,p119-129),离散系统的计算机模拟实例,范例 海港系统的卸载货物问题1(p110-111,p119-129)程序片段(船只到港时间均匀分布,船只卸货时间均匀分布)ShipBetweenTime(1)=unifrnd(15,145,1,1);%船只到港间隔时间随机化(均匀分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸货时间随机化(均匀分布)通用程序haibor.m可实现多次模拟，并将结果保存到H.txtdelete H.txt%清除历史数据harbor(100,15,145,45,90)load H.txt;Hmean=mean(H);%导入H并按列取平均值,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题2(p110-111,p119-129)程序片段(船只到港时间指数分布,船只卸货时间均匀分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港间隔时间随机化(指数分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸货时间随机化(均匀分布)通用程序haibor2.m可实现多次模拟，结果保存到H2.txtdelete H2.txt%清除历史数据harbor2(100,60,45,90)load H2.txt;Hmean2=mean(H2);%导入H2并按列取平均值,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3(p110-111,p119-129)程序片段(船只到港时间离散分布,船只卸货时间离散分布)1 编写船只到港间隔离散累积分布函数并作阶梯图：xs=15:10:145;for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;endpx=0.009,0.029,0.035,0.051,0.090,0.161,0.200,0.172,0.125,0.071,0.037,0.017,0.003;Fx=0;for i=1:length(px),Fx=Fx,sum(px(1:i);endplot(10,x,Fx,-rs);hold on;stairs(0,x-5,145,Fx,1);set(gca,xtick,0:5:145);set(gca,xgrid,on);axis tight;,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3(p110-111,p119-129)程序片段(船只到港时间离散分布,船只卸货时间离散分布)2 编写船只到港间隔离散累积分布反函数并作线性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,0,x,Fxi,linear);r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi);if xs(1)=xi(index(1)-1)=xs(length(xs)rnd=rnd,xi(index(1)-1);endend%以上程序已编写通用M函数文件 harborrnd(xs,px,n)%即给出n个满足离散分布(x,px)的船只到港间隔随机数,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3(p110-111,p119-129)程序片段(船只到港时间离散分布,船只卸货时间离散分布)3 编写船只卸货时间离散累积分布函数并作阶梯图：xs=45:5:90;for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;endpx=0.017,0.045,0.095,0.086,0.130,0.185,0.208,0.143,0.091;Fx=0;for i=1:length(px),Fx=Fx,sum(px(1:i);endplot(40,x,Fx,-rs);hold on;stairs(40,x-2.5,90,Fx,1);set(gca,xtick,40:2.5:90);set(gca,xgrid,on);axis tight;,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3(p110-111,p119-129)程序片段(船只到港时间离散分布,船只卸货时间离散分布)4 编写船只卸货时间离散累积分布反函数并作线性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,0,x,Fxi,linear);r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi);if xs(1)=xi(index(1)-1)=xs(length(xs)rnd=rnd,xi(index(1)-1);endend%以上程序已编写通用M函数文件 harborrnd(xs,px,n)%即给出n个满足离散分布(x,px)的船只卸货时间随机数,离散系统的计算机模拟实例,范例 海港系统的卸载货物问题3(p110-111,p119-129)程序片段(船只到港时间指数分布,船只卸货时间均匀分布)5 模拟船只到港间隔/卸货时间均为离散分布的海港系统ShipBetweenTime(1)=harborrnd(sbtxs,sbtpx,1);%船只到港间隔时间随机化(离散分布)ShipUnloadTime(1)=harborrnd(sutxs,sutpx,1);%船只卸货时间随机化(离散分布)通用程序haibor3.m可实现多次模拟，结果保存到H3.txtdelete H3.txt%清除历史数据load harbor.mat%载入数据harbor3(100,sbtxs,sbtpx,sutxs,sutpx)load H3.txt;Hmean3=mean(H3);%导入H3并按列取平均值,连续系统的计算机模拟实例,范例 某军导弹基地发现正北方向120km处有一艘敌舰以90km/h的速度向正东方向行驶，该基地即刻发射导弹进行拦击，导弹速率450km/h，制导系统确保在任一时刻导弹都能对准敌舰(1)试问导弹何时何处击中敌舰(2)如果敌舰即刻发现导弹，并以垂直导弹方向135km/h的速度逃逸，试问导弹何时何处击中敌舰,连续系统的计算机模拟实例,1 将随等距时间连续变化的状态变量轨迹x(t),y(t)用欧拉法离散化：,连续系统的计算机模拟实例,2 编写程序模拟导弹拦击敌舰过程x1=0;y1=0;x2=0;y2=120;t=0.001;v1=450;v2=90;dis=120;axis(0,40,0,140);grid on;set(gca,nextplot,add);for k=1:1000 x1=x1+v1*t*(v2*k*t-x1)/sqrt(v2*k*t-x1)2+(dis-y1)2);y1=y1+v1*t*(dis-y1)/sqrt(v2*k*t-x1)2+(dis-y1)2);x2=x2+v2*t;y2=y2;plot(x1,y1,ro,x2,y2,bs);frame(k)=getframe;if sqrt(x1-x2)2+(y1-y2)2)=0.1,break;endendT=k*t,x1,y1,