计算机模拟技术

1、计算机模拟技术课程名：计算机模拟技术计算机模拟是在科学研究中常采用的一种技术,特别是在科学试验环节，利用计算机模拟非常有效。所谓计算机模拟总是用计算机来模仿真实的事物，用一个模型（物理的-实物模拟；数学的-计算机模拟）来模拟真实的系统，对系统的内部结构、外界影响、功能、行为等进行实验，通过实验使系统达到优良的性能，从而获得良好的经济效益和社会效益。计算机模拟方面的研究始于六十年代，早期的研究主要用于国防和军事领域（如航空航天、武器研制、核试验等），以及自动控制等方面。 随着计算机应用的普及，应用范围也在扩大， 现在已遍及自然科学和社会科学的各个领域。在农业方面，我国从80年代开始进行作物生长发

2、育模拟模型和生产管理系统的研究，目前有一定基础的：在小麦方面有北农大、中科院；棉花方面有中国农业大学、中国棉花所；水稻方面有 江西农科院；在土壤水份、水资源及灌溉方面西北农业科技大学。目前影响较大的有比较成形的有江苏省 农科院。目前的主要成果有：我国主要农作物栽培模拟优化决策系统RCSODS（水稻）和WCSODS（小麦-江苏省农科院）、MCSODS（玉米-河南省农科院）、CCSODS（棉花-中国农业大学）等。计算机模拟特别适合于实验条件苛刻、环境恶劣（如真空、高温、高压、有毒有害的场所）、试验周期长，花费大的场合。农作物的生产系统就很适合于计算机模拟：农作物的生产受各种条件的影响，不同作物、不

3、同品种也 有差异。比如，要想提高一种作物的产量，就先要作试验，通过试验了解这种作物的特性：抗旱性、耐寒 性、对氮、磷、钾哪种肥更有效等。但农业的田间实验不能保证精度（除人为可控条件外，还有许多随机因素）、周期长（周期一年），耗费大。可通过计算机模拟来实现： 先建立这种作物生产系统的数学模型（依靠专业知识或试验数据。一般来说，诸如作物产量和农业环境的关系可用微分方程或其它方程来描述）,通过计算机模拟来找出这种作物的生长与农业环境相互作用的关系，以及各种条件之间的协迫情况。不仅可大大 节省实验经费、加快研究进度 （周期一年的实验结果几秒钟内即可得到），这种模拟软件的开发还可与农业生产管理系统，决策

4、系统相联系，实现对农作物生产的预测、分析、调控、设计的数字化和科学化。作为一门课程，不是研究某个特定系统的模拟问题，而是了解计 握基础知识，掌握建模及动态模拟的一般方法。第一章计算机模拟概述1.1计算机模拟技术研究对象在一个计算机模拟问题中，我们研究的对象是一个系统。系统:一些具有特定功能的、相互之间按一定规律联系着的实体的集合。如作物的生产系统可看作由作物、环境、技术、经济等要素构成的。各要素之间相互影响、相互联系，称为系统的 相关性；一个系统是一个整体，整体内的各个部分不能分割，各因素之间必须相互协调，不能在任何一个环节出问题，才能 使系统达到优良的状态，称为系统的完整性。目标计算机模拟的

5、目标是了解系统的各个实体之间的相互制约关系，从而使系统在预定的目标下达到最优和完善。 如在作物生产系统中，怎样控制、实施各水、肥、栽培技术等，从而使产量最高，以 获得最优的经济效益。方法模拟的方法是先建立系统与环境相互作用的数学模型，用数学模型来类比、模仿现实系统（一个数学模型就是从数学上表达系统各因素之间的数量关系，或各因素之间协调的规则；从整个模拟过程来看就是一个算法， 或一系列数据，这些数据综合描述一个系统过程或现象的重要行为 ），然后在数学模型 和对系统深刻了解的基础上，开发模拟软件，用影响系统目标的因素作为输入，通过计算机技术来表达系 统各因素作用的状态。 从数学的角度来看，模拟的过

6、程就是对数学模型求解的过程，并把系统过程演示出 来。基础知识可见，对一个系统进行计算机模拟，(1) 要对模拟对象有深刻的了解。如：一个公交车调试系统，要编制一个好的调度程序，必须先对现行系统作周密的调查，搞清哪些是影响调度程序的成分及实体，如现有车辆数、每车载客数、每趟车花费的时间、沿途乘客的密集程度、乘客的一般去向、乘客高峰期的人数等。只有经过周密的调查研究，才能 形成一个完整的模拟系统；作物生长模拟系统中，也要搞清影响作物生长的各因素，具备丰富的专业知识，否则不能建立起精确的模型 。在对一个系统的动态特性不完全清楚的情况下，有必要通过实验获取数据， 以用于数学模型的建立。(2) 要有数学知

7、识。一般研究的系统较复杂，不能用简单的函数或方程来描述，要综合使用各种数学 方法，才能使模型准确、可靠。(3) 计算机知识和编程技巧。软件应完整地实现系统的数学描述，输出应直观、形象，如三维可视化输出等。软件的开发是项目的重要工作，也直接影响模拟的效果。软件的编写可使用任何编程语言，如C、VB、Java等。专门的语言如 GPSS, GPSS是面向对象问题的离散事件的专用模拟语言，优其适用于排队 系统。1961年，旧M公司发表GPSS的第一个版本，后来又有其它公司的各种版本。标准的版本有52个模块，每个模块用特定的名称和图形来表示其功能。现在一般常用的语言都有模拟库(已编好的用于实现模拟功能的函

8、数)专门用于模拟软件的开发。1.2系统的分类可模拟的系统各种各样，不同类型的系统用不同的模型来描述。系统的分类方法很多，重要的分类方 法是按系统的状态是否随时间变化.来分：一个行为与时间有关的系统称为动态系统待研系统：静态系统：系统的行为与时间无关； 用静态模型来描述，一般为数学方程、逻辑表达式等。如，电路的布尔表达式；电路中电压与电流的关系；系统的稳态解公式等动态系统：连续系统：系统的状态随时间连续变化；常用微分方程来描述；方程对所有时间点有效。如，卫星运行轨道，作物生长量等确定性系统：系统的输出完全由其输入来描述，即系统输入与输出按某种规则对应集中参数模型：用常微分方程来描述，即方程中的导

9、数不是偏导数分布参数模型：用偏微分方程来描述，但一般用集中参数模型近拟表示随机系统：系统的输出是随机的，有规律的存在一族随机变量，且随机变量序列与时间有关(随机过程)。在确定性系统的模拟中使用随机变量的研究方法 称为家特卡罗方/离散系统：系统状态的变化只在离散的时间发生；动态方程只在离散点上有效； 一般为随机系统。如，库存问题；企业的管理系统等 离散时间系统：时间步长固定；常用差分方程来描述离散事件系统：用事件来表示系统在时间间隔内的变化；常用概率模型来描述1.3建立数学模型对一个系统，确定其类型后有助于选择合适的方法建模，建模的一般步骤可分两个大的阶段：(1) 实质内容模型阶段：首先对模拟对

10、象进行调查(了解系统，搜索模拟所需的信息 )、实验(参数估讨等统计推断方法，确定参数及参数的敏感性 )、分析(将信息分类、量化，确定描述系统的规则 )，尽可能全 面地掌握系统的基本特性、运动规律以及中间状态 (最好是有对系统有深入了解的专业人员参与 )，通过分 析和逻辑推理，揭示系统内的规律。(2) 形式数量阶段：在调查、实验、分析的基础上，进一步揭示系统内部的数量关系，并对其进行数 学处理，即对系统用数学形式来描述：用变量描述系统状态，用各种数学方程定量表示各变量之间的相互 联系，用递归方程描述系统状态的发展趋势。多数情况下，建立一个好的(能真实的描述系统， 有代表性，能准确的模拟系统的数量

11、信息，与实际系统较吻合)数学模型不易，优其是复杂多变的系统或系统本身的特性尚不完全清楚的情况(如对农作物开发新品种)。所以在数学模型建立后， 必须进行模拟验证，如与真实系统相差较大， 则要重新建模或修改模型。 所以一个好的数学模型必须经过多次模拟，不断修改、完善，才能得到。一个数学模型是描述系统行为的一个算法或一系列方程，工程中对系统建模应先建立系统的需求规格 说明，在制定模拟规划之前予以充分讨论，过程中需考虑以下因素：(1) 模型中需考虑哪些有效的因素：一个系统可能有不同的行为(如一个作物生产系统中有作物的长势、产量、质量、病虫害等)，但不一定所有这些因素都要建立在模型中。模型中需要考虑的因

12、素是：真正 能简化系统的建模；对系统易于建模、测试和维护；使用较少的计算资源；对研究的系统有直接作用的。(2) 建模细化到什么程度：根据系统的需求确定细化的程度，是只须建立一个简单的模型，还是要对 系统行为精确描述。在模型的准确性和花费之间求得平衡。(3) 与系统有相互作用的哪些外部环境 考虑在建模中：如在作物生长模型中，气候、水、肥等。(4) 在建模中采取什么技术：首先，是基于物理的方程，还是基于测试数据。如果基于物理的方程， 是用微分方程，还是差分方程，考虑不考虑随机因素等。这个问题常由专业人员根据系统的专业知识来决 定；如果基于试验数据，则建立经验方程。(5) 建模时必须获得什么样的数据

13、：如林木生长模型中，需要胸径、树高、材积等数据。另外还需考 虑这些数据的方便的输出形式，以及模型需要的计算资源，如模型需要占用的内存、磁盘空间、消耗的 CPU时间等。(6) 建模和测试模型需多长时间，多少人力、物力、财力：随着模型复杂性的增加，成本也会增加。 一般可从简单开始，随着应用逐步完善。(7) 如何验证和确认模型：必须确认建立的模型被正确实现，模型所描述的行为与真实系统匹配到可 接收的程度，才能有价值。那么以什么标准来衡量。现在已有验证、确认与签定(VV&A)技术来证明和验证模拟的精确性。以上问题应列在系统需求说明书中，并应用于最高层次：在实际的模拟问题中，可能将系统分解为子

14、系统、组件，在对各组件、子系统建模时仍须遵从以上规则。模拟步骤：(1) 明确系统(2) 建立模型(3) 模型变换(4) 软件设计开发(5) 测试检验(6) 评估、对结果的评价和分析一个模拟项目中各项工作的过程应该是一个迭代过程，如图1-1所示。下面通过实例来说明模拟过程。假设两机相距10公里以内可实施攻击,1.4应用举例1.4.1两物体追逐问题。设有一架歼击机追踪一架敌方轰炸机, 且须在12分钟内完成追击任务， 否则认为追击失败。 设两机初始位置如图 1-2所示。问题：对轰炸机的任 一条特定航线，模拟歼击机的追击过程。图1-2两机追踪模拟分析：在这个系统中，是对指定的轰炸机的一条航线而言，模拟

15、歼击机的追击情况：歼击机按什么航 线飞行，何时完成追击任务。能否歼灭敌机由在规定的时间内两机随时间变动的距离而定，所以实施功击 的时间、距离是要输出的数据。在时刻t两机距离由歼击机在 t时的位置决定，而t时的位置依赖于其速度和航向。为使模型简单，作如下简化：(1) 设两机在同一平面飞行。于是三维问题转化为二维问题。(2) 设歼击机的速度(必须考虑的因素)Vf是常数(20km/mim)。变速须用微分方程描述，而常速即可用 一般方程表达，求解更简单。(3) 设歼击机的航向(必须考虑的因素)在 t (设1分钟)改变一次，而在1分钟以内操作不变。这样曲 线运动即变为折线运动。(4) 轰炸机航向(航线)

16、可任意，比如现给定一条航线，如表 1-1所示。表1-1轰炸机航线t012345678901112Xb (t)809099108116125133141151160169179180Yb (t)50484541353227212225293033(5) 歼击机初始位置：Xf(0)=0 , Yf (0)=0(初始条件)建立数学模型：设在t时刻两机位置为(XF(t),YF(t)、(XB(t),YB(t),两机连线与水平线夹角为。，则1分钟后歼击机 的位置为：XF(t+1) = X F(t) +Vf Cos。(1.1)YF(t+1) = Y F(t) +VFSin 0(1.2)由(XF(t),YF(t

17、)计算(XF(t+1),YF(t+1)时涉及角0 ,从图中看出：Sin 0 = Y b (t) - Yf (t) / D(t)(1.3)Cos。= X b (t) - X f (t) / D(t)(1.4)而D(t) = SQRY B(t) - Y f (t) 2 + X b (t) - X f (t) 2(1.5)式(1.5) - (1.1)即是这个问题的数学模型。先算出两机之距离， 不断判断是否在12分钟内到达可攻击的距离之内。程序流程图如图1-3所示。此例是连续系统，确定性模型，用一组方程来描述。下例是一个离散事件系统，是随机系统，用概率 模型来描述。计算新位置计算D（t）t = t+1

18、结束图1-3追击模拟流程图VB程序见实例。1.4.2库存问题。商业部门为了合理利用有限的流动资金，每项商品都要在库存与销售之间求得平衡: 库存量太大会增加管理费用、积压资金；库存量太小又可能造成缺货，也会造成销售损失、信誉损失。所以，当库存量不满足某一时段的顾客需求时， 就要到厂家订货。 这就需要采取一种策略： 当库存量（比如布 匹）降到P匹布时（称为定货点）则向厂家订货，定货量为 Q匹（称为定货量）。假设现在有5种策略（方案）。从中选择一种使花费最少。表1-2库存策略策略PQ11251502125250315025041752505175300已知条件：（1）从发出订货单到收到货物需3天，即

19、第i大订货，第i+3天收到。（2）每匹布的保管费为 0.75元，缺货损失为1.8元/匹，订货费（包括手续费、采购差旅费及其他费用） 为750元。 需求量为一个0-99之间的均匀分布随机数。（4）原始库存量为115匹，并设第一天没发出订货。分析：库存问题是商业上的一个重要问题在数学上有专门的模型研究，存储论也是运筹学的一个重要 分支。库存问题用计算机模拟最合适，若通过销售来找最优方案，势必造成经济损失。这是典型的离散事 件系统，用概率模型来描述。这里已给定概率模型：X U 0,99八0分布函数为：F(x) = x / 99以下框图对每种策略模拟180天,即密度函数为：p(x) = 1/990&l

20、t; x< 99x<00 < x < 99x>99选出效益最好的方案，如图 1-5所示。图1-5库存问题模拟程序框图C程序见实例。);以上两个例子一个连续系统，状态随时间连续变化，用方程来描述；一个离散系统，到货和销售都按 一些离散的步骤进行，存在随机性，只能用概率模型来描述。但解决问题的过程、分析方法类似：先分析 系统，使对系统有充分的了解；建立数学模型，设定一些初始条件(如t = 0时的状态，开始时的库存量系统状态的变化对应一组方程或一组规则，随着时是的变化，系统状态改变，当一个周期结束时，收集模 拟过程的统计数据 (即问题的解)。如果程序设计的好，就会使模拟

21、非常逼真，就象一个真正的系统在演示 一样。从理论上讲，任何问题都可用计算机模拟。计算机模拟具有经济、安全可靠、周期短等优点。对任何 问题，只要建立起数学模型，改变参数值及变量值，就可模拟各种情况下的系统运行情况，从模拟输出的结果，可分析系统内各因素的权重及其制约关系，帮助决策者作出合理的决策，克服盲目性，使系统在实 际运行过程中取得最好的效益。1.4.3传染病传播问题。假设某一地区有一种传染病正在流行，那么政府、医务部门就要采取措施来 控制这种疾病的传播，要使采取的措施能够有效的控制传染速度.，就必须先搞清被传染的人数跟哪些因素 有关，被传染的人数是一个什么样的发展趋势，从而来有效的预测和控制

22、，下面建立描述被传染人数的数 学模型。传染病的传播涉及因素很多，不可能通过一两次简单的假设就能建好完善的数学模型。这里的作法是:先作出最简单的假设，看会得到什么样的结果，然后对不合理的地方再行修改，逐步得到较满意的模型。先讨论一个粗略的模型。模型I:假设：(1) 每个病人在单位时间内传染其他人的人数为一个常数后。(2) 一人得病后，经久不愈，该人在传染期内不会死亡。记时刻t的得病人数为y(t)，开始模拟时有y0个传染病人，则在 t时间内增加的病人人数为y(t+ t) - y(t) = k 0 y(t) t由导数定义得(在假设(1)、(2)下的数学模型)：-dy/dt = k 0 y(t)(1.

23、6)y(0) = vo初值问题(1.6)的解为：y(t) = y 0 ek0 t(分析：)这个结果表明，得病人数将按指数形式无限增加，当t8时，y(t)roo,显然与实际不符,说明上面的假设条件不合理。 事实上，一个地区的总人数大致可视为常数 (不考虑疾病传播期间出生的、迁出的、死亡的)，所以一个病人在单位时间内能传播的人数k0是在改变的：在初期，k0较大，随着病人的增多，健康人数的减少，被传染的机会也将减少，所以k0将变小。对上述模型进行改进：记 t时刻的健康人数人 S(t)，当总人数不变时，k0随S(t)减少而减少。模型II:假设：(1) 总人数为常数n,七时刻的健康人数为S(t)，得病人

24、数为Y(t),则(1.7)k(医学上称为传架理度.)(1.8)(1.9)(1.10)Y(t) + S(t) = n(2) 单位时间内一个病人传染的人数与当时的健康人数成正比，比例系数为(3) 同模型I的假设(2)。由假设(2), (1.6)式中的k。应为k S (t),即:/ dy(t) / dt = k S(t) Y(t)、y(0) = y 0将(1.7)式代入(1.8)式，得(上述假设下的数学模型)：:dy(t) / dt = k Y(t) (n- Y(t) y(0) = y 0初值问题(1.9)的解为：Y(t) = n / 1+(n / y 0-1) e k n t (分析:)由(1.1

25、0)式得:dy/dt = kn 2 (n / y") e * n t / 1 + ( n / y 0-1 ) e k n t 2(1.11)(1.10)式是被传染人数随时间变化的关系；(1.11)式是被传染病人的变化率与时间的关系，如图 1-6和图1-7所示。dy/dtyotl图1-6病人人数变化曲线图1-7病人变化曲线这个模型可预报疾病高峰到来的时间：令d y / dt = 0 ,得极大值点:ti = ln ( n / y o - 1 ) / k n(1.12)由(12)式可知，当传染强度k或人数n较大时，t1变化较小，表明传染高峰到来的快，这与实际情况相吻合。但由(1.10)式知

26、当t8时，y(t) r n,这意味着最终人人都被传染，这与实际不符，原因在于假设不合理。再改进：可将人员分为三类：第一类为传染者(y);第二类为易受传染者(s)，即这一类是非传染者，但能得病而成为传染者；第三类为除前两者之外的人(r),包括患病死去的、痊愈后具有长期免疫力的、被隔离的。用y(t)、s(t)、r(t)分别表示这三类人数。模型III :假设：(1)总人数为常数n,则：y(t) + s(t) + r(t) = n(1.13)同模型II的假设(2)(3)单位时间内病愈免疫的人数与当时病人的人数成正比，比例系数为m (恢复系数)由假设，有dr / dt = m y(t)由于引入了 r(t

27、)，贝U模型II的方程(1.8)应改为:(1.14)dy / dt = k S (t) y (t) -dr / dt(1.15)(1.15)式表示单位时间内病人人数的增加应等于被传染的人数减去病愈的人数。从(1.13) - (1.15)中消去dy,并设 S (0) = So, r (0) = ro 得dS/dt = - k S(t) y(t) dr/dt = m y(t) y(t) + s(t) + r(t) = n S(0) = S0 r (0) = r 0(1.16)模型(1.16)较好地描述了传染病传播问题。通过以上实例，对数学模型的建立就有了一个基本的思路，对计算机模拟技术、模拟的过程

28、、问题的 方法步骤也有了一个概括的了解。后面两章下分别对连续系统和离散系统讨论基本的模拟方法。 对例1.4.1中的飞机追击问题采用极坐标(r, 0 ),相应的方程为:dr/dt = Vb Cos 0 -Vfr - d 0 /dt = -Vb Sin 0其中r为两机距离，0为两机连线的角度。编程模拟两机追击情况。(2)球摆：如图1-8所示，设绳长为L,夹角为0 ,球 度。(0) = 0 ,初始偏角0 0 ,确定摆动周期。提示：影响摆动周期的因素： 摆球重力； 摆球质量； 摆球尺寸相对于绳很小，故可视为质点建模； 绳子质量忽略； 磨擦力忽略； 空气阻力不考虑。建模：影响摆球运动的重力(使摆球回到中

29、心位置的力 力与绳子平行，它不影响摆球运动，忽略。 F = -mg Sin 0 其中a是摆球的切线加速度，由于 L = 0 ”,故得运动方程:0 " = -L/g - Sin 00 (0) = 0 0 , 。 (0) = 0m,初始速 /A 0图1-8球摆)F的方向与绳子垂直；对绳子产生拉紧力的重 将牛顿定律 F = ma代入上式得：a = -g Sin 0 ,第二章连续系统的计算机模拟技术连续系统的状态随时间连续的动状变化，这种变化依赖相关的因素，且遵从一定的规律，所以一般来说可用数学方程描述(代数方程、微分方程，此外还有传递函数、结构图及状态方程等)。问题：对数学模型怎样求解，求

30、解过程也就是模拟的过程，即程序的算法。对前面的飞机追击的问题，其模型是一组代数方程，而且是显式.的，所以，作法是：假设 1分钟改变一次航向，每隔 1分钟计算一次系统的状态，而且 由模型知道了歼击机前一分钟的状态就可以算出后一分钟的位置。这样，通过迭代，就可求出问题的解。对一般的模拟问题，模型用微分方程表达(方程中除含有自变量外， 还有导数或偏导数 一分布参数型可 通过一些变换转换为集中参数型)，而对一个常微分方程(只有线性的或几种特殊的能求出通解及特解，即 解析解)一般来说求解析解是不可能的，(即使能求出解析解)在计算机上求解要用 数值积分法：将连续的系统离散化，把微分方程转化为差分方程，通过

31、迭代运算，求出问题的数值解。通过正确的控制步长、选择 合理的算法即能达到要求的精度，这就是连续系统模拟的主要技术。数值积分法种类很多，本章介绍几种简单常用的方法。2.1 欧拉(Euler)法设连续系统的数学模型为f (t,y)可近似的看成常数f(tn,yn),这样可用小矩形的(2.2)dy/dt = f (t,y)(2.1)L y(0) = y 0为了模拟系统状态 y随时间t的变化，需求解微分方 程(2.1)的数值解(不是解析解)，为此，把模拟周期分为 若干小区间，比如分为N个相等的小区间，如图 2-1所示。只须在每个离散的点上计算系统的状态。在区间(tn, tn+1)上求积分，得tn -1y

32、(t n+1) - y(t n) = tn f (t, y)dt积分的几何意义是小曲边梯形的面积。(tn, tn+1)之间的tn+1时的积分值为如果小区间取的足够小，则在区间 面积代替小曲边梯形的面积，于是得在y(tn+1 ) 2 yn + f (tn,yn) - h = y n+1其中h = T /N ,将(2.2)式写为差分方程形式，得yn+1 = yn + f (tn,yn) hn = 0,1,2,(2.3)这就是欧拉公式：它是一个递推的差分方程， 任何一个新的数值解 yn+1均基于前一个数值解 yn以及导 数值f(tn,yn)求得，只要给出初始条件 yo及步长h,就可算出y1,由y可算

33、出 住 如此迭代计算事,直到满足所需计算的范围。欧拉法也叫折线法，特点是方法简单，计算量小，但计算精度也底。扩展I:改进的欧拉法(梯形法)为了提高计算的精度，可对欧拉法进行改进：从几何意义上看到，用小梯形的面积来代替小矩形的面积，必能提高精度。这样(2.3)式即可写成：yn+1 = yn + h / 2 - f (tn,yn) + f (t n+1,y n+1)=yn + h / 2 - f n + f n+1(2.4)但(2.4)式是隐式公式，即公式右端含有yn+1,而这是未知的待求量，故梯形法不能自行启动运算，要借助于其它算法：如用欧拉法算出yPn+1作为梯形法中ycn+1的预估值，再进行

34、计算，这样，公式可写为：预估：yPn+1 = yn + f (t n,yn) - h校正：yCn+1 = yn + h / 2 - f (t n,y n) + f (t n+1,yPn+1) *(2.5)=yn + h / 2 - f n + f Pn+1(2.5)式也称为预估校正法，计算量稍大(需要附加的计算)，但精度要比欧拉法高，稳定性好。扩展II :多个输入的多阶系统在(2.1)式所表示的数学模型中，y是系统的状态，称之为状态罂量，它随时间t而连续变化。(2.1)式是 最简单的情况，只有一个状态变量，而它直接依赖时间而变化，在实际模拟中，状态变量可能不只一个， 而影响系统状态的因素更多，

35、这些因素是随时间变化的同，而系统状态的改变也正是由于这些因素随时间 变化而变化。如在林木的生长系统中，要考察的状态可能有树高、胸径、材积等；在作物的生长系统中， 考察的状态可能有株高、叶子的片数、叶子的宽度等，而影响系统状态的因素如水、肥、气候、土质、病 虫害等，而这些因素都随时间而变化，如水份会随时间而被蒸发，肥会随时间被作物吸收，气候等自然环 境更是随时间变化无常。所以与(2.1)式相应的数学模型应为：J dyk / dt = f (x1(t), X2(t), -, xm(t), y 1, y2,yq)k=1,2,q(2.6)I y k (0) = y k, o其中y k为系统状态变量，X

36、i为系统输入信号，这样的系统称为具有m个输入、q阶系统。(2.6)式是q个微分方程的方程组，q个初始条件。可以把(2.6)式表示成向量的形式，就和(2.1)式在形式上一致了：记 Y为q个状态变量和向量，X为m 个输入信号向量，即：Y = (y1(t), y2(t),-q yt) , X = (X1 (t), X2 (t), -，湖)则(2.6 )式写为:/ dY/dt = f (X(t), Y)(2.6)'L Y(0) = Yo对于(2.6)式所表达的数学模型，与(2.3)式相对应的欧拉公式为：J y k, n+1 = y k, n + f(X 1(t n), X 2 (t n),，X

37、 m(t n), y1, n, y2, n,，yq, n) ' h(2.7)k=1,2，,q或写成向量的形式：Yn+1 = Yn+f (X n , Y n) - h(2.7)'例2.1设两种物质A和B合到一起产生化学反应，生成新物质 C,假设1克A和1克B结合能产生 2克C,形成C的速率与 A和B的数量乘积成正比。同样C也可分解为 A和B, C分解的速率正比于C的数量。问题：在给定 A和B的数量后，模拟有多少 C物质产生出来，以及达到稳定的时间。解：在任何时刻，设 a, b, c分别是A, B, C的数量，则它们增加和减少的速度服从下列微分方程：'da / dt = k

38、2 C - k 1 a - b< db /dt = k2 C - k1 a b(2.8)、dc / dt = 2 k 1 a b - 2 k2 C其中k1和k2是比例常数(实际问题中k1和k2会随温度、压力等发生变化，但在模拟过程中为简化模型， 可 视为常数)，是模型中的参数，-k1 ab中的负号是因为 A和B是减少的过程。假设模拟从 t = to( 一般取0) 开始，使t以At时间间隔增加(由步长 t和模拟周期可定出所分时间段数 )，则相应的欧拉公式为(2.9)式。 给出常数k1和k2值以及A、B、C的初始数量，即得迭代公式：a n+1 = a n + (k2 Cn - k an ,

39、bn) , tb n+1 = b n + (k2 Cn - k 1 an ' bn) ' t<C n+1 = C n + (2k 1 an ,加-2k2 Cn)， t(2.9)a (0) = a 0 , b (0) = b o , c (0) = 0i k1 = k 1,0 , k2 = k 2, 0由(2.9)式，从 t = 0 开始，由 a 0、b0、c 0 可算出 a 1、b、c 1,又可算出 a 2、b2、C2, (a 2 = a(2 t)图2-2例2-1模拟流程图以At为间隔，进行N = T/次计算，就可算出周期T的系统状态，从而得出模拟结果。根据数学模型(2.

40、9)就可设计编写模拟程序，可选用任何编程语言，在一些流行的语言中，对常用的模 拟算法(比如后面要介绍的龙格-库塔(Runge-kutta)方法)都有相应的用于模拟计算的软件包。public class Rate(static double k1,k2,h,t;开始static double A=new double53;static double B=new double53;static double C=new double53;static int i;public static void main(String args)(A1=100.0;B1=50.0;C1=0.0;t=0;h=0

41、.1;k1=0.008;k2=0.002;for(i=1;i<53;i+)System.out.println(i+" "+t+", "+Ai+","+Bi+","+Ci);strut(i);static void strut(int i)Ai+1=Ai+(k2*Ci-k1*Ai*Bi)*h;Bi+1=Bi+(k2*Ci-k1*Ai*Bi)*h;Ci+1=Ci+2.0*(k1*Ai*Bi-k2*Ci)*h; t=t+h;欧拉法是最简单的数值积分法，在介绍更好的数值积分法之前，先介绍几个关于数值积分的基本概念。

42、1、单步法与多步法数值积分法是用递推公式求解，如果仅由前一时刻的数值yn就能算出后一时刻的数值yn+1，则称为单步西，反之，如果求yn+1时需用到yn, yn-1, yn-2,等多个值，则称为多步西.。如，欧拉法是单步法，而 改进和欧拉法是多步法。单步法由递推公式自身就能启动运算(由初值即可算出y1,由y1可算出V2,),所以它是能自启动的算法；多步法在开始的时候要先用其它的方法计算该时刻前面的函数值，所以不能自启动运算。一般来说，由于多步法更能充分利用多个时刻的信息，所以模拟速度快，精度高，但计算量要大一些，后面还要介绍多步法的算法。2、显式与隐式在递推公式中，计算yn+1时所用到的数据均已

43、算出的计算公式称为显式公式.，相反，在算式中隐含未知量yn+1称为隐式公式。如欧拉法中递推公式是显式的，而梯形法中递推公式是隐式的。在隐式公式中，必须先用显式公式估计一个值，再用隐式公式迭代，即预估-校正法。隐式与显式相比，有明显的高精度和稳定性。3、误差在数值解法的过程中，每一步计算都会产生误差，误差的来源有两个方面：一是计算误差(即计算机计 算本身的误差)，一个是公式误差(用差分方程代替微分方程)。分别称之为舍人供着.和做理谖卷。舍入误差：计算机的字长是有限的，数字不能表示的完全精确，实际上计算的结果是用有限精度的有 理数(如计算机中使用的浮点数)来近似无限精度的实数，所以在对动态系统模拟

44、的过程中，舍入误差是不 可避免的。通常舍入误差的大小与积分步长成反比，步长 h越小，计算次数多则舍入误差大，但不能随意 加大步长，否则将产生大的截断误差甚至影响系统稳定性。在给定运算序列的条件下，唯一可减少舍入误 差的方法是增加数字表示的精度，如将单精度浮点数(有7位数的精度)表示改为双精度数(有15位数的精度)。实际的动态系统模拟时，多采用双精度数据类型，尽管这样变量占用的存储空间大，处理时间稍长， 但对提高模拟的精度来讲是值得的。截断误差：是用差分方程代替微分方程产生的误差，即用数值解代替微分方程的精确解产生的误差，所以是公式本身的误差，一般用台劳级数来分析积分公式的精度：假设在t n时积

45、分(精确值)已经算出，则用台劳级数可求得t n+1时的精确解：2r(r)r+1y(t n+1)= y(t n + h) = y(t n) + y n)(+ h / 2! - y II (n) + , +h / r! - y ( t n) + O(h )(2.10)如果一个数值解法是用前r+1项来近似的计算y(t n+1),则后面的各项(记为)O(h r+1)是在这一步计算中引进的附加误差，称为这个算法的(局部)截断误差。误差O(h r+1)与h r+1同阶(h0)，即计算的精度保特了 r阶,此时称这个方法是r阶的精度。一个数值方程的阶数可视为衡量这一方法的精确度的重要标志,不同的数 值解法有不

46、同的精度。欧拉法只取(2.10)中的前两项作近似计算，所以是一阶精度的算法；梯形法相当于取(2.10)式中的前三项，故是二阶精度的方法。4、稳定性如果系统是稳定的(系统的状态随时间的推移逐步稳定在某个水平上)，则在模拟的迭代过程中，数值积分的解也应该是稳定的，但由于初始数据的误差及在迭代运算中的舍入误差会对后面的计算结果产生影 响。当步长h选择不合理时，可能使模拟结果不稳定。对于欧拉法，可用下面的检验方程(其中入为方程的 特征根)：j y' = y(2.11)-y(0) = y 0(注意其精确解为y = yoe")来检验步长对数值解稳定性的影响：对(2.11)表示的方程，欧拉

47、公式为：r y n+1 = y n + 入 y n h = (1+ 入 h)y n(2.12)L y(0) = y 0所以，要使数值解稳定，必须使：| 1+入h| < 1 ，解得|入h|<2或h<2T (T = 1/入是系统时间常数)。所以要保证欧拉方法计算的稳定性，步长h必须小于系统时间常数的2倍。在(2.11)中特征根入在一定数量级范围变动，现令 入=1 ,来作一个具体的模拟，此时 h应小于2,否 贝U将不稳定。取 h = 0.01, h = 0.05 , h = 1.0, h = 1.9 , h = 2.0 , h = 2.1六个不同的步长，y0= 1。模拟结果 为：h

48、 = 0.01 :表现出较好精度h = 0.05:虽然近似解接近精确解，但存在误差h = 1.0：解在一步后趋于零，并一直保持，虽然稳定，但明显精度不高h = 1.9:解振荡，幅度值逐衰减，并趋于稳定h = 2.0：解不衰减，等幅振荡h = 2.1 :解的振荡幅度值递增，表明系统不稳定，数值解发散数值解图如图2-3所示。2.2 龙格-库塔(Runge-Kutta)法R-K方法的基本思想是用台劳展开式的前几项来对微分方程求近似解。再以模型(2.1)为例：jy ' = f (t, y)(2.1)1 y (t 0) = y 0假设从t 0开始，以h增长，hi = t 0 + h, yi =

49、y (t 0 + h),在t 0附近展开成台劳级数，保留h2项，则有：yi = y 0 + f (t 0 , y 0) - h + h2 /2 ( Sf / 8 y - dy / dt + Sf / 8 t)(2.13)(此式括号中的导数是在(t 0 , y 0)处的导数值)为了求(2.13 )的解，假设(2.13 )的解可写成如下的形式：yi = y 0 + (b1 k1 + b2 k2 ) - h 其中k1 = f (t 0 , y 0)1(2.14)k 2 = f (t 0 +C2 h , y 0 + a1 k h)-(注意(2.13 )式是关于函数y的导数的算式，而(2.14)式中k1

50、和k2都是y在某点处的导数，相差的只是 常量级的系数不同，问题正是要定出这些系数，从而得到数值解表达式，为此)对k2式右端的函数在(t0 , y0)处展开台劳级数，保留h项，得：k 2 Q f (t 0 , y0) + (C2 - 8 f / 8 t + a1 k1 - S f / 8 y - dy / dt) - h把k、k2代入(2.14)中，得yi = y 0 + bif (t 0 , y 0) h + b2 h f (t 0 , y0)+ (C2 S f / 8 t + ai ki S f / 8 y dy / dt) h (2.15)将(2.15)式与(2.14)式比较，得关于系数的

51、方程：(2.16)b1 = b2,得一个解:, b1 + b2 = 1b2 C2 = 1/2-b2 a1 = 1/2方程组(2.16)中四个未知量，三个方程，故有无穷多解，求出一个解，比如令b1 = b2 = 1/2 a1 = C2 = 1代入(2.14),得如下一个公式：y1 = y 0+1/2 (K1 + K2 ) h其中k = f (t 0 , y 0).I k 2 = f (t 0 + h , y 0 + k1 h) -这是从t 0计算t 1时刻的公式，写成一般的形式，得如下递推公式:yn+1 = y n + h/2 (K1 + K2 )其中k1 = f (t n , y n)(2.1

52、7)I k 2 = f (t n + h , y n + k1 h) -模型(2.1)的数值解公式(2.17)即称为R-K公式，这种数值解法即称为 R-K方法。在(2.13)式中，只保留了 h2项，故公式(2.17)的精度是2阶的，公式(2.17)称为二阶R-K公式。目前在实际模拟问题中，四阶 R-K公式用的最为普遍。在推导四阶 R-K公式时，在台劳公式中保留 到h4项，推导过程与前面类似。一个四阶R-K公式可以是下面的形式：yn+1 = y n + h/6 (K1 + 2K2 +2K3 +K4)其中f k1 = f (t n , y n)Ik 2 =f (tn+h /2 , y n + k1

53、 h/2)>(2.18)k 3 =f (tn+h /2 , y n + k2 h/2)、k 4 =f (tn+h , y n + k3 h)问题：关于R-K方法还有下面一些问题需了解。(1) 多种形式：方程组(2.16)有无穷多解，我们取了一个解，得到了公式(2.17),也可以取其它的解，所以每一阶R-K公式都有多种形式。二阶R-K公式常用的除(2.17)式外，还有yn+1 = y n + K2 hk1 = f (t n , y n)' k 2 = f (t n + h /2 , y n + k1 h/2)(对应(2.16)的解b1 = 0 , b2 = 1 , C2 = a 1

54、 = 1/2 )。四阶R-K公式常用的除(2.18)式外，还有yn+1 = y n + h/8 (K1 + 3K2 +3K3 +K4) k1 = f (t n , y n)Y k 2 = f (t n + h /3 , y n + k1 h/3)k 3 = f (t n + h 2/3 , y n + kh/3+k2 h)'k 4 = f (t n + h , y n + hk1 -hk2 + hk3)(2) 精度：也可推导3阶、5阶或更高阶的R-K公式，但在一般工程中，四阶公式就完全能达到要求 的精度，而且四阶公式是最常用的，所以一般不使用更高阶的公式。另一个特殊情况：一阶 R-K公

55、式，只含有h的1次项，即：yn+1 = y n + h f (t n , y n),这就是欧拉公式, 所以欧拉公式也可看面一阶R-K公式，其精度最低。R-K方法的精度取决于步长h及求解方法。一般来说，为达到同样的精度，四阶方法的步长可以比二阶方法的步长大10倍，而四阶方法每步的计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小。 R-K方法可使用较大的步长也是其一个特点。(3) 单步法：不论几阶的R-K公式都是单步法，在计算 yn+i时只用到y n,可自启动，使用的存储量也小。另外，无论几阶的R-K公式，算式中均有两部分组成：一部分是上一步的结果 y n，第二部分是h=t n+n-1中对f(t,y)的积分，它是步长 h乘以各点斜率的加权平均