第22卷第1期海南大学学报自然科学版Vol . 22 No . 1 2004 年 3 月NATURAL SCIENCE JO URNAL OF HAINAN UNIVERSITY Mar . 2004 文章编号:1004 - 1729 (2004) 01 - 0089 - 07计算机模拟在数学建模中的应用欧宜贵,李志林,洪世煌(海南大学信息科学技术学院 , 海南海口 570228)摘要:阐述了计算机模拟在数学建模中的作用,给出了蒙特卡洛方法和离散系统模拟方法实现的具体过程,并通过具体的实例分析,说明计算机模拟方法在数学建模中的有效性.关键词:计算机模拟;数学建模;蒙特卡洛方法;离散系统; Matlab 6. 0中图分类号: O 141文献标识码: A1概述计算机科学技术的迅猛发展,给许多学科带来了巨大的影响.计算机不但使问题的求解变得更加方便、快捷和精确,而且使得解决实际问题的领域更加广泛.计算机适合于解决那些规模大、难以解析化以及不确定的数学模型.例如对于一些带随机因素的复杂系统,用分析方法建模常常需要作许多简化假设,与面临的实际问题可能相差甚远,以致解答根本无法应用,这时模拟几乎成为人们的唯一的选择.在历届的美国和中国大学生的数学建模(MCM)中,学生们经常用到计算机模拟方法去求解、检验等.计算机模拟(computer simulation)是建模过程中较为重要的一类方法(见文献[ 1 ]) .所谓计算机模拟,就是用计算机程序在计算机上模仿各种实际系统的运行过程,并通过计算了解系统随时间变化的行为或特性.它是在已经建立起的数学、逻辑模型之上,通过计算机实验,对一个系统按照一定的决策原则或作业规则,由一个状态变换为另一个状态的行为进行描述和分析.计算机模拟实质上是计算机建模,而计算机模型就是计算机方法和理论(如程序、流程图、算法等) ,它是架于计算机理论和实际问题之间的桥梁.它与数学建模的关系如图 1 :一般说来,在下列情况中,计算机模拟能有效地解决问题.图1计算机模拟流程图1) 难于用数学公式表示的系统 ,或者没有建立和求解数学模型的有效方法;收稿日期: 2003 - 09 - 02基金项目:海南大学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 次抽样, 得到样本值( x i ,φi ) ,i = 1 , 2 , ⋯, n , 统计出满足不等式: x i ≤lsin φi 的次数 m ( m <2下面使用 Matlab 语言(见文献[ 3 ]) 来编程 ,这需建立一个 function yy = simu ( n , L , a)m = 0 ;for k = 1∶nx = unifrnd (0 , a / 2) ;y = unifrnd (0 , pi ) ;if x < = 0. 5 3L 3 sin ( y)m = m + 1 ; elseendend p = m/ npi m = 1/ pn) . 然后进行计算机模拟.M - 文件:输入n , L , a,运行后,即可估算出概率p的估计值m/ n.若选取L∶a= 1∶2 ,我们还可得到π的近似值.Monte Carlo 方法适用范围很广泛 ,它既能求解确定性的问题 ,也能求解随机性的问题以及第1期欧宜贵等:计算机模拟在数学建模中的应用91科学研究中的理论问题.例如,利用Monte Carlo方法可近似地计算定积分,还可以近似求解非线性规划问题,请参看文献[ 4 ,5 ] .注 1 :对随机现象进行模拟,实质上要给出随机变量的模拟,也就是说利用计算机随机地产生一系列数值(称为随机数) ,它们的出现要服从一定的概率分布.目前,经常使用的是按照一定的算法产生的随机数.下面列举的是Matlab中各种分布下产生随机数的命令: 常见的分布函数Matlab 语句均匀分布U [0,1] 均匀分布U [ a , b ] 指数分布 E (λ)正态分布N (μ,σ) 二项分布 B ( N , p) Poisson 分布P (λ) R = rand ( m , n)R = unifrnd ( a , b , m , n) R = exprnd (λ, m , n)R = normrnd (μ,σ, m , n) 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模拟过程流程图第二,计算出按照2种不同方案经T天后企业的92海南大学学报自然科学版2004 年利润值;第三,比较大小,从中选出一个更优的方案.引入下列符号:D :每天需求量;Q1 :方案甲当天的库存量;Q2 :方案甲当天的库存量;S1 :方案甲前1 d的销售量;S21 :方案乙前1 d的销售量;S22 :方案乙前2 d的销售量;S3 :方案甲当天实际销售量;S4 :方案乙当天实际销售量;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 < TQ1 = S1 ; Q2 = ( S 21 + S 22) / 2 ;D = normrnd (135 ,22. 4) ;if D < Q1S3 = Q1 ;elseS3 = D ;endif D < Q2S4 = Q2 ;elseS4 = D ;endL 1 = 53 S 3 - 2. 53 Q1 ; L 2 = 53 S4 - 2. 53 Q2 ;TL 1 = TL 1 + L 1 ; TL 2 = TL 2 + L 2 ;k = k + 1 ;endS1 = 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= 2AB= 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轴,建立空间直角坐标系.由题中所给数据可求得旋转抛物面的方程是:60x=y2+z2.根据光路的几何原理和空间解析几何的知识,易推出结论.线光源发出的光线经抛物面反射后如果能到达B、C两点,则反射点应该在z= 0的抛物线60x=y2上,如图3所示.其中F是焦点.由题意可知 B (25015,1300) , C (25015,2600) ,F (15 , 0) .)ⅰ 能够到达B点的入射光线方程的求法设入射光线的斜率是k .在抛物线上任取一点P ( y02 , y0) , 则直线 B P 的斜率为图3 抛物面反射60。