火电厂仿真课件第5章 面向微分方程的数值积分法仿真

第五章面向微分方程的数值积分法仿真数值积分是数值分析的一个基本问题。也是复杂计算问题中的一个基本组成部分。数值积分往往用极简单的方法就能较好地得出对所求解的具体数值问题的解答。但数值积分的难点在于计算时间有时会过长，有时会出现数值不稳定现象。另外，数值积分的理论性较强。其理论和方法都已经比较成熟，计算精度也比较高。5.1仿真中研究数值积分法的意义5.2数值积分法仿真的基本原理5.3欧拉(Euler)法5.4龙格-库塔（Runge-Kutta）法5.5计算方法的选择5.6计算步长的选择5.7面向微分方程的仿真程序构成5.1仿真中研究数值积分法的意义由模拟仿真可知，一个n阶连续系统可以由n个积分器组成的模拟仿真模型予以描述。因此，利用数字计算机进行连续系统仿真时，从本质上讲，就是在数字计算机上构造出几个数字积分器，进行几次数值积分运算。在系统的仿真中，连续系统的数学模型一般都可用微分方程来表示，但它们通常是高阶的微分方程，由前所述，高阶的微分方程可化为若干个一阶微分方程组成的一个方程组。（5-2）或用状态方程表示：（5-1）所以，连续系统的仿真就是从给定的初始条件出发，对描述系统动态特性的常微分方程或常微分方程组进行求解，从而得到系统在一定输入作用下的变化过程。在求解这些微分方程(组)时，最常用、也是最有效的一种方法就是数值积分法。数值解的一种近似方法。数值积分法是求解一系列微分方程：（5-3）5.2数值积分法仿真的基本原理若对微分方程（5-3）两端同时取积分，可得当时，上式变为：（5-4）（5-5）将积分项拆成两项（5-6）故上式可写为：（5-7）——此式是方程（5-3）在tn+1时刻的精确解。

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

欧拉法是一种最简单的数值积分法，对于方程：在区间[t0，

tn+1]上求积分，得到：f(t,y)在[tn,tn+1]区间上包围的面积Qn若区间[tn，tn+1]足够小，则[tn，tn+1]上的f（t，y）可近似地看成常数f（tn，yn）

。故可用矩形面积Qn近似代替即：tntn+1f（tn，yn）y(tn)Qn于是有：（5-9）其中：——称为计算步长或步距将此式写成差分方程为：（5-10）——著名的欧拉公式（n=1，2，3……）5.3.2改进的欧拉法

如果用梯形面积而不是矩形面积来代替每一个小区间上的曲线积分，就可以提高计算精度，梯形法的计算公式为：（5-11）式中的右端含有待求量yn+1，因而它是隐函数形式。这种方法不能自行启动运算，需要依赖其它算法的帮助。每次计算都用欧拉法算出y(t

n+1

)的近似值，以此计算近似值，然后利用梯形公式（5-11）求出修正后的。即有：帮助方法：预估式校正式（5-12）——改进的欧拉公式5.3.3几个基本概念

简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn+1，这种方法称为单步法，它是一种自行启动的算法。如果求yn+1时需要tn,

tn-1,

tn-2……时刻yn,

yn-1,

yn-2……的值，这种方法为多步法（改进的欧拉法为两步法），它是一种不能自行启动的算法。1、单步法与多步法简单的欧拉法表达式的右端，计算所用的数据均已求出，这种公式称为显式公式。改进的欧拉法表达式的右端，有待求量，这种公式称为隐式公式。隐式公式不能自行启动，需要用预估-校正法。单步法和显式在计算上比多步法和隐式方便，但有时为了满足精度、稳定性等方面的要求，需要采用隐式算法。2、显式与隐式3、截断误差这里用泰勒级数为工具来分析数值积分公式的精度。假定yn是精确的，用泰勒级数表示处的精确解，即：显然，简单的欧拉法是从以上精确解中取前两项之和来近似计算，每一步由这种方法引入的误差称为局部截断误差，简称截断误差。简单的欧拉法的截断误差为：不同的数值方法有不同的截断误差。一般若截断误差为，则方法为r阶的。所以方法的阶数可以作为衡量方法精确度的一个重要标志。（5-13）由此可见，简单的欧拉法是1阶精度；改进的欧拉法由于采用了平均斜率，相当于取了泰勒级数的前3项，因此为2阶精度。

分析欧拉法截断误差的思想，同样也适用于其它数值积分方法。4、舍入误差由于计算机的字长有限，数字不能表示得完全精确，在计算过程中，不可避免地会产生舍入误差。舍入误差与计算步长成反比。如果计算步长小，计算次数多，则舍入误差就大。舍入误差除了与计算机字长有关以外，还与计算机所使用的数字系统、数的运算次序以及计算所用的子程序的精度等因素有关。5、数值解的稳定性问题采用数值积分法求解稳定的常微分方程，应该保持原系统的稳定特征。但是：（1）在计算机逐步计算时，初始数据的误差及计算过程的舍入误差对后面的计算结果将产生影响。（2）如果计算步长取的不合格，有可能使仿真出现不稳定的结果。下面我们简单讨论一下这个问题。差分方程的解与微分方程的解类似，可分为特解和通解两部分。与稳定性有关的是方程的通解，它取决于差分方程的特征根是否满足稳定性条件。显然，要使该差分方程是稳定的，必须使下式成立。即：表明：为使数字仿真稳定，对计算步长应有所限制。另外，稳定性还与系统的特性以及数值积分方法有关。上述分析欧拉法稳定性的思想，同样适用于其它数值积分方法。（5-14）（5-15）例如：为了考查欧拉法的稳定性，我们研究检验方程（TestEquation）：其中，为方程的特征根，对此有：5.4龙格-库塔（Runge-Kutta）法由前面的分析可知，将泰勒展开式多取几项以后截断，就能提高截断误差的阶数和计算精度。然而，直接采用泰勒展开方法要计算函数的高阶导数，运用起来不便。龙格-库塔方法的基本思想是：用几个点上的函数值的线性组合代替函数的各阶导数，然后按泰勒级数展开确定其中的系数，这样既可避免计算高阶导数，又可提高积分的精度。龙格-库塔法有多种形式，以下从实用角度直接给出公式的形式和相应的精度。5.4.1龙格-库塔方法的基本思想5.4.2二阶龙格-库塔方法2阶龙格-库塔方法的公式为：（5-16）上式表示的数值解是用的泰勒级数在2阶导数以后截断所求得的，因此称为2阶方法。故2阶龙格-库塔法与式改进的欧拉法相比，实质完全相同。所以改进的欧拉法实质上是2阶龙格-库塔法。截断误差为：实时仿真的2阶龙格-库塔方法:（5-17）其截断误差为：5.4.3四阶龙格-库塔方法4阶龙格-库塔法是一种最常用的方法。其经典表达式为：（5-18）其截断误差为：显然，4阶龙格-库塔法的计算量较大，但计算精度较高，在比较不同算法的计算精度时，常以它的计算结果作为标准。实时仿真的4阶龙格-库塔方法（5-19）以上公式都是标量形式，如果要换成向量形式，只要把式中的标量y，f，k换成向量Y，F，K即可。从理论上讲，可以构造任意阶数的龙格-库塔方法，但是，精度的阶数与计算函数值f

的次数之间并非等量增加的关系，见下表所列：由此可见，4阶经典龙格-库塔方法有其优越性，而4阶以上的龙格-库塔方法计算f所需的次数比阶数多，增加了计算量，从而限制了更高龙格-库塔方法的应用。对于大量的实际问题，4阶方法已可满足精度要求，所以得到了广泛的应用。我们仍采用检验方程进行讨论，对它利用泰勒级数展开得：5.4.4龙格-库塔公式的稳定区域对于，有，将它代入上式得：（5-20）（5-21）令：——将其代入上式得到该式的稳定条件为：（5-22）由此稳定条件，下表给出了各阶龙格-库塔公式的稳定区域。表5.2龙格-库塔公式的稳定区域在使用龙格-库塔公式时，选取的步长应使落在稳定区域内。否则，在计算时会产生很大的误差，从而得不到稳定的数值解。例如：用4阶龙格-库塔方法解：（1）用解析法求解：本例是稳定的；（2）用数值法求解：当h=0.1时，数值解是稳定的;当h=0.2时，数值解就不稳定了，这时因为：此数值在稳定区间以外，所以数值解不能收敛。这种对步长有限制的数值积分法称为条件稳定积分法。从4阶龙格-库塔方法的稳定条件：中可以看出，系统的特征根越大，需要的积分步长就越小，这一点可作为选择步长的依据。步长的大小，除了与数值积分方法的阶数有关外，还与方程本身的性质有关。除以上介绍的欧拉法、龙格-库塔法外，数值积分方法还有许多种，如亚当姆斯方法、吉尔方法等等。此处不一一介绍。5.4.5四阶龙格-库塔法仿真程序设计目前，应用于系统仿真的商用软件包诸多，但无论是开发者还是应用者，了解程序的设计思想、分析程序的结构与功能，完善和编写仿真程序都是必要的。因此，此处通过对经典龙格-库塔方法（即式（5-18））仿真程序的介绍，使读者了解和掌握编写仿真程序的基本技术。为了便于说明程序的实质，以下用具体的例子来讨论。假设有一个2阶系统：令：则：以下是用C语言编写的原理性程序：

程序中变量说明：i，n: 分别为循环控制和积分器个数变量。h，t，tmax：分别为步长、时间和仿真总时间变量。yn[10]:

存放tn时刻y值的临时变量数组。y[10]:

存放y值的变量数组。doty[10]:

存放的数组。k1[10]，k2[10]，k3[10]，k4[10]:分别为存放K1，K2，K3，K4的数组。程序清单：（“/*”和“*

/”之间的文字是非执行语句，作为注释使用。）#include<stdio.h>voiddiffeq();inti,n;floath,tmax,t,yn[10],y[10],doty[10],k1[10],k2[10],k3[10],k4[10];main() /*主函数*/{n=2; /*设置积分器个数*/y[1]=0;y[2]=0 /*设置积分器初值*/

h=0.01;tmax=20;t=0/*设置步长、仿真总时间和初始时间*/

whilt(t<=tmax) /*判断是否仿真时间到*/{printf(“t=%f”,t)； /*打印时间*/

printf(“y[1]=%fy[2]=%f\r\n”,y[1],y[2]);/*打印积分器输出*/for(i=1,i<=n,i++)

yn[i]=y[i]; /*保留tn时刻的初值*/dif(); /*计算*/for(i=1;i<=n;i++){k1[i]=doty[i]; /*计算K1*/y[i]=yn[i]+(h/2)*k1[i];}t=t+h/2dif(); /*计算*/for(i=1;i<=n;i++){k2[i]=doty[i]; /*计算K2*/y[i]=yn[i]+(h/2)*k2[i];}dif(); /*计算*/for(i=1;i<=n;i++){k3[i]=doty[i]; /*计算K3*/y[i]=yn[i]+(h/2)*k3[i];}t=t+h/2dif(); /*计算*/for(i=1;i<=n;i++){k4[i]=doty[i]; /*计算K4*/y[i]=yn[i]+(h/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);

/*计算tn+1时刻的y值*/}}}voiddif() /*计算变量导数的函数*/{floatu;u=1; /*计算u=1(t)*/

doty[1]=y[2];

doty[2]=4*u-2.828*y[2]-4*y[1];}下面是程序的主函数main()的处理流程图

编写以下实时仿真的4阶龙格-库塔法的仿真程序,对以下2阶系统进行数字仿真:练习:并通过仿真实验,验证使该算法稳定的步长h取值区间.5.5计算方法的选择数值积分方法很多，在实际使用时存在一个选择问题。对于一个具体问题，如何选择具有一定的难度，至今尚无一种确定的方法。一般来说，数值积分方法的选择应考虑的因素有：1、精度要求数值积分法的精度受截断误差舍入误差积累误差的影响，而这些误差与积分方法、步长、计算时间、计算机精度等有关。一般：（1）积分方法阶数越高，截断误差越小，精度越高；（2）步长越小，精度越高；（3）多步法的精度高于单步法；（4）隐式算法的精度高于显式算法。因此，当需要高精度时，可采用高阶、多步、较小步长、隐式算法。但是，步长的减小往往会增加迭代次数并增大舍入误差和积累误差。与此相反，在精度要求不高时，最好使用低阶方法。计算速度取决于计算步数、每步积分所需的时间。而每步的计算时间又与积分方法有关，它主要取决于计算导数的次数（4阶龙格-库塔方法每步要计算4次导数），导数的计算是最费时的；为了加快计算速度，在积分方法已定的条件下，应要保证精度的前提