第5章 数值积分法仿真

1、数值积分是数值分析的一个基本问题。数值积分是数值分析的一个基本问题。也是复杂计算问题中的一个基本组成部分。也是复杂计算问题中的一个基本组成部分。数值积分往往用极简单的方法就能较好地得数值积分往往用极简单的方法就能较好地得出对所求解的具体数值问题的解答。出对所求解的具体数值问题的解答。但数值积分的难点在于计算时间有时会过长，但数值积分的难点在于计算时间有时会过长，有时会出现数值不稳定现象。有时会出现数值不稳定现象。另外，数值积分的理论性较强。其理论和方另外，数值积分的理论性较强。其理论和方法都已经比较成熟，计算精度也比较高。法都已经比较成熟，计算精度也比较高。5.1 仿真中研究数值积分法的意义仿

2、真中研究数值积分法的意义数值解的一种近似方法。对于连续系统的高阶微分方程，可化为若干个一阶微分方程组成的方程组。数值积分法是求解微分方程：00)(),(ytyytfdtdyy 例如：下式所示的状态方程BuAxx可以化为一个一阶微分方程组ubxaxaxaxubxaxaxaxubxaxaxaxnnnnnnnnnnn2211222221212112121111所以，连续系统的仿真就是从给定的初始条件出发，对描述系统动态特性的常微分方程或常微分方程组进行求解，从而得到系统在一定输入作用下的变化过程。在求解这些微分方程时，最常用、也是最有效的一种方法就是数值积分法。5.2 数值积分法仿真的基本原理数值积

3、分法仿真的基本原理对微分方程（5-1）两端同时取积分，可得dyftytytt0),()()(0dyftytytt0),()()(0当 时，上式变为 ：1nttdtytftytynttn10),()()(01将积分项拆成两项 dtytfdtytftytynnnttttn10),(),()()(01)(nty故上式可写为： dtytftytynnttnn1),()()(1此式是方程（5-1）在tn+1时刻的精确解。 在数值解法中，希望用近似解 ：nnnQyy1代替准确解，其中 ： 1ny为)(1ntyny为)(ntynQ为1),(nnttdtytf的近似值 httnn1令：称为计算步长或步距 是从

4、常微分方程（是从常微分方程（5-1）出发建立的离散数学）出发建立的离散数学模型模型差分方程。差分方程。 由此可见，数值积分法就是在已知微分方程初值的情况下，求解该方程在一系列离散点 处的近似值，其特点是步进式根据初始值逐步递推地计算出以后各时刻的值。从式（5-8）可知，数值积分法的主要问题归结为如何对f（t，y）进行数值积分求出f（t，y）在区间tn， tn+1上定积分的近似值Qn。采用不同的方法求Qn，就出现了各种各样的数值积分方法。不同的数值积分将对求解的精度、速度和数值稳定性会产生不同的影响，这将在下述内容中具体介绍。 121,nntttt数值积分法种类繁多，在此从实用角度介绍几种基本的

5、方法 5.3 欧拉欧拉(Euler)法法5.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著名的欧拉公式著名的欧拉公式 5.3.2 改进的欧

6、拉法改进的欧拉法 如果用梯形面积而不是矩形面积来代替每一个小区间上的曲线积分，就可以提高计算精度，梯形法的计算公式为： ),(),(2111nnnnnnytfytfhyy（5-11）式中的右端含有待求量yn+1，因而它是隐函数形式。这种方法不能自行启动运算，需要依赖其它算法的帮助。 每次计算都用欧拉法算出y（t n+1 ）的近似值 ，以此计算近似值 ，然后利用梯形公式（5-11）求出修正后的 。即有： 帮助方法：帮助方法：Pny1),(111PnnPnytfyCny1),(1nnnpnythfyy),(),(2111PnnPnnnCnytfytfhyy预估式 校正式 （5-12）改进的欧拉公式

7、改进的欧拉公式 5.3.3 几个基本概念几个基本概念 简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn+1，这种方法称为单步法，它是一种自行启动的算法。如果求yn+1时需要tn , tn-1 , tn-2 时刻yn , yn-1 , yn-2 的值，这种方法为多步法（改进的欧拉法为两步法），它是一种不能自行启动的算法。1、单步法与多步法、单步法与多步法 简单的欧拉法表达式的右端，计算 所用的数据均已求出，这种公式称为显式公式。 改进的欧拉法表达式的右端，有待求量 ，这种公式称为隐式公式。 隐式公式不能自行启动，需要用预估-校正法。 单步法和显式在计算上比多步法和隐式方便，但有时为了满足精度

8、、稳定性等方面的要求，需要采用隐式算法。 2、显式与隐式、显式与隐式1ny1ny3、截断误差、截断误差 这里用泰勒级数为工具来分析数值积分公式的精度。假定yn是精确的，用泰勒级数表示 处的精确解，即： 1nt)(0)(!)(! 2)()()(1)()2(2)1(1rnrrnnnnhtyrhtyhthytyty显然，简单的欧拉法是从以上精确解中取前两项之和来近似计算，每一步由这种方法引入的误差称为局部截断误差，简称截断误差。简单的欧拉法的截断误差为： )(0)(211hytynn不同的数值方法有不同的截断误差。一般若截断误差为 ，则方法为r阶的。所以方法的阶数可以作为衡量方法精确度的一个重要标志

9、。)(01rh（5-13）由此可见，简单的欧拉法是1阶精度；改进的欧拉法由于采用了平均斜率，相当于取了泰勒级数的前3项，因此为2阶精度。 分析欧拉法截断误差的思想，同样也适用于其它数值积分方法。4、舍入误差、舍入误差 由于计算机的字长有限，数字不能表示得完全精确，在计算过程中，不可避免地会产生舍入误差。舍入误差与计算步长成反比。如果计算步长小，计算次数多，则舍入误差就大。舍入误差除了与计算机字长有关以外，还与计算机所使用的数字系统、数的运算次序以及计算所用的子程序的精度等因素有关。 5、数值解的稳定性问题、数值解的稳定性问题 采用数值积分法求解稳定的常微分方程，应该保持原系统的稳定特征。但是：

10、 （1）在计算机逐步计算时，初始数据的误差及计算过程的舍入误差对后面的计算结果将产生影响。 （2）如果计算步长取的不合格，有可能使仿真出现不稳定的结果。下面我们简单讨论一下这个问题。差分方程的解与微分方程的解类似，可分为特解和通解两部分。与稳定性有关的是方程的通解，它取决于差分方程的特征根是否满足稳定性条件。例如，为了考查欧拉法的稳定性，我们研究检验方程（Test Equation）：yy其中， 为方程的特征根，对此有： 显然，要使该差分方程是稳定的，必须使下式成立。 nnnnyhhyyy)1 (111 h即：2h表明：为使数字仿真稳定，对计算步长应有所限制。表明：为使数字仿真稳定，对计算步长

11、应有所限制。 另外，稳定性还与系统的特性以及数值积分方法有关。上述分析欧拉法稳定性的思想，同样适用于其它数值积分方法。 （5-14）（5-15）5.4 龙格龙格-库塔库塔（Runge-Kutta）法法由前面的分析可知，将泰勒展开式多取几项以后截断，就能提高截断误差的阶数和计算精度。然而，直接采用泰勒展开方法要计算函数的高阶导数，运用起来不便。龙格-库塔方法的基本思想是：用几个点上的函数值的线性组合代替函数的各阶导数，然后按泰勒级数展开确定其中的系数，这样既可避免计算高阶导数，又可提高积分的精度。龙格-库塔法有多种形式，以下从实用角度直接给出公式的形式和相应的精度。5.4 .1 龙格龙格-库塔方

12、法的基本思想库塔方法的基本思想5.4.2 二阶龙格二阶龙格-库塔方法库塔方法2阶龙格-库塔方法的公式为：),(),()(2/(121211hkyhtfkytfkkkhyynnnnnn（5-16）上式表示的数值解是用的泰勒级数在2阶导数以后截断所求得的，因此称为2阶方法。故2阶龙格-库塔法与式改进的欧拉法相比，实质完全相同。所以改进的欧拉法实质上是2阶龙格-库塔法。)(03h截断误差为：实时仿真实时仿真的的2阶龙格阶龙格-库塔方法库塔方法)2,2(),(12121khyhtfkytfkhkyynnnnnn（5-17）)(03h其截断误差为：5.4.3 四阶龙格四阶龙格-库塔方法库塔方法4阶龙格-

13、库塔法是一种最常用的方法。其经典表达式为：),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn（5-18）其截断误差为：)(05h显然， 4阶龙格-库塔法的计算量较大，但计算精度较高，在比较不同算法的计算精度时，常以它的计算结果作为标准。实时仿真实时仿真的的4阶龙格阶龙格-库塔方法库塔方法)2103,54()52,53()52,52()5,5(),()105515(2441521413121543211khkhyhtfkhkkhyhtfkkhyhtfkkhyhtfkytfkkkkkkhyynnnn

14、nnnnnnnn（5-19）以上公式都是标量形式，如果要换成向量形式，只要把式中的标量y，f，k换成向量Y，F，K即可。从理论上讲，可以构造任意阶数的龙格-库塔方法，但是，精度的阶数与计算函数值f 的次数之间并非等量增加的关系，见下表所列：表5.1 f 的计算次数与精度阶数的关系每一步计 算f的次数234567N8精度阶数234456n-2由此可见，4阶经典龙格-库塔方法有其优越性，而4阶以上的龙格-库塔方法计算f所需的次数比阶数多，增加了计算量，从而限制了更高龙格-库塔方法的应用。对于大量的实际问题，4阶方法已可满足精度要求，所以得到了广泛的应用。我们仍采用检验方程 进行讨论，对它利用泰勒级

15、数展开得：5.4.4 龙格龙格-库塔公式的稳定区域库塔公式的稳定区域yy)(0!1)(11rinriinnhyrhyy对于 ，有 ，将它代入上式得：yyyyii)()(0)(!1)(! 211 121rnrnhyhrhhy（5-20）（5-21）令： 将其代入上式得到该式的稳定条件稳定条件为：hh1!1!21121rhrhh（ 5-22 ）由此稳定条件，下表给出了各阶龙格-库塔公式的稳定区域。表表5.2 龙格龙格-库塔公式的稳定区域库塔公式的稳定区域r1h所在的实稳定区域1h1（2，0）22/12hh （2，0）36/2/132hhh（2.51，0）424/6/2/1432hhhh（2.78，

16、0）在使用龙格-库塔公式时，选取的步长应使落在稳定区域内。否则，在计算时会产生很大的误差，从而得不到稳定的数值解。例如：用4阶龙格-库塔方法解：1)0(,20yyy （1）用解析法求解：本例是稳定的；（2）用数值法求解： 当h=0.1 时，数值解是稳定的; 当h=0.2时，数值解就不稳定了，这时因为：4)20(2 . 0hh此数值在稳定区间以外，所以数值解不能收敛。这种对步长有限制的数值积分法称为条件稳定积分法稳定积分法。从4阶龙格-库塔方法的稳定条件：中可以看出，系统的特征根越大，需要的积分步长就越小，这一点可作为选择步长的依据。步长的大小，除了与数值积分方法的阶数有关外，还与方程本身的性质

17、有关。078. 2h除以上介绍的欧拉法、龙格-库塔法外，数值积分方法还有许多种，如亚当姆斯方法、吉尔方法等等。此处不一一介绍。5.5 计算方法的选择计算方法的选择数值积分方法很多，在实际使用时存在一个选择问题。对于一个具体问题，如何选择具有一定的难度，至今尚无一种确定的方法。一般来说，数值积分方法的选择应考虑的因素有：1、精度要求、精度要求数值积分法的精度受截断误差舍入误差积累误差的影响，而这些误差与积分方法、步长、计算时间、计算机精度等有关。一般：（1）积分方法的阶数越高，截断误差越小，精度越高；（2）步长越小，精度越高；（3）多步法的精度高于单步法；（4）隐式算法的精度高于显式算法。 因此

18、，当需要高精度时，可采用高阶、多步、较小步长、隐式算法。 但是，步长的减小往往会增加迭代次数并增大舍入误差和积累误差。与此相反，在精度要求不高时，最好使用低阶方法。计算速度取决于计算步数、每步积分所需的时间。而每步的计算时间又与积分方法有关，它主要取决于计算导数的次数（4阶龙格-库塔方法每步要计算4次导数），导数的计算是最费时的；为了加快计算速度，在积分方法已定的条件下，应要保证精度的前提下，尽量选择较大步长，以缩短积分时间。2、计算速度、计算速度数值解的稳定性必须保证，否则，计算结果将失去真实意义。从稳定性来看，不同的数值积分算法有不同的稳定性。应用时应控制步长h，使数值积分算法在稳定域内。

19、3、数值稳定性、数值稳定性步长的选择很重要，步长过大会增大截断误差，甚至出现数值不稳定现象，过小了又因增加了步数，而使舍入误差增大。所以，仿真的总误差与步长的关系不是单调函数关系，而是一个具有极值的函数。如图所示。5.6 计算步长的选择计算步长的选择误差总误差截断误差舍入误差步长由图可知，存在一个最佳步长。当积分方法确定以后，在选取步长时，需要考虑的一个重要因素就是仿真系统的动态响应特性。如果系统的动态响应快，导数变化激烈，则应取高阶方法和小步长进行计算。为了保证数值稳定性，步长应限制在最小时间常数（相当于最大特征值的倒数）的数量级上。用经验方法经验方法选取步长的两种方法：（1）由系统方程中最

20、小时间常数Tmin来决定：min)05. 02 . 0(Th （2）由系统开环频率特性的剪切频率 来决定：c1)205(ch5.7 面向微分方程的仿真程序构成面向微分方程的仿真程序构成一般的仿真程序的组成可用下图表示：图中各方块的功能为：主程序：主程序：实现对整个仿真计算的逻辑控制。输入或预置参数模块：输入或预置参数模块：输入系统的参数初值、计算步长、计算时间等参数。运行管理模块：运行管理模块：这是数字仿真程序的核心，对仿真计算进行时间控制，以保证计算机按要求进行计算及输出。计算模块：计算模块：根据被仿真的系统及所选的仿真方法编写的计算程序。输出及显示模块：输出及显示模块：将仿真结果以数据或图

21、表的形式输出给用户。4阶龙格-库塔法仿真程序设计 目前，应用于系统仿真的商用软件包诸多，但无论是开发者还是应用者，了解程序的设计思想、分析程序的结构与功能，完善和编写仿真程序都是必要的。因此，此处通过对经典龙格-库塔方法（即式（7-22）仿真程序的介绍，使读者了解和掌握编写仿真程序的基本技术。 为了便于说明程序的实质，我们用具体的例子来讨论。 假设有一个2阶系统： )( 1, 0)0(, 0)0(,44828. 222tuyyuydtdydtyddtdyyyy21,令：则：122214828. 24,yyudtdyydtdy以下是用C语言编写的原理性程序：程序中变量说明：i，n: 分别为循环控制和积分器个数变量。h，t，tmax： 分别为步长、时间和仿真总时间变量。yn10: 存放tn时刻y值的临时变量数组。y10: 存放y值的变量数组。doty10: 存放 的数组。k110，k210，k310，k410:分别为存放K1，K2，K3，K4的数组。)( fy 程序清单 ：（“/*”和“* /”之间的文字是非执行语句，作为注释使用。 ）#includevoid diffeq();int i,n;float h, tmax, t, yn10,y10,doty10,k110,k21