计算机模拟在数学建模中的应用.doc

井冈山大学大学生创新性实验计划项目计算机模拟在数学建模中的应用（2010.112011.11）指导老师：邓志云 副教授成 员：赵愈旭 吴德龙 刘 波 叶启斌71目录一、计算机模拟建模概述1（一）计算机模拟的定义与特点1（二）计算机模拟模型的建立方法2（三）计算机模拟的应用4二、系统数学模型7（一）系统数学模型7（二）连续系统的数学模型71、微分方程72、传递函数83、状态空间表达式84、系统结构图9（三）离散系统的数学模型91、差分方程92、状态空间表达式103、Z函数104、结构图表示11三、连续系统的计算机模拟12（一）数值积分法121、数值积分的几个基本概念122、欧拉（Euler）法133、龙格库塔(RungeKutta) 法154、变步长的龙格库塔法175、线性多步法18（二）连续系统模拟及模拟程序方法211、连续系统模拟简述212、模拟程序设计方法22（三）计算机模拟算法的MATLAB实现231、数值积分法解常微分方程的MATLAB实现232、其他求解问题的MATLAB实现23（四）典型连续系统模拟实例261、传染病的传播模型262、商品广告模型293、追逐问题模型32四、离散事件系统的计算机模拟37（一）离散事件系统的基本要素37（二）离散事件系统的模型371、到达模型372、服务模型373、排队模型38（三）典型离散事件系统的模拟381、排队系统模拟382、库存系统模拟47（四）随机数与随机变量的生成531、均匀分布随机数532、随机变量的生成533、随机数字和随机变量生成的MATLAB实现55五、模拟结果的校核、验证与确认58（一）VV&A的基本概念58（二）VV&A的一般原则59（三）VV&A技术与主要方法611、系统模型校核的常用方法612、系统模型验证的常用方法62参考文献70一、计算机模拟建模概述计算机模拟是介于运筹学、数理统计与计算机科学等学科之间的一门边缘学科。它利用计算机对系统或过程进行动态仿真，以安全和经济的方法获得动态运行的数量结果，为我们提供可靠的决策依据。因此，它已成为系统分析、战略研究、运筹规划和预测、决策的强有力工具。在历届的美国和中国大学生数学建模竞赛（MCM）中，学生们经常用到计算机模拟方法去求解、检验等。计算机模拟（computer simulation）是建模中较为重要的一类方法。本项目将讨论如何利用计算机技术对连续系统和离散系统进行模拟。由于计算机模拟以重复性运算见长，故它为研究模拟方法提供了极为合适的手段。计算机模拟是一种极为广义的数值计算方法。通过本项目的研究，我们将了解蒙特卡洛方法的思想；掌握对连续系统或离散系统进行模拟的方法；掌握由实际问题怎样建立计算机模拟模型以及应用MATLAB编程语言进行计算。（一）计算机模拟的定义与特点在实际问题中，大量问题很难用数学模型来描述，或有些问题虽建立起了数学模型，但由于模型中随机因素很多，难于用解析的方法求解，这时就需要借助于模拟的方法。例如对于一些带随机因素的复杂系统，用分析方法建模常常需要做许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用，这时模拟几乎成为人们的唯一选择。模拟又称为仿真。它的基本思想是建立一个试验的模型，这个模型包含有所研究的系统中的主要特点。通过这个实验模型的运行，获得所要研究系统的必要信息。一些较简单的问题可以用手工的方法求解，而比较复杂的问题则要借组计算机来进行模拟计算。计算机科学技术的迅猛发展，给许多学科带来了巨大的影响。计算机不但使问题的求解变得更加的方便、快捷和精确，而且使得解决实际问题的领域更加广泛。计算机适合于解决那些规模大、难以解析化以及不确定的数学模型。在一定假设条件下，利用数学运算模拟系统的运行，也可称为数学模拟，现代的数学模拟都是在计算机上进行的，因此称为计算机模拟。计算机模拟是一种人造的实验手段。通过模拟，我们能够对所研究的对象进行类似于物理、化学实验那样的实验，并根据实验数据估算研究对象各种期望的真实特征。它和在现实系统中进行实验的主要区别在于模拟实验依据的是被研究对象的模拟模型以及相应的人造环境，所借助的实验工具是计算机的数值计算与逻辑判断两大功能。由于计算机运行速度高，在很短的时间内即可完成成千上万次的“实物实验”。从而使问题得以轻松解决。与其他一些定量分析相比，计算机模拟技术还有以下几个方面的特点：（1）模拟时间的可伸缩性。使用模拟技术可以把几个月甚至几年的客观系统的活动在几分钟甚至几秒钟内模拟出来。另外，通过适当的安排，模拟也可以对真实系统中变化的细微情节“放宽时间”进行研究。（2）模拟运行的可控性。在模拟模型运行过程中，可以根据需要随时暂停模型的运行，并得到暂停运行时的各项数据，以取得系统发展中任一阶段的有关信息，并且不会影响最终的实验结果。（3）模拟实验的优化性。模拟可以在不同条件下多次重复一项试验。在运行当中，除决策变动外其他一切不变，从而得出相对的最佳决策。这种只改变相关因素而其他条件保持不变的实验优化，在比较分析系统的相关因素是能发挥极大的作用。当然，模拟技术也有着它的局限性，主要表现为：（1）模拟不是最优化技术，它只是针对各个不同的具体决策，通过反复实验比较得出一个较好的结论，但它不能保证是最优的。（2）模拟仅仅是一种评价性的技术，不能自己产生决策，产生方案。（3）在模拟实验运行中，通常要使用大量的随机数，这些随机抽样也会造成模拟的误差，这种误差在其他定量分析技术中一般是不存在的。（二）计算机模拟模型的建立方法模型是系统的一种描述，是为了研究目的而开发的对真实系统进行模仿的一种形式。模型应当是有效的，以便使用模型产生的各种判断，就像用实际系统本身实验所得来的判断一样。图1.1表示一个模型的概念。我们在形式上加以理想化以说明模型通常是真实系统的一种简化。系统和模型都用参数表示他们的特征和属性，真实系统和模型的输入应当基本上相同。然而，二者的输出却不一定完全一致。当系统和模型都可以认为是输出对输入的变换函数时，则一个有效模型的输出可用来预测或推断它代表的真实系统的输出。参数输入参数输入系统（或过程）输出对应对应推断模型输出图1.1 模型的图解概念模型通常具有以下几个特点：（1）它是客观事物的模仿或抽象；（2）它由与被分析问题有关的因素构成；（3）它体现了有关因素之间的联系。为了构造一个客观系统的模型，首先必须对所研究的系统进行观察，获得概念，形成认识。然后再将这种认识抽象出来，并用某种形式的信息表达出来。模型的一般构造过程如图1.2所示。观察认识抽象综合系统信息模型图1.2 模型的形成从信息论的观点来看，模型可以定义为一个信息的整体。在不同的模型中，信息的含义是不同的，例如在物理模型中，信息体现为模型的物理特性；在数学模型中，信息表现为解析式或数值方程的形式；而在模拟模型中，信息一般表现为逻辑流程图或计算机程序的形式。通常，对不同类型系统或过程建立模拟模型的途径也不相同。其中之一是先建立系统或过程的数学模型，然后对数学模型作出相应处置，产生能由计算机运行求出数值解的模拟模型。例如化工生产过程，电机调速控制等连续系统和铁路编组站驼峰解体车列的计算机采样控制等离散事件系统均属于这一范畴，相应数学模型分别为微分或差分状态方程，通常在将它们转化为模拟模型后，采用计算机对其进行模拟实验。建立这类数学模型的任务主要是确定模型的结构和参数。模型建立的方法主要有：1、演绎法建模利用人们已有的关于系统的知识采用演绎的方法来建立数学模型。演绎法是根据已知的原理推导出未知原理的一种推理方法。用这种方法建立模型是，是通过系统本身机理（物理、化学规律等）的分析确定模型的结构（也可能确定某些参数），从理论上推导出系统的数学模型。在这个过程中，要用到一些与系统有关的基本定律、系统本身的物理结构和参数。2、归纳法建模根据对一个已经存在的系统观察、测量所得到的大量输入、输出数据，推断出被研究系统的数学模型。这种利用归纳法得出的数学模型称为经验模型。用归纳法来获得数学模型的基本依据是，在一定的输入激励条件下，系统的动态行为必然体现在它的输入和输出数据之中。利用输入、输出数据所提供的信息，有可能借助于数学方法建立起该系统的数学模型。实际工作中，在建立数学模型时，往往是将演绎法和归纳法结合起来，尽量利用被研究系统的已知规律或原理，采用演绎法确定模型的数学结构。然后，利用实测数据将系统模型中的尚未知道的部分（如模型参数）辨识出来。但是，许多诸如社会经济系统、生态系统、交通运输系统等，本身的机理是很复杂的，人们对这些系统的机理目前还不是了解的很清楚。面对这种情况，往往不得不主要依靠对系统的实测数据，利用归纳法推断出系统的数学模型。系统的数学模型可以按照基本的数学描述分类，如表1.1所示。表1.1 系统数学模型的分类模型分类静态系统模型动态系统模型连续系统模型离散系统模型集中参数分布参数时间离散离散事件数学描述代数方程微分方程传递函数状态方程偏微分方程差分方程脉冲传递方程离散状态方程概率分布排队论应用举例系统稳态解工程动力学系统动力学热传导场计算机数据采样系统简单的随机服务系统将数学模型转换成计算机模拟模型通常有三种途径。第一种是通过直接编程由计算机来求问题的解析解，大多数有代数方程、差分方程、概率分布、排队论描述的数学模型属于这一范畴。第二种是对数学模型进行某种数学处理后再编程由计算机来求问题的解析解，如由微分方程描述的数学模型，只有在对微分方程进行离散化处理后才能编程由计算机来求解。第三种是在对数学模型所作描述的基础上，通过对问题的进一步分析，增补一系列求解的规则与逻辑条件，以流程图的形式建立起相应的模拟模型，然后编程由计算机来求问题的解。当数学模型仅仅是作为对问题的一种描述而无法求解或者是有无穷多个解是，往往有必要通过最后途径得到一些可行解。模拟分为静态模拟和动态模拟。数值积分中的模特卡罗方法就是典型的静态模拟。动态模拟可分为连续系统模拟和离散事件系统模拟，连续系统模拟研究系统的状态随时间连续变化的情形，其模型一般是微分方程模型。建模时首先确定系统的连续状态变量，然后将它在时间上进行离散化处理，并由此模拟系统的运行状态。离散事件系统模拟讨论的是系统状态只在一些离散时间点上，由于随机事件驱动而发生变化，其模型一般用流程图或网络来表示。为了模拟系统必须设置一个模拟时钟将时间从一个时刻向另一个时刻推进，并且可随时反映系统时间的当前值。模拟时间推进方式有两种：下次事件推进法和均匀间隔时间推进法。模拟离散系统常用下次事件推进法，其过程是：置模拟时间初值为0，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件并加入到队列中去；跳到下一事件，计算系统的状态，重复这一过程直到满足某个终止条件为止。要建立一个有效的离散事件系统模拟模型，有必要根据研究目的，并通过对系统的分析从以下几个方面对系统进行抽象：（1）确定系统的范围与组成；（2）定义系统的参变量集合；（3）确定事件类型及其发生时点；（4）设置模拟时间推进结构；（5）确定各类事件与系统状态变量的数学逻辑关系；（6）确定约束条件与目标函数。（三）计算机模拟的应用计算机模拟运用计算机语言编写程序来模拟现实世界。尽管计算机程序只有很少的基本语句，但其所能完成的工作似乎很少受到限制，它可以进行数值计算，可以表示逻辑关系，可以表示活动、事件、进程和过程，可以表示模糊量，可以表示随机量。特别是对于采用数学模型方法解决不了的问题，计算机模拟可以凭借经验数据，依据数量关系，逻辑关系直接模仿客观对象，求得问题的解。因此，计算机模拟在工业、农业、商业、交通运输及军事等各领域的规划、调度、设计、决策及人员培训等方面发挥着越来越重要的作用。计算机模拟应用得如此广泛，主要在于它在解决如下问题上有较大优势：1、无法实施的一类问题在现实世界中有许多问题无法付诸实施来解决，例如预测问题就属于此范畴。假若要预测未来五年或十年的经济计划执行情况，我们不可能让国民经济去实际运行一段时间以取得有关经济指标。但可以构造一个经济模拟模型，利用尽可能收集到的数据，根据不同设想对其进行模拟实验，从而得到各种预测的经济计划指标。又如在研究某地区的抗震能力时，对于某些设施可能遭受的破坏程度，可以在若干假定条件下，应用计算机模拟来进行估计。2、大量方案的比较和选优当某个计划的执行存在多种方案时，例如要设计一个铁路车站，由于线路配置和线路参数的变化，会存在大量的备选方案，若想把全部方案都计算出来进行比较选优，其工作量之大是难以想象的，这一类问题适宜于用计算机模拟来解决。企业的财务计划和作业计划的模拟，以及电子线路设计的模拟等都属于这种情况。在国外，已经研究出了许多相应的计算机模拟软件用来对某些系统的设计进行比较选优，以获得经济合理的方案。3、不易为人们所了解的大型复杂系统有一些复杂系统，如生成系统、交通运输系统、航空航天系统等，其运行情况和使用效果难以用一般的理论分析或数学解析的方法来进行。为此，可以把整个系统分成若干个子系统并形成模型，然后将这些子模型装配成一个总模型，利用模拟实验来检验和判明系统的内部性能和内在联系，作为系统最佳运行的依据。4、有危险的现象对于某些有危险的训练，例如飞行员培训，往往需要利用模拟飞行器在地面上做充分的训练后才能进行实际飞行训练，以保证训练的安全。另一种情况是在研制新产品或者进行新的开拓性投资时，也存在着大量的风险。一旦遇到挫折或失败，会遭到巨大的经济损失，甚至对企业也带来致命的打击。对于这一类问题，要事先进行各方面的调研和分析，依据各种因素的逻辑关系构造模拟模型，并且设法变换有关的数据进行模拟实验，然后制定出对策。在工业发达国家，应用计算机模拟进行投资风险性分析，已经得到了迅速的发展和应用。5、无法重复的现象对于大型工程项目的建设，例如，新建一个海港、一条铁路或一个机场，一旦建成后发现有问题，想再改建或重建，需要花费大量的人力、物力和财力。为此，再设计如何确定设计参数，安排施工进度和步骤，预测和估计交通设施的运行情况，以及完工后对其他方面的影响等等。都要求事先给出答案。显然，把要设计的交通系统构造成模型，利用历史数据的分析结果。对将要新建的交通设施进行反复进行模拟实验，这是一种技术上可行的方法。这样一来，无法重现的现象可以再计算机中反复的重演，为我们提供有关的实验数据。6、开销过高的实物实验某些实物实验，例如新研制的武器系统图，倘若要通过实际实验来对其作战效果和性能指标进行评价和预测，往往开销过于昂贵。军事演习也是如此，一次军事演习仅仅消耗的动力和弹药就相当可观，如果采用计算机模拟实验来带起，就能节约大量的人力和财力。随着计算机模拟智能化的逐步发展，人们进一步研究了怎样用计算机模拟来模仿管理中的决策过程，以便以计算机支持决策。一般说来，当对一个问题进行决策时有许多方法和工具可以利用。图1.3描述了选择决策方法和工具，展开决策的一般性程序。由图可见，计算机模拟技术在管理决策过程中占有相当重要的地位。N问题的提出一般的系统分析是否的有现成的属性模型确定合适的解法、调用相应的标准程序构造新模型系统分析：找出实体、属性、活动、目标函数和可控变量能否新建一个解析模型能否建立一个计算机模拟模型能否用列举法求解可否用计算机进行决定模拟决策构造模拟模型编写和调试程序构造模型编制程序运算直观分析决策采用其他实验等方法选取解答改变参数进行模拟实验确定模型参数、取值范围和精度，设计输入数据，对问题求解验证结果的正确性结果可否接受？决策NNNNNYYYYYY图1.3 决策方法的选择二、系统数学模型系统模拟技术所研究的是将实际物理系统抽象为数学模型，然后将数学模型转化为可以在计算机上运行的计算机模拟模型，最后利用计算机对其可识别的模拟模型进行数值计算的过程。（一）系统数学模型系统数学模型是相对于系统真实现象和真实关系的静、动态特性和有用信息的简化、结合和抽象，是系统的本质描述，可以提供所研究系统的信息。它以数学公式，图表，图线等形式表示，其中采用数学方法来描述实际系统本质特性的方程（组）或不等式（组），称为系统的数学模型。建立系统的数学模型，一般是根据系统实际结构、参数以及计算精度的要求，抓住主要因素，略去一些次要因素，使系统的数学模型既能准确地反映系统的本质，又能简化分析计算的工作。比较常见的方法是机理法和实验法。机理法是根据系统及元部件中各变量之间所遵循的物理、化学等客观定律，列出系统各变量的数学表达式，然后建立系统的数学模型；实验法是采用某些特定的检测仪器与检测方法，在现场对控制信号加入特定的信号（比如单位脉冲信号，单位阶跃信号等），对系统的响应输出进行测量和分析，得到相关的实验数据，从而建立系统的数学模型。系统数学模型可以根据模型的时间集合进行分类为连续系统时间模型和离散系统时间模型。连续系统时间模型的时间用实数表示，即系统的状态可以在任意时刻点获得。通常系统的状态取值连续，可以用常微分方程或偏微分方程表示；离散系统时间模型的时间用整数表示，即系统的状态只能在离散的时刻点上获得。值得注意的是，离散时间点只是表示时间的序次，而不表示具体的时刻点。通常系统的状态仅在特定的时刻点取值变化，可以用差分方程组描述系统。通过建模，得到系统数学模型后，就面临着系统数学模型的求解、分析、校验和验证、修正和应用等后处理。一般，计算机模拟过程是一个需要反复的过程，也只有通过不断的校验、验证，才能找到符合实际情况的高质量数学模型，找到计算准确高效的模型算法。（二）连续系统的数学模型1、微分方程设连续系统的输入为u(t)，输出为y(t)，他们之间的关系可以用系统的高阶微分方程表示，即 (2.1)式（2.1）中，ai(i=0,1,，n-1)，bj (j=0,1,m)为微分方程的常系数。对应的初始条件为： (2.2)通常，式（2.1）右边u(t)的导数项最高阶次为mn，此时式（2.1）、式（2.2）就称为n阶线性系统的数学模型。2、传递函数在线性定常系统中，常用传递函数来描述系统的输入与输出的关系，传递函数也可以扩展到一定的非线性系统中。线性定常系统的传递函数，是指初始条件为零时，系统输出量（响应函数）的Laplace变换与输入量（激励函数）的Laplace变换之比。对方程（2.1）两边取Laplace变换，并假设y(t)和u(t)及其各阶导数的初值均为零，即 （2.3）则得到 （2.4）式（2.4）中，Y(s)=Ly(t),U(s)=Lu(t)，故系统的传递函数G(s)为 （2.5）为了确定复杂的系统的传递函数，需要消去那些组成系统各环节的中间量，才能获得整个系统的输入和输出之间的关系。串联系统传递函数等于组成该系统各环节传递函数之积；并联系统传递函数等于组成该系统各环节的传递函数之和。3、状态空间表达式状态空间的分析方法是贯穿于现代控制理论中重要的分析工具和方法，是建立在状态概念之上，可用于多输入、多输出系统，系统可能是线性的或非线性的，定常的或时变的。系统的状态空间表达式通常由下述的状态方程和输出方程组成。 （2.6） （2.7）式（2.6）由n个一阶微分方程组成，称为状态方程；式（2.7）由l个线性代数方程组成，称为输出方程。其中，X=x1x2xnT为系统的n维状态向量，y=y1y2ylT为系统的l维输出向量，u=u1u2umT为系统的m维输入向量（也称为控制向量）。矩阵Ann为系统的状态矩阵，由控制对象的参数决定；矩阵Bnm为系统的输入矩阵；矩阵Cln为系统的输出矩阵；矩阵Dlm为系统的传输矩阵，若该系统的传递函数为真分式，则矩阵Dlm为零矩阵。4、系统结构图系统结构图（又称为方块图）是系统中每个元件或环节的功能和信号流向的图解表示，各方块表明了系统中各种元件间的相互关系。这种表示方式直观，能够清楚地表明实际系统中的信号流动情况，优于抽象的数学表达式。对于单输入单输出的线性系统，结构图与传递函数之间容易转换；对于多输入多输出或具有非线性环节的系统，也可通过面向结构图的模拟方法得到系统的动态特性。系统动态结构图的组成符号主要有以下4种：信号线。表示系统中信号的流通方向，并标明信号对应的变量。引出点。表示信号从该点取出，从同一信号线上取出信号，其大小、性质完全相同。比较点。表示两个或两个以上的信号在该点进行叠加。方框。表示输入、输出信号之间的动态传递关系。只要依据信号的流向，将各元件的方块连结起来，就能够容易地组成整个系统的结构图，并且通过结构图，还可以评价每一个元件对系统性能的影响。结构图比物理系统本身更容易体现系统的函数功能。结构图还包含了与系统动态特性有关的系统，但是它不包括与系统物理结构有关的信息。因此，许多完全不同和根本无关的系统，可以用同一个结构图来表示。下图就是线性系统状态空间表达式（2.6）（2.7）的模拟结构图。Bu(t)CX +BU +DU+y(t)XAXAX+CD图2.1 线性系统状态空间的模拟结构图（三）离散系统的数学模型1、差分方程离散时间系统的输入与输出都是离散的时间序列，常用差分方程来描述。（2.8）式（2.8）中，输出变量的初始条件为y(0),y(T),Y(n-1)TT为序列的时间间隔，有时T并不真正代表时间点或时刻的含义，而是表示序列的位置点或位置权，通常式（2.8）也写为 （2.9）即 （2.10）描述离散时间系统的差分方程中，有延时、相乘、相加三种基本单元，如图2.2所示。可利用这三种单元的组合描述系统特性。他们的关系式分别为：延迟器：y(n+1)=u(n)加法器：y(n)=u1(n)+u2(n)乘法器：y(n)=u1(n)u2(n)Dy(n)u(n)(a)(b)u1(n)y(n)u2(n)(c)u1(n)y(n)u2(n)图2.2 离散系统中的基本单元(a)延时器；(b)加法器；(c)乘法器2、状态空间表达式状态方程为： （2.11）输出方程为： （2.12）系统的初始状态向量为X(0)，式（2.11）和式（2.12）中，X(k)状态列向量；u(k)输入列向量；Y(k)输出列向量；系统状态转移矩阵；m输入矩阵；C输出矩阵；D传输矩阵。3、Z函数有离散时间系统时域分析可知，系统由激励引起零状态（系统初始条件为零）响应为激励u(n)与系统单位样值响应h(n)的卷积，即y(n)=u(n)h(n)。对其进行Z变换得，式中，G(z)为h(n)的Z变换，类似于连续系统的传递函数G(s)。设式（2.10）表示的系统初始条件为零，即y(k)=u(k)=0 (k0)对其两边取Z变换，得 （2.13）即系统的传递函数为 （2.14）4、结构图表示离散系统结构图表示和连续系统类似，只要将每个环节的传递函数改成Z函数即可。三、连续系统的计算机模拟（一）数值积分法数值积分法不仅方法种类多，而且有较强的理论性，本节由浅入深地介绍几种常见的数值解法。主要为欧拉法、单步法中的四阶龙格-库塔法以及线性多步法。利用数值积分法对连续系统进行数字模拟的基本原理和具体方法，并给出数值积分法中几个常用的算法以及实现这些算法的计算程序，为连续系统模拟提供理论基础。1、数值积分的几个基本概念（1）单步法与多步法数值积分法都用递推公式求解，而递推公式是步进式的，当由前一时刻的数值yn就能计算出后一时刻的数值yn+1时，这种方法称为单步法，它是一种能自启动的算法, 欧拉法是单步法。反之，如果求yn+1时要用到 yn, yn-1, yn-2 ,等多个值，则称为多步法，由于多步法在一开始时要用其它方法计算该时刻前面的函数值，所以它不能自启动，以上所讲的改进的欧拉法就是一个多步法的例子。（2）显式与隐式在递推公式中，计算yn+1时所用到的数据均已算出的计算公式称为显示公式。相反，在算式中隐含有末知量yn+1 的则称为隐含公式。这需用另一个显示公式估计一个值，再用隐式公式进行迭代运算，这就是预估校正法，这种方法不能自启动。用单步法与显示公式在计算上比用多步法和隐式公式方便。但有时要考虑精度与稳定性时，则必须采用隐式公式运算。（3）截断误差一般使用台劳级数来分析数值积分公式的精度。为简化分析，假定前一步得的结果yn是准确的，则用台劳级数求得tn+1处的精确解为 (3.9)与前面的递推公式比较, 欧拉法是从以上精确解中只取前两 项之和来近似计算yn+1的，由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又称为局部离散误差。欧拉法的局部截断误差为不同的数值解法，其局部截断误差也不同。一般若截断误差为O(h2)，则方法称为r阶的，所以方法的阶数可以作为衡量方法精确度的一个重要标志。欧拉法是一阶精度的数值解法，而改进的欧拉法(梯形法)相当于取台劳级数的前三项，故其解为二阶精度。（4）舍入误差由于计算机的字长有限，数字不能表示得完全精确，在计算过程中不可避免地会产生舍入误差。舍入误差与步长 h 成反比，如计算步长小，计算次数多则舍入误差大。产生舍入误差的因素较多，除了与计算机字长有关外，还与计算机所使用的数字系统，数的运算次序及计算f(t,y)所用的程序的精确度等因素有关。（5）数值解的稳定性以上对连续系统的模拟用到了差分方程，把初始值代入递推公式进行迭代运算，如果原系统是稳定的，由数值积分法求得的解也应该是稳定的。但是，在计算机逐次计算时，初始数据的误差及计算过程中的舍入误差对后面的计算结果会产生影响。而且计算步长选择得不合理时，有可能使模拟结果不稳定。以下我们简要地讨论这一问题。差分方程的解与微分方程的解类似，均可分为特解与通解两部分。与稳定性有关的为通解部分，它取决于该差分方程的特征根是否满足稳定性要求。例如，为了考虑欧拉法的稳定性，可用检验方程y=y。其中为方程的特征根。对此方程有要使该微分方程稳定，须使下式成立|1+h|1有此可得|h|2或h8精度阶数234456n-2可见精度的阶数与计算函数值f的次数之间的关系并非等量增加。（4）若在展开成泰勒级数时，只取h这一项，而将h2以及h2以上的高次项略去，可得这是欧拉公式，因此欧拉公式也可看成是一阶龙格库塔公式，它的截断误差正比h2，其精度最低。 如果将h取得很小，能否用欧拉公式得到很高的精度呢？从理论上讲是可以的，但是实际上由于计算机字长有限，在计算中有舍入误差。步长h越小，计算的步数越多，舍入误差的积累会越来越大，故用欧拉法时不能用很小的h。四阶经典龙格库塔格式的MATLAB程序为：function x,y=marunge4(dyfun,xspan,y0,h)%用途：4阶经典龙格库塔格式解微分方程y=f(x,y),y(x0)=y0%格式：x,y=marunge4(dyfun,xspan,y0,h)%dyfun为函数f(x,y)，xspan为求解区间x0,xn，y0为初值，h为步长，x返回节点，y返回数值解format shortx=xspan(1):h:xspan(2);y(1)=y0;for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(n); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x;y=y;例：取h=0.1，求解下列初值问题：并与精确解y(x)=2ex-x-1进行比较。解：在MATLAB命令窗口中执行： clear; dyfun=inline(x+y); x,y=marunge4(dyfun,0,0.5,1,0.1);x;yans = 0 0.1000 0.2000 0.3000 0.4000 0.5000 1.0000 1.1103 1.2428 1.3997 1.5836 1.7974 y1=2.*exp(x)-x-1;y1ans = 1.0000 1.1103 1.2428 1.3997 1.5836 1.7974 (y1-y)ans = 1.0e-005 * 0 0.0169 0.0375 0.0621 0.0915 0.12644、变步长的龙格库塔法以上提到，步长h的选择十分重要，h太大不能达到预期的精度要求，h太小则增加了模拟时间，且有可能影响计算精度，要克服这一问题，就要采用变步长方法，使计算量尽可能的小，且精度也合乎要求。步长的自动控制是根据每一步的误差的大小来实现的。为了得到每一步的局部误差的估计值，可以采用两种不同阶次的递推公式（一般采用低一阶的RK公式，同时计算yn+1并取两者之差），要使计算量最少，可以选择系数，要求两个公式中的Ki相同，使中间结果对两种RK公式都适用，则这两个公式的计算结果之差就作为误差估计。1957年Merson（默森）给出了一个四阶变步长的方法，称为龙格库塔默森法。其四阶公式为 (3.11)计算估计误差的三阶公式如下其绝对误差为这一算法是四阶精度，三阶误差，故称为RKM34法，目前使用较多。其稳定域和一般RK4法相近，缺点是计算量稍大，每步计算5次f值，除用绝对误差外，也可用相对误差，最大相对误差RE为 在每一步计算后，先计算误差，根据误差的大小来控制步长，控制策略如下： （1）当误差ERR大于预先设定的最大允许误差Emax时，则缩小步长，一般将步长减半，并以新步长重新计算以后再比较。 （2）当误差ERR小于预先定的最小允许误差Emin时，步长增加一倍，以新步长计算下一步并计算误差。 （3）如步长小于某一下限hmin时不再减半，以免增加模拟时间，舍入误差过大。 （4）如步长大于某上限hmin（指打印间隔），则不再加大步长。龙格库塔默森法虽然增加了计算量，但在微机日益普及的今天，这一方法是值得广泛采用的。5、线性多步法 上几节我们讨论了单步法，即在计算yn+1值时，只要知道前一步的yn和fn的值就可以进行计算，而每往后计算一步都可以使用前面已经算出的结果。多步法就不同了，在计算yn+1时，不仅要用到前一步yn和fn的值, 还要使用yn+1, yn+2及fn-1, fn-2的值，一般而言，充分利用前面多步信息来计算yn+1，可期望加快模拟速度，也会使精度得到提高，这显然比单步法要好一些，这就是构成多步法计算的基本思想。线性多步法中以亚当姆斯（Adams）法使用得最为普遍，下面介绍这种方法。欧拉法在计算函数y(t)时，在tntn-1中积分，用矩形公式求得。欧拉法计算误差比较大，为了提高计算精度，将矩形改为梯形, 即在求积分时使用下式 (3.12)则整个数值积分公式变成 (3.13)在(3.13)式中，为求yn+1值时，要用到yn+1本身，而这是一个未知数，故此式称为二阶隐式亚当姆斯公式，也称为梯形公式。Adams法的几何意义可用下图表示F(t,y)fn+1fnf1f0=f(t0,y)f(t)t0t1tntn+1th图3.2 用梯形公式进行数值积分在求yn+1，时由于要用到yn+1本身，那就要用迭代法求解，即先猜出一个初值，然后用下式求出再用求出 (3.14)如此重复下去，直到与的差很小为止，此时可以认为，其截断误差为h2。还有二阶亚当姆斯显式公式，其可表示如下 （3.15）（3.16） 表2-1到出了各阶亚当姆斯法计算公式中的系数值，这是一个多步公式，为了计算yn+1,不仅要知道yn,还要知道yn-1,yn-2,因此这些值都要保存，阶数越搞，要用到的值越多，存储量也就大了。多步法的缺点是不能自启动，即开始时要用单步法，然后才用多步法。表3.1亚当姆斯系数表名 称B1B0B1B2B3一阶显式01000二阶显式03/2-1/200三阶显式023/12-16/125/120四阶显式055/24-59/2437/24-9/24一阶隐式10000二阶隐式1/21/2000三阶隐式5/128/12-1/1200四阶隐式19/2419/24-5/241/240在隐式多步法中，一般先用显式多步法计算处值，然后用隐式多步法作一次校正计算，这样使计算过程得到简化。例如，用一阶显式计算初值，即 （3.17） （3.18） 以上计算公式与二阶龙格库塔法公式相同。这种在多步法公式中先用显式公式估计一个值的方法 叫预估校正法，它的优点是在精度相同（四阶），步长h相同时，预估校正法，每步计算两次右函数，而龙格库塔法要计算四次右函数，因此这种方法应用较为普遍。四阶亚当姆斯预估校正法公式为：预估 校正 亚当姆斯的显示公式与隐式公式的特点如下：（1）相同阶数的隐式公式的系数值比显式公式小（一阶例外），因而隐式公式比显式公式精确得多（2）隐式公式的稳定性亦比显式的好。下表给出了Adams法的稳定区。对于同样阶次，隐式的稳定域比显式的要大随着阶数的增加，Adams法的稳定域逐步减小 阶数 1 2 3 4 显式 2 1 6/11 -3/10 隐式 - - -6 -3 （3）从计算上看，显式比隐示计算量小。隐示公式需要计算fn+1,通常需要用Adams显式公式对它提供一个首次近似值yn+1。这种将显式公式和隐示公式联合起来使用以改进稳定性和精度的方法就是上面介绍的预估校正法。亚当姆斯法每步只需要计算一次f 函数值，因而计算量较小，这是线性多步法的共同优点。但它需要知道多个出发值才能计算，而微分方程只能提供一个初值，显然多步法不象单步法那样可以自启动，它必须用其他方法先获得所求时刻以前多步的解。为了保证出发值的计算精度，一般采用比较高阶的单步法，以较小步长计算。此外，多步法改变步长不能象单步法那样方便。四阶Adams预测校正公式的MATLAB程序：function x,y=maadams4(dyfun,xspan,y0,h)%用途：4阶Adams预报-校正格式解微分方程y=f(x,y),y(x0)=y0%格式：x,y=maadams4(dyfun,xspan,y0,h)%dyfun为函数f(x,y)，xspan为求解区间x0,xn，y0为初值，h为步长，x返回节点，y返回数值解format shortx=xspan(1):h:xspan(2);xx,yy=marunge4(dyfun,x(1),x(4),y0,h);y(1)=yy(1);y(2)=yy(2);y(3)=yy(3);y(4)=yy(4);for n=4:(length(x)-1) p=y(n)+h/24*(55*feval(dyfun,x(n),y(n)-59*feval(dyfun,x(n-1),y(n-1)+37*feval(dyfun,x(n-2),y(n-2)-9*feval(dyfun,x(n-3),y(n-3); y(n+1)=y(n)+h/24*(feval(dyfun,x(n-2),y(n-2)-5*feval(dyfun,x(n-1