计算机模拟在数学建模中的应用
第 22 卷第 1 期
海南大学学报自然科学版
Vol . 22 No . 1
文章编号:1004 - 1729 (2004) 01 - 0089 - 07
计算机模拟在数学建模中的应用
欧宜贵 , 李志林 , 洪世煌
(海南大学信息科学技术学院 , 海南海口 570228)
摘 要: 阐述了计算机模拟在数学建模中的作用 , 给出了蒙特卡洛方法和离散系统模拟方法实 现的具体过程 , 并通过具体的实例分析 , 说明计算机模拟方法在数学建模中的有效性. 关键词: 计算机模拟 ; 数学建模; 蒙特卡洛方法; 离散系统; Matlab 6. 0 中图分类号 : O 141
文献标识码: A
1 概 述
计算机科学技术的迅猛发展 , 给许多学科带来了巨大的影响. 计算机不但使问题的求解变
得更加方便、快捷和精确 , 而且使得解决实际问题的领域更加广泛. 计算机适合于解决那些规模大、难以解析化以及不确定的数学模型. 例如对于一些带随机因素的复杂系统 , 用分析方法建模 常常需要作许多简化假设 , 与面临的实际问题可能相差甚远 , 以致解答根本无法应用 , 这时模拟几乎成为人们的唯一的选择. 在历届的美国和中国大学生的数学建模(MCM) 中 , 学生们经常用到计算机模拟方法去求解、检验等. 计算机模拟(computer simulation) 是建模过程中较为重要的一
类方法(见文献[ 1 ]) .
所谓计算机模拟 , 就是用计算机程序在计算机
上模仿各种实际系统的运行过程 , 并通过计算了解
系统随时间变化的行为或特性. 它是在已经建立起
的数学、逻辑模型之上 , 通过计算机实验 , 对一个系统按照一定的决策原则或作业规则 , 由一个状态变换为另一个状态的行为进行描述和分析.
计算机模拟实质上是计算机建模 , 而计算机模
型就是计算机方法和理论 ( 如程序、流程图、算法 等) ,它是架于计算机理论和实际问题之间的桥梁. 它与数学建模的关系如图 1 :
一般说来 , 在下列情况中 , 计算机模拟能有效
地解决问题.
图 1 计算机模拟流程图
收稿日期: 2003 - 09 - 02
1) 难于用数学公式表示的系统 ,或者没有建立和求解数学模型的有效方法;
基金项目: 海南大学 2002～2003 年度教学研究项目“数学建模教育及对学生综合能力培养
的研究与实践”资助 作者简介: 欧宜贵(1965 - ) ,男 , 湖北钟祥人 , 海南大学信息科学技术学院副教授 , 博士.
90
海南大学学报自然科学版
2004 年
2) 虽然可以用解析的方法解决问题 ,但数学的分析与计算过于复杂 ,此时计算机模拟可能
提供简单可行的求解方法;
3) 希望能在较短的时间内观察到系统发展的全过程 ,以估计某些参数对系统行为的影响; 4) 难以在实际环境中进行实验和观察时 ,计算机模拟是唯一可行的方法 ,例如太空飞行的
研究;
5) 需要对系统或过程进行长期运行比较 ,从大量方案中寻找最优方案.
计算机模拟是系统随时间变化而变化的动态写照 , 因此 , 在通常情况下 , 模拟是按时间来划分的. 目前 , 计算机模拟大致可分成静态模拟( static simulation) 和动态模拟( dynamic simulation) .
数值积分中的蒙特卡洛(Monte Carlo) 方法是典型的静态模拟; 动态模拟又分为连续系统模拟和
离散系统模拟. 下面将主要讨论数学建模竞赛活动中经常用到的Monte Carlo 方法和离散系统的
模拟方法. 实际上 , 对连续系统的模拟 , 是将连续状态变量在时间上进行离散化处理 , 并由此模拟系统的运行状态.
2 Mo nte Carlo 方法
Monte Carlo 方法是计算机模拟的基础 ,其历史源于 1777 年法国科学家蒲丰提出的一种计
算圆周率π的方法 ———随机投针法 , 即著名的蒲丰投针问题(见文献[ 2 ]) .
Monte Carlo 方法的基本思想是首先建立一个概率模型 ,使所求问题的解正好是该模型的参
数或其他有关的特征量. 然后通过模拟一统计 , 即多次随机抽样实验 , 统计出某事件发生的百分
比. 只要实验次数很大 , 该百分比便近似于事件发生的概率. 这实际上就是概率的统计定义.
Monte Carlo 方法属于试验数学的一个分支.
例如 , 为了对蒲丰投针问题进行模拟 , 我们先要建立如下的概率模型:
设“X ”是一随机变量, 它服从区间[ 0 , a/ 2 ]上的均匀分布. 同理, φ是服从区间[ 0 , π]上的均匀分布. 按照某种抽样法, 产生随机变量的可能值, 例如进行 n 次抽样, 得到样本值( xi , φi ) ,
i = 1 , 2 , ⋯, n , 统计出满足不等式: xi ≤ i 的次数 m ( m
l
2
x = unifrnd (0 , a / 2) ; y = unifrnd (0 , pi ) ;
下面使用 Matlab 语言(见文献[ 3 ]) 来编程 , 这需建立一个 function yy = simu ( n , L , a)
m = 0 ;
for k = 1∶n
if x
3 m = m + 1 ; 3
else
end
end p = m/ n
pi m = 1/ p
n) . 然后进行计算机模拟.
M - 文件:
输入 n , L , a , 运行后 , 即可估算出概率 p 的估计值 m/ n . 若选取 L ∶a = 1∶2 , 我们还可得到π的近似值.
Monte Carlo 方法适用范围很广泛 ,它既能求解确定性的问题 ,也能求解随机性的问题以及
第 1 期
欧宜贵等:计算机模拟在数学建模中的应用
91
科学研究中的理论问题. 例如 , 利用Monte Carlo 方法可近似地计算定积分 , 还可以近似求解非线性规划问题 , 请参看文献[ 4 ,5 ] .
注 1 :对随机现象进行模拟 , 实质上要给出随机变量的模拟 , 也就是说利用计算机随机地产生一系列数值(称为随机数) ,它们的出现要服从一定的概率分布. 目前 , 经常使用的是按照一定
的算法产生的随机数. 下面列举的是 Matlab 中各种分布下产生随机数的命令:
常见的分布函数
均匀分布 U [ 0 , 1 ]
均匀分布 U [ a , b ]
指数分布 E (λ)
正态分布 N (μ, σ)
Matlab 语句
R = rand ( m , n)
R = unifrnd ( a , b , m , n) R = exprnd (λ, m , n) R = normrnd (μ, σ, m , n)
二项分布 B ( N , p)
Poisson 分布 P (λ)
R = binornd ( N , p , m , n)
R = poissrnd (λ, m , n)
以上语句均产生 m ×n 的矩阵.
3 离散系统的模拟
顾名思义 , 这是对离散系统进行模拟. 离散系统
(discrete system) 是指系统状态只在有限的时间点或可 显然状态量的变化只是在离散的随机时间点上发生.
数的时间点上有随机事件发生的系统. 例如排队系统 ,
假设离散系统状态的变化是在一个时间点上瞬间完成 的.
下面讨论某企业生产的库存系统的计算机模拟方 法 , 这是排队系统的一个实际应用的典型例子.
3 . 1 问题的提出 在销售部门、工厂等领域中都存在
库存问题 , 库存太多造成浪费以及资金积压 , 库存少
了不能满足需求也造成损失. 部门的工作人员需决定
何时进货 , 进多少 , 使得所花费的平均费用最少 , 而收
益最大 , 这就是库存问题.
某企业生产易变质的产品. 当天生产的产品必须
售出 , 否则就会变质. 该产品单位成本为 2. 5 元 , 单位
产品售价为 5 元. 企业为避免存货过多而造成损失 , 拟
从以下 2 种库存方案中选出一个较优的方案:
方案甲:按前 1 d 的销售量作为当天的库存量;
方案乙:按前 2 d 的平均销售量作为当天的库存
量.
假定市场对该产品的每天需求量是一个随机变 量 , 但从以往的统计分析得知它服从正态分布: N
(135 , 22. 4) .
3 . 2 模型的建立 计算机模拟的基本思路:
第一 , 获得市场对该产品的需求量的数据;
第二 , 计算出按照 2 种不同方案经 T 天后企业的
图 2 模拟过程流程图
92
海南大学学报自然科学版
2004 年
利润值;
第三 , 比较大小 , 从中选出一个更优的方案. 引入下列符号:
D :每天需求量;
Q 1 :方案甲当天的库存量; Q 2 :方案甲当天的库存量; S 1 :方案甲前 1 d 的销售量; S 21 :方案乙前 1 d 的销售量; S 22 :方案乙前 2 d 的销售量; S 3 :方案甲当天实际销售量; S 4 :方案乙当天实际销售量; L 1 :方案甲当天的利润; L 2 :方案乙当天的利润; TL 1 :方案甲累计总利润; TL 2 :方案甲累计总利润; T :预定模拟天数.
模拟过程的流程图如图 2.
3 . 3 模型的求解 利用 Matlab 编程来实现这一过程 ,这需要建立如下的 M - 文件:
function[ TL 1 , TL 2 ] = kucun ( T , S1 , S21 , S 22)
TL 1 = 0 ; TL 2 = 0 ; k = 1 ; while k
Q 1 = S1 ; Q2 = ( S 21 + S 22) / 2 ; D = normrnd (135 ,22. 4) ;
if D
S 3 = Q1 ; else
S 3 = D ; end
if D
S 4 = Q2 ; else
S 4 = D ; end
L 1 = 5 3 S 3 - 2. 5 3 Q1 ; L 2 = 5 3 S4 - 2. 5 3 Q2 ;
TL 1 = TL 1 + L 1 ; TL 2 = TL 2 + L 2 ; k = k + 1 ;
end
S 1 = S 3 ; S22 = S21 ; S21 = S4 ;
给出初值 , 反复运行上述程序 , 通过比较即可得出哪一个方案的优劣.
至于排队系统的其他实际应用例子 , 如订票系统、加工制造系统、交通控制系统、计算机系统等等 , 请参看文献[ 6 ] .
第 1 期
欧宜贵等:计算机模拟在数学建模中的应用
93
此外 , 对于更复杂系统的数值模拟 , 如飞机起飞调度问题、公交车的调度问题、北京 SARS
疫情过程分析等等 , 本文在此不作进一步的讨论 , 请参看文献[ 7 ,8 ] .
4 建模实例
2002 年全国大学生数学建模竞赛的“车灯线光源的优化设计”问题是一道从实际问题提
炼、简化而来的数学问题. 由于理论上的困难 , 很难得到满足设计要求的最优长度的线光源. 本
文利用计算机模拟技术 , 可求得满足设计要求的近似最优线光源 , 从而体现了数学建模中计算
机模拟的重要性.
4 . 1 问题重述 安装在汽车头部的车灯的形状为一旋转抛物面 , 车灯的对称轴水平地指向正
前方 , 其开口半径 36 mm ,深度 21. 6 mm. 经过车灯的焦点 , 在与对称轴相垂直的水平方向 , 对称地放置一定长度的均匀分布的线光源. 要求在某一设计规范标准下确定线光源的长度.
该设计规范在简化后可描述如下:在焦点 F 正前方 25 m 处的 A 点放置一测试屏 , 屏与 FA
垂直 , 用以测试车灯的反射光. 在屏上过 A 点引出一条与地面相平行的直线 , 在该直线 A 点的同侧取 B 点和 C 点 , 使 AC = 2 AB = 2. 6 m. 要求 C 点的光强度不小于某一额定值(可取为 1 个单位) , B 点的光强度不小于该额定值的 2 倍(只须考虑一次反射) . 在满足该设计规范的条件下 ,
计算线光源长度 , 使线光源的功率最小.
4 . 2 问题分析 由于线光源是均匀分布的 , 要使线光源功率最小 , 其长度也应该较小. 但若线光源的长度太小 , 有可能出现 C 点的光强度小于额定值; 若线光源的长度过大 , 虽然能同时满足 B 、C 两点光强度的要求 , 但线光源的功率也增大了. 我们的目的就是在 B 、C 两点光强度满足题目要求的情况下 , 求出最优的线光源长度. 又由于到达屏上某一点的光线数目与该点的光
强度成正比 , 因此 , 可以将题中条件转化为:到达 C 点的光线数目不小于某一额定值 , 到达 B 点
的光线数目不小于该额定值的 2 倍.
然后 , 在抛物线上任取一点 , 并利用光路的可逆性 , 分别求出能够到达 B 点和 C 点的入射光线方程. 若入射光线与线光源所在直线的交点的纵坐标的绝对值不大于线光源长度的一半 ,
即与线光源有交点 , 则表示该光线经反射后能够到达屏上的 B 点或 C 点. 这可通过计算机模拟
来实现.
4 . 3 模型的基本假设 1) 线光源看成是无数个点光源叠加而成;2) 不考虑光在抛物面上的折射 ,
并且光在传播过程中 , 其强度不受空气的影响;3) 不考虑车灯前配置镜面对反射光方向的影响.
4 . 4 模型的建立及求解 以抛物面的顶点为原点 O ,对称轴为 x 轴 ,过原点 O 且与线光源平行的直线为 y 轴 , 过顶点且与 xy 平面垂直的直线为 z 轴 , 建立空间直角坐标系. 由题中所给数
据可求得旋转抛物面的方程是:60 x = y + z .
根据光路的几何原理和空间解析几何的知识 , 易推出结论.
22
线光源发出的光线经抛物面反射后如果能到达 B 、C 两点 , 则反射点应该在 z = 0 的
2
抛物线 60 x = y 上 , 如图 3 所示. 其中 F 是焦点.
由题意可知 B ( 25015 , 1300 ) , C ( 25015 , 2600 ) ,
F (15 , 0) . )
ⅰ 能够到达 B 点的入射光线方程的求法
设入射光线的斜率是 k . 在抛物线上任取一点
2
P ( y , y0) , 则直线 B P 的斜率为
图 3 抛物面反射
60