计算机仿真1

计计 算算 机机 仿仿 真真 1.1.系统仿真概论系统仿真概论2.时间步长法时间步长法3.事件表法事件表法4.Monte Carlo 方法方法计算机仿真是利用计算机对一个实际系统的计算机仿真是利用计算机对一个实际系统的结构和行为进行动态演示，以评价或预测该结构和行为进行动态演示，以评价或预测该系统的行为效果。它是解决较复杂的实际问系统的行为效果。它是解决较复杂的实际问题的一条有效途径。题的一条有效途径。计算机仿真在航空、机电、冶金等工业部门计算机仿真在航空、机电、冶金等工业部门及社会经济、交通运输、生态系统等方面有及社会经济、交通运输、生态系统等方面有着广泛的应用，已成为分析、研究和设计各着广泛的应用，已成为分析、研究和设计各种系统的重要手段。种系统的重要手段。一一.什么是计算机仿真什么是计算机仿真绪 论二二.为什么要进行计算机仿真为什么要进行计算机仿真1)1)在一个实际系统还没有建立起来之前，要对系统在一个实际系统还没有建立起来之前，要对系统的行为或结果进行分析研究时，计算机仿真是一的行为或结果进行分析研究时，计算机仿真是一种行之有效的方法。种行之有效的方法。2)2)在有些真实系统上做实验会影响系统的正常运行，在有些真实系统上做实验会影响系统的正常运行，如在生产中任意改变工艺参数可能会导致废品，如在生产中任意改变工艺参数可能会导致废品，在经济活动中随意将一个决策付诸行动可能会引在经济活动中随意将一个决策付诸行动可能会引起经济混乱。起经济混乱。3)3)当人是系统的一部分时，他的行为往往会影响实当人是系统的一部分时，他的行为往往会影响实验的效果，这时最好对系统进行仿真研究。验的效果，这时最好对系统进行仿真研究。4)4)在实际系统上做实验时，很难保证每次操作的在实际系统上做实验时，很难保证每次操作的条件相同，因而对实验结果的好坏很难作出正确条件相同，因而对实验结果的好坏很难作出正确的判断。的判断。5)5)在实际系统上做实验可能会花费很长的时间，在实际系统上做实验可能会花费很长的时间，费用太高或者有危险，使得实验不易进行。费用太高或者有危险，使得实验不易进行。6)6)有些系统一旦建立起来之后就无法复原。例如有些系统一旦建立起来之后就无法复原。例如要投资建立一个大型企业，要分析它建成之后的要投资建立一个大型企业，要分析它建成之后的经济效益和社会效益，不能用建立起来试试看的经济效益和社会效益，不能用建立起来试试看的办法，因为建成后就无法回到原来的状态了。办法，因为建成后就无法回到原来的状态了。1)1)难以用数学公式描述的系统，或者没有求解的难以用数学公式描述的系统，或者没有求解的有效数学方法；有效数学方法；2)2)虽然可以用解析方法解决问题，但数学的分析虽然可以用解析方法解决问题，但数学的分析与计算过于复杂，这时计算机仿真可能提供简单与计算过于复杂，这时计算机仿真可能提供简单可行的求解方法；可行的求解方法；3)3)希望能在较短的时间内观察到系统发展的全过希望能在较短的时间内观察到系统发展的全过程，以估计参数对系统行为的影响；程，以估计参数对系统行为的影响；4)4)难以在实际环境中进行实验和观察时，计算机难以在实际环境中进行实验和观察时，计算机仿真是唯一可行的方法，例如太空飞行的研究；仿真是唯一可行的方法，例如太空飞行的研究；5)5)需要对系统或过程进行长期运行的比较，从大需要对系统或过程进行长期运行的比较，从大量方案中寻找最优方案。量方案中寻找最优方案。三三.适合用计算机仿真解决的问题适合用计算机仿真解决的问题1.1.系统仿真概论系统仿真概论1.1 1.1 系统系统1)1)系统系统 一些具有特定的功能、相互之间以一定的规一些具有特定的功能、相互之间以一定的规律联系的对象所组成的总体。律联系的对象所组成的总体。2)2)系统边界系统边界 每个系统都存在于一定的环境中每个系统都存在于一定的环境中,系统与环境之系统与环境之间的分界称为边界间的分界称为边界.边界确定了系统的范围，边界以外对系统的作边界确定了系统的范围，边界以外对系统的作用称为用称为系统的输入系统的输入，系统对边界以外环境的作用，系统对边界以外环境的作用称为称为系统的输出系统的输出.3)3)系统的三要素系统的三要素实体实体：系统内的对象，系统的组成元素：系统内的对象，系统的组成元素.实体确定了系实体确定了系 统的构成，也确定了系统边界。统的构成，也确定了系统边界。属性属性：系统内实体的有效特征，它可以是文字型、数：系统内实体的有效特征，它可以是文字型、数 字型或逻辑型字型或逻辑型.属性也称为描述变量，描述每一属性也称为描述变量，描述每一 实体的特征。实体的特征。活动活动：系统内部发生的任何变化过程。：系统内部发生的任何变化过程。活动定义了系活动定义了系 统内部实体之间的相互作用，从而确定了系统统内部实体之间的相互作用，从而确定了系统 内部发生的变化过程。内部发生的变化过程。4)系统状态系统状态 系统的状态是指在某一时间点上，实体及其属性值的集合。在任意给定的时间，对系统所有实体、属性和活动的情况，都用系统状态加以描述。随机型活动通常用概率分布加以描述，如电随机型活动通常用概率分布加以描述，如电路系统中的噪声、电话系统中通话时间的长路系统中的噪声、电话系统中通话时间的长短等，属于随机系统。短等，属于随机系统。若系统中的某项活若系统中的某项活动结果受随机变化动结果受随机变化的影响，系统输入的影响，系统输入可能得到各种各样可能得到各种各样的结果，则称为随的结果，则称为随机系统。机系统。5)5)系统类型系统类型若执行系统中某项活动的结果，其输出完全可以用输入加以描述，则为确定型系统确定型系统确定型系统随机系统随机系统按活动是否具有随机性分按活动是否具有随机性分5)5)系统类型系统类型连续型系统连续型系统离散型系统离散型系统按状态是否随时间连续变化分按状态是否随时间连续变化分系统S=时间基；输入集；内部状态集；状态转移函数；输出集；输出函数 时间基是描述系统变化的时间坐标；1.2 1.2 模型模型模型模型模型实际系统本质的抽象与简化实际系统本质的抽象与简化 物理模型物理模型物理模型采用一定的比例按真实系统的采用一定的比例按真实系统的“样子样子”制作制作数学模型数学模型数学模型用数学表达式形式来描述系统的内在规律用数学表达式形式来描述系统的内在规律系统模型水平系统模型水平系统模型水平模型描述的详细程度模型描述的详细程度物理仿真物理仿真：在没有计算机以前，仿真都是利用实物或者它的模型来进行研究的，又称物理仿真。优点：直接、形象、易信缺点：模型受限、易破坏、难以重用。数学仿真：数学仿真：对实际系统进行抽象，并将其特性用数学关系加以描述而得到系统的数学模型，对数学模型进行实验的过程称为数学仿真。优点：方便、灵活、经济缺点：受限于系统建模技术，即系统数学模型不易建立。计算机技术的发展为数学仿真创造了环境。1.3 仿真仿真计算机仿真计算机仿真在研究系统过程中，根据相似原理，利用计算机来逼真模仿研究对象。研究对象可以是真实的系统，也可以是设想中的系统。计算机仿真是将研究对象进行数学描述，建模编程，且在计算机中运行实现它不怕破坏、易修改、可重用。计算机仿真可以用于研制产品或设计系统的全过程中，包括方案论证、技术指标确定、设计分析、生产制造、试验测试、维护训练、故障处理等各个阶段。仿真算法仿真算法 仿真算法是将系统模型转换成仿真模型的一类算法，在数字仿真模型中起核心和关键作用。仿真结束输出结果问题的分析项目研究计划仿真步骤及流程仿真步骤及流程1 1系系系系统统统统分分分分析析析析2 2模模模模型型型型构构构构造造造造3 3运运运运行行行行与与与与改改改改进进进进4 4输输输输出出出出结结结结果果果果运行与分析运行与分析进一步运行进一步运行实实 验验 设设 计计是是是是否否否否建立模型收集数据编编 制制 程程 序序程程 序序 验验 证证模模 型型 确确 认认否否否否否否否否否否否否2 时间步长法2.1 2.1 仿真时钟仿真时钟 在进行系统仿真时，可以把整个仿真过程分为许在进行系统仿真时，可以把整个仿真过程分为许多相等的时间间隔，时间间隔的长度可根据实际多相等的时间间隔，时间间隔的长度可根据实际问题分别取作秒、分、小时、天、周、月、年等。问题分别取作秒、分、小时、天、周、月、年等。程序中按此间隔前进的虚拟时钟就是程序中按此间隔前进的虚拟时钟就是仿真的时钟仿真的时钟，此时间间隔称为此时间间隔称为仿真时钟的步长仿真时钟的步长。选取系统的一个初始状态作为仿真时钟的零点，选取系统的一个初始状态作为仿真时钟的零点，仿真时钟每步进一次，就对系统的所有实体、属仿真时钟每步进一次，就对系统的所有实体、属性和活动进行一次前面的扫描考察，按照预定的性和活动进行一次前面的扫描考察，按照预定的计划和目标进行分析、计算和记录系统状态的变计划和目标进行分析、计算和记录系统状态的变化，这个过程一直进行到仿真时钟结束为止。此化，这个过程一直进行到仿真时钟结束为止。此即为时间步长仿真法。即为时间步长仿真法。2.2 2.2 流程图流程图仿真结束初初 始始 化化时间步进时间步进1考察实体或活动考察实体或活动发生事件或变化发生事件或变化处理事件子程序处理事件子程序更新状态数据更新状态数据是是仿真是否完毕仿真是否完毕否否输出结果输出结果是是否否2.3 2.3 时钟步长法特点时钟步长法特点1)1)仿真时钟以等步长前进；仿真时钟以等步长前进；2)2)在一个步长内，认为系统的状态不变；在一个步长内，认为系统的状态不变；3)3)时钟每步进一次，就要对系统进行一次时钟每步进一次，就要对系统进行一次全面的考察，即使状态无变化也要扫描。全面的考察，即使状态无变化也要扫描。例例1 1 池水含盐问题池水含盐问题水：2000 m3盐：2 kg 注入6 m6 m3 3/分钟分钟0.5kg/m0.5kg/m3 3流出4 m4 m3 3/分钟分钟要使池中盐水的浓要使池中盐水的浓度达到度达到0.2kg/m0.2kg/m3 3，需经过多长时间？需经过多长时间？1)1)系统分析系统分析实体：盐水实体：盐水属性：体积，含盐量，含盐率属性：体积，含盐量，含盐率活动：注入，流出活动：注入，流出2.4 2.4 实例分析实例分析在这个问题中，系统的状态随时间连续变化，在这个问题中，系统的状态随时间连续变化，要对这样一个连续系统进行仿真，必需在一系要对这样一个连续系统进行仿真，必需在一系列离散时间点上进行考察，一般取等间隔时间列离散时间点上进行考察，一般取等间隔时间点，设间隔为点，设间隔为 t t（即仿真时钟步长）。（即仿真时钟步长）。取取T=0T=0作为系统仿真的初始时刻，当池水中盐作为系统仿真的初始时刻，当池水中盐的浓度达到的浓度达到0.2kg/m0.2kg/m3时仿真结束。时仿真结束。注水速度：注水速度：6,6,排水速度：排水速度：4 4每分钟水的体积增加每分钟水的体积增加,6,64 42 2注入水的含盐率：注入水的含盐率：0.5,0.5,每分钟向池内注入盐每分钟向池内注入盐 60.5 60.53 3最终含盐率：最终含盐率：SF=0.2SF=0.2t t时刻时刻水的体积水的体积 VT VT水的含盐量水的含盐量 ST ST水的含盐率水的含盐率 SR SR每分钟向池外流出盐每分钟向池外流出盐 4SR 4SR每分钟池内盐增加每分钟池内盐增加 3 3 4SR 4SR2)2)模型构造模型构造收集、整理、分析数据收集、整理、分析数据当前时刻当前时刻t t水的体积水的体积 VT VT水的含盐量水的含盐量 ST ST水的含盐率水的含盐率 SR SR下一时刻下一时刻 t+t+t t 水的体积水的体积:VT+6:VT+6 t t 水的含盐量水的含盐量:ST+3-4 SR :ST+3-4 SR t t 水的含盐率水的含盐率:水的含盐量水的含盐量/水的体积水的体积 系统的状态用系统的状态用 VT VT，ST ST，SR SR 表示表示分析系统的状态更新规则分析系统的状态更新规则初初 始始 化化:VT=2000,ST=2,SR=0.001,t=0,t=1,T=0 t t+t,T T+t计算池水体积计算池水体积VT，含盐量，含盐量ST，含盐率，含盐率SR否否SR SF?T=10?否否是是打印打印t，VT,ST,SR下一个观察周期开始，下一个观察周期开始，T 0是是结束打印打印t，VT,ST,SR仿仿真真流流程程在这个问题中，系统的状态随时间连续变化，在这个问题中，系统的状态随时间连续变化，系统离散化后，仿真结果与离散化的时间间隔系统离散化后，仿真结果与离散化的时间间隔即仿真时钟步长有关，显然步长越小，结果越即仿真时钟步长有关，显然步长越小，结果越精确。但步长越小，计算量越大。精确。但步长越小，计算量越大。本例中可取本例中可取 t=5 min t=5 min t=1 mint=1 min t=0.1 mint=0.1 min等进行计算，并与解析结果比较。等进行计算，并与解析结果比较。3)3)运行与改进运行与改进4)4)设计格式，输出仿真结果设计格式，输出仿真结果t tVTVTSTSTSRSR10102020202031.6936531.693650.015690.0156920202040204060.8100860.810080.029090.02909808021402140224.5663224.56630.103970.10397909021602160250.2118250.21180.114780.1147817017023402340441.2655441.26550.188570.1885718018023602360463.5802463.58020.196430.1964318518523702370474.6208474.62080.200260.20026 t=1 mint=1 min的结果的结果思考写出写出MATLABMATLAB程序。分别设置程序。分别设置 t=5 min t=5 min，t=1 mint=1 min，t=0.1 mint=0.1 min 比较仿真结果。比较仿真结果。？追击问题的计算机模拟追击问题的计算机模拟练习如图，正方形如图，正方形ABCDABCD的四的四个顶点各有一个人，在某个顶点各有一个人，在某一时刻，四人同时出发以一时刻，四人同时出发以匀速匀速v v走向顺时针的下一走向顺时针的下一个人，如果他们的方向始个人，如果他们的方向始终保持对准目标，则最终终保持对准目标，则最终将按螺旋状曲线汇合于中将按螺旋状曲线汇合于中心点心点A A，1)1)求出每个人的轨迹；求出每个人的轨迹；2)2)若四个人的速度不全相若四个人的速度不全相等，结果如何？等，结果如何？A AA AB BC CD D例例2 2 库存问题库存问题方案编号方案编号订货点订货点P P订货量订货量QQ方案方案1 1125125150150方案方案2 2125125250250方案方案3 3150150250250方案方案4 4175175250250方案方案5 5175175300300某自行车商店的仓库管理员采用一种简单的订货策某自行车商店的仓库管理员采用一种简单的订货策略略,当库存量降到当库存量降到P P辆时就向厂家订货，每次订辆时就向厂家订货，每次订QQ辆，辆，如果某一天的需求量超过了库存量，商店就有销售如果某一天的需求量超过了库存量，商店就有销售损失和信誉损失，损失和信誉损失，但如果库存过多，但如果库存过多，将导致资金积压和将导致资金积压和保管费增加，若现保管费增加，若现有以下五种库存策有以下五种库存策略，请帮仓库管理略，请帮仓库管理员员确定花费最少的确定花费最少的策略策略。注意注意;P Q;P Q其它已知条件其它已知条件1)1)从发出订货单到收到货物隔从发出订货单到收到货物隔3 3天；天；2)2)每辆自行车每天的保管费为每辆自行车每天的保管费为0.750.75元；元；3)3)每辆自行车每天的缺货费为每辆自行车每天的缺货费为1.81.8元；元；4)4)每次的订货费为每次的订货费为7575元；元；5)5)原始库存为原始库存为115115辆，并假设第一天没有发出订货单；辆，并假设第一天没有发出订货单；6)6)自行车每天的需求量服从自行车每天的需求量服从0 09999之间的均匀分布。之间的均匀分布。1)1)系统分析系统分析实体实体1 1：自行车：自行车 属性：保管费，缺货费，库存量，需求量属性：保管费，缺货费，库存量，需求量活动：进货，销售活动：进货，销售实体实体2 2：订单：订单属性：从发出订货单到收到货物间隔属性：从发出订货单到收到货物间隔,订货费订货费活动：发出订单，收到订货活动：发出订单，收到订货这是一个随机离散系统的仿真问题。这是一个随机离散系统的仿真问题。2)2)模型建立模型建立输入参数输入参数 订货点订货点 P P；订货量；订货量 Q Q 注意注意 P QP 0)是是QL1，s=1产生顾客离开时刻产生顾客离开时刻更新更新事件表更新有关变量更新有关变量更新更新事件表出纳员闲出纳员闲置置S=0否否上述仿真过程做了如下假设上述仿真过程做了如下假设先到先服务；先到先服务；顾客不会因排队时间太长而离开；顾客不会因排队时间太长而离开；注意：假设注意：假设2 2可修改。可修改。3)3)程序运行、分析程序运行、分析4)4)设计格式，输出仿真结果设计格式，输出仿真结果为了了解仿真过程，需要输出一些中间过程变量及参为了了解仿真过程，需要输出一些中间过程变量及参量的值，可按事件处理的顺序来列表：量的值，可按事件处理的顺序来列表：输出参数输出参数 顾客的平均等待时间顾客的平均等待时间CMCMCT/NCT/N最大队长最大队长QM QM 在程序中记录在程序中记录出纳员的工作强度出纳员的工作强度SRSRN N个顾客的总服务时间个顾客的总服务时间/仿真结束时间仿真结束时间序号序号事件类型事件类型时时 刻刻顾客编号顾客编号出纳员状态出纳员状态队长队长1 11 12.52.51 11 10 02 21 15 52 21 11 1思考思考编写出编写出MATLABMATLAB程序。程序。如何修改假设如何修改假设2?2?，写出修改后的程序。，写出修改后的程序。？4 Monte Carlo 方法4.1 Monte Carlo 4.1 Monte Carlo 方法方法 Monte Carlo Monte Carlo 是世界著名赌城摩洛哥的蒙特卡是世界著名赌城摩洛哥的蒙特卡罗，在二次世界大战中，美国军方将一项绝密研罗，在二次世界大战中，美国军方将一项绝密研究计划的代号命名为究计划的代号命名为Monte Carlo Monte Carlo，其目的是研，其目的是研究铀裂变过程中链式反应的能量计算问题。究铀裂变过程中链式反应的能量计算问题。Monte Carlo Monte Carlo 方法是随机模型的计算机仿真方法，方法是随机模型的计算机仿真方法，在用传统方法难以解决的问题中，有很大一部分在用传统方法难以解决的问题中，有很大一部分可以用概率模型描述，由于这类问题含有不确定可以用概率模型描述，由于这类问题含有不确定的随机因素，难以用定量分析法得到解析结果，的随机因素，难以用定量分析法得到解析结果，在这种情况下，在这种情况下，Monte Carlo Monte Carlo 方法是非常有效。方法是非常有效。Monte Carlo Monte Carlo 方法的历史可追溯到方法的历史可追溯到17771777年法国科学家年法国科学家浦丰浦丰(Buffon)(Buffon)提出的一种计算圆周率的随机实验方法提出的一种计算圆周率的随机实验方法随机投针法随机投针法通过计算可知针与平行线通过计算可知针与平行线相交的概率为：相交的概率为：p p=2=2 l l/d d其中其中其中其中d d 为平行线间距离为平行线间距离l l 为针的长度，为针的长度，l l d d将针投将针投n n次，若有次，若有m m次与平次与平行线相交，则可用行线相交，则可用m/nm/n作作为为p p的近似值，由此可计的近似值，由此可计算出算出 的近似值为的近似值为2nl/md2nl/mdMonte Carlo Monte Carlo 方法的步骤方法的步骤4.2 Monte Carlo 4.2 Monte Carlo 方法的步骤和数学原理方法的步骤和数学原理问题分析建模问题分析建模问题分析建模问题分析建模建立问题的概率模型建立问题的概率模型建立问题的概率模型建立问题的概率模型抽抽抽抽 样：按照假设的分布，产生随机数样：按照假设的分布，产生随机数样：按照假设的分布，产生随机数样：按照假设的分布，产生随机数计算有关结果计算有关结果计算有关结果计算有关结果抽样计算结束？抽样计算结束？抽样计算结束？抽样计算结束？否否依照多次抽样试验的结果依照多次抽样试验的结果依照多次抽样试验的结果依照多次抽样试验的结果估计有关统计参数估计有关统计参数估计有关统计参数估计有关统计参数是是算出模型解算出模型解算出模型解算出模型解的近似值的近似值的近似值的近似值结果满意？结果满意？结果满意？结果满意？是是结束结束否否Monte Carlo Monte Carlo 方法的数学原理方法的数学原理1)1)解的稳定性问题解的稳定性问题：设：设 x x1,x x2 x xn是服从同一分布的随机是服从同一分布的随机变量，且有有限的数学期望变量，且有有限的数学期望 和方差和方差 2，根据大数定律，根据大数定律，x x1,x x2 x xn的算术平均当的算术平均当nn时以概率时以概率1 1收敛到收敛到，即，即可知，可知，Monte Carlo Monte Carlo 方法收敛的速度是方法收敛的速度是和一般的数值方法相比较慢。和一般的数值方法相比较慢。2)2)收敛速度收敛速度由中心极限定理由中心极限定理4.3 4.3 随机数的生成方法随机数的生成方法用用Monte Carlo Monte Carlo 方法解题时，需要根据随机变方法解题时，需要根据随机变量的分布生成随机变量的若干个具体的取值，量的分布生成随机变量的若干个具体的取值，叫做抽样。随机变量的抽样方法很多，不同的叫做抽样。随机变量的抽样方法很多，不同的分布采用的抽样方法不尽相同。分布采用的抽样方法不尽相同。0,1 0,1 区间上的均匀分布随机变量的抽样是其区间上的均匀分布随机变量的抽样是其中最基本的问题，在计算机上，其他分布的随中最基本的问题，在计算机上，其他分布的随机数都是在此基础上产生的。机数都是在此基础上产生的。目前计算机上常用的高级语言都有产生均匀分目前计算机上常用的高级语言都有产生均匀分布随机数的系统函数，我们可以直接使用而不布随机数的系统函数，我们可以直接使用而不必关心其实现原理。必关心其实现原理。仿真对随机数的要求仿真对随机数的要求1 1）准确性：即生成的随机数准确的服从要求的分布）准确性：即生成的随机数准确的服从要求的分布2 2）快速性：有些随机仿真问题中，需要生成几万甚）快速性：有些随机仿真问题中，需要生成几万甚至几十万个随机数，这样就要求生成随机数的速度要至几十万个随机数，这样就要求生成随机数的速度要快。这一速度极大的影响仿真的速度。快。这一速度极大的影响仿真的速度。下面介绍几种从均匀分布随机数产生非均匀分布随机下面介绍几种从均匀分布随机数产生非均匀分布随机数的方法。用数的方法。用 r,r1,r2,r3,r,r1,r2,r3,代表独立的代表独立的0,10,1均匀均匀分布的随机数。分布的随机数。设连续型随机变量设连续型随机变量X X的分布函数为的分布函数为F F(x x),),若若F F(x x)存存在反函数在反函数F F1 1(x x)，则随机数，则随机数 x x=F=F1 1(r r)是是X X的一的一个样本。个样本。例例 指数分布随机数的生成指数分布随机数的生成F F(x x)1 1e e-x x,x0 x0 ，F F1 1(x x)=-ln(1-=-ln(1-x x)/)/抽样方法为：抽样方法为：x x=-ln(1-=-ln(1-r r)/)/1 1）分布函数抽样法）分布函数抽样法生成均匀分布随机变数生成均匀分布随机变数 r r (0,1)(0,1)若若 F Fk1 r r F Fk 则随机数则随机数 x x=x xk 是是X X的一个样本。的一个样本。X Xx x1 1x x2 2x xk kx xk+1k+1P Pp p1 1p p2 2p pk kp pk+1k+1F FF F1 1 F F2 2F Fk kF Fk+1k+1设离散型随机变量设离散型随机变量X X的分布率为的分布率为其中其中F Fk=p=p1+p+p2+p+pk 显然有显然有 F F1 F F2.F Fn F Fn+1.1.1令令 F F0=0=0 x x1 1x x2 2x xk k-1-1x xk k0 0 F F1 1 F F2 2F Fk k1 1F Fk kr r有些随机变量的分布函数很难求出其反函数，有些随机变量的分布函数很难求出其反函数，F F1 1(x x)，若能推倒出其近似表达式，若能推倒出其近似表达式G G 1 1(x x),),则随则随机数机数 x x=G=G1 1(r r)是是X X的一个近似样本。的一个近似样本。例例 标准正态分布随机数的生成标准正态分布随机数的生成可用以下两种近似方法可用以下两种近似方法2 2）近似抽样法）近似抽样法标准正态分布的分布函数为标准正态分布的分布函数为方法方法1 1对反函数对反函数F F1 1(x x)作逼近作逼近可得抽样方法为可得抽样方法为其中其中方法方法2 2由中心极限定理可知：当由中心极限定理可知：当n n较大时较大时近似服从标准正态分布，可得抽样方法为近似服从标准正态分布，可得抽样方法为实际中常取实际中常取n=12n=12，得，得4.4 4.4 实例分析实例分析实例分析实例分析计算二重积分计算二重积分1)1)建立待计算积分的概率模型建立待计算积分的概率模型。设二维随机变量设二维随机变量 r=(X,Y)r=(X,Y)服从区域服从区域 上的均匀分布上的均匀分布,则随机则随机变量变量u(X,Y)u(X,Y)的数学期望为的数学期望为E u(X,Y)=E u(X,Y)=例例1 1 用用Monte Carlo Monte Carlo 方法计算积分方法计算积分2)2)抽样计算抽样计算 生成二维均匀分布随机数生成二维均匀分布随机数 r1,r2,r3,rr1,r2,r3,rn,计算计算u(u(ri ri)i=1,2,ni=1,2,nI I的近似值为的近似值为1)1)建立概率模型建立概率模型。当粒子在某内点当粒子在某内点B B时，统计与其相连的临时，统计与其相连的临点数点数k,k,并将临点依次用并将临点依次用1 1到到k k编号，然后从编号，然后从1 1到到k k中等概率地任选一个数，则粒子移到对中等概率地任选一个数，则粒子移到对应此编号的点。若此点为边界点，则游动应此编号的点。若此点为边界点，则游动结束，记录边界吸收的能量，然后开始下结束，记录边界吸收的能量，然后开始下一轮游动。一轮游动。例例2 2 二维随机游动问题二维随机游动问题 A A i i i 1,101,10组成，粒子在任何一组成，粒子在任何一个内点以等概率向其临近点移动，若个内点以等概率向其临近点移动，若到达边界点到达边界点A Ai则被吸收。边界获得能则被吸收。边界获得能量量f f i,求由求由Q Q点出发的粒子被边界吸收点出发的粒子被边界吸收的能量。的能量。如图所示的网格，由内点如图所示的网格，由内点 Q,P Q,P i,i i=1,2,3,4,5 =1,2,3,4,5 及边界点及边界点2)2)抽样计算抽样计算在第在第mm轮游动中，当粒子游动到某内点轮游动中，当粒子游动到某内点B B时，生成均匀分布时，生成均匀分布随机数随机数 r r，若与，若与B B其相连的临点数为其相连的临点数为k k，则将区间，则将区间(0,1)(0,1)分成分成k k个小区间，若个小区间，若r r落在第落在第i i个小区间，则粒子游动到个小区间，则粒子游动到B B的第的第i i个临个临点，如此继续，直至粒子被某边界点点，如此继续，直至粒子被某边界点A A i吸收，则本轮游动吸收，则本轮游动结束，记录边界吸收的能量结束，记录边界吸收的能量E Em=m=E Ei，然后开始下一轮游动。，然后开始下一轮游动。则由则由Q Q点出发的粒子被边界吸收的能量的点出发的粒子被边界吸收的能量的近似值为近似值为