第八章基于MATLAB的科学计算—常微分方程数值解法

科学计算理论、方法及其基于MATLAB的实现与分析常微分方程（组）数值解法1常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的NEWTON第二定律（）002XTXTFDM和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题（）0YTFADT的解为（）DFETYTAAT0但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题（）0,YTTTFDF其中（）TYTTN21（）YTFTFTFN,21常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。1常微分方程（组）的CAUCH问题数值解法概论假设要求在点（时刻），处初值问题（）的解KTN,21的（近似）值,如果已求得时刻的值或它的近似值0,YTKTKTY如时刻的值，那么将式（）的两端在区间上积分K0T1,KT10DTYTFTYTDTYKKTKKT11,1可得11DTYTFTYTKTKK1,1或12DTYTFYTYKTTKK1,11显然,为了利用式（11）或（12）求得的精确值（近似值），1K1KY必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解,可以说,不同的式值积分方法将给出不同的CAUCH问题的数值解法。12最简单的数值解法EULER方法假设要求在点（时刻），处初值问KHTK0NTF0NK,21题（）的解的近似值。首先对式（）的两端积分，得TY（13）DTYFTYTDTKKTKT11,1对于式（13）的右边，如果用积分下限处的函数值代替KTKTYF,被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面积），则有（14）KTKKTYHFDTYFYTK,11进而得到下式给出的递推算法EULER方法（15）01,21,YTNTFKK例用EULER方法解如下初值问题，取，30H10SIN23YTTDT解由（15）得1010,2SIN633YKTYKKK结果如下OPENEULER_METHODM05152537891124M如果取，其结果如下图所示1HEULER_METHOD5553781123413改进的EULER方法对于（15）的右边，如果被积函数用积分限和处的函数KT1KT值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有（16）11,2,1KKTKKTYFTYFHDTYFYTK进而得到下式给出的递推算法（17）011,2,2YTNKYTFYTFHKKK通常算法（17）比EULER方法（15）的精度高，但是，按算法（17）求时要解（非线性）方程（组），这是算法（17）不如EULER方法1KY的方面，为了）尽可能地保持算法（17）精度高的优点；）尽可能地利用EULER方法计算简单的长处；人们采取了如下的称之为改进的EULER方法的折衷方案预测（18）KKYTHFY,01修正（19）NKYTFKKK,21,2010YT例EULER方法与改进的EULER方法的比较下图是当时比较的结果3HOPENIMPROVED_EULER_METHODM051525381234D14EULER方法和改进的EULER方法的误差分析由TAYLOR公式（19）22111,HOTYFTYTTYTKKKKK说明EULER方法的截断误差是，类似地，由（20）3211,HOTYTFTYTKKKK（21）32111,THTFTTKKKK以及22HOTYTHTYTTYKKKKK11让式（20）的两端减式（21）的两端，可得（23）31321,24,HOTYFTYFHTYKKKKKK从上述推导EULER方法、改进的EULER方法的过程以及例、例容易看出，改进的EULER方法EULER方法的精度高，其原因在于1在推导EULER方法时，我们是用待求解函数在一点处的变化TY率代替在区间上的平均变化率KTYF,TY,1KT（24）TYHFDYFKT,12而在推导改进的EULER方法时，我们是用待求解函数在两点TY处变化率的平均值代替在区间上1,21KKTYFTYFTY,1K的平均变化率；显然，通常比更接近于在1,21KKTYFTYFKTYF,TY区间上的平均变化率。由此启发人们适当地选取区间上,1KT,1K函数若干点处的变化率，用它们加权平均值代替在区间YTY上的平均变化率，近似解的精度应更高。,1KT下面将要介绍的RUNGEKUTTA法就是基于上述想法得到的。2RUNGEKUTTA法RUNGEKUTTA法是按选取区间上函数变化率的,1KTTYYTF,个数的多少和截断误差的阶数来区分的一系列方法，如MHO二阶的RUNGEKUTTA法（改进的EULER方法）（25）21121,KHYTFKKK三阶的RUNGEKUTTA法（26）3211312146,KHYTFKYTFKKK四阶的RUNGEKUTTA法1古典形式（27）432114231216,KKHYTFHYTFKKKKK2GILL公式具有减小舍入误差的优点（28）432114213216,2,KKHYHTFKYHTFTFKKKKKKK4RUNGEKUTTA法的一般形式29NIIKKKKAHYTHY11,其中称为增量函数INCREMENTFUNCTION以及TK,（30）1,11,121223121,NNNKNKNKKKHKQHKQYHPTFKTFYTF需要特别指出的是在确定（29）、（30）中的参数，和时，IAIPIJQ应该使（29）的右端和适当阶的（如阶）TAYLOR展式一致，这样，至少对于底阶的RK法来说，参与加权的斜率个数与方法的阶数是一N致的。例如一阶的RK法即EULER法，局部的截断误差（TRUNCATIONERROR）是，所以，当微分方程的阶是一次函数时是精确的；二阶2HORK法即改进的EULER法，局部的截断误差是，所以，当微分方3HO程的阶是二次函数时是精确的；类似地，四阶RK法局部的截断误差是。5H但是，一般五阶以及五阶以上的RK法的局部截断误差不再具有上述的特点，因此用“性价比”来衡量，四阶RK得到了广泛的应用。下面是BUTCHERSFIFTHORDERRK方法（1964）（28）654311722790KKHYK（29）543216453242132781277,69,8,HKHKHKYTFKYHTFHKKTFKKKKKK4线性多步法与常微分方程数值解法的分类一般的常微分方程初值问题的数值解法都是以递推的形式给出的，即是递推算法，根据递推算法，在计算时，已经得到了前面各时刻1KY的近似值，前面介绍的各种数值解法都有一个共同点IYK,10在计算时，只用到了前一个时刻（当前时刻）的“信息”预测1KYKTKY未来某时刻系统的“状态”，这样的数值解法称为单步法，对T1KTY于一个动态过程在时刻的状态而言，不仅前一个时刻的TY1K1KT信息对它有影响，前若干个时刻的信息通常对它也有影响，显然，KY单步法的缺点是没有充分利用已得到的信息。基于上述考虑，适当取前若干个时刻的信息，IYMK,1,并用的线性组合代替在区间上的平均变化率所给出IYTF,TY,1KT的方法称为线性多步法。41线性多步法从二步法谈起在改进的EULER方法中预测（48）KKYTHFY,01较正（49）NKYTFKKK,21,2010YT如果改进预测的手段，如利用数值微分的研究成果503112,2HOTYFTYTHXFXFKKKKK那么，由于（48）的局部截断误差是，精度高于EULER公式的局3部截断误差，所以，下列预测校正方法2HOPREDICTOR（51）KKYTHFY,210CORRECTOR（52）NKTFKKK,21,010YT尽管步长增大了，精度却提高了。公式（49）、（50）给出的数值方法称为非自始的（NONSELFSTARTING）二步法，其原因在于利用该方法求解微分方程初值问题时，自身不能提供所需的初值。1Y43线性多步法的进一步扩展在积分公式（60）DTYTFYDTYTFTTDKNKNKNTKKTNY111,1中，用个节点，的次插值多项IIITFT,1,KKN式代替被积函数，得到的就是一般的线性多步法。TPNY如当时，1（61）KKYTHF,21如当时，2N（62）KKKYTFYTFY,2311特别地，用四个节点，的三次插值IIIT,12,3多项式代替被积函数，得到TP3YF（63）3211975924KKKKKFFFHY四个节点，的三次插值多项式代替，IIITFT,TP3又得到（64）21115924KKKKFFFHY其中误差分别为、，可见，公式（64）的精度57070Y比公式（63）的精度高。ADAMS预报校正公式（65）2111359249375NNNNFFFHY其中，NTF,11,NTF说明公式（65）自身不能提供最初所需的三个初值，需要RUNGEKUTTA法的配合。44一般线性多步方法一般的步线性多步方法可以写成L110LIIKLIKKFY（824）1011LKLKKLKLKFFH这里，当时，上述公式就是显式公式，时为单步,PPYXF法1、数值积分构造方法数值积分方法是构造线性多步法的一种途径ADAMS方法就是利用插值多项式进行积分得出来的又如，由SIMPSON求积公式有，901435111KKKKKKYHXFFXFHXY其中，对应的数值公式XH（825）111KKKFFY称为SIMPSON方法SIMPSON方法是四阶方法，它的局部截断误差为9015KKYHT2、TAYLOR展开构造方法构造线性多步方法的另一重要途径是利用TAYLOR展开方法例如确定公式（824）中的个参数，方法是将12L11010,LL写成，并将所有项在进行TAYLOR展开，比较两边相同IKFIKXYKXH幂次前的系数，若,2111,100MJIJIJLILISLI，（826）110LIMILIMM则公式（824）是阶方法，局部截断误差主项为，11KMKXYHCT其中1101LIMILIMMC一般应用的线性多步方法都是大于或等于一阶的方法，这样的方法亦称为与微分方程（83）相容的方法阶方法的构造是选择参数P，使（826）中的为11010,LLM例6推导最高阶的二步线性多步方法解1）二步显式线性多步方法为10101KKKKFHYY共四个待确定参数由公式（826），令,26,2,11110110解得又因为，则得阶数最高2,4,5,4101004的二步显式线性多步方法为，2454111KKKFHYY其局部截断误差主项为，是三阶方法61KKXT2）二步隐式线性多步方法为101101KKKKKFFHYY共五个待确定参数由公式（826），令,2,1110110,64,261111解得则得SIMPSON方法3/,3/4,001111KKKFFHY【注】两种多步法的构造公式中，数值积分方法是有局限的，它只对能将微分方程（83）转化为等价的积分方程的情形才能适用，而用TAYLOR展开则可构造任意多步法公式3、几个常用的线性多步法的预测校正公式（1）ILNE预测校正法ILNESPREDICTORCORRECTOR预测公式，2234101KKKF