1.1.1模拟技术在许多数学方法中，一般都要用到解析论证和数值计算的技巧。

但是，许多现实的系统是很复杂的，其中的随机性因素往往难以用数学公式表示出来，也就很难使用数学推导或数值计算机的手段来分析系统、预测系统的性能。

因此，产生了对系统进行模拟的技术。

在使用计算机对系统进行模拟之前，一般要清楚模拟的一般步骤和方法。

后面将从较小的模拟实例，对模拟技术进行简要的介绍。

模拟过程的一般过程为：（1）分析问题，收集资料。

需要搞清楚问题要达到的目标，根据问题的性质收集有关随机性因素的资料。

这里用得较多的知识为概率统计方面。

在这个阶段，还应当估计一下待建立的模拟系统的规模和条件，说明哪些是可以控制的变量，哪些是不可控制的变量。

（2）建立模拟模型，编制模拟程序。

按照一般的建模方法，对问题进行适当的假设。

也就是说，模拟模型未必要将被模拟系统的每个细节全部考虑。

模拟模型的优劣将通过与实际系统有关资料的比较来评价。

如果一个“粗糙”的模拟模型已经比较符合实际系统的情况，也就没有必要建立费时、复杂的模型。

当然，如果开始建立的模型比较简单，与实际系统相差较大，那么可以在建立了简单模型后，逐步加入一些原先没有考虑的因素，直到模型达到预定的要求为止。

编写模拟程序之前，要现画出程序框图或写出算法步骤。

然后选择合适的计算机语言，编写模拟程序。

模拟实现的工具较多，如数学软件类：Matlab、Mathematica、Maple、MathCAD等，还有其它的高级语言工具，如：Visual C++、C++ Builder、Delphi、Borland C++3.1等。

（3）运行模拟程序，计算结果。

为了减小模拟结果的随机性偏差，一般要多次运行模拟程序，还有就是增加模拟模型的时段次数。

（4）分析模拟结果，并检验。

模拟结果一般说来反映的是统计特性，结果的合理性、有效性，都需要结合实际的系统来分析，检验。

以便提出合理的对策、方案。

以上步骤是一个反复的过程，在时间和步骤上是彼此交错的。

比如模型的修改和改进，都需要重新编写和改动模拟程序。

模拟结果的不合理，则要求检查模型，并修改模拟程序。

1.1.2模拟时间利用计算机进行模拟时，有两种控制模拟时间的方法。

一种是固定时间增量法，另一种是可变时间增量法，又叫面向事件法。

固定时间增量法，是选用一段合适的时间作单位，然后每隔一个单位时间就计算一次有关参数的值，到达预定的模拟时间后，模拟程序结束。

在编写这种程序时，一般可以建立一个“模拟时钟”变量。

程序的主体框架一般时个大的循环，循环变量，则为模拟时间；在每个循环体内，就是对每个时段作处理。

例如，有些排队论模型，可能就是以每隔一段时间（一天或者一个月）进行处理。

采用可变时间增量法编写的模拟程序也有一个“模拟时钟”变量，但它是在一个事件发生时，“模拟时钟”才向前推进。

需要注意的是，该模拟方法每一步经过的时间是可变的，而且会自动寻找下一个最早使系统状态发生变化的事件。

整个模拟直到“模拟时钟”到达指定的时间长度为止。

可以参考有关离散系统仿真的内容。

1.1.3 模拟语言运用计算机进行模拟，还要选用适当的程序设计的语言。

当然，象C/C++、Pascal 、Fortran 这样的高级语言是可以用于模拟的实现，特别是在面向对象这样的程序设计思想的引入，使得采用这样的语言实现模拟要更方便些。

但是，用这样的语言编写的程序一般都很长，而且编写复杂，调试费时。

因此人们研究初了许多专门用于模拟的语言，如GPSS 、SIMULA 等。

这里本书不讨论其他的模拟语言，而主要讲解如何使用Matlab 语言编写模拟程序。

Matlab 不仅数值计算功能强大，而且由于其语言的简洁和高级，编写的代码少，而且容易调试，实现模拟模型很快。

1.1.4 随机数的模拟计算机模拟主要用于模拟复杂过程或现象的一些方法。

采用计算机模拟一个实验或一个过程，那么用不同的数据重复计算机模拟就能得出统计学结论。

使用这种研究方法得到的结论可能在数学上不很精确，但其精确性对于我们了解所模拟的过程已经足够了。

考虑落在单位区间（0，1）中的一个实数序列。

简单的说，如果这些数是杂乱地分布于整个区间中，且其排列次序似乎也无章法可循，那么这一序列就称为是随机地。

下列序列的数就不是随机的：（1） 该序列的数是单调增加的或单调减小的；（2） 后一个数是关于前一个数的连续函数，如)(1-=i i x f x ；1.1.5 随机数的产生大多数计算机系统都有随机数生成器，MATLAB 也有自己的随机数生成器，但是这些系统产生的随机数一般称为伪随机数，它们不是真正随机的。

有许多产生随机数的算法，下面介绍一种算法还是比较令人满意的。

产生均匀分布于开区间(0,1)中的随机数 ,,21x x 。

算法如下： 取一个整数0l ，使得121310-<<l ， 对 ,3,2,1=i ，计算：)12,7mod(3115-=-i i l l1231-=i i l x这里的i l 是范围在12131-<<i l 中的整数。

初始整数0l 称为这个序列的种子。

1和21474836471231=-这个梅森(Mersenne)质数之间的任何一个整数都可取为种子。

下面就用Matlab 编写实现该算法的程序，函数名为mrand ： function r=mrand global LL=mod(16807*L,2147483647); r=L*4.6566128752459e-10; 调用说明：如果要调用该函数mrand ，需要在程序中添加如下2条命令：global L L=3第一句：global L 声明变量L 为全局变量 第二句：给L 赋予初值下面集中随机数也是常用的，都可以借助上面的函数mrand 实现。

fix:取整函数，向靠近0的方向取整（1） 要产生在(a,b)上均匀分布的随机数x ，则x=(b-a)*mrand+a（2） 产生集合{0,1,2,…,n}中的随机数I ，则I=fix((n+1)*mrand)（3） 产生从j 到k(j ≤k)的随机整数X ，则X=fix((k-j+1)*mrand)+j注意：一般说来直接采用Matlab 中的rand 就可以产生（0，1）之间均匀分布的随机数。

下面通过自定义的Matlab 函数用于模拟随机变量，当然也可以改写为非函数实现。

1.1.6 模拟均匀分布随机变量的函数function r=rnd_u(a,b) %产生在[a,b]间均匀分布的随机数 r=a+(b-a)*rand; return1.1.7 模拟指数分布随机变量的函数function r=rnd_beta(lamada) %模拟指数分布%lamad表示指数分布的参数r = -log(rand)/lamada;return1.1.8模拟正态分布随机变量的函数function r=rnd_normal(arg_mean,arg_segema)%arg_mean 均值%arg_segema 标准差r = arg_mean + arg_segema*randn;%不是randreturn1.2 蒙特卡罗模拟法蒙特卡罗（Monte Carlo）方法，也称为随机模拟（random simulation）。

基本思想：为了解决数学、物理、工程技术等方面的问题，首先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。

1.2.1模拟寻求近似圆周率平行四边形ABCD的面积为422=圆的面积ππ=s=21现在模拟产生在正方形ABCD中均匀分布的点n个，如果这n个点中有m个点在该圆内，则圆的面积与正方形ABCD的面积之比可近似为m/n；即nm nm 44≈⇒≈ππ通过模拟求圆周率的程序如下：n=yesinput('请输入产生点的个数：',500);m=0;for i=1:nif (-1+2*rand)^2+(-1+2*rand)^2<=1 m=m+1; end endmypi=4*m/n注意：以上运行结果如果n 得输入都相同，但结果都一般都会完全相同，这是由于产生的点是随机的，自然结果也就不同了。

1.2.2 用蒙特卡罗法估算定积分例：求定积分31102=⎰dx x 。

31是该定积分的精确解。

思想：对于⎰ba dx x f )(，如果0)(≥x f ，则可以通过模拟的方法计算其定积分。

构造一个矩形包含了曲边梯形，),max(b a d >产生n （足够大）个在举行区域内的点，如果落在由函数)(x f 构成的曲边梯形内的点为m 个，则所求定积分为d a b nm dx x f ba)()(-≈⎰求解程序： n=10^6; a=0; b=1;d=max(a,b)+1; m=0;for i=1:n,x=a+rand*(b-a); y=d*rand; if y<=x^2, m=m+1; end ends=m/n*d*(b-a)运行上面的程序，求得结果为≈⎰dx x 1020.33321600000000这个结果非常接近于真实值。

下面编写一个通用的程序，前提是被积函数0)(≥x f 通用函数：function s=cmmmonitixing(funname,a,b,d,n) if nargin~=5, n=10000;%缺省 endif n<10000,n=10000;%至少产生1万个点 endif nargin<4d=max(a,b)+1;%缺省 endm=0;%计数器for i=1:n,x=a+rand*(b-a); y=d*rand;if y<=feval(funname,x), m=m+1; end ends=m/n*d*(b-a)%计算近似定积分 为了计算上面的定积分，按如下步骤： 第一步：编写被积函数的M 文件 function y=myfunx2(x) y=x^2;第二步：调用上面的通用程序 输入：cmmmonitixing('myfunx2',0,1)ans =0.33680000000000cmmmonitixing('myfunx2',0,1,3,100000)ans =0.33429000000000通过上面的调用可以看出，当n 取得越大，所计算得到的定积分更接近于真实值。

1.2.3 用蒙特卡罗法估计体积有这样一个问题：位于锥面S1:222y x z +=上方和球面S2:1)1(222=-++z y x 内部的区域的体积。

首先分析，这块区域包含于以11,11≤≤-≤≤-y x 和20≤≤z 为界的盒子内，这个盒子的体积为8。

理论上，所求区域的体积为球面S2面积的一般加上圆锥的体积，设为V .1415926.33113142123≈=⋅+⎪⎪⎭⎫⎝⎛=πππV 因此，需要生成在这个盒子内生成随机点，并将落在这块区域内的随机点个数于生成的随机点的总数之比乘以8（盒子的体积），即为所求区域的面积。