蒙特卡罗方法与计算机模拟

1、蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟 蒙特卡罗与计算机模拟内容：内容：计算机模拟计算机模拟(或称仿真或称仿真)是一种广义数值计算是一种广义数值计算 方法，适合解决一些规模大、难以解析化方法，适合解决一些规模大、难以解析化 以及受随机因素影响的不确定数学模型以及受随机因素影响的不确定数学模型目的：目的：了解蒙特卡罗方法的基本思想，掌握利用了解蒙特卡罗方法的基本思想，掌握利用 Matlab对离散对离散/连续系统进行模拟的方法连续系统进行模拟的方法要求：要求：掌握掌握Matlab随机数函数，处理应用问题随机数函数，处理应用问题了解蒙特卡罗方法的起源和基本思想了解蒙特卡罗方法的起源和基本思

2、想掌握离散系统和连续系统计算机模拟实例掌握离散系统和连续系统计算机模拟实例掌握随机函数掌握随机函数 rand unifrnd normrnd exprnd了解了解Matlab仿真模块仿真模块 Simulink 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒙特卡罗方法的起源和基本思想 蒙特卡罗方法蒙特卡罗方法(Monte Carlo method)，或称计，或称计算机随机模拟方法，是一种基于算机随机模拟方法，是一种基于“随机数随机数”的计的计算方法。源于美国在第算方法。源于美国在第二二次世界大战研制原子弹次世界大战研制原子弹的的“曼哈顿计划曼哈顿计划”，该计划的主持人之一数学家该计划的主持人

3、之一数学家冯冯诺伊曼诺伊曼用驰名世界的用驰名世界的赌城赌城摩纳哥的摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色来命名这种方法，为它蒙上了一层神秘色彩。彩。 蒙特卡罗蒙特卡罗方法的方法的基本思想基本思想很早以前就被人们所很早以前就被人们所发现和利用。早在发现和利用。早在17世纪，人们就知道世纪，人们就知道用事件发用事件发生的生的“频率频率”来决定事件的来决定事件的“概率概率”。19世纪人世纪人们用们用蒲丰蒲丰投针投针的方法来的方法来计算计算圆周率圆周率，上，上世纪世纪40年年代电子计算机的出现，特别是近年来代电子计算机的出现，特别是近年来高速电子计高速电子计算机算机的出现，

4、使得用数学方法在计算机上大量、的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。快速地模拟这样的试验成为可能。蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验： 法国科学家蒲丰法国科学家蒲丰(Buffon)在在1777年提年提出的蒲丰投针实验是早期几何概率一个出的蒲丰投针实验是早期几何概率一个非常著名的例子。非常著名的例子。蒲丰投针实验蒲丰投针实验的重要的重要性并非是为了求得比其它方法更精确的性并非是为了求得比其它方法更精确的值，而是它值，而是它开创了开创了使用使用随机数处理确定性数学问随机数处理确定性数学问题题的先河，是用的

5、先河，是用偶然性方法去解决确定性计算偶然性方法去解决确定性计算的的前导，前导，由此由此可以可以领略到从领略到从“概率土壤概率土壤”上开出的上开出的一朵瑰丽的鲜花一朵瑰丽的鲜花蒙特卡罗蒙特卡罗方法方法(MC) 蒲丰投针实验蒲丰投针实验可归结为下面的数学问题可归结为下面的数学问题：平面：平面上画有距离为上画有距离为a的一些平行线，向平面上任意投一的一些平行线，向平面上任意投一根长为根长为l(la)的针，假设针落在任意位置的可能性的针，假设针落在任意位置的可能性相同，试求相同，试求针与平行线相交的概率针与平行线相交的概率P(从而求从而求)蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒲丰投针实验近

6、似计算圆周率蒲丰投针实验蒲丰投针实验： 如右图所示，以如右图所示，以M表示针落下后的中点，表示针落下后的中点，以以x表示表示M到最近一条到最近一条平行线的距离，以平行线的距离，以表示针与此线的交角：表示针与此线的交角：针落地的所有可能结果满足：针落地的所有可能结果满足：其样本空间视作矩形区域其样本空间视作矩形区域, 面积是面积是:针与平行线相交的条件：针与平行线相交的条件：它是样本空间它是样本空间子集子集A，面积是：，面积是：syms l phi; int(l/2*sin(phi),phi,0,pi); %ans=l因此，针与平行线相交的概率为：因此，针与平行线相交的概率为：从而有：从而有：

7、特别当特别当 时时02,0 xa( )2Sa 0sin2,0 xl0( )sin2S Al dl( )/ ( ) 2 /p S A Sl a 2l ap2al1 p蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验的计算机模拟：的计算机模拟：format long; %设置设置15位显示精度位显示精度a=1; l=0.6; %两平行线间的宽度两平行线间的宽度和和针长针长figure; axis(0,pi,0,a/2); %初始化绘图板初始化绘图板set(gca,nextplot,add); %初始化绘图方式为叠加初始化绘图方式为叠加counter

8、=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.); counter=counter+1; %统计针与线相交的次数统计针与线相交的次数frame(counter)=getframe; %描点并取帧描点并取帧 endendfren=counter/n; pihat=2*l/(a*fren) %用用频率频率

9、近似计算近似计算figure(2)movie(frame,1) %播放帧动画播放帧动画1次次蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒲丰投针实验近似计算圆周率蒲丰投针实验蒲丰投针实验的计算机模拟：的计算机模拟： 意大利数学家拉泽里尼得到了准确到意大利数学家拉泽里尼得到了准确到6位小数位小数的的值，不过他的实验因为值，不过他的实验因为太准确太准确而受到了质疑而受到了质疑蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟蒲丰投针实验计算圆周率蒙特卡罗投点法蒙特卡罗投点法是是蒲丰投针实验蒲丰投针实验的推广：的推广： 在一个边长为在一个边长为a的正方形内随机投点，该点落的正方形内随机投点，该点落

10、在此正方形的内切圆中的概率应为该内切圆与正在此正方形的内切圆中的概率应为该内切圆与正方形的面积比值，即方形的面积比值，即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);22/2:/4aaxyo(a/2,a/2)蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟常见分布的随机数产生语句 蒙特卡罗方法的蒙特卡罗方法的关键步骤关键步骤在于在于随机数随机数的产生，的产生，

11、计算机产生的随机数都不是计算机产生的随机数都不是真正的随机数真正的随机数(由算法由算法确定的缘故确定的缘故)，如果，如果伪随机数伪随机数能够通过一系列统计能够通过一系列统计检验，我们也可以将其检验，我们也可以将其当作真正的随机数当作真正的随机数使用：使用： 蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟MATLAB随机数的“重置”问题 Matlab的随机数是的随机数是伪随机数伪随机数，但在一定的信度之下可，但在一定的信度之下可以看作真正的随机数。问题是以看作真正的随机数。问题是rand函数函数产生的随机数从一产生的随机数从一个随机数序列中取出来个随机数序列中取出来，而每次启动，而每次启动Ma

12、tlab时，时，rand的状的状态都会被重置态都会被重置(相当于把序列的指针移到了随机数序列的开相当于把序列的指针移到了随机数序列的开始始)，换言之第一次启动，换言之第一次启动Matlab调用的第调用的第n次次rand函数与下函数与下一次启动调用的第一次启动调用的第n个个rand函数函数产生相同的数值产生相同的数值。 如果想如果想打乱这种状态打乱这种状态，可以为，可以为rand指定一个指定一个与当前时与当前时间相关间相关的初始状态，而不用默认状态：的初始状态，而不用默认状态：rand(state,sum(100*clock);或者或者 rand(state,sum(100*clock)*ran

13、d);若使每次运行程序产生的值是相同的rand(seed,0.1);rand(1)非常见分布的随机数的产生非常见分布的随机数的产生1.逆变换方法逆变换方法1 1 ( ) (0,1) ( ) F xUuFu由定理，要产生来自的随机数，只要先产生来自随机数，然后计算 即可。具体步骤如下：(1) (0,1)U上生成均匀分布的随机数 。-1(2) , ( ) ( ) XXUFFx计算则为来自 分布的随机数.蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟非常见分布的随机数的产生举例说明：举例说明：1 连续型随机变量连续型随机变量(以指数分布为例以指数分布为例)：e01-e0( )( )( )0000t

14、xxtxf tf t dt F xtx 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)蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模

15、拟非常见分布的随机数的产生2 离散型随机变量离散型随机变量(以离散分布为例以离散分布为例)：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)2.Acceptance-Rejection 方法 最早由最早由 Von Neumann提出，现在已经提出，现在已经广泛应用于各种随机数的生成。广泛应

16、用于各种随机数的生成。离散系统的计算机模拟实例范例范例 报童售报问题报童售报问题:一报童每天清晨从邮局订购报一报童每天清晨从邮局订购报纸后零售，每份报纸进价纸后零售，每份报纸进价0.35元，售价元，售价0.5元，邮元，邮局要求最低订购数量为局要求最低订购数量为60份，根据过去经验一个份，根据过去经验一个报童一天平均售出报纸报童一天平均售出报纸120份份(且满足泊松分布且满足泊松分布)，未售出的报纸只要没有破损可退给邮局，试求报未售出的报纸只要没有破损可退给邮局，试求报童每天清晨订购多少份报纸可获最大利润？童每天清晨订购多少份报纸可获最大利润？1 数学建模数学建模()(,)() ()outout

17、ininoutininoutinoutinoutinnpnpnnProfit nnnppnnPin=0.35元元 Pout=0.5元元Nin60,200区间区间Nout=poissrnd(120)随机数随机数离散系统的计算机模拟实例2 计算机模拟计算机模拟当订购当订购 122 份报纸时，可获得最大利润份报纸时，可获得最大利润 17.3584 元元蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物问题海港系统的卸载货物问题蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题

18、1程序片段程序片段(船只到港时间船只到港时间均匀分布均匀分布,船只卸货时间船只卸货时间均匀分布均匀分布)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.tx

19、t;Hmean=mean(H); %导入导入H并按列取平均值并按列取平均值蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题2程序片段程序片段(船只到港时间船只到港时间指数指数分布分布,船只卸货时间船只卸货时间均匀分布均匀分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港间隔时间随机化船只到港间隔时间随机化(指数分布指数分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸货时间随机化船只卸货时间随机化(均匀分布均匀分布)通用程序通用程序haibor2

20、.m可实现多次模拟，结果保存到可实现多次模拟，结果保存到H2.txtdelete H2.txt %清除历史数据清除历史数据harbor2(100,60,45,90)load H2.txt;Hmean2=mean(H2); %导入导入H2并按列取平均值并按列取平均值蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题3程序片段程序片段(船只到港时间船只到港时间离散分布离散分布,船只卸货时间船只卸货时间离散分布离散分布)1 编写编写船只到港间隔船只到港间隔离散累积分布函数并作阶梯图：离散累积分布函数并作阶梯图：xs=15:1

21、0: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;蒙特卡罗方法蒙特卡罗方法与计

22、算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题3 程序片段程序片段(船只到港时间船只到港时间离散分布离散分布,船只卸货时间船只卸货时间离散分布离散分布)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(i

23、ndex(1)-1); endend%以上程序已编写通用以上程序已编写通用M函数文件函数文件 harborrnd(xs,px,n)%即给出即给出n个满足离散分布个满足离散分布(x,px)的的船只到港间隔船只到港间隔随机数随机数蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题3程序片段程序片段(船只到港时间船只到港时间离散分布离散分布,船只卸货时间船只卸货时间离散分布离散分布)3 编写编写船只卸货时间船只卸货时间离散累积分布函数并作阶梯图：离散累积分布函数并作阶梯图：xs=45:5:90; for i=1:length

24、(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;4042.54547.55052.55557.56062.56567.57072.57577.58082.58

25、587.59000.20.40.60.81蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题3程序片段程序片段(船只到港时间船只到港时间离散分布离散分布,船只卸货时间船只卸货时间离散分布离散分布)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(inde

26、x(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); endend%以上程序已编写通用以上程序已编写通用M函数文件函数文件 harborrnd(xs,px,n)%即给出即给出n个满足离散分布个满足离散分布(x,px)的的船只卸货时间船只卸货时间随机数随机数蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟离散系统的计算机模拟实例范例范例 海港系统的卸载货物海港系统的卸载货物问题问题3程序片段程序片段(船只到港时间船只到港时间指数指数分布分布,船只卸货时间船只卸货时间均匀分布均匀分布)5 模拟模拟船只到港间隔船只到港间隔 / 卸货时间卸货时间均为离散分布的海港

27、系统均为离散分布的海港系统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)l

28、oad H3.txt;Hmean3=mean(H3); %导入导入H3并按列取平均值并按列取平均值蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例范例范例 某军导弹基地发现正北方向某军导弹基地发现正北方向120km处有一艘处有一艘敌舰以敌舰以90km/h的速度向正东方向行驶，该基地即的速度向正东方向行驶，该基地即刻发射导弹进行拦击，导弹速率刻发射导弹进行拦击，导弹速率450km/h，制导系，制导系统确保在任一时刻导弹都能对准敌舰统确保在任一时刻导弹都能对准敌舰(问题问题1) 试问导弹何时何处击中敌舰试问导弹何时何处击中敌舰o(x,y)120-y90t-x22120tan

29、90()()450dyydxtxdydxdtdtsyms x y t dydt dxdt; solve(dydt/dxdt=(120-y)/(90*t-x),dydt2+dxdt2=4502,dydt,dxdt);ans=ans.dxdt,ans.dydt; dxdt=ans(2,1); dydt=ans(2,2); pretty(dxdt); pretty(dydt);蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例1 将随等距时间将随等距时间连续变化的状态变量轨迹连续变化的状态变量轨迹x(t),y(t)用用欧拉法欧拉法离散化：离散化：12212290450(90)(

30、120)120450(90)(120)kkkkkkkkkkkxxxkxyyyykxy1111(,)()kkkkkkkkkkkkkkyyf tyttdxxxdtdyyydt向前欧拉法蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例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)/sq

31、rt(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%微分方程求解和计算机模拟过程已整合进微分方程求解和计算机模拟过程已整合进daodan1.m蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例3 导弹拦击敌舰过程模拟动画导弹拦击敌舰过程模拟动画(AVI视

32、频视频)warning offmovie2avi(frame,daodan1,compression,Indeo5,quality,100,fps,6); %将模拟过程编码成将模拟过程编码成AVI视频视频导弹于导弹于0.2770小时在小时在(25.0018千米千米, 119.9306千米千米)击中敌舰击中敌舰蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例(问题问题2) 如果敌舰即刻发现导弹，并以垂直导弹方如果敌舰即刻发现导弹，并以垂直导弹方向向135km/h速度逃逸，试问导弹何时何处击中敌舰速度逃逸，试问导弹何时何处击中敌舰1 建立模型建立模型121121221221

33、22112222tantan()()450()()135dyyydxxxdxyydyxxdydxdtdtdydxdtdto(x1,y1)x2-x1y2-y1(x2,y2)蒙特卡罗方法蒙特卡罗方法与计算机模拟与计算机模拟连续系统的计算机模拟实例2 求解并模拟模型求解并模拟模型 (下为片段下为片段,完整见完整见daodan2.m)syms x1 y1 x2 y2 t dy1dt dx1dt dy2dt dx2dt; solve(dy1dt/dx1dt=(y2-y1)/(x2-x1),dy2dt/dx2dt=(x2-x1)/(y2-y1),dy1dt2+dx1dt2=4502,dy2dt2+dx2dt2=1352,dy1dt,dx1dt,dy2dt,dx2dt);for k=1:1