火电厂仿真(课件)ppt.ppt

火电厂仿真太原理工大学电气与动力工程学院热能系太原理工大学电气与动力工程学院热能系3.提高机组的安全性、效率和使用寿命 由于培训和研究开发试验不在现场进行，可以避免现场的人为失误、减少停机次数、提高设备的安全性、提高机组使用效率及寿命。 按真实系统的物理性质构造系统的物理模型，按真实系统的物理性质构造系统的物理模型，再现系统的一些特性。再现系统的一些特性。例如例如:新型锅炉生产前，各部件按一定比例缩小，新型锅炉生产前，各部件按一定比例缩小，构成一个几何相似的物理模型。构成一个几何相似的物理模型。 对此分析研究可发现该锅炉在结构和形态上对此分析研究可发现该锅炉在结构和形态上是否合理，安装、操作和检修是否便利。是否合理，安装、操作和检修是否便利。 这种模型仅反映了其结构特性，而未能反映这种模型仅反映了其结构特性，而未能反映锅炉内部传热学、热力学、流体力学的特性。锅炉内部传热学、热力学、流体力学的特性。一、计算机仿真的概念第二节计算机仿真技术简介第二节计算机仿真技术简介 按真实系统的数学关系构造系统的数按真实系统的数学关系构造系统的数学模型，即将实际系统的运动规律用数学学模型，即将实际系统的运动规律用数学形式表达出来，并在数学模型上进行试验形式表达出来，并在数学模型上进行试验，再现系统的某些特性。，再现系统的某些特性。数学模型能精确反映系统内部的各种静数学模型能精确反映系统内部的各种静态和动态特性，如锅炉运行中的燃料化学反态和动态特性，如锅炉运行中的燃料化学反应、传热过程、能量储存与释放、工质循环应、传热过程、能量储存与释放、工质循环流动的特性。流动的特性。数学仿真一般是利用计算机对系统数学数学仿真一般是利用计算机对系统数学模型进行运算和试验的。因此，利用计算机模型进行运算和试验的。因此，利用计算机实现的系统数学仿真也称为：实现的系统数学仿真也称为：计算机仿真计算机仿真。 将真实系统的一部分用数学模型将真实系统的一部分用数学模型描述，并放到计算机上运行，而另一描述，并放到计算机上运行，而另一部分则构造其物理模型或直接采用实部分则构造其物理模型或直接采用实物，然后将它们连接成系统，并在此物，然后将它们连接成系统，并在此系统上进行试验，再现系统的某些特系统上进行试验，再现系统的某些特性。性。这种仿真又称：这种仿真又称：半物理仿真。半物理仿真。4. 计算机仿真计算机仿真定义：定义：将一个能近似描述实际系统的数学模型进将一个能近似描述实际系统的数学模型进行第二次模化，转变为一个仿真模型，并将其放行第二次模化，转变为一个仿真模型，并将其放到计算机上进行运行，以再现系统某些特性的过到计算机上进行运行，以再现系统某些特性的过程。程。第一次模化：研究实际系统与数学模型之间的关系；第一次模化：研究实际系统与数学模型之间的关系；第二次模化：研究数学模型与计算机之间的关系。第二次模化：研究数学模型与计算机之间的关系。计算机仿真是本课程讨论的主题。计算机仿真是本课程讨论的主题。5.计算机仿真系统计算机仿真系统用于系统仿真的一套软、硬设备。用于系统仿真的一套软、硬设备。其主体是计算机。其主体是计算机。按计算机的类型分，有以下仿真系统：按计算机的类型分，有以下仿真系统：1 1、模拟仿真系统、模拟仿真系统2 2、数字仿真系统、数字仿真系统3 3、数模混合仿真系统、数模混合仿真系统用模拟计算机组成的仿真系统用模拟计算机组成的仿真系统。优点：优点：运算是运算是 “ “并行并行”的、的、 “ “连续连续”的的。因此，运算速度快，且更接近连。因此，运算速度快，且更接近连续系统。续系统。缺点：缺点：运算精度低、线路复杂、对采样和逻运算精度低、线路复杂、对采样和逻辑系统的仿真比较困难、仿真的自动辑系统的仿真比较困难、仿真的自动化程度低（依靠排题板接线）。化程度低（依靠排题板接线）。用数字计算机组成的仿真系统用数字计算机组成的仿真系统。优点：优点：被仿真的系统包含在一组程序中，仿真自被仿真的系统包含在一组程序中，仿真自动化程度高，使用方便；运算精度高；易实现逻动化程度高，使用方便；运算精度高；易实现逻辑处理和非线性环节；程序和参数修改容易。辑处理和非线性环节；程序和参数修改容易。缺点：缺点：运算过程是运算过程是“串行串行”的，运算速度相对较的，运算速度相对较低、实时仿真和寻优计算等不如模拟仿真系统快低、实时仿真和寻优计算等不如模拟仿真系统快。 用微型数字计算机阵列组成的仿真系统，用微型数字计算机阵列组成的仿真系统， 称为称为 用模拟计算机和数字计算机组成用模拟计算机和数字计算机组成的仿真系统。的仿真系统。由于模拟计算机和数字计由于模拟计算机和数字计算机的优缺点是互补的。因此，该系统达算机的优缺点是互补的。因此，该系统达到了扬长避短的目的。到了扬长避短的目的。该系统适用于：该系统适用于：（1 1）要求与实物连接，又有许多复杂函数需计算的实时）要求与实物连接，又有许多复杂函数需计算的实时仿真；仿真；（2 2）需要进行反复迭代计算（如统计分析、参数寻优）需要进行反复迭代计算（如统计分析、参数寻优）的仿真；的仿真；（3 3）对计算机控制系统的仿真；）对计算机控制系统的仿真；三、火电厂仿真技术的应用三、火电厂仿真技术的应用 火电厂计算机仿真技术是一门综火电厂计算机仿真技术是一门综合技术。它是以数学、物理、化学合技术。它是以数学、物理、化学、传热学、热力学、流体力学、控、传热学、热力学、流体力学、控制理论、计算机技术、热能动力、制理论、计算机技术、热能动力、电工学、热工仪表及电气仪表等多电工学、热工仪表及电气仪表等多学科专业的理论为基础，以计算机学科专业的理论为基础，以计算机和各种物理效应设备（它再现真实和各种物理效应设备（它再现真实环境）为工具，对真实火电厂发电环境）为工具，对真实火电厂发电机组及其运行进行仿真的技术。机组及其运行进行仿真的技术。 火力发电机组的设备庞大，系统复杂。火力发电机组的设备庞大，系统复杂。因此：因此：对设备的研究、对系统的试验、对参对设备的研究、对系统的试验、对参数的校正、对运行安全和经济的分析、对运数的校正、对运行安全和经济的分析、对运行值班员的培训等，行值班员的培训等，在实际发电机组上直接在实际发电机组上直接进行，或是很困难或根本不可能。利用仿真进行，或是很困难或根本不可能。利用仿真技术的特点、在模型上进行试验研究，在仿技术的特点、在模型上进行试验研究，在仿真机上培训运行值班员，则上述问题可迎刃真机上培训运行值班员，则上述问题可迎刃而解。而解。 等之分。等之分。数学模型：数学模型：是系统的数学描述，是系是系统的数学描述，是系统研究的基础，是计算机仿真的依据。统研究的基础，是计算机仿真的依据。实物模型：实物模型：根据相似性建立的形象模型，具有实体性属物理模拟试验技术范畴。（不介绍）数学模型：数学模型：用符号和数学方程式表示的系统模型，具有抽象性。（本课程内容）系系 统统 模模 型型按描述的状态分按描述的方式分按系统的性质分按求解的方法分按获取的方法分静态模型动态模型连续模型离散模型线性模型非线性模型分析求解模型数字求解模型理论模型黑箱模型反映变量间变量间的相互函数关系反映变量与时间变量与时间的函数关系是时间的连续函数 是时间的离散函数可以用线性微分方程描述不能用线性微分方程描述利于解析法求解 利于计算机求解 理论推导所得的模型试验研究所得的模型 注：（1）一般，仿真所要建立的是系统的动态模型，静态模型包含在动态模型之中。（2）电厂的系统仿真，往往综合采用连续离散线性非线性理论黑箱等模型，以保证模型精度。线性特性汽机功率调节级压力抽气量汽机功率(凝汽式机组) 非线性特性炉膛传热温度流量压差转速油压 可以应用叠加原理不能应用叠加原理必须有效地协调线性与非线性之间的相互关系必要时,非线性特性线性化 工质状态工质状态 （温度） 工况工况环境环境条件条件 （过热器管长） 时间时间建模时多采建模时多采用用分区集总分区集总方法方法，即将，即将三维空间三维空间的分布参数简化的分布参数简化为一维空间。否为一维空间。否则无法求解。则无法求解。汽机甩负荷转速烟温主汽温时间常数小响应快燃料汽压减温水量主汽温时间常数大响应慢建模时应处理好模型的简洁与精确之间的关系对响应速度差异很大的对象，一般分别建对响应速度差异很大的对象，一般分别建立动态数学模型，可采用不同的时间步长立动态数学模型，可采用不同的时间步长进行仿真计算。进行仿真计算。分解集合把实际系统分成若干子系统，分别建立其模型，然后综合。把若干子系统视为一个整体，统一建立其模型。模型精细复杂模型简洁粗糙因此，在建立复杂系统的动态模型时，应善于对系统进行合理的分解，即按不同的化学、物理过程（如燃烧、传热、流动、做功等），把系统分成许多子系统，这些子系统并不都是简单地按设备的类型来划分的，而根据系统的属性来划分更为合理。在分解系统时，正确处理各子系统间的关系是十分重要的，因此要正确地选择系统中的状态变量，充分保证信息的传递如实际过程那样进行。为满足实时运行要求，仿真机应采用速度快、为满足实时运行要求，仿真机应采用速度快、容量大的计算机。容量大的计算机。 （7 7）模型计算涉及的数据量大）模型计算涉及的数据量大数学模型运算涉及的数学模型运算涉及的I/OI/O变量、中间变量、常变量、中间变量、常数等数量巨大。例如：数等数量巨大。例如：一台采用常规仪表盘一台采用常规仪表盘/ /台的台的300MW300MW燃煤发电机组燃煤发电机组的全范围仿真机有：的全范围仿真机有：50005000多个多个I/OI/O变量；变量； 10000 10000多多个中间变量；个中间变量； 建模是综合知识的体现，它涉及的知识广泛。建模是综合知识的体现，它涉及的知识广泛。如：如：理论知识理论知识物理学化学数学热力学传热学物理学化学数学热力学传热学专业知识专业知识锅炉原理汽轮机原理电机原理机械原理、控锅炉原理汽轮机原理电机原理机械原理、控制理论制理论实践知识实践知识结构知识运行知识试验技术结构知识运行知识试验技术系统数学模型建立在每一个方框基础上系统数学模型建立在每一个方框基础上分解分解模型。模型。准确的模型是取得正确仿真结果的基础 系统简化 与 的正确与否密切相关 假设条件 物理规律依赖于对系统内部 化学规律 的熟悉程度和概括能力 运作机理 是指与热平衡相关的系统。例如主燃料系统、风烟系统、炉膛、主蒸汽系统、汽轮机和抽汽系统、发电机、凝汽器、凝结水系统、给水系统、循环水系统、冷却塔和和加热器疏水系统等。ucdtudcdtudcyadtdyadtydadtydnnnnnnnnnnn12211101111自由项强迫项(2-1)式中：y=y（t）系统的输出量 u=u（t）系统的输入量 ai， cI系数dtdp pny+a1pn-1y+an-1py+any=c0pn-1u+c1pn-2u+cn-1unjnjjjnjjnupcypa0101或：或：（其中a0=1）njjjnnjjjnpapcuy0101(2-2)()()()()()()()()()()()(2211222221212112121111tuytaytaytadtdytuytaytaytadtdytuytaytaytadtdynnnnnnnnnnnaij随时间而变时变系统aij不随时间变定常系统Ui(t)全部恒等于零齐次线性微分方程组Ui(t)有一个不为零非齐次线性微分方程组式中：yi 局部输出量； ui 局部输入量； aij, 系数（2-3）)(,1122tudtyddtyddtdyytfdtydnnnn令 y1 y2 y3 yn （2-4）一种表达式为则一阶微分方程组的另)(,211132221tuyyytfdtdyydtdyydtdyydtdynnnnnn式中：y1，y2，yn为自变量t的n个未知函数。求解式（2-4）相当于求解式（2-5）(2-5)()()()()()()(12110111sUcsUscsUscsYassYasYsasYsnnnnnnn（2-6）即： sn 替换 dn/dtn s为拉普拉斯算子 Y(s)替换y(t) Y(s)为输出量y(t)的拉普拉斯变换 U(s)替换u(t) U(s)为输入量u(t)的拉普拉斯变换 此式与式（2-2）比较可知，在初值为零的情况下，用算子P表示的式子与用传递函数G(s)表示的式子在形式完全相同。nnnnnnnasasascscscsUsYsG11112110)()()(（2-7）将上式整理可得系统传递函数G（s）的规范式为：nnnnnnnnajajajcjcjcjcjUjYjG)()()()()()()()()(111122110（2-8）上述模型只描述了系统输入与输出之间的关系未描述系统内部的情况称为：外部模型仿 真在计算机上对系统的数学模型进行试验必须在计算机上复现（实现）系统即复现系统输入量复现系统输出量复现系统的内部变量又称：状态变量仿真要求将系统辨识或其它方法建立的外部模型转化为：内部模型内部模型状态空间模型状态空间模型状态方程状态方程控制理论中称之为：实现问题实现问题传递函数G(s)U1（s）U2（s）UL（s）Y1（s）Y2（s）Ym（s）当引入状态变量后，框图变成了：状态变量xi UYlnmU、Y、X均为时间的函数。luuuU21（2-9）myyyY21（2-10） m个输出变量yI，输出变量的集合可用输出向量Y来表示；nxxxX21（2-11） n个状态变量xI，状态变量的集合可用状态向量X来表示：),(),()(00ttUtXftX),(),()(00ttUtXgtY（2-12）（2-13）系统的状态方程系 统 的 输 出 方程式中：X（t0）初始状态向量 U（t0,t）初始时刻t0到t的系统输入向量 f 表示一个单值函数 g 也表示一个单值函数若用微分方程表示，则其表示式为：)(),()(tUtXFtX)(),()(tUtXGtY（2-14）（2-15）若用线性微分方程表示，则其表示式如下：)()()()()(tUtBtXtAtX)()()()()(tUtDtXtCtY（2-16）（2-17）系统数学模型的状态方程规范式在状态空间中： A（t）(nn阶)系统矩阵（系数矩阵、状态矩阵）B（t）(nl阶)控制矩阵（驱动矩阵）C（t）(mn阶)输出矩阵D（t）(ml阶)传递矩阵（传输矩阵）当A(t)、B(t)、C(t)、D(t)随时间而变化时，系统称为变系数系统（时变系统）；当A(t)、B(t)、C(t)、D(t)不随时间而变时，可用A、B、C、D来表示，系统则称为常系数系统（定常系统）。D(t)U(t)项表示系统的输入和输出通过D(t)直接联系的部分。在普通的控制系统中，通常D(t)=0，则式（2-17）变成：（2-18）)()()(tXtCtYa、微分方程的强迫项无无导数项时的状态方程（规范式）)(11111tucyadtdyadtdadtydnnnnnnn在该情况下，式（2-1）中的微分方程为：(2-19）其中：nnnnnnndtydxdtydxxdtydxxdtdyxxyx/1112223121今引进n个状态变量一阶微分方程组nnndtydx )()(112112111222111tucxaxaxaxatucyadtdyadtydadtydannnnnnnnnnnn把上述一阶微分方程组写成矩阵形式，可得：ucxxxaaaaxxxXnnnnnn1211212100101000010（2-20）XxxxYn0, 0, 0, 1 0, 0, 0, 1 21(2-21)121101000010aaaaAnnn100ncB0 , 0 , 0 , 1CCXYBuAXX（222）（223）显然，状态变量的初值为：)0()0()0()0()0()0(121nnyyyxxx系统. 例例 ：已知描述某系统的微分方程为：已知描述某系统的微分方程为：uydtdydtyddtyd54322233且：且：3)0(,2)0(, 1)0(yyy 试求该系统的试求该系统的状状态方程，并将给出的初态方程，并将给出的初始条件用状态变量的初值表示。始条件用状态变量的初值表示。解：令22321dtydyxdtdyyxyx 则有uxxxxxxxx523432133221将上式写成矩阵形式，uydtdydtyddtyd54322233 321001xxxy状态变量的初值为：3)0()0(2)0()0(1)0()0(321yxyxyx 即可得状态方程：uxxxxxx500234100010321321规范式（规范式（）b、微分方程强迫项有有导数项时的状态方程（规范式、）ucdtudcdtudcyadtdyadtydadtydnnnnnnnnnnn12211101111dtdp pny+a1pn-1y+an-1py+any=c0pn-1u+c1pn-2u+cn-1unjnjjjnjjnupcypa0101或：（其中a0=1）njjjnnjjjnpapcuy0101njjjnxpau0引进n个状态变量：xpdtxdxxxpdtxdxxpxdtdxxxxxnnnnn111122223121即： pjx=xj+1 (j=0, 1, 2, , n-1)（2-24）则式（2-24）可写为：xpaxpapxaxaxpaunnnnnjjjn022101001njnjjnxpaxau101njjjnnnuxaxxp（*）将式（2-24）代入式（2-2）可得：njnjnjjjnjjnjjnxpapcypa01001（*）由式（*）和（*）可得：（2-25）BuAXX1010111njnjjjnjjnCXxcxpcy规范式（）121101000010aaaaAnnn其中：100B,0, 2, 1cccCnn)0()0()0()0()0()0(121nnxxxxxx规范式（）的应用：是系统的数学模型可直接写成状态方程形式，且初始条件也是按各状态变量的初值给出时，计算将十分方便。假定给出的微分方程组为：njnjjjnjjnupcypa00（2-30）ucupcupcyapyaypaypannnnnnn1101110即u和y的导数阶次相等若令：nnnjjjnnnpxucyapucyapucyapucyap)()()()(1111100则有：ucyapxnnn若将式（2-30）两边同时取不定积分一次并令：122112001)()()(nnnnnpxucyapucyapucyap则有：ucyaxpxnnnn111并令：233113002)()()(nnnnnpxucyapucyapucyap则有：ucyaxpxnnnn2212同理有ucyaxpxjjjj1ucyapxxjjjj1或：则有：ucyapxxjjjj111)()()()(1122112001ucyaucyapucyapucyapxjjjjjjj即：（j=0，1，2，n）（2-31）（#）又因为在j=0时，由式（#）有：00010ucyaxx 则在a0=1时，有：求初值用求初值用ucxy01uaccxaucucxaxucyaxxuaccxxaucucxaxucyaxxuaccxxaucucxaxucyaxxuaccxxaucucxaxucyaxxnnnnnnnnnnnnnnnnnnnn)()()()()()()()(0120111101111011111202312201232232101211101121121此时由式（#）可得： 由上述讨论，可得式（2-30）在a0=1时的状态方程和输出方程如下：BuAXuaccaccaccaccxxxxaaaaxxxxXnnnnnnnnnn0101202101121121121000100010001DuCXucxxxxynn01210, 0, 0, 1 （2-32）（2-33）规范式（规范式（）uBAXuccccxxxxaaaaxxxxXnnnnnnnn121121121121000100010001CXxxxxynn1210,0,0, 1 （2-32）（2-33）规范式（规范式（）jiijiijijucyax1)(1)(1)0()0()0(其中： ( j=1，2，n )udtduydtdydtyddtyd654322233式中：u为单位阶跃函数，且初始条件为：0)0(u3)0(,2)0(, 1)0(yyy 试求该系统的状态方程和输出方程，并将给出的初始条件用状态变量的初值表示。解：根据题意，有：n=3， a0=1， a1=2， a2=3， a3=4，c0=0， c1=0， c2=5， c3=6可由规范式（），得系统的状态方程和输出方程：uxxxxxxX6500041030123213213210, 0, 0, 1 xxxy其中，各状态变量由式（2-31）可求得：uydtdydtyducyaucyapucyapxydtdyucyaucyapxyucyax532)()()(2)()(222211002311002001jiijiijijucyax1)(1)(1)0()0()0(可求得：100343)0()0()0()0()0()0()0()0()0()0()0(422)0()0()0()0()0()0()0(1)0()0()0()0(2212211003111002001ucyayayucyaucyaucyaxyayucyaucyaxyucyax 应该指出：一个系统可以引入不同组合的状态变量。那么，所设的状态变量不同，获得的状态方程也会不同。即一个外部模型，可能有很多不同的内部模型（状态方程）一个系统的状态方程形式并非是唯一的。)()()()(1122112001ucyaucyapucyapucyapxjjjjjjj（j=0，1，2，n）（2-31）由式（由式（2-31）可知，所设的各状态变量与可知，所设的各状态变量与u的导数项有关，因此而得的状态方程的导数项有关，因此而得的状态方程:BuAXuaccaccaccaccxxxxaaaaxxxxXnnnnnnnnnn0101202101121121121000100010001DuCXucxxxxynn01210,0,0, 1（2-32）（2-33）规范式（规范式（）其解其解可能可能不是唯一的。不是唯一的。ucupcupcyapyaypaypnnnnnnn110111（2-34）如取以下如取以下n个变量作为一组状态变量：个变量作为一组状态变量：uxuuuuyxuxuuuyxuxuuyxuyxnnnnnnnn1112211012221031110201 （3-35）且式中，且式中， 0， 1， 2， n分别为：分别为：a0=1011221103122133021122011100nnnnnnaaaacaaacaacacc（2-36）则可保证状态方程解的唯一性。根据上面一组状态变量，对于式（2-34）所描述的系统，可得出如下规范式：DuCXYBuAXX（2-37）（2-38）式中：00121121121121,0, 0, 0, 1 ,100001000010,cDCBxxxxXaaaaAxxxxXnnnnnnnnn状态方程规规范范式式（）1对于以上各规范式：若A A，B B，C C，D D随时间变化称为：变系数系统变系数系统或时变系统时变系统若A A，B B，C C，D D不随时间变称为：常系数系统常系数系统或定常系统定常系统2一个系统的数学模型并非唯一的。表示一个系统的模型，可用不同的数学方法。最常用的方法为传递函数，其次是微分方程，状态方程。3采用微分方程微分方程进行连续系统的仿真时：微分方程一阶微分方程组代数方程采用变量分离法一般把转换成转换成 （a）应根据初始值是否为零，选择相应的规范式。初始值为零时可选用规范式（）或规范式（）。 （b）若高阶微分方程的初始条件是输入u，输出y及其导数的初值，应变换成状态方程的初值。 （c）若微分方程的强迫项含有导数项，可采用规范式（）或规范式（），但此时可能没有唯一的解，若要获得唯一的解，则应采用规范式（）。5、对已确定的系统，仿真所采用的数学方法可以不同，但研究结果应与系统的实际情况一致。数值积分是数值分析的一个基本问题。数值积分是数值分析的一个基本问题。也是复杂计算问题中的一个基本组成部分。也是复杂计算问题中的一个基本组成部分。数值积分往往用极简单的方法就能较好地得数值积分往往用极简单的方法就能较好地得出对所求解的具体数值问题的解答。出对所求解的具体数值问题的解答。但数值积分的难点在于计算时间有时会过长但数值积分的难点在于计算时间有时会过长，有时会出现数值不稳定现象。，有时会出现数值不稳定现象。另外，数值积分的理论性较强。其理论和方另外，数值积分的理论性较强。其理论和方法都已经比较成熟，计算精度也比较高。法都已经比较成熟，计算精度也比较高。3.1 仿真中研究数值积分法的意义仿真中研究数值积分法的意义数值解的一种近似方法。对于连续系统的高阶微分方程，可化为若干个一阶微分方程组成的方程组。数值积分法是求解微分方程：00)(),(ytyytfdtdyy 例如：下式所示的状态方程BuAxx可以化为一个一阶微分方程组ubxaxaxaxubxaxaxaxubxaxaxaxnnnnnnnnnnn2211222221212112121111所以，连续系统的仿真就是从给定的初始条件出发，对描述系统动态特性的常微分方程或常微分方程组进行求解，从而得到系统在一定输入作用下的变化过程。在求解这些微分方程时，最常用、也是最有效的一种方法就是数值积分法。3.2 数值积分法仿真的基本原理数值积分法仿真的基本原理对微分方程（3-1）两端同时取积分，可得dyftytytt0),()()(0dyftytytt0),()()(0当 时，上式变为 ：1nttdtytftytynttn10),()()(01将积分项拆成两项 dtytfdtytftytynnnttttn10),(),()()(01)(nty故上式可写为： dtytftytynnttnn1),()()(1此式是方程（3-1）在tn+1时刻的精确解。 在数值解法中，希望用近似解 ：nnnQyy1代替准确解，其中 ： 1ny为)(1ntyny为)(ntynQ为1),(nnttdtytf的近似值 httnn1令：称为计算步长或步距 是从常微分方程（是从常微分方程（3-1）出发建）出发建立的离散数学模型立的离散数学模型差分方程。差分方程。 dtytftytynnttnn1),()()(1nnnQyy11ny为)(1ntyny为)(ntynQ为1),(nnttdtytfnnnQyy11ny为)(1ntyny为)(ntynQ为1),(nnttdtytf由此可见，数值积分法就是在已知微分方程初值的情况下，求解该方程在一系列离散点 处的近似值，其特点是步进式根据初始值逐步递推地计算出以后各时刻的值。从式（3-8）可知，数值积分法的主要问题归结为如何对f（t，y）进行数值积分求出f（t，y）在区间tn， tn+1上定积分的近似值Qn。采用不同的方法求Qn，就出现了各种各样的数值积分方法。不同的数值积分将对求解的精度、速度和数值稳定性会产生不同的影响，这将在下述内容中具体介绍。 数值积分法种类繁多，在此从实用角度介绍几种基本的方法 3.3 欧拉欧拉(Euler)法法3.3.1 简单欧拉法简单欧拉法 欧拉法是一种最简单的数值积分法，对于方程：00)(),(ytyytfdtdyy 在区间tn， tn+1上求积分，得到： dtytftytynnttnn1),()()(1若区间tn， tn+1足够小，则tn， tn+1上的f（t，y）可近似地看成常数f（tn，yn） 。故可用矩形面积近似代替),(1ytfnntt即：tntn+1f（tn，yn）于是有：11),()(nnnnnyhytfyty将此式写成差分方程为： , 3 , 2 , 1 , 0),(1nhytfyynnnn著名的欧拉公式著名的欧拉公式 3.3.2 改进的欧拉法改进的欧拉法 如果用梯形面积而不是矩形面积来代替每一个小区间上的曲线积分，就可以提高计算精度，梯形法的计算公式为： ),(),(2111nnnnnnytfytfhyy（3-11）式中的右端含有待求量yn+1，因而它是隐函数形式。这种方法不能自行启动运算，需要依赖其它算法的帮助。 每次计算都用欧拉法算出y（t n+1 ）的近似值 ，以此计算近似值 ，然后利用梯形公式（3-11）求出修正后的 。即有： 帮助方法：帮助方法：Pny1),(111PnnPnytfyCny1),(1nnnpnythfyy),(),(2111PnnPnnnCnytfytfhyy预估式 校正式 （3-12）改进的欧拉公式改进的欧拉公式 3.3.3 几个基本概念几个基本概念 简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn+1，这种方法称为单步法，它是一种自行启动的算法。如果求yn+1时需要tn , tn-1 , tn-2 时刻yn , yn-1 , yn-2 的值，这种方法为多步法（改进的欧拉法为两步法），它是一种不能自行启动的算法。1、单步法与多步法、单步法与多步法 简单的欧拉法表达式的右端，计算 所用的数据均已求出，这种公式称为显式公式。 改进的欧拉法表达式的右端，有待求量 ，这种公式称为隐式公式。 隐式公式不能自行启动，需要用预估-校正法。 单步法和显式在计算上比多步法和隐式方便，但有时为了满足精度、稳定性等方面的要求，需要采用隐式算法。 2、显式与隐式、显式与隐式3、截断误差、截断误差 这里用泰勒级数为工具来分析数值积分公式的精度假定yn是精确的，用泰勒级数表示 处的精确解，即： 1nt)(0)(!)(! 2)()()(1)()2(2)1(1rnrrnnnnhtyrhtyhthytyty显然，简单的欧拉法是从以上精确解中取前两项之和来近似计算，每一步由这种方法引入的误差称为局部截断误差，简称截断误差。简单的欧拉法的截断误差为： )(0)(211hytynn不同的数值方法有不同的截断误差。一般若截断误差为 ，则方法为r阶的。所以方法的阶数可以作为衡量方法精确度的一个重要标志。)(01rh4、舍入误差、舍入误差 由于计算机的字长有限，数字不能表示得完全精确，在计算过程中不可避免地会产生舍入误差。舍入误差与计算步长成反比。如果计算步长小，计算次数多，则舍入误差就大。舍入误差除了与计算机字长有关以外，还与计算机所使用的数字系统、数的运算次序以及计算所用的子程序的精度等因素有关。 采用数值积分法求解稳定的常微分方程，应该保持原系统的稳定特征。但是：（1）在计算机逐步计算时，初始数据的误差及计算过程的舍入误差对后面的计算结果将产生影响。（2）如果计算步长取的不合格，有可能使仿真出现不稳定的结果。下面我们简单讨论一下这个问题。 差分方程的解与微分方程的解类似，可分为特解和通解两部分。与稳定性有关的是方程的通解，它取决于差分方程的特征根是否满足稳定性条件。例如，为了考查欧拉法的稳定性，我们研究检验方程（Test Equation）：yy其中， 为方程的特征根 nnnnyhhyyy)1 (111 h即：2h表明：为使数字仿真稳定，对计算步长应有所限制表明：为使数字仿真稳定，对计算步长应有所限制。 另外，稳定性还与系统的特性以及数值积分方法有关。上述分析欧拉法稳定性的思想，同样适用于其它数值积分方法。 （3-14）（3-15）3.4 龙格龙格-库塔库塔（Runge-Kutta）法法由前面的分析可知，将泰勒展开式多取几项以后截断，就能提高截断误差的阶数和计算精度。然而，直接采用泰勒展开方法要计算函数的高阶导数，运用起来不便。龙格-库塔方法的基本思想是：用几个点上的函数值的线性组合代替函数的各阶导数，然后按泰勒级数展开确定其中的系数，这样既可避免计算高阶导数，又可提高积分的精度。龙格-库塔法有多种形式，以下从实用角度直接给出公式的形式和相应的精度。3.4 .1 龙格龙格-库塔方法的基本思想库塔方法的基本思想3.4.2 二阶龙格二阶龙格-库塔方法库塔方法2阶龙格-库塔方法的公式为：),(),()(2/(121211hkyhtfkytfkkkhyynnnnnn（3-16）上式表示的数值解是用的泰勒级数在2阶导数以后截断所求得的，因此称为2阶方法。故2阶龙格-库塔法与式改进的欧拉法相比，实质完全相同。所以改进的欧拉法实质上是2阶龙格-库塔法。)(03h截断误差为：实时仿真实时仿真的的2阶龙格阶龙格-库塔方法库塔方法)2,2(),(12121khyhtfkytfkhkyynnnnnn（3-17）)(03h其截断误差为：3.4.3 四阶龙格四阶龙格-库塔方法库塔方法4阶龙格-库塔法是一种最常用的方法。其经典表达式为：),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn（3-18）其截断误差为：)(05h显然， 4阶龙格-库塔法的计算量较大，但计算精度较高，在比较不同算法的计算精度时，常以它的计算结果作为标准。实时仿真实时仿真的的4阶龙格阶龙格-库塔方法库塔方法)2103,54()52,53()52,52()5,5(),()105515(2441521413121543211khkhyhtfkhkkhyhtfkkhyhtfkkhyhtfkytfkkkkkkhyynnnnnnnnnnnn（3-19）以上公式都是标量形式，如果要换成向量形式，只要把式中的标量y，f，k换成向量Y，F，K即可。从理论上讲，可以构造任意阶数的龙格-库塔方法，但是，精度的阶数与计算函数值f 的次数之间并非等量增加的关系，见下表所列：表 3.1f 的计算次数与精度阶数的关系每一步计算 f的次数234567N8精度阶数234456n-2 由此可见，4阶经典龙格-库塔方法有其优越性，而4阶以上的龙格-库塔方法计算f所需的次数比阶数多，增加了计算量，从而限制了更高龙格-库塔方法的应用。对于大量的实际问题，4阶方法已可满足精度要求，所以得到了广泛的应用。我们仍采用检验方程 进行讨论，对它利用泰勒级数展开得：3.4.4 龙格龙格-库塔公式的稳定区域库塔公式的稳定区域yy)(0!1)(11rinriinnhyrhyy对于 ，有 ，将它代入上式得：yyyyii)()(0)(!1)(! 211 121rnrnhyhrhhy（3-20）（3-21）令： 将其代入上式得到该式的稳定条件稳定条件为：hh1!1!21121rhrhh（ 3-22 ）由此稳定条件，下表给出了各阶龙格-库塔公式的稳定区域。表表3.2 龙格龙格-库塔公式的稳定区域库塔公式的稳定区域r1h所在的实稳定区域1h1（2，0）22/12hh （2，0）36/2/132hhh（2.51，0）424/6/2/1432hhhh（2.78，0）在使用龙格-库塔公式时，选取的步长应使落在稳定区域内。否则，在计算时会产生很大的误差，从而得不到稳定的数值解。例如：用4阶龙格-库塔方法解：1)0(,20yyy （1）用解析法求解：本例是稳定的；（2）用数值法求解： 当h=0.1 时，数值解是稳定的; 当h=0.2时，数值解就不稳定了，这时因为：4)20(2 . 0hh此数值在稳定区间以外，所以数值解不能收敛。这种对步长有限制的数值积分法称为条件 稳定积分法稳定积分法。从4阶龙格-库塔方法的稳定条件：中可以看出，系统的特征根越大，需要的积分步长就越小，这一点可作为选择步长的依据。步长的大小，除了与数值积分方法的阶数有关外，还与方程本身的性质有关。078. 2h除以上介绍的欧拉法、龙格-库塔法外，数值积分方法还有许多种，如亚当姆斯方法、吉尔方法等等。此处不一一介绍。3.5 计算方法的选择计算方法的选择数值积分方法很多，在实际使用时存在一个选择问题。对于一个具体问题，如何选择具有一定的难度，至今尚无一种确定的方法。一般来说，数值积分方法的选择应考虑的因素有：1、精度要求、精度要求数值积分法的精度受截断误差舍入误差积累误差的影响，而这些误差与积分方法、步长、计算时间、计算机精度等有关。一般：（1）积分方法的阶数越高，截断误差越小，精度越高；（2）步长越小，精度越高；（3