计算机模拟
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
的体积为 2 2 1 4 。
创建p_jifen.m如下：
clear; N = 1000000; MATLAB 程序 x = unifrnd(-1,1,N,1); p_jifen.m y = unifrnd(-1,1,N,1); z = rand(N,1); v1 = (x.^2+y.^2)<1; v2 = (x.^2+z.^2)<1; Nv = sum(v1&v2); I = Nv/N*4
执行后得到 I 2.6681 。 N 值越大越精确。
（三）生成随机数的函数
在 Matlab 中可产生满足各种分布的随机数： 1．unidrnd(N,m,n) 产生服从 1,2,…,N 的等概 率的离散均匀分布的 mn 阶的随机数矩阵； unidrnd(N) 产生相应的随机数。
它的概率分布为 P k 1 N , k 0,1,2,..., N 。 2． binornd (k, p, m, n) 生成参数为 k,p 的 mn 阶的二项分布随机数矩阵；binornd (k, p)生成 k,p 的相应的随机数。 它的概率分布为 P 1 p, P 1 1 p 。
连续系统动态模拟：研究系统的状态随时间 连续变化的情形，其模型一般是微分方程模 型。建模时首先确定系统的连续状态变量，然 后将它在时间上进行离散化处理, 并由此模拟 系统的运行状态。
离散事件系统模拟：离散事件系统模拟讨论 的是系统状态只在一些离散时间点上，由于随 机事件驱动而发生变化，其模型一般用流程图 或网络来表示。
1. 符号假设
[1] 观察所对目标的指示正确与否？
模拟试验有 2 种结果：每种结果概率都是 1/2。 与投掷一枚硬币方式相仿，出现硬币正面时为指示正 确，反之为不正确。
[2] 指示正确时，我方的射击效果如何？
模拟试验有 3 种结果：毁伤一门火炮的可能性为 1/3(即 2/6)，毁伤两门的可能性为 1/6，没能毁伤敌火炮 的可能性为 1/2(即 3/6)。 与投掷骰子的方式相仿，若出现是１、２、３点，则 认为没击中敌人；若出现是４、５点，则认为毁伤敌人 一门火炮；若出现是６点，则认为毁伤敌人两门火炮。
6. 结果比较
理论计算和模拟结果的比较
分类 项目 模 理 拟 论 无效射击 ０．６５ ０．７５ 有效射击 ０．３５ ０．２５ 平均值 ０．５ ０．３３
E1 0 13 4 3 1 2 =0.5 20 20 20
E=7/20=0.35
5. 理论计算
确 0 观察所对目标指示不正 设： j 1 观察所对目标指示正确 A0：射中敌方火炮的事件；A1：射中敌方一门火炮的事件； A2：射中敌方两门火炮的事件． 则由全概率公式： E = P(A0) = P(j=0)P(A0∣j=0) + P(j=1)P(A0∣j=1) 1 1 1 = 0 0.25 2 2 2 P(A1) = P(j=0)P(A1∣j=0) + P(j=1)P(A1∣j=1) 1 1 1 1 = 0 2 2 3 6 P(A2) = P(j=0)P(A2∣j=0) + P(j=1)P(A2∣j=1) 1 1 1 1 = 0 2 2 6 12 1 1 E1 = 1 2 0.33 6 12
7．poissrnd , m, n 产生 m n 阶的服从参数 为 的帕松分布的随机数矩阵； poissrnd 产生 服从参数为 的帕松分布的随机数。 若离散型随机变量 X 的所有可能取值为 0,1,2,…,且取各个值的概率为
P( X k )
k e
若顾客到达某商店的间隔时间服从参数为 0.1 的指数分 布，则指两个顾客到达商店的平均间隔时间是 10 个单位 时间，即平均 10 个单位时间到达 1 个顾客。 若该商店在单位时间内到达的顾客数服从参数为 0.1 的 帕松分布，则指一个单位时间内平均到达 0.1 个顾客。
k e
k!
, k 0,1, 2, 。
计算机数学模拟
（一）模拟(Simulation)概念
1. 模拟：利用物理的、数学的模型来类比、 模仿现实系统及其演变过程，以寻求过程规律 的一种方法。 2. 思路： 建立一个包含所研究系统的主要特 点试验模型。通过对这个试验模型的运行，获 得所要研究系统的必要信息。
物理模拟
形式：对实际系统及其过程用功能相似的实 物系统去模仿。例如，军事演习、船艇实验、 沙盘作业等。 优点：直观。 缺点：花费较大、周期较长，在模型上改变 系统结构和系数都较困难。且许多系统无法进 行物理模拟，如社会经济系统、生态系统等。
用蒙特卡洛方法计算二重积分 计算二重积分 I f x, y dxdy ，f x, y 0 ，
D
D 它的几何意义是以 f x, y 为曲顶、 x, y D 。
为底的曲顶柱体 V 的体积。 ⑴假设 V 包含在体积已知的立体 内。 ⑵在 内产生一个服从均匀分布的随机数 ，则 P 落在V内 V的体积/的体积。
1. 问题分析 需要模拟以下两件事：
[1] 观察所对目标的指示正确与否？
模拟试验有 2 种结果：每种结果概率都是 1/2。 与投掷一枚硬币方式相仿，出现硬币正面时为指示正 确，反之为不正确。
[2] 指示正确时，我方的射击效果如何？
模拟试验有 3 种结果：毁伤一门火炮的可能性为 1/3(即 2/6)，毁伤两门的可能性为 1/6，没能毁伤敌火炮 的可能性为 1/2(即 3/6)。 与投掷骰子的方式相仿，若出现是１、２、３点，则 认为没击中敌人；若出现是４、５点，则认为毁伤敌人 一门火炮；若出现是６点，则认为毁伤敌人两门火炮。
⑶产生在 内的 N 个服从均匀分布的随机 数，若 Nc 个落在 V 内，则 I 的体积 NV / N 。
算例 I
x2 y 2 1
1 x2 dxdy
曲顶柱体 V ： x2 y 2 1， x 2 z 2 1 ， z 0 。 包含 V 的立方体 ： x 1, y 1 ， 0 z 1 。
数学模拟-计算机模拟
在一定的假设条件下，运用数学运算模拟系 统的运行，称为数学模拟。现代数学模拟都在 计算机上进行，故也称计算机模拟。 优点：计算机模拟可反复进行，改变系统的 结构和系数都比较容易。 实际的一些带随机因素的复杂系统，用分析 方法建模常需要作许多简化假设，使模型与实 际问题可能相差甚远，导致解答根本应用。这 时，计算机模拟几乎成唯一选择。
2. 符号：
i ：模拟的打击次数；
k1 ：没击中敌火炮的射击总数； k2 ：击中敌一门火炮的射击总数； k3 ：击中敌两门火炮的射击总数。
E ：有效射击比率；
E1 ：20 次射击平均每次毁伤敌人的火炮数。
通过模拟（下页） ，从模拟结果可算出：
7 13 4 3 E 0.35 ， E1 0 1 2 0.5 。 20 20 20 20
静态模拟和动态模拟
静态模拟(static simulation)：通过模拟方法 只要求求得问题的结果的模拟方法。 典型实例是用蒙特卡罗（Monte Carlo）方法 求数值积分的值。 动态模拟(dynamic simulation)：通过模拟方 法研究系统的状态随时间连续变化的情形。 动态模拟又可分为连续系统模拟和离散事件 系统模拟。
Hale Waihona Puke Baidu
（二）蒙特卡洛方法
蒙特卡洛(Monte Carlo)方法是计算机模拟 的基础。 它得名于赌城——摩纳哥的蒙特卡洛。
基本思路： ⑴建立一个概率模型，使所求问题的解正 好是该模型的参数或特征量。 ⑵通过模拟一统计试验，统计某事件发生 的比率，当试验次数很大时比率就近似事件发 生的概率。 ⑶利用所建概率模型，求要估计的参数。
3. 计算框图：
初始化:i=0,k1=0,k2=0,k3=0 i=i+1 Y
1，2，3
硬币正面?
N
骰子点数?
4，5 6
k1=k1+1
k2=k2+1
k3=k3+1 Y
k1=k1+1
i＜20? N
E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20 停止
4. 模拟结果
试验 序号 １ ２ ３ ４ ５ ６ ７ ８ ９ １０ 投硬币 结 果 正 正 反 正 正 反 正 正 反 反 ∨ ∨ ∨ ∨ 指示 正确 ∨ ∨ ∨ ∨ ∨ ∨ ３ ６ ∨ ∨ １ ２ 指 示 不正确 掷骰子 结 果 ４ ４ ∨ ∨ ∨ ∨ ∨ ∨ 消灭敌人火炮数 ０ １ ２ ∨ ∨
2 2
。
当研究对象视为大量相互独立的随机变量之和，且其 中每一种变量对总和的影响都很小时，可认为该对象服 从正态分布。 例如，机械加工得到的零件尺寸的偏差、射击命中点 与目标的偏差、各种测量误差、人的身高、体重等，都 可近似看成服从正态分布。
6. exprnd , m, n 产生 m n 阶期望值为 的 指数分布的随机数矩阵；exprnd 产生期望值 为 的指数分布的随机数。 若连续型随机变量 X 的概率密度函数为
5. normrnd , , m, n 产生 m n 阶的服从均值
、方差 的正态分布的随机数矩阵；
normrnd , 产生服从均值 、 方差 的正 态分布的随机数。
正态分布 N ,
2
的密度函数 f x
1 e 2
x 2
（四）实例：炮击效果的计算
在我方前沿防守地域，敌方以一个炮排(含两 门火炮)为单位对我方干扰和破坏。为躲避我方 打击，敌方的阵地进行伪装并变换射击地点。 经观察发现，我方指挥对敌方目标的指示 50 ％准确， 在指示正确时而我方火力有 1/3 的射击 效果能毁伤敌人一门火炮， 1/6 的射击效果能 有 全部消灭敌人。 现希望能用某种方式把我方将要对敌人实施 的 20 次打击结果显现出来，确定有效射击的比 率及毁伤敌方火炮的平均值。
试验 序号 １１ １２ １３ １４ １５ １６ １７ １８ １９ ２０
投硬币 结 果 正 反 正 反 正 正 正 正 反 正
指示 正确 ∨
指 示 不正确
掷骰子 结 果 ２
消灭敌人火炮数 ０ １ ２ ∨ ∨
∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ６ ６ ４ ２ ４ ３
∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨
从以上模拟结果可计算出：
当只知道一个随机变量取值在 a, b 内，但不知道（也 没理由假设）它在何处取值的概率大，在何处取值的概 率小，就只好用 U a, b 来模拟。
4． rand (m, n) 生成 mn 阶的 0,1 上服从均匀 分布的随机数矩阵； rand 生成在 0,1 上服从均匀 分布的随机数。
k!
, k 0,1, 2, ， 其中 0 为常数， 则
则称 X 服从参数为 的帕松分布。 参数为 的帕松分布的期望值为 。
帕松分布在排队系统、产品检验、天文、物理等领域 有广泛应用。
指数分布与帕松分布的关系
如相继两个事件出现的间隔时间服从参数为 的指数分布，则在单位时间间隔内事件出现的次 数服从参数为 的泊松分布。即单位时间内该事 件出现 k 次的概率为 P( X k ) 反之亦然。
e x , x 0, f x 其中 0 为常数， x 0. 0,
则称 X 服从参数为 的指数分布。 参数为 的指数分布的期望值为 1 。
排队服务系统中，顾客到达率为常数时的到达间隔、 故障率为常数时零件的寿命都服从指数分布。指数分布 在排队论、可靠性分析中有广泛应用。
P k C p 1 p
k n k nk
, k 0,1,..., n 。
3．unifrnd(a,b,m, n) 产生 mn 阶在 a, b 上服 从均匀分布 U a, b 的随机数矩阵；unifrnd(a,b) 产生在 a, b 上服从均匀分布 U a, b 的随机数。
创建p_jifen.m如下：
clear; N = 1000000; MATLAB 程序 x = unifrnd(-1,1,N,1); p_jifen.m y = unifrnd(-1,1,N,1); z = rand(N,1); v1 = (x.^2+y.^2)<1; v2 = (x.^2+z.^2)<1; Nv = sum(v1&v2); I = Nv/N*4
执行后得到 I 2.6681 。 N 值越大越精确。
（三）生成随机数的函数
在 Matlab 中可产生满足各种分布的随机数： 1．unidrnd(N,m,n) 产生服从 1,2,…,N 的等概 率的离散均匀分布的 mn 阶的随机数矩阵； unidrnd(N) 产生相应的随机数。
它的概率分布为 P k 1 N , k 0,1,2,..., N 。 2． binornd (k, p, m, n) 生成参数为 k,p 的 mn 阶的二项分布随机数矩阵；binornd (k, p)生成 k,p 的相应的随机数。 它的概率分布为 P 1 p, P 1 1 p 。
连续系统动态模拟：研究系统的状态随时间 连续变化的情形，其模型一般是微分方程模 型。建模时首先确定系统的连续状态变量，然 后将它在时间上进行离散化处理, 并由此模拟 系统的运行状态。
离散事件系统模拟：离散事件系统模拟讨论 的是系统状态只在一些离散时间点上，由于随 机事件驱动而发生变化，其模型一般用流程图 或网络来表示。
1. 符号假设
[1] 观察所对目标的指示正确与否？
模拟试验有 2 种结果：每种结果概率都是 1/2。 与投掷一枚硬币方式相仿，出现硬币正面时为指示正 确，反之为不正确。
[2] 指示正确时，我方的射击效果如何？
模拟试验有 3 种结果：毁伤一门火炮的可能性为 1/3(即 2/6)，毁伤两门的可能性为 1/6，没能毁伤敌火炮 的可能性为 1/2(即 3/6)。 与投掷骰子的方式相仿，若出现是１、２、３点，则 认为没击中敌人；若出现是４、５点，则认为毁伤敌人 一门火炮；若出现是６点，则认为毁伤敌人两门火炮。
6. 结果比较
理论计算和模拟结果的比较
分类 项目 模 理 拟 论 无效射击 ０．６５ ０．７５ 有效射击 ０．３５ ０．２５ 平均值 ０．５ ０．３３
E1 0 13 4 3 1 2 =0.5 20 20 20
E=7/20=0.35
5. 理论计算
确 0 观察所对目标指示不正 设： j 1 观察所对目标指示正确 A0：射中敌方火炮的事件；A1：射中敌方一门火炮的事件； A2：射中敌方两门火炮的事件． 则由全概率公式： E = P(A0) = P(j=0)P(A0∣j=0) + P(j=1)P(A0∣j=1) 1 1 1 = 0 0.25 2 2 2 P(A1) = P(j=0)P(A1∣j=0) + P(j=1)P(A1∣j=1) 1 1 1 1 = 0 2 2 3 6 P(A2) = P(j=0)P(A2∣j=0) + P(j=1)P(A2∣j=1) 1 1 1 1 = 0 2 2 6 12 1 1 E1 = 1 2 0.33 6 12
7．poissrnd , m, n 产生 m n 阶的服从参数 为 的帕松分布的随机数矩阵； poissrnd 产生 服从参数为 的帕松分布的随机数。 若离散型随机变量 X 的所有可能取值为 0,1,2,…,且取各个值的概率为
P( X k )
k e
若顾客到达某商店的间隔时间服从参数为 0.1 的指数分 布，则指两个顾客到达商店的平均间隔时间是 10 个单位 时间，即平均 10 个单位时间到达 1 个顾客。 若该商店在单位时间内到达的顾客数服从参数为 0.1 的 帕松分布，则指一个单位时间内平均到达 0.1 个顾客。
k e
k!
, k 0,1, 2, 。
计算机数学模拟
（一）模拟(Simulation)概念
1. 模拟：利用物理的、数学的模型来类比、 模仿现实系统及其演变过程，以寻求过程规律 的一种方法。 2. 思路： 建立一个包含所研究系统的主要特 点试验模型。通过对这个试验模型的运行，获 得所要研究系统的必要信息。
物理模拟
形式：对实际系统及其过程用功能相似的实 物系统去模仿。例如，军事演习、船艇实验、 沙盘作业等。 优点：直观。 缺点：花费较大、周期较长，在模型上改变 系统结构和系数都较困难。且许多系统无法进 行物理模拟，如社会经济系统、生态系统等。
用蒙特卡洛方法计算二重积分 计算二重积分 I f x, y dxdy ，f x, y 0 ，
D
D 它的几何意义是以 f x, y 为曲顶、 x, y D 。
为底的曲顶柱体 V 的体积。 ⑴假设 V 包含在体积已知的立体 内。 ⑵在 内产生一个服从均匀分布的随机数 ，则 P 落在V内 V的体积/的体积。
1. 问题分析 需要模拟以下两件事：
[1] 观察所对目标的指示正确与否？
模拟试验有 2 种结果：每种结果概率都是 1/2。 与投掷一枚硬币方式相仿，出现硬币正面时为指示正 确，反之为不正确。
[2] 指示正确时，我方的射击效果如何？
模拟试验有 3 种结果：毁伤一门火炮的可能性为 1/3(即 2/6)，毁伤两门的可能性为 1/6，没能毁伤敌火炮 的可能性为 1/2(即 3/6)。 与投掷骰子的方式相仿，若出现是１、２、３点，则 认为没击中敌人；若出现是４、５点，则认为毁伤敌人 一门火炮；若出现是６点，则认为毁伤敌人两门火炮。
⑶产生在 内的 N 个服从均匀分布的随机 数，若 Nc 个落在 V 内，则 I 的体积 NV / N 。
算例 I
x2 y 2 1
1 x2 dxdy
曲顶柱体 V ： x2 y 2 1， x 2 z 2 1 ， z 0 。 包含 V 的立方体 ： x 1, y 1 ， 0 z 1 。
数学模拟-计算机模拟
在一定的假设条件下，运用数学运算模拟系 统的运行，称为数学模拟。现代数学模拟都在 计算机上进行，故也称计算机模拟。 优点：计算机模拟可反复进行，改变系统的 结构和系数都比较容易。 实际的一些带随机因素的复杂系统，用分析 方法建模常需要作许多简化假设，使模型与实 际问题可能相差甚远，导致解答根本应用。这 时，计算机模拟几乎成唯一选择。
2. 符号：
i ：模拟的打击次数；
k1 ：没击中敌火炮的射击总数； k2 ：击中敌一门火炮的射击总数； k3 ：击中敌两门火炮的射击总数。
E ：有效射击比率；
E1 ：20 次射击平均每次毁伤敌人的火炮数。
通过模拟（下页） ，从模拟结果可算出：
7 13 4 3 E 0.35 ， E1 0 1 2 0.5 。 20 20 20 20
静态模拟和动态模拟
静态模拟(static simulation)：通过模拟方法 只要求求得问题的结果的模拟方法。 典型实例是用蒙特卡罗（Monte Carlo）方法 求数值积分的值。 动态模拟(dynamic simulation)：通过模拟方 法研究系统的状态随时间连续变化的情形。 动态模拟又可分为连续系统模拟和离散事件 系统模拟。
Hale Waihona Puke Baidu
（二）蒙特卡洛方法
蒙特卡洛(Monte Carlo)方法是计算机模拟 的基础。 它得名于赌城——摩纳哥的蒙特卡洛。
基本思路： ⑴建立一个概率模型，使所求问题的解正 好是该模型的参数或特征量。 ⑵通过模拟一统计试验，统计某事件发生 的比率，当试验次数很大时比率就近似事件发 生的概率。 ⑶利用所建概率模型，求要估计的参数。
3. 计算框图：
初始化:i=0,k1=0,k2=0,k3=0 i=i+1 Y
1，2，3
硬币正面?
N
骰子点数?
4，5 6
k1=k1+1
k2=k2+1
k3=k3+1 Y
k1=k1+1
i＜20? N
E=(k2+k3)/20 E1=0*k1/20+1*k2/20+2*k3/20 停止
4. 模拟结果
试验 序号 １ ２ ３ ４ ５ ６ ７ ８ ９ １０ 投硬币 结 果 正 正 反 正 正 反 正 正 反 反 ∨ ∨ ∨ ∨ 指示 正确 ∨ ∨ ∨ ∨ ∨ ∨ ３ ６ ∨ ∨ １ ２ 指 示 不正确 掷骰子 结 果 ４ ４ ∨ ∨ ∨ ∨ ∨ ∨ 消灭敌人火炮数 ０ １ ２ ∨ ∨
2 2
。
当研究对象视为大量相互独立的随机变量之和，且其 中每一种变量对总和的影响都很小时，可认为该对象服 从正态分布。 例如，机械加工得到的零件尺寸的偏差、射击命中点 与目标的偏差、各种测量误差、人的身高、体重等，都 可近似看成服从正态分布。
6. exprnd , m, n 产生 m n 阶期望值为 的 指数分布的随机数矩阵；exprnd 产生期望值 为 的指数分布的随机数。 若连续型随机变量 X 的概率密度函数为
5. normrnd , , m, n 产生 m n 阶的服从均值
、方差 的正态分布的随机数矩阵；
normrnd , 产生服从均值 、 方差 的正 态分布的随机数。
正态分布 N ,
2
的密度函数 f x
1 e 2
x 2
（四）实例：炮击效果的计算
在我方前沿防守地域，敌方以一个炮排(含两 门火炮)为单位对我方干扰和破坏。为躲避我方 打击，敌方的阵地进行伪装并变换射击地点。 经观察发现，我方指挥对敌方目标的指示 50 ％准确， 在指示正确时而我方火力有 1/3 的射击 效果能毁伤敌人一门火炮， 1/6 的射击效果能 有 全部消灭敌人。 现希望能用某种方式把我方将要对敌人实施 的 20 次打击结果显现出来，确定有效射击的比 率及毁伤敌方火炮的平均值。
试验 序号 １１ １２ １３ １４ １５ １６ １７ １８ １９ ２０
投硬币 结 果 正 反 正 反 正 正 正 正 反 正
指示 正确 ∨
指 示 不正确
掷骰子 结 果 ２
消灭敌人火炮数 ０ １ ２ ∨ ∨
∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨ ６ ６ ４ ２ ４ ３
∨ ∨ ∨ ∨ ∨ ∨ ∨ ∨
从以上模拟结果可计算出：
当只知道一个随机变量取值在 a, b 内，但不知道（也 没理由假设）它在何处取值的概率大，在何处取值的概 率小，就只好用 U a, b 来模拟。
4． rand (m, n) 生成 mn 阶的 0,1 上服从均匀 分布的随机数矩阵； rand 生成在 0,1 上服从均匀 分布的随机数。
k!
, k 0,1, 2, ， 其中 0 为常数， 则
则称 X 服从参数为 的帕松分布。 参数为 的帕松分布的期望值为 。
帕松分布在排队系统、产品检验、天文、物理等领域 有广泛应用。
指数分布与帕松分布的关系
如相继两个事件出现的间隔时间服从参数为 的指数分布，则在单位时间间隔内事件出现的次 数服从参数为 的泊松分布。即单位时间内该事 件出现 k 次的概率为 P( X k ) 反之亦然。
e x , x 0, f x 其中 0 为常数， x 0. 0,
则称 X 服从参数为 的指数分布。 参数为 的指数分布的期望值为 1 。
排队服务系统中，顾客到达率为常数时的到达间隔、 故障率为常数时零件的寿命都服从指数分布。指数分布 在排队论、可靠性分析中有广泛应用。
P k C p 1 p
k n k nk
, k 0,1,..., n 。
3．unifrnd(a,b,m, n) 产生 mn 阶在 a, b 上服 从均匀分布 U a, b 的随机数矩阵；unifrnd(a,b) 产生在 a, b 上服从均匀分布 U a, b 的随机数。