《蒙特卡罗方式和计算机模拟资料》由会员分享，可在线阅读，更多相关《蒙特卡罗方式和计算机模拟资料（41页珍藏版）》请在金锄头文库上搜索。
1、蒙特卡罗与计算机模拟 内容 计算机模拟 或称仿真 是一种广义数值计算方法 适合解决一些规模大 难以解析化以及受随机因素影响的不确定数学模型目的 了解蒙特卡罗方法的基本思想 掌握利用Matlab对离散 连续系统进行模拟的方法要求 掌握Matlab随机数函数 处理应用问题了解蒙特卡罗方法的起源和基本思想掌握离散系统和连续系统计算机模拟实例掌握随机函数randunifrndnormrndexprnd了解Matlab仿真模块Simulink 蒙特卡罗方法的起源和基本思想 蒙特卡罗方法 MonteCarlomethod 或称计算机随机模拟方法 是一种基于 随机数 的计算方法 源于美国在第二次世界大战研制
2、原子弹的 曼哈顿计划 该计划的主持人之一数学家冯 诺伊曼用驰名世界的赌城 摩纳哥的MonteCarlo 来命名这种方法 为它蒙上了一层神秘色彩 蒙特卡罗方法的基本思想很早以前就被人们所发现和利用 早在17世纪 人们就知道用事件发生的 频率 来决定事件的 概率 19世纪人们用蒲丰投针的方法来计算圆周率 上世纪40年代电子计算机的出现 特别是近年来高速电子计算机的出现 使得用数学方法在计算机上大量 快速地模拟这样的试验成为可能 蒲丰投针实验近似计算圆周率 蒲丰投针实验 法国科学家蒲丰 Buffon 在1777年提出的蒲丰投针实验是早期几何概率一个非常著名的例子 蒲丰投针实验的重要性并非是为了求得比
3、其它方法更精确的 值 而是它开创了使用随机数处理确定性数学问题的先河 是用偶然性方法去解决确定性计算的前导 由此可以领略到从 概率土壤 上开出的一朵瑰丽的鲜花 蒙特卡罗方法 MC 蒲丰投针实验可归结为下面的数学问题 平面上画有距离为a的一些平行线 向平面上任意投一根长为l l a 的针 假设针落在任意位置的可能性相同 试求针与平行线相交的概率P 从而求 蒲丰投针实验近似计算圆周率 蒲丰投针实验 如右图所示 以M表示针落下后的中点 以x表示M到最近一条平行线的距离 以 表示针与此线的交角 针落地的所有可能结果满足 其样本空间视作矩形区域 面积是 针与平行线相交的条件 它是样本空间 子集A 面积是
4、 symslphi int l 2 sin phi phi 0 pi ans l因此 针与平行线相交的概率为 从而有 特别当时 蒲丰投针实验近似计算圆周率 蒲丰投针实验的计算机模拟 formatlong 设置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 样本空间 fori 1 nifx i l sin phi i 2 满足此条
5、件表示针与线的相交plot phi i x i r counter counter 1 统计针与线相交的次数frame counter getframe 描点并取帧endendfren counter n pihat 2 l a fren 用频率近似计算 figure 2 movie frame 1 播放帧动画1次 蒲丰投针实验近似计算圆周率 蒲丰投针实验的计算机模拟 意大利数学家拉泽里尼得到了准确到6位小数的 值 不过他的实验因为太准确而受到了质疑 蒲丰投针实验计算圆周率 蒙特卡罗投点法是蒲丰投针实验的推广 在一个边长为a的正方形内随机投点 该点落在此正方形的内切圆中的概率应为该内切圆与正方
6、形的面积比值 即n 10000 a 2 m 0 fori 1 nx 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随机数的 重置 问题 Matlab的随机数是伪随机数 但在一定的信度之下可以看作真正的随机数 问题是rand函数产生的随机数从一个随机数序列中取出来 而
7、每次启动Matlab时 rand的状态都会被重置 相当于把序列的指针移到了随机数序列的开始 换言之第一次启动Matlab调用的第n次rand函数与下一次启动调用的第n个rand函数产生相同的数值 如果想打乱这种状态 可以为rand指定一个与当前时间相关的初始状态 而不用默认状态 rand state sum 100 clock 或者rand state sum 100 clock rand 若使每次运行程序产生的值是相同的rand seed 0 1 rand 1 非常见分布的随机数的产生 1 逆变换方法 非常见分布的随机数的产生 举例说明 1 连续型随机变量 以指数分布为例 symstxlam
8、bda Fx int lambda exp lambda t t 0 x 分布函数symsr 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 离散型随机变量 以离散分布为例 x 2 4 6 8 px 0 1 0 4 0 3 0 2 以下为程序片段Fx 0 forn 1 length px Fx Fx sum px 1 n endr ra
9、nd index find r Fx x index 1 1 已编写通用离散分布随机数产生程序scatrnd x px n 2 Acceptance Rejection方法最早由VonNeumann提出 现在已经广泛应用于各种随机数的生成 离散系统的计算机模拟实例 范例报童售报问题 一报童每天清晨从邮局订购报纸后零售 每份报纸进价0 35元 售价0 5元 邮局要求最低订购数量为60份 根据过去经验一个报童一天平均售出报纸120份 且满足泊松分布 未售出的报纸只要没有破损可退给邮局 试求报童每天清晨订购多少份报纸可获最大利润 1 数学建模 Pin 0 35元Pout 0 5元Nin 60 200
10、 区间Nout poissrnd 120 随机数 离散系统的计算机模拟实例 2 计算机模拟当订购122份报纸时 可获得最大利润17 3584元 离散系统的计算机模拟实例 范例海港系统的卸载货物问题 离散系统的计算机模拟实例 范例海港系统的卸载货物问题1程序片段 船只到港时间均匀分布 船只卸货时间均匀分布 ShipBetweenTime 1 unifrnd 15 145 1 1 船只到港间隔时间随机化 均匀分布 ShipUnloadTime 1 unifrnd 45 90 1 1 船只卸货时间随机化 均匀分布 通用程序haibor m可实现多次模拟 并将结果保存到H txtdeleteH txt
11、 清除历史数据harbor 100 15 145 45 90 loadH txt Hmean mean H 导入H并按列取平均值 离散系统的计算机模拟实例 范例海港系统的卸载货物问题2程序片段 船只到港时间指数分布 船只卸货时间均匀分布 ShipBetweenTime 1 exprnd 60 1 1 船只到港间隔时间随机化 指数分布 ShipUnloadTime 1 unifrnd 45 90 1 1 船只卸货时间随机化 均匀分布 通用程序haibor2 m可实现多次模拟 结果保存到H2 txtdeleteH2 txt 清除历史数据harbor2 100 60 45 90 loadH2 txt
12、 Hmean2 mean H2 导入H2并按列取平均值 离散系统的计算机模拟实例 范例海港系统的卸载货物问题3程序片段 船只到港时间离散分布 船只卸货时间离散分布 1 编写船只到港间隔离散累积分布函数并作阶梯图 xs 15 10 145 fori 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 fori 1 length px Fx Fx sum px 1 i endplot 10 x Fx rs ho
13、ldon stairs 0 x 5 145 Fx 1 set gca xtick 0 5 145 set gca xgrid on axistight 离散系统的计算机模拟实例 范例海港系统的卸载货物问题3程序片段 船只到港时间离散分布 船只卸货时间离散分布 2 编写船只到港间隔离散累积分布反函数并作线性插值 Fxi 0 0 001 1 eps xi interp1 Fx 0 x Fxi linear r rand 1 n rnd fori 1 nindex find r i Fxi ifxs 1 xi index 1 1 xs length xs rnd rnd xi index 1 1 e
14、ndend 以上程序已编写通用M函数文件harborrnd xs px n 即给出n个满足离散分布 x px 的船只到港间隔随机数 离散系统的计算机模拟实例 范例海港系统的卸载货物问题3程序片段 船只到港时间离散分布 船只卸货时间离散分布 3 编写船只卸货时间离散累积分布函数并作阶梯图 xs 45 5 90 fori 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 fori 1 length px Fx Fx sum px 1 i endplot 4
15、0 x Fx rs holdon stairs 40 x 2 5 90 Fx 1 set gca xtick 40 2 5 90 set gca xgrid on axistight 离散系统的计算机模拟实例 范例海港系统的卸载货物问题3程序片段 船只到港时间离散分布 船只卸货时间离散分布 4 编写船只卸货时间离散累积分布反函数并作线性插值 Fxi 0 0 001 1 eps xi interp1 Fx 0 x Fxi linear r rand 1 n rnd fori 1 nindex find r i Fxi ifxs 1 xi index 1 1 xs length xs rnd rn
16、d xi index 1 1 endend 以上程序已编写通用M函数文件harborrnd xs px n 即给出n个满足离散分布 x px 的船只卸货时间随机数 离散系统的计算机模拟实例 范例海港系统的卸载货物问题3程序片段 船只到港时间指数分布 船只卸货时间均匀分布 5 模拟船只到港间隔 卸货时间均为离散分布的海港系统ShipBetweenTime 1 harborrnd sbtxs sbtpx 1 船只到港间隔时间随机化 离散分布 ShipUnloadTime 1 harborrnd sutxs sutpx 1 船只卸货时间随机化 离散分布 通用程序haibor3 m可实现多次模拟 结果保存到H3 txtdeleteH3 txt 清除历史数据loadharbor mat 载入数据harbor3 100 sbtxs sbtpx sutxs sutpx loadH3 txt Hmean3 mean H3 导入H3并按列取平均值 连续系统的计算机模拟实例 范例某军导弹基地发现正北方向120km处有一艘敌舰以90km h的速度向正东方向行驶 该基地即刻发射导弹进行拦击 导弹速率450km
