SS03_数值积分法仿真

西南交通大学电气工程学院电子信息系[e-mail]系统仿真SystemSimulation第1页第三章数值积分法仿真第2页Overview数值积分方法的原理是什么？病态系统的特点和仿真算法选取？算法的稳定性分析？第3页第一节数字仿真原理在连续系统的仿真中，数值积分法可分为两大类：单步法：以龙格-库塔法为代表多步法：以Adams法为代表数值积分法的要素：基本特性：稳定性空间特性：精度时间特性：速度第4页数值积分基本原理连续系统的仿真，主要是对一阶微分方程（组）的求解110),()(),()()(),(01mmmttmttmdttyftydttyftytytyfdtdymmmmttmmmQyytyhtyfdttyfQtyymm11)(),(),()(,1为常数：在一个积分步长内近似当时间间隔足够小可见仿真关键是对Qm准确，快速的求解第5页步长：将时间t离散t(k)(k=1,2,…n)，相邻两点的距离为步长，即h=t(k+1)-t(k)步进法：数值积分法求近似解根据初始值y0，按照离散的时间序列步进求解。t0t1t2t3…tny0y1y2y3…tn计算格式：由y(k)计算出y(k+1)（k=0,1,…,n)的递推公式。数值积分基本名词第6页数值积分的基本性能数值积分算法的性能包含：定性特征：稳定性时间粒度：计算速度空间粒度：计算精度不同的数值积分方法的具有不同的稳定性。同一个模型采用不同的积分算法和不同的积分步长h，稳定性不同。第7页计算速度和计算精度各种数值积分方法的差分方程是对原微分方程的近似逼近，并且因为计算机的字长有限，存在明显的截断误差。这些误差都和计算步距h密切相关，所以计算步距是影响计算精度、速度和稳定性的重要因素。h取得较大，计算时间少，截断误差大；h取得较小，截断误差就会减小，但在给定时间范围内，计算次数必然增加，使误差积累增加。第8页截断误差、累计舍入误差与步长h截断误差、累计舍入误差与步长h关系如图。图中可知，两种误差对步距的要求是矛盾的，但两者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难的。一般只能根据经验确定一个附近的合理步长区，通常将步长h限制在系统的最小时间常数数量级上。第9页引理：泰勒级数:如果f(x)在x0点处任意阶可导,则在该邻域内的n阶泰勒公式为：第二节单步法)()(!)()(000)(xRxxnxfxfnnNnn称为拉格朗日余项其中：10)1()()!1()()(NNNxxNfxR单步数值积分法的核心就是泰勒级数的近似。第10页2.1一阶欧拉法对于一阶微分方程00)(),,(ytytyfdtdyhtyfyxRhdtdyyyhn),()]([000001项（欧拉递推）只保留)()(!)()(000)(00xRxxnxfxfxxnnNnn附近开始，根据泰勒公式，从起点),(1mmmmtyfKhKyy故一般的一阶欧拉法递推形式为：第11页一阶欧拉法图示2h截断误差数项略去了2阶以上高阶导阶导数项，一数式中的欧拉公式只取到泰勒级)(tytmy1myhmt),(tyfdtdy),(1mmmmtyfKhKyy1mtthmt),(1mmmmtyfKhKyy1mtdtdy第12页2.22阶龙格-库塔对于一阶微分方程00)(),,(ytytyfdtdy20000022201200000]),([21),()]([,htftyfyfhtyfyxRhdtydhdtdyyyhtttyyttyyn！保留附近展开泰勒级数在tfyftyfdtdyhytt,)),((,,,00000未知项为：及已知项为：对于上式，在当前时刻第13页2阶龙格-库塔),(),,()(][21),(101202001221102000100hbthKbyfKtyfKhKaKayhtftyyfhtyfyyttyy其中假设上式可写成hyfKbtfbtyfhbthKbyfKtK][),(),(1210010120202级数附近进一步展开为泰勒在将}][),({),(12100200101hyfKbtfbtyfhatyfhayy带入上式整理可得：第14页2阶龙格-库塔2000112100200101][21),(][),({),(:htftyyfhtyfyyhyfKbtfbtyfhatyfhayy！对比1,121,2121211212122122121bbaababaaaaa补充约束条件）（从而：21012KKhyy第15页2阶龙格-库塔)],(),([2)],(),([22),(),(1m*11m211121hKyhtfytfhyytfytfhyKKhyyhthKyfKtyfKmmmmmmmmmmmmmm）（故一般的二阶龙格-库塔法递推形式为：2阶龙格-库塔只取到泰勒级数展开式中y的二阶导数项，略去了三阶以上高阶导数项。其截断误差正比于步长h3为纪念提出该方法的德国数学家C.Runge和M.W.Kutta，称这种计算方法为二阶龙格-库塔法。第16页2阶龙格-库塔图示）（2111212),(),(KKhyyhKyhtfKytfKmmmmmm)(tytmy1myhmt1mtthmt),(1mmmmtyfKhKyy1mtdtdy）（2121KK1K),(112mmytfK预报校正第17页比较)(tytmy1myhmt),(tyfdtdy),(1mmmmtyfKhKyy1mtthmt),(1mmmmtyfKhKyy1mtdtdy）（2111212),(),(KKhyyhKyhtfKytfKmmmmmm)(tytmy1myhmt1mtthmt),(1mmmmtyfKhKyy1mtdtdy）（2121KK1K),(112mmytfK预报校正110),()(),()()(),(01mmmttmttmdttyftydttyftytytyfdtdy第18页高阶龙格-库塔(RK-4)),()2,2()2,2(),(226342312143211hthKyfKhthKyfKhthKyfKtyfKKKKKhyymmmmmmmmmm）（一般在计算精度要求较高的情况下，多使用四阶龙格-库塔法。其计算公式为,其截断误差正比于步长h5第19页高阶龙格-库塔(RK-4)),()2,2()2,2(),(226342312143211hthKyfKhthKyfKhthKyfKtyfKKKKKhyymmmmmmmmmm）（份。点各取，份，各取，精度高，所以加权时点步长小、由于第和，式中取四点导数的加权是在各估值处的导数，14122h32iK第20页单步法的特点以上介绍的几种数值积分公式，有一个共同的特点，由于采用了泰勒级数展开，在本次计算中，仅仅用到前一步的计算结果，而不需要利用更前面各步的结果。这类计算方法称为单步法。单步法运算有下列优点：(1)需要存储的数据量少，占用的存储空间少。(2)给定初值，就可启动递推公式进行运算(自启动计算能力)(3)容易实现变步长第21页第三节变步长龙格-库塔法步长控制是实现高精度的仿真算法的手段之一。实现步长控制涉及：局部误差估计步长控制策略步长控制策略积分算法E（允许误差）+-误差估计h=?-e第22页3.1误差估计通常设法寻找一个低一阶的龙格-库塔公式，两者的结果之差可以设为误差。为减少计算量，Ki通常要求公用。Runge-Kutta-Merson法(RK34)),2)34(()2,8)3(()3,6)(()3,3(),(3415314213121hthKKKyfKhthKKyfKhthKKyfKhthKyfKtyfKmmmmmmmmmm）（误差）（级公式：三阶）（四阶五级公式：543143115411-89-26:129-36ˆ446KKKKhEKKKhyyKKKhyymmmmm）（计算结果取：541146KKKhyymm第23页Runge-Kutta-Fehlberg(RK45)计算公式为5阶6级，误差估计低阶公式为4阶五级，具有四阶误差估计和五阶精度，称为RK45法。RK45被公认为对非病态系统仿真的最有效的方法之一。61*111111161611)(-ˆ,0:),(1*iiiimmmijijiijjijkikiiiiiimmKcchyyEbaaKbhyhatfKCKChyy估计误差为：第24页Runge-Kutta-Fehlberg(RK45)iaibijcic*i1016/13525/21621/41/40033/83/329/326656/128251408/2565412/131932/2197-7200/21977296/219728561/564302197/410451439/216-83680/513-847/4104-9/50-1/561/2-8/272-3544/25661859/4104-11/402/550第25页RKF-12）（误差）（低阶公式：）（计算公式：各点导数：31112113211213121-512ˆ:255256ˆ510512),256)255(()2,2(),(KKhyyEKKhyyKKKhyyhthKKyfKhthKyfKtyfKmmmmmmmmmmmmm第26页RKS-34（1978，Shamping））（误差步的刚好为（注：）（低阶公式：）（计算公式：各点导数：543211115*5*4321143211321421312143-3-332ˆ:)14915332ˆ338),)(()32,3)3(()3,3(),(KKKKKhyyEKmKKKKKKhyyKKKKhyyhthKKKyfKhthKKyfKhthKyfKtyfKmmmmmmmmmmmmmmm第27页3.2步长控制步长控制策略一般分为：1）加倍-减半法2）最优步长法第28页步长控制：加倍-减半法加倍-减半法)1(:mmmyEe定义每步的局部误差为mmmmmmmmmhhhhhhtheneeeif221111minmaxminmax步长控制策略：第29页步长控制：最优步长法)1()()(:),(KmKmmmKmyhtehhEktf足够小，当步长某处变量偏导数的组合在积分区间可表示为阶积分算法，误差估计对于mKmKmmmmKmmmKmmmmhetyhyhthyhtee101010111110)()1()1()()1()(,)1最优步长：足够小下一步误差一步的步长本次积分成功，确定下若第30页1.2.2步长控制重新计算本步。的步长作为本次步长本次积分失败，计算处若,)20