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

《常微分方程初值问题数值解法》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法（162页珍藏版）》请在装配图网上搜索。

1、会计学1常微分方程初值问题数值解法常微分方程初值问题数值解法2 定理定理1 1 设 在区域 上连续，关于 满足利普希茨条件，则对任意 ，常微分方程(1.1),(1.2)式当 时存在唯一的连续可微解 .)(xyfR,),(ybxayxDyR,00ybax,bax 关于解对扰动的敏感性，有以下结论. 定理定理2 2 设 在区域 （如定理1所定义）上连续，且关于 满足利普希茨条件，设初值问题fDysxyyxfy)(),(0的解为 ，则),(sxy.),(),(21210ssesxysxyxxL第1页/共161页3 这个定理表明解对初值的敏感性，它与右端函数 有关，当 的Lips.常数 比较小时，解对

2、初值和右端函数相对不敏感，可视为好条件.若 较大则可认为坏条件，即为病态问题. 如果右端函数可导，由中值定理有 若假定 在域 内有界，设 ，则LfDfL.,),(),(),(212121之间在yyyyyxfyxyyxyyyxf),(Lyyxf),(.),(),(2121yyLyxyyxy第2页/共161页4 求解常微分方程的解析方法只能用来求解一些特殊类型的方程，实际中归结出来的微分方程主要靠数值解法. 它表明 满足利普希茨条件，且 的大小反映了右端函数 关于 变化的快慢，刻画了初值问题(1.1)和(1.2)式是否是好条件.fLfy第3页/共161页5 所谓数值解法数值解法，就是寻求解 在一系

3、列离散节点 )(xy121nnxxxx上的近似值 . ,121nnyyyy 相邻两个节点的间距 称为步长步长. nnnxxh1 如不特别说明，总是假定 为定数，),2, 1(ihhi这时节点为 . ),2, 1 ,0(0inhxxn 初值问题(1.1),(1.2)的数值解法的基本特点是采取“步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 第4页/共161页6 描述这类算法，只要给出用已知信息 ,21nnnyyy计算 的递推公式. 1ny 一类是计算 时只用到前一点的值 ，称为单步法单步法. 1nyny 另一类是用到 前面 点的值 ，1nyk11,knnnyyy称为 步法步法. k

4、 其次，要研究公式的局部截断误差和阶，数值解 与精确解 的误差估计及收敛性，还有递推公式的计算稳定性等问题. ny)(nxy 首先对方程 离散化，建立求数值解的递推公式.),(yxfy 第5页/共161页7 9.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法 积分曲线上一点 的切线斜率等于函数 的值. ),(yx),(yxf 如果按函数 在 平面上建立一个方向场，那么，),(yxfxy积分曲线上每一点的切线方向均与方向场在该点的方向相一致. 在 平面上，微分方程 的解 称作它的积分曲线积分曲线. xy)(xyy ),(yxfy 第6页/共161页8 基于上述几何解释，从初始点 出发，),(000

5、yxP先依方向场在该点的方向推进到 上一点 ，然后再1xx 1P从 依方向场的方向推进到 上一点 ，循此前进1P2xx2P做出一条折线 （图9-1).210PPP图9-1第7页/共161页9 一般地，设已做出该折线的顶点 ，过 依方向场的方向再推进到 ，显然两个顶点 的坐标有关系 ),(nnnyxPnP),(111nnnyxP1,nnPP),(11nnnnnnyxfxxyy即 ).,(1nnnnyxhfyy（2.1）这就是著名的欧拉欧拉（Euler）方法方法. 欧拉法实际上是对常微分方程中的导数用均差近似，即)(,()()()(1nnnnnxyxfxyhxyxy直接得到的.第8页/共161页1

6、0若初值 已知，则依公式（2.1）可逐步算出0y),(0001yxhfyy),(1112yxhfyy第9页/共161页11 例例1 1 求解初值问题 . 1)0(),10(2yxyxyy（2.2） 解解 欧拉公式的具体形式为 ).2(1nnnnnyxyhyy取步长 ，计算结果见表9-1. 1.0h 初值问题（2.2）的解为 ，按这个解析式子算出的准确值 同近似值 一起列在表9-1中，两者相比较可以看出欧拉方法的精度很差. xy21)(nxyny7848.10.14351.15.07178.19.03582.14.06498.18.02774.13.05803.17.01918.12.05090

7、.16.01000.11.019nnnnyxyx计算结果对比表 第10页/共161页127321.17848.10.14142.14351.15.06733.17178.19.03416.13582.14.06125.16498.18.02649.12774.13.05492.15803.17.01832.11918.12.04832.15090.16.00954.11000.11.0)()(19nnnnnnxyyxxyyx计算结果对比表 还可以通过几何直观来考察欧拉方法的精度. 假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图9-2）.)(nnxyy nP

8、)(xyy 1nnPP)(xyy nP第11页/共161页13图9-2 从图形上看，这样定出的顶点 显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的. 1nP 为了分析计算公式的精度，通常可用泰勒展开将 )(1nxy在 处展开，则有 nx第12页/共161页14)()(1hxyxynn在 的前提下，)(nnxyy )()(,(),(nnnnnxyxyxfyxf)(2)(211nnnyhyxy 称为此方法的局部截断误差. ).,()(2)()(12 nnnnnnxxyhhxyxy于是可得欧拉法（2.1）的误差 （2.3）),(22nxyh .)(,()()(11nnxxnndttytfxyx

9、y（2.4）如果对方程 从 到 积分，得 nx1nx),(yxfy ).,(1nnnnyxhfyy（2.1）第13页/共161页15右端积分用左矩形公式 近似.)(,(nnxyxfh 再以 代替ny),(nxy 如果在（2.4）中右端积分用右矩形公式 )(,(11nnxyxfh),(111nnnnyxhfyy（2.5）称为后退的欧拉法后退的欧拉法. 代替 也得到（2.1），)(1nxy1ny局部截断误差也是（2.3）. 近似，则得另一个公式 它也可以通过利用均差近似导数 ，即)(,()()()(11111nnnnnnnxyxfxyxxxyxy)(1nxy直接得到. 第14页/共161页16 欧

10、拉公式是关于 的一个直接的计算公式，这类公式称作是显式的显式的；1ny 后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的隐式的. 1ny1ny第15页/共161页17 隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化. 设用欧拉公式 ),()0(1nnnnyxfhyy给出迭代初值 ，用它代入（2.5）式的右端，使之转化为显式，直接计算得 )0(1ny),()0(11)1(1nnnnyxfhyy然后再用 代入（2.5）式，又有 )1(1ny).,()1(11)2(1nnnnyxfhyy第16页/共161页18如此反复进行，得 )., 1 ,0(),()(11)

11、1(1kyxfhyyknnnkn（2.6）由于 对 满足利普希茨条件(1.3). ),(yxfy 由(2.6)减(2.5)得 ),(),(11)(111)1(1nnknnnknyxfyxfhyy.1)(1nknyyhL由此可知，只要 迭代法(2.6)就收敛到解 .1hL1ny 从积分公式可以看到后退欧拉方法的公式误差与欧拉法是相似的. 第17页/共161页19 若用梯形求积公式近似等式（2.4）右端的积分并分别用 代替 则可得到比欧拉法精度高的计算公式1,nnyy),(),(1nnxyxy),(),(2111nnnnnnyxfyxfhyy（2.7）称为梯形方法梯形方法. 梯形方法是隐式单步法，

12、可用迭代法求解. 1)(,(nnxxdttytf.)(,()()(11nnxxnndttytfxyxy（2.4）第18页/共161页20);,()0(1nnnnyxfhyy 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得 ),(),(2)(1111)1(11knnnnknnyxfyxfhyy),(),(2)(11)1(1knnnnnknyxfyxfhyy（2.8）).,2, 1 ,0(k 同后退的欧拉方法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为 ),(),(2111nnnnnnyxfyxfhyy（2.7）第19页/共161页21如果选取 充分小，使得 h, 12hL则

13、当 时有 ，k1)(1nknyy 这说明迭代过程(2.8)是收敛的. 于是有 ,2)(11)1(11knnknnyyhLyy式中 为 关于 的利普希茨常数. ),(yxfyL第20页/共161页22 梯形方法虽然提高了精度，但其算法复杂. 在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.),(yxf 为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法. 具体地，先用欧拉公式求得一个初步的近似值 , 1ny 而迭代又要反复进行若干次，计算量很大，而且往往难以预测. 称之为预测值预测值，第21页/共161页23 这样建立的预测-校正系统通常称为改进的欧拉

14、公式：改进的欧拉公式： 预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正值校正值. .1ny1ny预测),(1nnnnyxfhyy校正).,(),(2111nnnnnnyxfyxfhyy（2.9）也可以表为下列平均化形式 ),(nnnpyxfhyy),(1pnncyxfhyy).(211cpnyyy),(),(2111nnnnnnyxfyxfhyy（2.7）),(),(2)(11)1(1knnnnnknyxfyxfhyy（2.8）第22页/共161页24 例例2 2 用改进的欧拉方法求解初值问题（2.2）. 解解 这里 改进的欧拉公式为

15、),2(nnnnpyxyhyy. 1)0(),10(2yxyxyy（2.2）),2(),(yxyyxf),2(1pnpncyxyhyy).(211cpnyyy第23页/共161页25仍取 ，计算结果见表9-2. 1.0h 同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度. 7321.17379.10.14142.14164.15.06733.16782.19.03416.13434.14.06125.16153.18.02649.12662.13.05492.15525.17.01832.11841.12.04832.14860.16.00954.10959.11.0)()(nnnnnn

16、xyyxxyyx计算结果对比2表9第24页/共161页26 初值问题(1.1)，(1.2)的单步法可用一般形式表示为 ),(11hyyxhyynnnnn（2.10）其中多元函数 与 有关，),(yxf 当 含有 时，方法是隐式的，若不含 则为显式方法，1ny1ny),(1hyxhyynnnn（2.11） 称为增量函数，),(hyx).,(),(yxfhyx所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.)(),(00yxyyxfy（1.1）（1.2）).,(1nnnnyxhfyy（2.1）)(2)(211nnnyhyxy （2.3）),(22nxyh 第2

17、5页/共161页27 对一般显式单步法则可如下定义. 定义定义1 1 设 是初值问题(1.1)，(1.2)的准确解，称)(xyy ),(,()()(11hxyxhxyxyTnnnnn（2.12）为显式单步法（2.11）的局部截断误差局部截断误差. 之所以称为局部的，是假设在 前各步没有误差. 1nTnx当 时，计算一步，则有 )(nnxyy ),()()(111hyxhyxyyxynnnnnn.),(,()()(11nnnnnThxyxhxyxy.)(),(00yxyyxfy（1.1）（1.2）),(1hyxhyynnnn（2.11）第26页/共161页28在前一步精确的情况下用公式（2.11

18、）计算产生的公式误差. 根据定义，欧拉法的局部截断误差 )(,()()(11nnnnnxyxhfxyxyT即为（2.3）的结果. 这里 称为局部截断误差主项. )(22nxyh )(21hOTn局部截断误差可理解为用方法（2.11）计算一步的误差，即)()()(nnnxyhxyhxy),()(232hOxyhn 显然),(1hyxhyynnnn（2.11）)(2)(211nnnyhyxy （2.3）),(22nxyh 第27页/共161页29 定义定义2 2 设 是初值问题(1.1)，(1.2)的准确解，若存在最大整数 使显式单步法(2.11)的局部截断误差满足 )(xyp),(),()()(

19、11pnhOhyxhxyhxyT（2.13）则称方法（2.11）具有 阶精度阶精度. p 若将（2.13）展开式写成 ),()(,(211ppnnnhOhxyxT则 称为局部截断误差主项局部截断误差主项. 1)(,(pnnhxyx 以上定义对隐式单步法（2.10）也是适用的. .)(),(00yxyyxfy（1.1）（1.2）),(1hyxhyynnnn（2.11）),(11hyyxhyynnnnn（2.10）第28页/共161页30 对后退欧拉法（2.5）其局部截断误差为 )(,()()(1111nnnnnxyxhfxyxyT这里 ，是1阶方法，局部截断误差主项为 . 1p)(22nxyh

20、)()()()()(2)(232hOxyhxyhhOxyhxyhnnnn ).()(232hOxyhn ),(111nnnnyxhfyy（2.5）第29页/共161页31 对梯形法（2.7）有 )()(2)()(111nnnnnxyxyhxyxyT 所以梯形方法是二阶的，其局部误差主项为)(123nxyh )()(2)()()(2)(!3)(2)(4232hOxyhxyhxyxyhxyhxyhxyhnnnnnnn ).()(1243hOxyhn ),(),(2111nnnnnnyxfyxfhyy（2.7）第30页/共161页32第31页/共161页33 上节给出了显式单步法的表达式),(1hy

21、xhyynnnn其局部截断误差为),(),()()(11pnhOhyxhxyhxyT对欧拉法 ，即方法为 阶，)(21hOTn1p).,(,(),(21nnnnnnnnyxhfyhxfyxfhyy（3.1） 若用改进欧拉法，它可表为 第32页/共161页34此时增量函数 ).,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx（3.2）与欧拉法的 相比，增加了计算一个右函数 的值，可望 . ),(),(nnnnyxfhyxf2p 若要使得到的公式阶数 更大， 就必须包含更多的 值. pf,)(,()()(11nnxxnndxxyxfxyxy（3.3） 从方程 等价的积分形式（2

22、.4），即 ),(yxfy 第33页/共161页35若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，必然要增加求积节点. 为此可将（3.3）的右端用求积公式表示为1.)(,()(,(1nnxxriininihxyhxfchdxxyxf点数 越多，精度越高，r上式右端相当于增量函数 ，),(hyx为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为 ),(1hyxhyynnnn（3.4）其中 ,)(,()()(11nnxxnndxxyxfxyxy（3.3）第34页/共161页36,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK（3

23、.5）,2ri这里 均为常数. ijiic, (3.4)和(3.5)称为 级显式龙格龙格- -库塔库塔（Runge-Kutta）法法，r简称R-K方法. 当 时，就是欧拉法，),(),(, 1nnnnyxfhyxr此时方法的阶为 . 1p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，2r),(1hyxhyynnnn（3.4）),(1hyxhyynnnn).,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx).,(1nnnnyxhfyy第35页/共161页37 下面将证明阶 . 2p 要使公式(3.4),(3.5)具有更高的阶 ,就要增加点数 . pr 下面就 推导

24、R-K方法. 2r,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK),(1hyxhyynnnn（3.4）（3.5）,2ri第36页/共161页38 对 的R-K方法，计算公式如下 2r).,(),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6）这里 均为待定常数.21221,cc 希望适当选取这些系数，使公式阶数 尽量高. p 根据局部截断误差定义，（3.6）的局部截断误差为 )()(11nnnxyxyT（3.7）),(),(21221nnnnnhfyhxfcyxfch第37页/共161页39这里 . ),

25、(),(nnnnnyxffxyy 为得到 的阶 ，要将上式各项在 处做泰勒展开，由于 是二元函数，故要用到二元泰勒展开，各项展开式为1nTp),(nnyx),(yxf),(!32)(4321hOyhyhyhyxynnnnn 其中 ,),(nnnnfyxfy,),(),()(,(nnnynnxnnnfyxfyxfxyxfdxdy ),(),(2),(2nnyynnnxynnnxxnyxffyxffyxfy （3.8）);,(),(),(nnynnnxnnyyxffyxfyxf第38页/共161页40),(212nnnfhyhxf将以上结果代入局部截断误差公式则有 )(),(),(),(),(23

26、2122121hOhfyxfhyxffcfchfyxfyxfhfhTnnnynnxnnnnnynnxnn).(),(),(2212hOfhyxfhyxffnnnynnxn).(),()21(),()21()1(3221222221hOhfyxfchyxfchfccnnnynnxn要使公式（3.6）具有 阶，必须使 2p).,(),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6）第39页/共161页41,0121cc即 . 1,21,212121222cccc非线性方程组（3.9）的解是不唯一的. 令 ，则得 02 ac.21,12121aac这样得到的公

27、式称为二阶R-K方法，如取 ，则2/1a.1,2/121221cc这就是改进欧拉法（3.1）.,02122c（3.9）.021212c).,(,(),(21nnnnnnnnyxhfyhxfyxfhyy第40页/共161页42若取 ，则 得计算公式 . 1a, 2/1, 1, 021221cc,21hKyynn称为中点公式中点公式，相当于数值积分的中矩形公式. （3.10）也可表示为 ).,(2,2(1nnnnnnyxfhyhxfhyy),(1nnyxfK （3.10）).2,2(12KhyhxfKnn 的R-K公式(3.6)的局部误差不可能提高到 . 2r)(4hO).,(2,2(1nnnnn

28、nyxfhyhxhfyy).,(),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6）第41页/共161页43 把 多展开一项，从（3.8）的 看到展开式中 的项是不能通过选择参数消掉的.2Kny ffffyxy 实际上要使 的项为零，需增加3个方程，要确定4个参数 ，这是不可能的. 3h21221,cc 故 的显式R-K方法的阶只能是 ，而不能得到三阶公式. 2r2p第42页/共161页44 要得到三阶显式R-K方法，必须 . 3r).,(),(),(),(232131321212213322111hKhKyhxfKhKyhxfKyxfKKcKcKchy

29、ynnnnnnnn（3.11）其中 及 均为待定参数.321,ccc32313212,.)()(33221111KcKcKchxyxyTnnn此时（3.4）,（3.5）的公式表示为 公式（3.11）的局部截断误差为 ,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK),(1hyxhyynnnn（3.4）（3.5）,2ri第43页/共161页45只要将 按二元函数泰勒展开，使 ，32,KK)(41hOTn可得待定参数满足方程 .61,31,21, 13223233222332232313212321cccccccc（3.12）第44页/共161页46

30、这是8个未知数6个方程的方程组，解也不是唯一的. 所以这是一簇公式. 满足条件（3.12）的公式（3.11）统称为三阶R-K公式. 一个常见的公式为 ).2,(),2,2(),(),4(62131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn此公式称为库塔库塔三阶方法. 第45页/共161页47 继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个： 可以证明其截断误差为 . )(5hO 四阶龙格-库塔方法的每一步需要计算四次函数值 ，f).,(),2,2(),2,2(),(),22(6342312143211hKy

31、hxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13）第46页/共161页48 例例3 3 设取步长 ，从 到 用四阶龙格-库塔方法求解初值问题 . 1)0(),10(2yxyxyy 解解 这里 ，经典的四阶龙格-库塔公式为yxyyxf2),(),22(643211KKKKhyynn,21nnnyxyK2.0h0 x1x第47页/共161页49,222223KhyhxKhyKnnn,222112KhyhxKhyKnnn.)(2334hKyhxhKyKnnn 表9-3列出了计算结果，同时了出了相应的精确解.7321.17321.10.16125.16125.1

32、8.04832.14833.16.03416.13417.14.01832.11832.12.0)(39nnnxyyx计算结果表 比较例3和例2的计算结果，显然以龙格-库塔方法的精度为高. 第48页/共161页50 但由于这里放大了步长 ，所以表9-3和表9-2所耗费的计算量几乎相同. )2.0(h 龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质. 反之，如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法. 四阶龙格-库塔方法每一步要4次计算函数 ,其计算量f比每一步只要2次计算函数 的二阶龙格-库塔方法-f改进的欧拉方法大

33、一倍.第49页/共161页51 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题. 在选择步长时，需要考虑两个问题： 1 怎样衡量和检验计算结果的精度？ 第50页/共161页52 2 如何依据所获得的精度处理步长？ 考察经典的四阶龙格-库塔公式 ).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13）

34、从节点 出发，先以 为步长求出一个近似值 ，nxh)(1hny第51页/共161页53,)(5)(11chyxyhnn（3.14）然后将步长折半，即取 为步长从 跨两步到 ，再求得一个近似值 ，每跨一步的截断误差是 ，因此有2hnx1nx)2(1hny52 hc,22)(5)2(11hcyxyhnn（3.15）比较（3.14）式和（3.15）式我们看到，步长折半后，由于公式的局部截断误差为 ，故有 )(5hO误差大约减少到 ，161第52页/共161页54.161)()()(11)2(11hnnhnnyxyyxy由此易得下列事后估计式 .151)()(1)2(1)2(11hnhnhnnyyyx

35、y这样，可以通过检查步长，折半前后两次计算结果的偏差 )(1)2(1hnhnyy即有来判定所选的步长是否合适. 具体地说，将区分以下两种情况处理： 第53页/共161页55 1. 对于给定的精度 ，如果 ，反复将步长折半进行计算，直至 为止. 这时取最终得到的 作为结果； )2(1hny 2. 如果 ，反复将步长加倍，直到 为止， 这种通过加倍或折半处理步长的方法称为变步长方法变步长方法.这时再将步长折半一次，就得到所要的结果. 表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的. 第54页/共161页56 9.4.1 收敛性与相容性收敛性与相容性 数值解法的基本思想是通过某

36、种离散化手段将微分方程转化为差分方程，如单步法(2.11)，即 ),(1hyxhyynnnn（4.1）它在 处的解为 ，而初值问题（1.1）,（1.2）在 处的精确解为 ，nxnynx)(nxy 记 称为整体截断误差. nnnyxye)(.)(),(00yxyyxfy（1.1）（1.2）第55页/共161页57 收敛性就是讨论当 固定且 时nxx 00nxxhn0ne的问题. 定义定义3 3 若一种数值方法对于固定的 ,当 时有 ，其中 是(1.1),(1.2)的准确解，则称该方法是收敛收敛的. nhxxn00h)(nnxyy )(xy 显然数值方法收敛是指 .0)(nnnyxye 对单步法（

37、4.1）有下述收敛性定理： .)(),(00yxyyxfy（1.1）（1.2）),(1hyxhyynnnn（4.1）第56页/共161页58 定理定理3 3 假设单步法（4.1）具有 阶精度，且增量函数 关于 满足利普希茨条件p),(hyxy,),(),(yyLhyxhyx（4.2）又设初值 是准确的，即 ，则其整体截断误差整体截断误差 0y)(00 xyy ).()(pnnhOyxy（4.3） 证明证明 设以 表示取 用公式（4.1）求得的结果，即1ny)(nnxyy ),),(,()(1hxyxhxyynnnn（4.4）则 为局部截断误差，11)(nnyxy),(1hyxhyynnnn（4

38、.1）),(1hyxhyynnnn（4.1）第57页/共161页59由于所给方法具有 阶精度，按定义2，存在定数 ，使pC.)(111pnnChyxy又由式（4.4）与（4.1），得 nnnnyxyyy)(11利用假设条件（4.2），有 . ),(),(,(hyxhxyxhnnnn,)()1(11nnnnyxyhLyy从而有 111111)()(nnnnnnyxyyyyxy.)()1(1pnnChyxyhL,),(),(yyLhyxhyx（4.2）),(1hyxhyynnnn（4.1）),),(,()(1hxyxhxyynnnn（4.4）第58页/共161页60即对整体截断误差 成立下列递推关

39、系式 nnnyxye)(,)1(11pnnChehLe（4.5）反复递推，可得 .1)1()1(0npnnhLLChehLe（4.6）再注意到当 时 Tnhxxn0,)()1(TLnhLneehL最终得下列估计式 ).1(0TLpTLneLCheee（4.7）第59页/共161页61由此可以断定，如果初值是准确的，即 ，则（4.3）式成立. 00e 依据这一定理，判断单步法（4.1）的收敛性，归结为验证增量函数 能否满足利普希茨条件（4.2）. 对于欧拉方法，由于其增量函数 就是 ，故当 关于 满足利普希茨条件时它是收敛的. ),(yxf),(yxfy 再考察改进的欧拉方法，其增量函数).,(

40、,(),(21),(nnnnnnnnyxhfyhxfyxfhyx给出，这时有).()(pnnhOyxy（4.3）,),(),(yyLhyxhyx（4.2）),(1hyxhyynnnn（4.1）第60页/共161页62),(),(21),(),(yxfyxfhyxhyx假设 关于 满足利普希茨条件，记利普希茨常数为 ,),(yxfyL.)21(),(),(yyLhLhyxhyx设 为定数），00(hhh 上式表明 关于 的利普希茨常数y),(,(),(,(yxfhyhxfyxfhyhxf则由上式推得 , )21(0LhLL因此改进的欧拉方法也是收敛的. 第61页/共161页63 类似地，也可验证

41、其他龙格-库塔方法的收敛性. 定理3表明 时单步法收敛，并且当 是初值问题（1.1），（1.2）的解，（4.1）具有 阶精度时，1p)(xyp),(,()()(1hxyxhxyhxyTn有展开式 )0),(,()0),(,(2)()(2 hxyxxyxhhxyhxyx).()0),(,()(2hOxyxxyh所以 的充要条件是 ，1p0)0),(,()(xyxxy（4.1）),(1hyxhyynnnn.)(),(00yxyyxfy（1.1）（1.2）第62页/共161页64而 ，于是可给出如下定义： )(,()(xyxfxy 定义定义4 4 若单步法（4.1）的增量函数 满足 ),()0,(y

42、xfyx则称单步法（4.1）与初值问题（1.1），（1.2）相容相容. 相容性是指数值方法逼近微分方程(1.1)，即微分方程(1.1)离散化得到的数值方法当 时可得到0h).,()(yxfxy第63页/共161页65 定理定理4 4 阶方法(4.1)与初值问题(1.1)，(1.2)相容的充分必要条件是p.1p 由定理3可知单步法（4.1）收敛的充分必要条件是(4.1)是相容的. 以上讨论表明 阶方法（4.1）当 时与(1.1),(1.2)相容，反之相容方法至少是1阶的. p1p第64页/共161页66 定义定义5 5 若一种数值方法在节点值 上大小为 的扰动，于以后各节点值 上产生的偏差均不超

43、过 ，则称该方法是稳定稳定的. 以欧拉法为例考察计算稳定性. 例例4 4 考察初值问题 .1)0(,100yyy其准确解 是一个按指数曲线衰减得很快的函数，xxy100e)(ny)(nmym如图9-3所示. 第65页/共161页67.)1001(1nnyhy若取 ，则欧拉公式的具体形式为 025.0h,5.11nnyy计算结果列于表9-4的第2列. 可以看到，欧拉方法的解 ny（图9-3中用号标出）在准确值 的上下波动，计算过程明显地不稳定. )(nxy图9-3 用欧拉法解方程 得 yy1000625.5100.0375.3075.025.2050.05.1025.049表欧拉方法节点计算结果

44、对比第66页/共161页68再考察后退的欧拉方法，取 时计算公式为 025.0h.5.311nnyy计算结果列于表9-4的第3列（图9-3中标以号），这时计算过程是稳定的. 0067.00625.5100.00233.0375.3075.00816.025.2050.02857.05.1025.04-9后退欧拉方法欧拉方法节点计算结果对比表 但若取 则计算过程稳定. nnyyh5.0,005.01第67页/共161页69 这表明稳定性不但与方法有关，也与步长 的大小有关，当然也与方程中的 有关. h),(yxf 为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，, yy（4.

45、8）其中 为复数. 例如在 的邻域，可展开为 ),(yx),(yxfy 模型方程为 对一般方程可以通过局部线性化化为这种形式.,)(,()(,(),(yyyxfxxyxfyxfyx第68页/共161页70略去高阶项，再做变换即可得到 的形式. uu 对于 个方程的方程组，也可线性化为 ，mAyy 这里 为 的雅可比矩阵 . Amm)(jiyf 若 有 个特征值 ，则 还可能是复数，Amm,21i 为保证微分方程本身的稳定性，还应假定 . 0)Re( 先研究欧拉方法的稳定性. 模型方程 的欧拉公式为 yy所以，为了使模型方程结果能推广到方程组，方程中 应为复数. yy第69页/共161页71.)

46、1(1nnyhy（4.9）设在节点值 上有一扰动值 ，它的传播使节点值 产生大小为 的扰动值，nyn1ny1n 假设用 按欧拉公式得出 的计算过程不再有新的误差，nnnyy*11*1nnnyy则扰动值满足 .)1(1nnh可见扰动值满足原来的差分方程（4.9）. 如果差分方程的解是不增长的，即有 ,1nnyy则它就是稳定的. 第70页/共161页72即图9-4 显然，为要保证差分方程（4.9）的解不增长，只要选取 充分小，h.11h（4.10）在 的复平面上，这是以 为圆心，1为半径的单位圆内部（图9-4）. h)0, 1( 这个圆域称为欧拉法的绝对稳定域，一般情形可由下面定义. 定义定义6

47、6 单步法（4.1）用于解模型方程（4.8），若得到的解 ，满足 ，则称方法（4.1）是绝对稳定绝对稳定的. nnyhEy)(11)(hE.)1(1nnyhy（4.9）,1nnyy使),(1hyxhyynnnn（4.1）, yy（4.8）第71页/共161页73 在 的平面上，使 的变量围成的区域，称为绝对稳定域绝对稳定域，h1)(hE 对欧拉法，,1)(hhE.11h给出，绝对稳定区间为 .02h它与实轴的交称为绝对稳定区间绝对稳定区间. 其绝对稳定域由 在例5中 ，01002,100h02.00 h即 为绝对稳定区间. 当取 时005.0h 例4中取 故它是不稳定的，025.0h它是稳定的

48、. 第72页/共161页74,2)(121nnyhhy故 .2)(1)(2hhhE绝对稳定域由 得到.12)(12hh 绝对稳定区间为 ，即 . 02h/20 h 类似可得三阶及四阶的R-K方法的 分别为 )(hE,!3)(!2)(1)(32hhhhE 用二阶R-K方法解模型方程 可得到 yy第73页/共161页75.!4)(!3)(!2)(1)(432hhhhhE由 可得到相应的绝对稳定域. 1)(hE 当 为实数时则得绝对稳定区间.分别为 三阶显式R-K方法： 即051.2h./51.20 h 四阶显式R-K方法： 即078.2h./78.20 h 图9-5给出了R-K方法 到 的绝对稳定

49、域. 从以上讨论可知显式的R-K方法的绝对稳定域均为有限域，都对步长 有限制. hh 如果 不在所给的绝对稳定区间内，方法就不稳定. 1p4p图9-5第74页/共161页76 例例5 5 分别取 ，及 用经典的四阶R-K方法（3.13）计算. , 1)0(),10(20yxyy1.0h, 2.0h 解解本例,20分别为 及 .24h前者在绝对稳定区间内，后者则不在.经典的四阶R-K方法为).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13）第75页/共161页770 .31250

50、 .6250 .1250 .2598. 42 . 01017. 01015. 01014. 01012. 01093. 01 . 00 . 18 . 06 . 04 . 02 . 05-943211hhxn计算结果表 以上结果看到，如果步长 不满足绝对稳定条件，误差增长很快. h用四阶R-K方法计算其误差见下表： 这里, 1,000yx第76页/共161页78 对隐式单步法，可以同样讨论方法的绝对稳定性，,111nnyhy例如对后退欧拉法，用它解模型方程可得故 .11)(hhE由 可得绝对稳定域为 .hhE11)(11h 它是以 为圆心，1为半径的单位圆外部，故绝对稳定区间为 . )0, 1(

51、0h第77页/共161页79 当 时，则 ，即对任何步长均为稳定的. 0 h0,21211nnyhhy故 .2/12/1)(hhhE对 有 ，0)Re(12/12/1)(hhhE故绝对稳定域为 的左半平面，h 对隐式梯形法，解模型方程 得 yy第78页/共161页80绝对稳定区间为 ，即 时梯形法均是稳定的. 0h h0 定义定义7 7 如果数值方法的绝对稳定域包含了 那么称此方法是A-稳定的. 隐式欧拉法与梯形方法的绝对稳定域均为 在具体计算中步长 的选择只需考虑计算精度及迭代收敛性要求而不必考虑稳定性，具有这种特点的方法特别重要. ,0)Re(hh,0)Re(hhh 由定义知A-稳定方法对

52、步长 没有限制.h第79页/共161页81 在逐步推进的求解过程中，计算 之前事实上已经求出了一系列的近似值 ，1nynyyy,10 如果充分利用前面多步的信息来预测 ，则可以期望会获得较高的精度. 1ny 这就是构造所谓线性多步法的基本思想. 本节主要介绍基于泰勒展开的构造方法. 构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程 两端积分后利用插值求积公式得到.),(yxfy 第80页/共161页82 一般的线性多步法公式可表示为 ,010kiinikiiniknfhyy（5.1）其中 为 的近似，iny)(inxy,),(0ihxxyxffinininin 计算时

53、需先给出前面 个近似值 , k110,kyyy再由（5.1）逐次求出 . ,1kkyy 如果计算 时，除了使用 的值，还用到kny1kny)2,2, 1 ,0(kiyin的值，则称此方法为线性多步法线性多步法. ii,为常数.及 不全为零，00则称为线性线性 步法步法. k若第81页/共161页83这时 可直接由（5.1）算出；kny 如果 ，称（5.1）为显式显式 步法步法，0kk 如果 ，则（5.1）称为隐式隐式 步法步法，0kk求解时与梯形法（2.7）相同，要用迭代法方可算出 .kny,010kiinikiiniknfhyy（5.1）),(),(2111nnnnnnyxfyxfhyy（2

54、.7）第82页/共161页84 定义定义8 8 设 是初值问题(1.1)，(1.2)的准确解，线性多步法（5.1）在 上的局部截断误差为 )(xyknx);(hxyLTnkn（5.2）. )()()(1010kiinikiiniknxyhxyxy （5.1）中系数 及 可根据方法的局部截断误差及阶确定，其定义为： ii,010kiinikiiniknfhyy（5.1）.)(),(00yxyyxfy（1.1）（1.2）若 ，则称方法（5.1）是 阶的阶的，)(1pknhOTp1p则称方法（5.1）与方程（1.1）是相容的相容的.第83页/共161页85 由定义8，对 在 处做泰勒展开，由于 kn

55、Tnx,)(!3)()(!2)()()()(32 nnnnnxyihxyihxyihxyihxy.)(!2)()()()(2 nnnnxyihxyihxyihxy代入（5.2）得 ,)()()()()(2210 npppnnnknxyhcxyhcxyhcxycT（5.3）（5.2）. )()()(1010kiinikiiniknknxyhxyxyT第84页/共161页86其中 ),(11100kc),()1(2101211kkkkc2)!1(1)1(2(!11211121kqqkqqqpkqkkqc., 3,2q（5.4） 若在公式（5.1）中选择系数 及 ，使它满足 ii.0,0110ppc

56、ccc由定义可知此时所构造的多步法是 阶的. p,010kiinikiiniknfhyy（5.1）第85页/共161页87).()(2)1(11pnpppknhOxyhcT（5.5）称右端第一项为局部截断误差主项局部截断误差主项， 称为误差常数误差常数. 1pc 根据相容性定义， , 即 ，1p010 cc, 1110k故方法（5.1）与微分方程（1.1）相容的充分必要条件是（5.6）成立. 且（5.6）.011kikiikii由（5.4）得 ,010kiinikiiniknfhyy（5.1）),(11100kc),()1(2101211kkkkc第86页/共161页88 当 时，若 ，则由（

57、5.6）可求得 1k01.1, 100此时公式（5.1）为 ,1nnnfhyy即为欧拉法. 从（5.4）可求得 ，故方法为1阶精度，02/12c),()(21321hOxyhTnn 这和第2节给出的定义及结果是一致的. 且局部截断误差为 , 1110k（5.6）.011kikiikii2)!1(1)1(2(!11211121kqqkqqqpkqkkqc第87页/共161页89 对 , 若 ，方法为隐式公式.1k01为了确定系数 ，可由 解得100,0210ccc.2/1, 1100于是得到公式 ),(211nnnffhyy即为梯形法. 由（5.4）可求得 , 故 ，所以梯形法是二阶方法，12/

58、13c2p其局部截断误差主项是 .12/)(3nxyh 这与第9.2节中的讨论也是一致的. 2)!1(1)1(2(!11211121kqqkqqqpkqkkqc第88页/共161页90 对 的多步法公式都可利用（5.4）确定系数 ，并由（5.5）给出局部截断误差. 2kii,).()(2)1(11pnpppknhOxyhcT（5.5）第89页/共161页91 考虑形如 kiiniknknfhyy01（5.7）的 步法，称为阿当姆斯阿当姆斯（Adams)方法方法. k 为显式方法， 为隐式方法，通常称为阿当姆斯显式与隐式公式，也称Adams-Bashforth公式与Adams-Monlton公式

59、. 0k0k 这类公式可直接由对方程 两端从 到 积分求得. 1knxknx),(yxfy 第90页/共161页92 也可以利用(5.4)由 推出，01pcc对比 kiiniknknfhyy01与 ,010kiinikiiniknfhyy可知此时系数 . 1,01210kk 显然 成立，下面只需确定系数 ，00ck,10可令 ，求得 .011kcck,10 若 ，则令 来求得 . 0k01kcc110,k2)!1(1)1(2(!11211121kqqkqqqpkqkkqc第91页/共161页93 以 为例，由 ，根据 3k04321cccc, 13210,5)32(2321,19)94(332

60、1.65)278(4321 若 ，则由前三个方程解得03,1223,1216,125210得到 的阿当姆斯显式公式是3k).51623(121223nnnnnfffhyy（5.8）第92页/共161页94由（5.4）求得 ，所以（5.8）是三阶方法，08/34c局部截断误差是).()(835)4(43hOxyhTnn 若 ，则可解得 03.83,2419,245,2413210于是得 的阿当姆斯隐式公式为 3k),5199(2412323nnnnnnffffhyy（5.9）它是四阶方法，局部截断误差是 ).51623(121223nnnnnfffhyy（5.8）第93页/共161页95).()

61、(720196)5(53hOxyhTnn（5.10） 类似的方法可求得阿当姆斯显式方法和隐式方法的公式，表9-5及表9-6分别列出了 时的阿当姆斯显式公式与阿当姆斯隐式公式，其中 为步数， 为方法的阶， 为误差常数.4, 3,2, 1kkp1pc第94页/共161页96720251)9375955(244483)51623(1233125)3(22221116-912334122311211nnnnnnnnnnnnnnnnnnpffffhyyfffhyyffhyyfhyycpk公式阿当姆斯显式公式表第95页/共161页971603)19106264646251(7205472019)5199(

62、2443241)85(1232121)(2217123434123231212111nnnnnnnnnnnnnnnnnnnnnnpfffffhyyffffhyyfffhyyffhyycpk公式阿当姆斯隐式公式表9第96页/共161页98 例例6 6 用四阶阿当姆斯显式和隐式方法解初值问题 .1)0(, 1yxyy取步长 . 1.0h 解解本题 . nnhxxyfnnnn1.0, 1)9375955(2412334nnnnnnffffhyy从四阶阿当姆斯隐式公式得到 从四阶阿当姆斯显式公式得到 ).24. 324. 09 . 07 . 39 . 55 .18(241123nyyyynnnn第97

63、页/共161页99)5199(2412323nnnnnnffffhyy由此可直接解出 而不用迭代，得到 3ny).324. 01 . 05 . 01 .22(9 .241123nyyyynnnn计算结果见表9-8，其中显式方法中的初值 及隐式方法中的 均用准确解 得到，3210,yyyy210,yyyxexyx)().324. 01 . 05 . 01 .229 . 0(241123nyyyynnnn对一般方程，可用四阶R-K方法计算初始近似. (表略，见教材. )第98页/共161页100 从以上例子看到同阶的阿当姆斯方法，隐式方法要比显式方法误差小. 这可以从两种方法的局部截断误差主项的系

64、数大小得到解释.)()1(11npppxyhc这里 分别为 及 . 1pc720/251720/19第99页/共161页101 考虑另一个 的显式公式 4k),(01122334nnnnnnffffhyy其中 为待定常数.3210, 由（5.4）可知 ，再令 得到 00c04321cccc,43210 可根据使公式的阶尽可能高这一条件来确定其数值. ,16)32(2321,64)94(3321.256)278(43212)!1(1)1(2(!11211121kqqkqqqpkqkkqc),(11100kc第100页/共161页102解此方程组得 .0,38,34,380123于是得到四步显式公

65、式 ),22(341234nnnnnfffhyy（5.11）称为米尔尼米尔尼(Milne）方法方法. 由于 ，故方法为4阶，其局部截断误差为 45/145c).()(45146)5(54hOxyhTnn（5.12） 米尔尼方法也可以通过方程 两端积分 ),(yxfy 第101页/共161页1034)(,()()(4nnxxnndxxyxfxyxy得到. .)(,()()(22nnxxnndxxyxfxyxy右端积分通过辛普森求积公式就有 ).4(3212nnnnnfffhyy（5.13）称为辛普森方法辛普森方法. ).()(906)5(52hOxyhTnn（5.14） 它是隐式二步四阶方法，其

66、局部截断误差为 若将方程 从 到 积分，可得 nx2nx),(yxfy 第102页/共161页104 辛普森公式是二步方法中阶数最高的，但它的稳定性差，为了改善稳定性，考察另一类三步法公式 ),(332211221103nnnnnnnfffhyyyy其中系数 及 为常数. 210,321, 若希望导出的公式是四阶的，则系数中至少有一个自由参数. 若取 ，则可得到辛普森公式. 11 若取 ，仍利用泰勒展开，由（5.4），令01,043210ccccc第103页/共161页105.81)278(416,27)94(38,9)32(24, 32, 1321232123212321220解此方程组得 .83,86,83,89,8132120于是有 ),2(83)9(8112323nnnnnnfffhyyy（5.15）则可得到 第104页/共161页106称为汉明汉明(Hamming)方法方法. 由于 ，故方法是四阶的，且局部截断误差 40/15c).()(406)5(53hOxyhTnn（5.16）第105页/共161页107 对于隐式的线性多步法，计算时要进行迭代，计算量较大. 为了避免进行