第5章面向微分方程的数值积分法仿真数值积分是数值分析的一个基本问题。也是复杂计算问题中的一个基本组成部分。数值积分往往用极简单的方法就能较好地得出对所求解的具体数值问题的解答。但数值积分的难点在于计算时间有时会过长,有时会出现数值不稳定现象。另外,数值积分的理论性较强。其理论和方法都已经比较成熟,计算精度也比较高。5.1仿真中研究数值积分法的意义数值解的一种近似方法。对于连续系统的高阶微分方程,可化为若干个一阶微分方程组成的方程组。数值积分法是求解微分方程:00)(),(ytyytfdtdyy例如:下式所示的状态方程BuAxx可以化为一个一阶微分方程组(5-1)(5-2)ubxaxaxaxubxaxaxaxubxaxaxaxnnnnnnnnnnn2211222221212112121111所以,连续系统的仿真就是从给定的初始条件出发,对描述系统动态特性的常微分方程或常微分方程组进行求解,从而得到系统在一定输入作用下的变化过程。在求解这些微分方程时,最常用、也是最有效的一种方法就是数值积分法。(5-3)5.2数值积分法仿真的基本原理对微分方程(5-1)两端同时取积分,可得dyftytytt0),()()(0dyftytytt0),()()(0当时,上式变为:1ntt(5-4)dtytftytynttn10),()()(01(5-5)将积分项拆成两项dtytfdtytftytynnnttttn10),(),()()(01(5-6))(nty故上式可写为:dtytftytynnttnn1),()()(1(5-7)——此式是方程(5-1)在tn+1时刻的精确解。在数值解法中,希望用近似解:nnnQyy1代替准确解,其中:(5-8)1ny为)(1ntyny为)(ntynQ为1),(nnttdtytf的近似值httnn1令:——称为计算步长或步距式(5-8)是从常微分方程(5-1)出发建立的离散数学模型——差分方程。由此可见,数值积分法就是在已知微分方程初值的情况下,求解该方程在一系列离散点处的近似值,其特点是步进式——根据初始值逐步递推地计算出以后各时刻的值。从式(5-8)可知,数值积分法的主要问题归结为如何对f(t,y)进行数值积分——求出f(t,y)在区间[tn,tn+1]上定积分的近似值Qn。采用不同的方法求Qn,就出现了各种各样的数值积分方法。不同的数值积分将对求解的精度、速度和数值稳定性会产生不同的影响,这将在下述内容中具体介绍。121,,,,nntttt数值积分法种类繁多,在此从实用角度介绍几种基本的方法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(5-9)将此式写成差分方程为:,3,2,1,0),(1nhytfyynnnn(5-10)——著名的欧拉公式5.3.2改进的欧拉法如果用梯形面积而不是矩形面积来代替每一个小区间上的曲线积分,就可以提高计算精度,梯形法的计算公式为:)],(),([2111nnnnnnytfytfhyy(5-11)式中的右端含有待求量yn+1,因而它是隐函数形式。这种方法不能自行启动运算,需要依赖其它算法的帮助。每次计算都用欧拉法算出y(tn+1)的近似值,以此计算近似值,然后利用梯形公式(5-11)求出修正后的。即有:帮助方法:Pny1),(111PnnPnytfyCny1),(1nnnpnythfyy)],(),([2111PnnPnnnCnytfytfhyy预估式校正式(5-12)——改进的欧拉公式5.3.3几个基本概念简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn+1,这种方法称为单步法,它是一种自行启动的算法。如果求yn+1时需要tn,tn-1,tn-2……时刻yn,yn-1,yn-2……的值,这种方法为多步法(改进的欧拉法为两步法),它是一种不能自行启动的算法。1、单步法与多步法简单的欧拉法表达式的右端,计算所用的数据均已求出,这种公式称为显式公式。改进的欧拉法表达式的右端,有待求量,这种公式称为隐式公式。隐式公式不能自行启动,需要用预估-校正法。单步法和显式在计算上比多步法和隐式方便,但有时为了满足精度、稳定性等方面的要求,需要采用隐式算法。2、显式与隐式1ny1ny3、截断误差这里用泰勒级数为工具来分析数值积分公式的精度。假定yn是精确的,用泰勒级数表示处的精确解,即:1nt)(0)(!)(!2)()()(1)()2(2)1(1rnrrnnnnhtyrhtyhthytyty显然,简单的欧拉法是从以上精确解中取前两项之和来近似计算,每一步由这种方法引入的误差称为局部截断误差,简称截断误差。简单的欧拉法的截断误差为:)(0)(211hytynn不同的数值方法有不同的截断误差。一般若截断误差为,则方法为r阶的。所以方法的阶数可以作为衡量方法精确度的一个重要标志。)(01rh(5-13)由此可见,简单的欧拉法是1阶精度;改进的欧拉法由于采用了平均斜率,相当于取了泰勒级数的前3项,因此为2阶精度。分析欧拉法截断误差的思想,同样也适用于其它数值积分方法。4、舍入误差由于计算机的字长有限,数字不能表示得完全精确,在计算过程中,不可避免地会产生舍入误差。舍入误差与计算步长成反比。如果计算步长小,计算次数多,则舍入误差就大。舍入误差除了与计算机字长有关以外,还与计算机所使用的数字系统、数的运算次序以及计算所用的子程序的精度等因素有关。5、数值解的稳定性问题采用数值积分法求解稳定的常微分方程,应该保持原系统的稳定特征。但是:(1)在计算机逐步计算时,初始数据的误差及计算过程的舍入误差对后面的计算结果将产生影响。(2)如果计算步长取的不合格,有可能使仿真出现不稳定的结果。下面我们简单讨论一下这个问题。差分方程的解与微分方程的解类似,可分为特解和通解两部分。与稳定性有关的是方程的通解,它取决于差分方程的特征根是否满足稳定性条件。例如,为了考查欧拉法的稳定性,我们研究检验方程(TestEquation):yy其中,为方程的特征根,对此有:显然,要使该差分方程是稳定的,必须使下式成立。nnnnyhhyyy)1(111h即:2h表明:为使数字仿真稳定,对计算步长应有所限制。另外,稳定性还与系统的特性以及数值积分方法有关。上述分析欧拉法稳定性的思想,同样适用于其它数值积分方法。(5-14)(5-15)5.4龙格-库塔(Runge-Kutta)法由前面的分析可知,将泰勒展开式多取几项以后截断,就能提高截断误差的阶数和计算精度。然而,直接采用泰勒展开方法要计算函数的高阶导数,运用起来不便。龙格-库塔方法的基本思想是:用几个点上的函数值的线性组合代替函数的各阶导数,然后按泰勒级数展开确定其中的系数,这样既可避免计算高阶导数,又可提高积分的精度。龙格-库塔法有多种形式,以下从实用角度直接给出公式的形式和相应的精度。5.4.1龙格-库塔方法的基本思想5.4.2二阶龙格-库塔方法2阶龙格-库塔方法的公式为:),(),())(2/(121211hkyhtfkytfkkkhyynnnnnn(5-16)上式表示的数值解是用的泰勒级数在2阶导数以后截断所求得的,因此称为2阶方法。故2阶龙格-库塔法与式改进的欧拉法相比,实质完全相同。所以改进的欧拉法实质上是2阶龙格-库塔法。)(03h截断误差为:实时仿真的2阶龙格-库塔方法)2,2(),(12121khyhtfkytfkhkyynnnnnn(5-17))(03h其截断误差为:5.4.3四阶龙格-库塔方法4阶龙格-库塔法是一种最常用的方法。其经典表达式为:),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyynnnnnnnnnn(5-18)其截断误差为:)(05h显然,4阶龙格-库塔法的计算量较大,但计算精度较高,在比较不同算法的计算精度时,常以它的计算结果作为标准。实时仿真的4阶龙格-库塔方法)2103,54()52,53()52,52()5,5(),()105515(2441521413121543211khkhyhtfkhkkhyhtfkkhyhtfkkhyhtfkytfkkkkkkhyynnnnnnnnnnnn(5-19)以上公式都是标量形式,如果要换成向量形式,只要把式中的标量y,f,k换成向量Y,F,K即可。从理论上讲,可以构造任意阶数的龙格-库塔方法,但是,精度的阶数与计算函数值f的次数之间并非等量增加的关系,见下表所列:表5.1f的计算次数与精度阶数的关系每一步计算f的次数234567N≥8精度阶数234456n-2由此可见,4阶经典龙格-库塔方法有其优越性,而4阶以上的龙格-库塔方法计算f所需的次数比阶数多,增加了计算量,从而限制了更高龙格-库塔方法的应用。对于大量的实际问题,4阶方法已可满足精度要求,所以得到了广泛的应用。我们仍采用检验方程进行讨论,对它利用泰勒级数展开得: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,0)在使用龙格-库塔公式时,选取的步长应使落在稳定区域内。否则,在计算时会产生很大的误差,从而得不到稳定的数值解。例如:用4阶龙格-库塔方法解:1)0(,20yyy(1)用解析法求解:本例是稳定的;(2)用数值法求解:当h=0.1时,数值解是稳定的;当h=0.2时,数值解就不稳定了,这时因为:4)20(2.0hh此数值在稳定区间以外,所以数值解不能收敛。这种对步长有限制的数值积分法称为条件稳定积分法。从4阶龙格-库塔方法的稳定条件:中可以看出,系统的特征根越大,需要的积分步长就越小,这一点可作为选择步长的依据。步长的大小,除了与数值积分方法的阶数有关外,还与方程本身的性质有关。078.2h除以上介绍的欧拉法、龙格-库塔法外,数值积分方法还有许多种,如亚当姆斯方法、吉尔方法等等。此处不一一介绍。5.5计算方法的选择数值积分方法很多,在实际使用时存在一个选择问题。对于一个具体问题,如何选择具有一定的难度,至今尚无一种确定的方法。一般来说,数值积分方法的选择应考虑的因素有:1、精度要求数值积分法的精度受截断误差舍入误差积累误差的影响,而这些误差与积分方法、步长、计算时间、计算机精度等有关。一般:(1)积分方法的阶数越高,截断误差越小,精度越高;(2)步长越小,精度越高