第章常微分方程初值问题的数值解法

1第7章常微分方程初值问题的数值解法71引言科学研究和工程技术中的许多问题在数学上往往归结为微分方程的求解问题。为了确定微分方程的解，一般要加上定解条件，根据不同的情况，这些定解条件主要有初始条件INITIALCONDITION和边界条件BOUNDARYCONDITION只含初始条件作为定解条件的微分方程求解问题称为初值问题INITIALVALUEPROBLEM例如天文学中研究星体运动,空间技术中研究物体飞行等，都需要求解常微分方程初值问题INITIALVALUEPROBLEMFORORDINARYDIFFERENTIALEQUATIONS只含边界条件作为定解条件的微分方程求解问题称为边值问题BOUNDARYVALUEPROBLEM除特殊情形外，微分方程一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，而且实际问题往往只要求在某一时刻解的函数值为了解决这个问题，有两种方法可以逼近原方程的解。第一种方法是将原微分方程化简为可以准确求解的微分方程，然后使用化简后的方程的解近似原方程的解；第二种方法是将求原微分方程的解析解转化为求原方程的数值解，这是实际中最常用的方法。本章将介绍求解常微分方程初值问题的常用的数值方法。第8章将介绍常微分方程的边值问题的常用的数值方法。为简明起见，本章主要介绍形如711,YTFATBA的初值问题的数值解法在介绍这些方法之前，需要了解常微分方程的一些相关定义和结果定义711函数称为在集合上关于变量满足LIPSCHITZ条件,FTY2DRYLIPSCHITZCONDITION，简称LIP条件，如果存在常数，使得0L7121212,FTFTY对所有都成立常数称为LIPSCHITZ常数LIPSCHITZCONSTANT，简称LIP12,TYD常数定义712如果对所有，都有12,TYTD，71312TY其中，那么称集合为凸集CONVEXSET01R注如果没有特别说明，本章中的集合为，2,TATBY显然它是凸集。定理711设函数在上连续如果,FTY,DTYATBY在上关于变量满足LIPSCHITZ条件，那么初值问题711对有唯一解,FTYDT2YT定理712设函数定义在凸集上如果存在常数，使得,FTY2DR0L714,YFTL对一切成立，那么在上关于变量满足LIPSCHITZ条件，且为,TYD,FTYLIPSCHITZ常数注因为LIPSCHITZ条件一般不易验证，所以往往用条件714代替LIPSCHITZ条件712。由定理711立得定理713设函数在上连续如果存在常,FTY,DTYATBY数，使得对一切成立，那么初值问题711对有唯一0L,YFTLATB解YT理论上，初始值是准确的，微分方程的计算过程也是准确的；而实际情况并非如此事实上，初始数据可能有误差，计算过程中的数字的舍入也会产生误差。对初值问题711，若初值有误差，即实际得到的初值是，再考虑到计算过程中的数字的舍入也会00产生误差，则得到的初值问题变为7150,ZTFZTATBA这种误差的扰动在计算过程中是否会增长很快，以至影响结果，这就是初值问题的适定性问题定义713设初值问题711和715的均有唯一解，分别是和，且YTZT连续。若对任意，存在常数，当时，对任意T00K0,，恒有，则称初值问题711是适定的WELLPOSED。,ABYTZ定理716设函数在上连续如果,FTY,DTYATBY在上关于变量满足LIPSCHITZ条件，那么初值问题711是适定的。,FTD注初值问题只有是适定的才有意义，因此本章考察的所有初值问题都是适定的。求微分方程数值解的主要问题1如何将原微分方程离散化,并建立求其数值解的递推公式2如何求递推公式的局部截断误差,数值解与精确解的误差估计；NYNT3研究递推公式的稳定性与收敛性372EULER方法及改进的EULER方法721EULER格式与梯形格式考虑一阶常微分方程的初值问题721,YTFATBA设，其中为等距节点，步长011NATTB0,1,NTHNHB在上对两边积分，得1,0,NT,YFT72211DNNTTYT7211EULER格式用左矩形求积公式计算式722右端积分项，得,1,D,NTNFYTHFTY代入式722右端，得7231,NNNTTFT用的近似值代入式723右端，记所得结果为，得NYTNY1Y72410,0,NNHFTYN并称式724为求解初值问题721的EULER方法EULERSMETHOD或EULER格式EULERSCHEME，称为差分方程DIFFERENCEEQUATION1,NNYFTY注EULER方法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它的精度不高，很少被采用，但是它反映了微分方程的数值解法的基本思想和特征若式722右边的积分由数值积分的右矩形公式近似，并用近似值替代，近NYNT似值替代，则可得到1NY1NT72510,0,1,NNYHFTYN并称式725为后退的EULER方法BACKWARDEULERSMETHOD或后退的EULER格式BACKWARDEULERSSCHEME是差分方程。11,NNYFTY注在平面上，微分方程的解称为积分曲线，积分曲线上一XOYT点的切线斜率等于函数在点的值如果在,T,TFX中每一点都画一条以在点的值为斜,DATBYTY,FTY,T欧拉LEONHARDPAULEULER1707年4月15日1783年9月18日是瑞士杰出的数学家，自然科学家他在数学分析，图论，力学，光学和天文学等很多领域都做出了卓越的贡献，被公认为有史以来最伟大的数学家之一4率并指向增加方向的有向线段即在上作出了一个由方程确定的方向场），TD,YFT那么从几何上看，微分方程的解就是位于此方向场中的曲线，它在所,YFTYT经过的每一点都与方向场的该点的方向相一致从初始点出发，过这点的积分曲线为，斜率为设在0,PAT00,YFA附近可用过点的切线近似表示，切线方程为当TYT00,FT时，的近似值为，并记为，这就是得到时计算1101,YFAT11的近似公式T100,FYTA当时，的近似值为，并记为于是就得到当时计2T2XY21,YT2Y2T算的近似公式2112,FTYT重复上面方法，一般可得当的计算的近似公式NTN11,NYFTT如果，则上面公式就是式724将连起来，就1,NHTNNP,10得到一条折线，所以EULER方法又称为折线法POLYGONMETHOD,见图721由公式724知，已知便可算出；已知，便可算出；如此继续下去，这种0Y1Y12Y只用前一步的值便可计算出的递推公式称为单步法ONESTEPMETHODKK图7217212梯形格式用梯形求积公式计算式722右端积分项，得,11,D,2NTNNHFYTFTYFTY代入式722右端，得72611,NNNNTTFTFT在式726中，将用近似替代，所得结果记为，得YY5727110,0,1,2,NNNHYFTYFTYN并称727为求解初值问题721的梯形方法TRAPEZOIDALMETHOD,或梯形格式TRAPEZOIDALSCHEME是差分方程11,2NNNHYFTYFTY722改进的EULER方法因为梯形求积公式的精度比矩形求积公式的精度高，所以求解初值问题721的梯形方法比EULER方法精度高；但梯形方法是隐式方法，不便于计算。为了克服这个困难，先用EULER格式724求得一个初步的近似值，称之为预测值；然后用公式726作一次1NY迭代得，即将校正一次，得如下方法1NY1N定义721对，称预测校正系统PREDICTORCORRECTORSYSTEM0,N72811,2NNNYHFTYFTY测预校正为改进的EULER格式MODIFIEDEULERSCHEME,称这种方法为求解初值问题721的改进的EULER方法MODIFIEDEULERMETHOD注改进的EULER格式728还可以表示为下列平均化形式729112,PNNCPNPCYHFTY例721取步长，分别用EULER方法及改进的EULER方法求解初值问题0H2,1,1TYTE解这个初值问题的准确解为根据题设知2TT2,TFTYEEULER方法的计算格式为210NTNNTYYE改进的EULER方法的计算格式为1212102,NNNTNTTNNTYYEYE或6121120,NNTPNTCPNPCNTYYE初始值均为，将计算结果列于表72100Y表721EULER方法改进的EULER方法准确值NNTNYNYNYT01000111027182820342377803459199212068475560858314508666425313127697831592749616072151414209354772598298226203596515318744513936444139676663616462081785678907157209615717646639647909209279638735818880911971072446711079362479191174799651423744181432308151020153982357185788825186830971从表721可以看出,EULER方法的计算结果与准确值差距较大,最多只有1位有效数字；而改进的EULER方法却至少有2位有效数字,这表明改进的EULER方法的精度比EULER方法高。723局部截断误差与阶无论是EULER方法，还是改进的EULER方法，得到都是近似解，因此首先要考虑的是这些数值解的误差为使讨论简单起见，不妨假定在计算时用到前面一步的值是准确值1NY，即NYTNYT定义722在的前提下，称误差为局部截断误差LOCALNYT1NTTRUNCATIONERROR若一种数值方法的局部截断误差，721011PNTOH则称这种方法的精度为阶的ACCURACYOFPTHORDERP定理722EULER方法是具有一阶精度的数值方法即其局部截断误差为2OH证设，由721知NYT,NNYTFTY将在处作TAYLOR展开1T721121,NNNNNNHTHTTTT由EULER方法724得72121,NNNNNNYFTYTFTYTHYT7将式7211与式7212相减得7213221NHYTYOH于是,即由定义821知EULER方法是具有一阶精度的数值方法证毕12P仿定理722的证明，并利用二元TAYLOR展开式得定理723改进的EULER方法是具有2阶精度的数值方法即其局部截断误差为3OH注从实际情况来看，在计算时用到前面一步的值不可能是准确值，而应1NYNYT该是的前面每一步都有误差，这样得到的与准确值的误差称为整体截断误1NY1N1NYT差作为误差估计，最根本的是估计整体截断误差而估计整体截断误差并不容易，好在整体截断误差与局部截断误差之间通常有一定的联系见定理844整体截断误差局部截断误差72141OH即单步法的整体截断误差比局部截断误差低一阶因此要构造高精度的计算方法，必须设法提高该方法的局部截断误差中的一般地，一个方法的整体截断误差阶越1P高，该方法的精度也越高由式7214知EULER方法的整体截断误差为,改进的EULER方法的整体截断误H差为2OH73RUNGEKUTTA方法731RUNGEKUTTA方法的基本思想对初值问题731,YTFATBA由72节知EULER方法每一步只计算一次的值，其整体截断误差为，它的计,FTYOH算格式可以写成73211,2NNYHKKFT而改进的EULER方法，每迭代一步要计算两个函数值，其整体截断误差为，它的计2H算格式可改写成73311221,NNHYKKFTY上述2种情况启发我们，能否通过在每步迭代中多计算几次函数的值，然后对这些,FTY8值作线性组合，使相应的数值方法的整体截断误差的阶数升高，从而提高该方法的精度,而库塔MARTINWILHELMKUTTA1867年11月3日1944年12月25日是德国数学家。这正是龙格库塔RUNGEKUTTA方法的基本思想由72节可知因为整体截断误差与局部截断误差的关系，所以要构造高精度的计算方法，只要设法提高该方法的局部截断误差即可；因此下面利用局部截断误差讨论RUNGEKUTTA方法。732二阶RUNGEKUTTA方法对初值问题731，以计算两次函数值为例，设一般计算格式为73411221,NNYHCKKFTQY其中。适当选择参数的值，使得在的假设下，局部截断误12C1,CNTY差的阶尽可能高为此,将在点作TAYLOR1NYT211KCHN,NTY展开，因为1,NNKFTYT的展开式为2K121122,1,NNNNNNTYNNFTQHYFFYQHKTQHKFTYTFTYFFFTOQHO考虑到，所以12C1212NCKTHCK7353,NNYTTQY再将在作TAYLOR展开1NYT736231NNNNTTTTO若要使局部截断误差，则通过比较式835和836知，参数3YOH必须满足12,CQ,7371221,CCQ这是3个未知量，2个方程式的方程组，因此解不惟一定义731满足条件737的数值方法734称为二阶RUNGEKUTTA方法RUNGEKUTTAMETHODOFORDERTWO简记为二阶RK方法。注二阶RUNGEKUTTA方法的局部截断误差都是3OH91若取,则式734化为121CQ73812121,NNHYKKFTY这就是改进的EULER方法7282若取,则式834化为120,CQ7391221,NNHYKKFTY并称之为中点方法MIDPOINTMETHOD3若取,则式734化为123,4CQ73101122143,NNHYKKFTY并称之为HEUN方法HEUNMETHOD例731取步长，分别用改进的EULER方法、中点方法及HEUN方法求解初值01H问题2,1,10TYTE解根据题设知，0H2,TFT改进的EULER方法的计算格式为1122012210,NNNTTNTYKKEYTE中点方法的计算格式为12012210,NNNTTNTYKKEYTEHEUN方法的计算格式为10112202321034,NNNTTNTYKKEYTE上述计算格式的初始值均为，将计算结果列于表731表731改进的EULER方法中点方法HEUN方法准确值NTNYNYNYNYT010000111034237780340944403413921034591992120858314508549096085597430866642531315927496158670781588600116072151414259829822588807725917859262035965153936444139225250392690203967666361656789071565938695665538157209615717790920927882694778910676796387358181072446711068931241070043671079362479191423744181419171261420621221432308151020185788825185203146185389208186830971从表731可以看出,三种RK方法均至少有2位有效数字。733三阶及四阶RUNGEKUTTA方法为了提高方法的精度,考虑每步计算三次函数值根据两次计算函数值的做法，,YXF很自然地取的形式为1NY731111232132,NNYHCKCKKFTPYFTQDH其中。适当选择参数,，使截断误差的阶123C1C2311NYT尽可能高由于在推导公式时只考虑局部截断误差，故设类似前面二阶RKNYT公式的推导方法，将在处作TAYLOR展开，然后再将在处作TAYLOR1NY,NT1NT展开，只要两个展开式的前四项相同，便有，而要两个展开式的前四410NYTH项相同，参数必须满足7312123122332,16CQDPPCPCD，这是7个未知数5个方程的方程组，解不是惟一的，可以得到很多公式定义732满足条件7312的数值方法7311称为三阶RUNGEKUTTA方法RUNGEKUTTAMETHODOFORDERTWO简记为三阶RK方法。注三阶RUNGEKUTTA方法的局部截断误差都是4HO11若在7312中取1231241,662CCPQD则得到一种常用的三阶RUNGEKUTTA公式RUNGERKUTTAFORMULAOFORDERTHREE7313112321324,NNNHYKKFTYFTHKH类似地，如果每步计算4次函数的值，那么可以导出局部截断误差为的5HO四阶RUNGEKUTTA公式详细推导过程这里略去，本节只给出一种最常用的四阶RK公式7314112342132436,NNNNHHYKKKFTYFTKK公式7314称为四阶经典RUNGEKUTTA方法CLASSICALRUNGEKUTTAMETHODOFORDERFOUR三阶、四阶RUNGEKUTTA公式都是显式单步公式例732取步长，分别用三阶RUNGEKUTTA方法7313和四阶经典RUNGE01HKUTTA方法7314求解初值问题2,1,10TYTE解根据题设知，0H2,TFT三阶RUNGEKUTTA方法8313的计算格式为01211232210130614,NNNNTTNTNTTYKKEYTEKK四阶经典RK方法7314的计算格式为120120121123422132430106,NNNNTTNTNTNTTYKKKEYTEKKT初始值均为，将计算结果列于表73200Y表732三阶RK方法四阶经典RK方法准确值NTNYNYNYT01000111034561110345910303459199212085944908666217086664253131606044716071813160721514142618628026203113262035965153965280539676019396766636165717823457208793572096157177959879779637718796387358181078866611079350181079362479191431704291432293571432308151020186758568186829266186830971从表732可以看出,三阶RK方法至少有3位有效数字；四阶经典RK方法至少有5位有效数字。例733试用EULER方法、改进的EULER方法及四阶经典RUNGEKUTTA方法在不同步长下计算初值问题2,1,10TYTE在12,14,18,20处的近似值，并比较它们的数值结果。解对上述三种方法，每执行一步所需计算的次数分别为1,2,42,TFTYE为了公平起见，上述三种方法的步长之比应为因此，在用EULER方法、改进的421EULER方法及四阶经典RK方法计算12,14,18,20处的近似值时，它们的步长应分别取为005,01,02，以使三种方法的计算量大致相等EULER方法的计算格式为2105NTNTYYE改进的ELUER方法的计算格式为131122012210,NNNTTNTYKKEYTE四阶经典RK方法的计算格式为0202112342213243006,NNNNTTNTNTNTTYKKKEYTEKKT上述计算格式的初始值均为，将计算结果列于表733010Y表733EULER方法步长H005改进的EULER方法步长H01四阶经典RK方法步长H02准确解NTNYNYNYNYT120769696008583145086637910866642514234022362598298226197405262035961651373786567890715719895357209615189743489410724467110792017610793624720169490133185788825186808524186830971从表733可以看出，在计算量大致相等的情况下，EULER方法计算的结果只有1位有效数字，改进的EULER方法计算的结果有2位有效数字，而四阶经典RK方法计算的结果却有4位有效数字，这与理论分析是一致的例721和例733的计算结果说明，在解决实际问题时，选择恰当的算法是非常必要的。需要指出的是RUNGEKUTTA方法是基于TAYLOR展开法，因而要求方程的解具有足够的光滑性如果解的光滑性差，使用四阶RUNGEKUTTA方法求得数值解的精度，可能不如改进的EULER方法精度高因此，在实际计算时，要根据具体问题的特性，选择合适的算法。734变步长的RUNGEKUTTA方法上面导出的RUNGEKUTTA方法都是定步长的，单从每一步来看，步长H越小，局部截断误差也越少但随着步长的减小，在一定范围内要进行的步数就会增加，而步数增加不仅增加计算量，还有可能导致舍入误差的积累过大由于RUNGEKUTTA方法是单步法，每一步计算步长都是独立的，所以步长的选择具有较大的灵活性因此根据实际问题的具体情况合理选择每一步的步长是非常有意义的正如第7章建立变步长的SIMPSON求积公式那样，我们来建立变步长的RUNGEKUTAA公式以四阶经典RUNGEKUTTA公式为例进行说明从初始点出发，先选一个步长，利用0TH14公式7314求出的近似值记为,由于公式的局部截断误差为,所以有1HY5HO,731551KTCH其中为常数然后将步长进行折半,即取步长为,由基点出发,到算一C20T012T步,由到再算一步,将求得的近似值记为,因此有012TT1Y7316521HYTC由式7315和7316得216HTY一般地,从出发,按上述做法可得到NT21,HNTY或写成731722115HHNNTY式7317是事后估计式,记731821HNY对给定的步长H，若是预先指定的精度，说明步长太大，应折半进行计算，直H至为止，这时取作为近似值若，则将步长加倍，直到为止，21HNY这时再将步长折半一次，就把这次所得结果作为近似值变步长的RUNGEKUTTA方法的计算步骤如下设误差上限为，误差最小下限为，步长最大值为，从出发进行计算,步长上下HNT为H1用步长和RUNGEKUTTA公式计算，用步长计算两步得，并计算；1HNY221HNY2若，说明步长过大，应将折半，返回步骤1重新计算；上3若，说明步长过小，在下一步将放大，但不超过下HH这种通过加倍和减半处理步长的RUNGEKUTTA方法就称为变步长RUNGEKUTTA方法RUNGEKUTTAVARINGSTEPSIZEMETHOD从表面上看，为了选择步长，每一步的计算量增加了，但从总体考虑还是合算的74单步法的相容性、收敛性与稳定性常微分方程的初值问题741,YTFATBA15的数值解法的基本思想是通过某种离散化手段将微分方程转化为差分方程，很自然地，我们必须回答以下三个问题1、这种转化是否合理即，相容性问题。2、由数值解法得到的数值解是否收敛到值问题741的解析解即，收敛性问题。3、数值解法是否稳定即，稳定性问题。本节仅介绍单步法的相容性、收敛性和稳定性，而多步法的的相容性、收敛性和稳定性的证明比较复杂，这里就不介绍了，有兴趣的读者请参阅文献1下面均设，其中为等距节点，011NATTB0,1,NTHN步长HB751相容性定义741若解初值问题741数值方法的计算格式可以写成74210,0,1,NNYHTY则称该方法为显式单步法EXPLICITONESTEPMETHOD其中称为增量函数,TYHINCREMENTFUNCTION注增量函数是依赖于微分方程中的的一个泛函例如EULER方法,YFT,FTY的增量函数是；改进的EULER方法的增量函数是,TYHFT1,2FTFTHFT定义742设求解处初值问题741的单步法为742，若关于连续，且对YH任意固定的，都有,TAB74300LIMLI,HHYTTT则称数值方法742与初值问题741是相容；或数值方法742与初值问题741是相容的CONSISTENT。因为，所以有0LI,HYTTYFT0LI,0HTYT，744式744称为相容性条件CONSISTENTCONDITION。注1若关于连续，则式743与式744式等价的；因此，若数值方法,TYH742的增量函数满足相容性条件744，则称数值方法742与初值问题741相容。注2相容性保证了，当时，数值方法的差分方程收敛于原初值问题的微分方程，0说明数值解法通过某种离散化手段将微分方程转化为差分方程求解原初值问题是合理的。注3容易验证，EULER方法和改进的EULER方法都与初值问题841相容。752收敛性定义743对单步法742，若对任意固定的，都有，则称方,NTAB0LIMNHYT16法742是收敛的CONVERGENT定义744若在计算除外时,用到的不是，而是近似值；即每一1NY0NYTNY步计算除局部截断误差外，还考虑前一步不准确而引起的误差，则称这种误差为整体截断误差GLOBALTRUNCATIONERROR整体截断误差包含各步的局部误差，也包括了以前各步的局部误差在逐步计算中的积累因此作为误差估计，最根本的是估计整体截断误差但估计整体截断误差并不容易，好在整体截断误差与局部截断误差之间通常有一定的联系见定理744,因此要构造高精度的计算方法，只要设法提高该方法的局部截断误差中的阶数即可P定理741若A单步法742具有阶精度；即它的局部截断误差为P，74511PNYTOH其中1,NNNYTHTB增量函数在区域内关于满足LIPSCHITZ,DTATBYY条件；即存在，对于，有0L12,YR746212,THTYLC初始值是准确的,00E则1当时，方法742是收敛的P2方法742的整体截断误差为，74711PNNEYTOH其中1,NNYHT证由条件A知，存在常数，使得C,11PNYT且由条件B得1,NNNNNNNNYTHTYHTYTTLTY因此111,NNNPYTTYHLYCH即有递推关系74811PNNEE利用递推关系748反复递推，得17121011NNNNPNPEHLEHLHLCHL注意到，当是某个固定的常数时，0NTH，1NNHLHE于是有7491LPNC因此当时，由式749知，方法742是收敛的，且747成立证毕0P注1在定理741的条件下，单步法742的整体截断误差比局部截断误差低一阶；即整体截断误差局部截断误差1OH注2要判断某个单步法是否收敛，只要判断增量函数是否关于满足LIPSCHITZ条件Y推论1EULER方法是收敛的证明因为EULER方法满足；且由71节的假定知EULER方法的增量函数0P关于满足LIPSCHITZ条件，故由定理744知EULER方法是收敛的证,TYHFTY毕推论2改进的EULER方法也是收敛的证明因为改进的EULER方法的增量函数是，1,2TYHFTYFTHYFT由71节的假定知存在，对于，有；0L12R1212YLY于是1211221212221,TYHTFFYFTHYFTFTHYFTLYYYLHY2212,HH其中是某个固定的常数,所以改进的EULER方法的增量函数关于满足HYLIPSCHITZ条件又因为改进的EULER方法满足；所以由定理741知改进的EULER方法也是20P18收敛的证毕注相容性保证了当步长时，数值方法的差分方程收敛到初值问题的微分方程；H收敛性保证了当步长时，数值方法的数值解收敛到初值问题的解析解。下面的定理阐述了单步法的收敛性与相容性的关系定理742若单步法742的增量函数在区域,TYH上关于满足LIPSCHITZ条件；即存在，对于,DTYATBYY0L，有，则当精度时，方法742收敛12R1212,HTLP的充分必要条件是相容性条件744成立。证明参见文献13。743稳定性在理论上，讨论方法的收敛性时，总是假定初始值是准确的，数值方法本身的计算也是准确的，而实际情况并非如此事实上，初始数据可能有误差，计算过程中的数字的舍入也会产生误差这种误差的扰动在计算过程中是否会增长很快，以至影响结果，这就是数值方法的稳定性问题定义745若一个数值方法在节点处的值有扰动，在以后的各点处NTNYMTN的值产生的偏差均不超过，则称这种方法是稳定的STABLEMY在选择数值方法时，我们总是希望选择稳定的数值方法定理743对初值问题741，若增量函数是连续的，在集合,TYH,0DTYATBH上关于满足LIPSCHITZ条件；即存在，对于，有YL12,R，741012,THTY其中是常数,则方法742是稳定的0H当考察某数值方法的稳定性时，为简明起见，通常将此方法用于简单的试验方程以检验该方法的稳定性试验TESTEQUATION方程为7411YTT其中是常数，可以是实数，也可以是复数定义746将单步法742用于试验方程7411，设得到的迭代格式为，其中步长是固定的若增量函数满足，则称此方法是1NNYGHYHGH1绝对稳定的ABSOLUTESTABLE。记，则在平面上，使的变量围成的区域成为绝对稳定域REGIONOFABSOLUTESTABLE，它与实轴的交集称为绝对稳定区间INTERVALOFABSOLUTESTABLE下面分别将EULER方法，后退的EULER方法，改进的EULER方法和梯形方法分别应用于19试验方程7411，并记，确定上述方法的绝对稳定区域H1将EULER方法应用于试验方程749，得1,NNYFTY,74121NHY其中若有误差而变为，相应地有误差而变为，则GHNN11NY74131NY将式7412与式7413相减,并记，得NE11,NNNEHYHE因此74141NE由此可见，若要求误差不增长，则比值必须满足1NH，即7415因为可以是复数，所以在的复平面上，式7415表示EULER方法的绝对稳XIY定区域是半径为1，圆心为的圆域内部，不包括圆周如图741,0图741绝对稳定区间是此时显然要求，即RE2,0HRE0这说明当步长满足其中时，EULER方法02HH是稳定的因为稳定性对步长有限制，因此称这样的稳定为条件稳定CONDITIONALSTABLE2将后退的EULER方法应用于试验方程7411，并记11,NNYFTY,得即,741611NNH1NNH其中若有误差而变为，相应地有误差而变为，则1GHYY1NY2074171NNYYH将式7416与式7417相减,并记，得E,111NNNNEYYEHH由此得74181NE由此可见，若要求误差不增长，则比值必须满足1NH或,7419H1即74201因此在复平面上，式7420表示后退的EULER方法的绝对稳定区域是半径为1，圆心为的圆域的外部，不包括圆周如图7421,0图742当时，在复平面上，式7420表示后退的EULER方法的绝对稳定区域是左RE0半平面，绝对稳定区间是，即这说明后退的EULER方ER,0HH法对任意步长都是稳定的，此时称后退的EULER方法是无条件稳定的UNCONDITIONALHSTABLE3将改进的EULER方法一种二阶RUNGEKUTTA方法1,2NNNNHYFTYFTHYFT应用于试验方程7411，得,742121NN其中仿上述过程，得2GHH2121NHE若要求误差不增长，则比值必须满足,即21N211H或74222于是改进的EULER方法的绝对稳定区域由确定，由得22RE改进的EULER方法的绝对稳定区间此时显然要求，即RE,0H0这说明当步长满足其中时，改进的02REHEEULER方法是稳定的，且是条件稳定的4将梯形方法11,2NNNHYFTYFTHY应用于试验方程7411，得,742312NN其中仿上述过程，得12HG12NEH若要求误差不增长，则比值必须满足1NE即74242H12当时，在复平面上，式7424表示隐式梯形方法的绝对稳定区域是左半平面，RE0绝对稳定区间是，即这说明隐式梯形方法对任意步长ER,0H都是稳定的，即隐式梯形方法也是无条件稳定的。H从前面的讨论可以看出，隐式方法的稳定区域如后退EULER方法和隐式梯形方法分别比显式方法如EULER方法和改进的EULER方法的稳定区域大，而且在条件下后退RE0EULER方法和隐式梯形方法都是无条件稳定的，而EULER方法和改进的EULER方法却是条件稳定的。一般地，显式RUNGEKUTTA方法皆具有有限的绝对稳定区域，即对步长都有限制，H因此它们都是条件稳定的。而隐式方法在某些情况下是无条件稳定的。2275线性多步法73节介绍的RUNGEKUTTA法是单步法，即计算时，只用到前面一步的信息1NYNY为了提高精度，三阶及三阶以上的RUNGEKUTTA法用的都是内部的点处的函数值1,NT信息，这样做的不足是计算量增加；而且这些新增加的函数值信息不保留到以后步骤的运算，在算法上似乎有点浪费这就启发我们如果能充分利用已知信息,来1NY2计算，那么不但有可能提高精度，而且大大减少了计算量；这就是构造所谓线性多步1NY法的基本思想又因为误差随着的增大而趋于增大，所以当求时，使用NYT这些以前更精确的数据信息的方法也是合理的定义751求解初值问题751,TFATBYA的K步线性多步法LINEARKSTEPMULTISTEPMETHOD的一般公式是1011NNNKY75201111,NNKNKHFTYFTYFTYFTY其中和是与无关的常数,，IK,IHBAN,2,NTAN当时，式752中不含，称方法752为显式的EXPLICIT或OPEN；101,NFTY当时，式752中含，称方法752为隐式的IMPLICIT或CLOSED11,NFT为清晰简明起见，本节主要介绍最常用的线性多步法四阶ADAMS方法；本节最后还简介了其他一些较常用的线性多步法。751ADAMS方法定义752对,称3,41NN753112301235937,9,NNNNNHYFTYFTYFTYFTY为四阶显式ADAMS方法EXPLICITADAMSMETHODOFORDERFOUR对,称,3N75411120299,5,4,NNNNNHYFTYFTYFTYFTY为四阶隐式ADAMS方法IMPLICITADAMSMETHODOFORDERFOUR推导式753和式754将区间等分,ABN23,011NATTB其中为等距节点，步长0,1NTHNHA阿达姆斯JOHNCOUCHADAMS1819年6月5日1892年1月21日英国数学家和天文学家他最著名的成就是仅用数学理论和方法准确地预测了海王星的存在和位置。在每个小区间上，对两边积分，得1,0,NT,YFT7551,DNTNYTFS对运用插值型求积公式，便可得到相应的线性多步法公式1,DNTFSY将式755右边积分的被积函数用插值多项式来近似代替这里，我们取三次插值多项式插值节点除,外，通常还要在内再取两点作为插值节点但区间1NT1,NT内的点处的函数值是未知的,而却是已知的，因此可取插值节1NT,32Y点,则三次插值多项式为23NT123233111112223,NNNNNNNNNTTTTTPFTFTYTTFT23331,NYTTTFT12323112123331,6,26NNNNNNNTTTFYTTTFYTHHT插值余项为7565131231,4NNNNYRTTTTTT其中于是有1451D,TFY7573,FYTPTR将式755右端的积分中的用代替，令，然后积分，SS01NSTUH并用的近似值替代，得IYTIY,21,ITN1110023123,D23D6,3,65,9,7,9,4NNNNNNNNHFTUUFTYUHYHFTFTYFTYFTY即得式753若插值节点取为作三次插值多项式211,NNTT2412311121,NNNNNNTTTPTFTYTTT1112222,NNNNFTYTTT11112332112,61,6NNNNNNNTTTFYTTTTFTYHH其中,ITI插值余项为,75852311221,4NNNNNYRTTTTTTT其中于是有252D,TFY7593,FYTPTR将式755右端的积分中的用代替，令，然后积分，再分别将SSNSTUH,用近似值,表示，得2NYT1NTNYT12N1Y100112112D2D6,269,9,5,4NNNNNNHFUUFTYUHTYHFFTYFTYFTY即得式754定理751若，则四阶显式ADAMS方法753的局部截断4,C,YTFTAB误差为751055111312,70NNNNETHYOHTT四阶隐式ADAMS方法754的局部截断误差为751155112219,NNNNYTTT证明由余项表示式756得，四阶显式ADAMS方法的局部截断误差为,75121511123D24NTNNNNNETYSTSTTS其中3T25因为，所以在上连续,亦在上连续由介5C,YTAB5YTAB3,NTAB值定理知在上存在最大值记为和最小值记为，于是有3NTMM记51MM,152123DNTNNNTYSTSTTS则有111231123DDNNTTNNNNNNNSTSTTTTTST令，上式化为0UH15510,NMUMHUU即或551223030NHT1532NMT由介值定理知存在，使得NTT,155111235D2NTNNNYYSTSTTS即1551231D0NTNNNHSTSTTSY代入式7512，得751351131,70NNNEHYTT注这里不能用积分第二中值定理证明，因为未必是连续的5由余项表示式758得，四阶隐式ADAMS方法的局部截断误差为15112112D4NTNNNNNEYTYSTTSTS由介值定理，仿式7513的证明，得512219,70NNNEHTT证毕742初始值的计算四阶显式ADAM方法753计算时需要知道,而微分方程1NYNY12N3Y751只提供了一个初始值,还有3个初值需要用其他方01法来计算四阶隐式ADAM方法754计算时需要知道,而微分方程751只提1NYNY12N供了一个初始值，还有2个初值需要用其他方法来计算0Y2,因此，线性多步法不是能自行启动的方法，它需要其他方法提供初始值26因为微分方程的初始值在定解时起着重要作用，所以补充的初始值对精度要求较高对于四阶ADAMS公式（显式的或隐式的），其局部截断误差为，因此用于确定补充5HO初始值的其他方法的局部截断误差至少应不低于，否则由于初始值不准确，会影响5到最后的结果对四阶ADAMS公式，最常用的补充初始值的方法是用四阶RUNGEKUTTA方法，计算四阶ADAMS方法中所缺少的初值还可以用TAYLOR展开法，如将在处XY0T作TAYLOR展开20001,2YTYTYT取它的前若干项，使余项满足精度要求753ADMAS预测校正方法四阶显式ADAMS方法，在补充缺少的初值后就可计算了每步只算一个函数值，局部截断误差可达到，这是单步法难以达到的而四阶隐式ADAMS方法，每步需要解一5HO个非线性方程与四阶显式ADAMS方法相比，在计算上增加了许多困难但是同阶的隐式ADAMS方法的精度及稳定性比同阶的显式ADAMS方法好得多在实际应用中，为了保留隐式ADAMS方法的优点，一般将显式公式与隐式公式联合使用，用显式公式提供预测值，用隐式公式加以校正这样，既不用求解非线性方程，又能达到较高的精度这种方法称为预测校正方法定义754对，称3,41NN0112311259,37,9,2,5,0,124NNNNKKNHYFTYFTYFTYFTYK7514为四阶ADAMS预测校正迭代系统FORTHORDERADAMSPREDICTORCORRECTORITERATIVESYSTEM式7514表明四阶ADAMS预测校正迭代系统就是先用四阶显式ADAMS公式给出较好的近似值，再用四阶隐式ADAMS公式进行迭代可以证明，在一定条件下,1KYNK于是，可在迭代几次后，取作为的近似值1N式7514用四阶隐式ADAMS公式进行迭代，计算量比较大在实际应用中，一般希望用隐式公式迭代时，迭代的次数不要太多，如果能不进行迭代而能达到要求最好为此，仿照改进的EULER方法，我们