1、2022-12-61常微分方程初值问题的数值解法常微分方程初值问题的数值解法8.2 几种简单的单步法 8.3 Runge Kutta 公式 8.4 单步法的收敛性、相容性与稳定性8.5 线性多步法第八章第八章8.1 引言数值算例2022-12-62本章着重讨论一阶常微分方程初值问题 00)(yxyyxfdxdyxy）（），（在区间a,b上的数值解法。8.1 引言一、问题 这些问题多数情况下求不出解析解，只能用近似的方法求解。常用的近似方法有两类。一类称为近似解析法，如级数解法，逐次逼近法等。另一类称为数值解法，它可以给出解在一些离散点上的近似值。2022-12-63若f（x,y）在区域D=上连
2、续，且关于y满足 李普希兹（Lipschitz）条件,即存在常数L,使Rybxa，2121),(),(yyLyxfyxf21,yy对D 内任两个 均成立，其中L是与x，y无 关的常数，则上面的初值问题存在唯一解，且解是连续可微的。二、解的存在唯一性Remark：在f(x,y)对y可微的情况下，若偏导数有界，则可取 ，此时有之间。与介于21212121y,)(),(),(),(yyyLyyyxfyxfyxfyyxfLDyx),(max),(此时Lipschitz连续条件显然成立。这是验证该条件的最简便的方法。2022-12-64解的适定性是指解的存在唯一性以及数值稳定性。此处主要是指解对于右端项
3、以及初值扰动的适应性。关于适定性有如下的结论：三、解的适定性定理：若f（x,y）在区域D=上满足Lipschitz连续条件，则初值问题是适定的。Rybxa，四、等价的积分方程若y(x)是初值问题的解，对方程两边同时积分，利用初始条件可得：xxdttytfxyxy0)(,()()(0该方程为与初值问题同解的积分方程，我们可以从积分方程出发去构造初值问题的求解公式。2022-12-65 常微分方程初值问题的数值解是求上述初值问题的解y(x)在区间a，b中的点列），（nihxxiii101上的近似值 。称为步长，iyih一般情况下我们取等步长，记为h。五、数值解法)(nxy初值问题的解析解（理论解）
4、用 表示,数值解法的精确解用 表示，并记fn=f(xn,yn),而 。ny)(,()(nnnxyxfxy求初值问题的数值解一般是逐步进行的，即计算出yn之后计算yn+1。2022-12-66数值解法一般分为：（1）单步法：在计算 时，只用到 ,和 ，即前一步的值。nx1nyny1nx （2）多步法：计算 时，除用到 ,和 以外，还要用 和 ,即前 k步的值。1nypny1nxnxnypnx)0;,2,1(kkp数值解法的分类单步法和多步法都有显式和隐式方法之分。显式和隐式的单步法可以分别写成：），（hyxhyynnnn,1），（hyyxhyynnnnn,11对多步法来说，显式和隐式方法具有相同
5、的意义。2022-12-67六、离散化方法建立常微分方程初值问题数值解法的过程，就是通过一定的离散化方法，将对连续性问题的求解转化为求解常微分方程在有限个离散节点上解的近似值的过程，这个过程通常称为数值离散数值离散常用的数值离散化方法有：n差商代替微商法nTaylor展开法 n数值积分法 2022-12-68一、显式欧拉（一、显式欧拉（Euler）公式公式 设节点为 ，欧拉方法的计算公式），（nnnhxxn,2100），）（，（2101nyxhfyynnnn8.2 几种简单的单步法这是一种最简单的显式单步方法，该方法可以通过不同的途径获得。),2,1,0)(,(1nyxfhyynnnn1、差商
6、法用两点差商公式 代替导数 ,再用 表示 的近似值,则得到nnnnxxxyxy11)()(ny)(nxy)(nxy2022-12-69假设在 附近把y（x）展成Taylor级数nx)(2)()()(2nnnnxyhxyhxyhxy取h的线性部分,并用 表示 的近似值,得ny)(nxy),2,1,0)(,(1nyxhfyynnnn显式欧拉（显式欧拉（Euler）方公式方公式（续）2、Taylor展开法3、数值积分法1)(,()()(1nnxxnndttytfxyxy对微分方程两端从xn到xn+1积分，得等价的积分方程对右端的积分部分采用左矩形公式近似即得Euler公式。Remark：Taylor
7、展开法与数值积分法是构造微分方程数值解的两类主要的方法。2022-12-610欧拉（欧拉（Euler）方法的几何意义方法的几何意义Euler方法有明显的几何意义。如右图所示，一阶常微分方程初值问题的解曲线y(x)过点P0(x0,y0)。从P0出发以f(x0,y0)为斜率作一直线段，与x=x1相交于点P1(x1,y1)，显然有y1=y0+hf(x0,y0)。同理，再由P1出发以f(x1,y1)为斜率作一直线段，与x=x2相交于点P2(x2,y2)，显然有y2=y1+hf(x1,y1)。这样一直做下去，得到一条折线P0P1P2，作为y(x)的近似曲线。因此，显示Euler公式又称为Euler折线法
8、。2022-12-611二、隐式二、隐式Euler公式公式dttytfxyxynnxxnn），（）（）1)(1中将积分用右矩形公式),(,)(2)(2baabfbfabdxxfba）（）（）（代入，有,2,1,0),(,()()(111nxyxhfxyxynnnn上式是一个隐式的单步方法，称为隐式欧拉公式或后退的欧拉公式。利用此公式，每一步都要把上式作为yn+1的一个方程来求解。从数值积分的误差分析，很难期望隐式欧拉法比显式欧拉法更精确。若在等价积分方程 从而得到,2,1,0),(111nyxhfyynnnn2022-12-612隐式隐式Euler公式（续）公式（续）通常情况下，隐式欧拉公式很
9、难直接求出yn+1的值，故常用迭代法求解。在实际计算时，该公式通常与显式Euler公式结合使用，并由显式Euler公式的结果作为迭代的初始值，从而有如下数值格式1,2,1),(),()(11)1(1)0(1Nnyxhfyyyxhfyysnnnsnnnnn对 循环计算，若 (为给定的误差限)，则取 作为yn+1的近似值,2,1,0s)(1)1(1snsnyy)1(1sny2022-12-613隐式隐式Euler公式（续）公式（续）由于f(x,y)关于y满足Lipschitz条件，故有1)0(111)(111)(111)1(1)(),(),(nnsnsnnnsnnnsnyyhLyyhLyxfyxf
10、hyy故当 时该迭代法收敛到隐式Euler公式的解yn+1，其中L为Lipschitz常数1hL2022-12-614三、梯形公式三、梯形公式为了得到更精确的方法，在等价的积分方程中用梯形公式),(),(1223bafabbfafabdxxfba）（）（）（）（,2,1,0),(),(2111nyxfyxfhyynnnnnn近似积分项，再分别用yn，yn+1代替y(xn)和y(xn1)，即可得到梯形公式：该方法也是一种隐式的单步方法。对该方法，从n0开始计算，每步都要求解yn+1的一个方程。一般来说，这是一非线性方程，可迭代计算如下：2022-12-615,2,1,0,),(),(2),(11
11、 1101syxfyxfhyyyxhfyysnnnnnsnnnnn梯形公式（续）梯形公式（续）使用上式时，先用第一式计算出yn+1的近似值 ，再用第二式反复进行迭代，得到数列 ，用 来控制是否继续进行迭代，其中为允许误差。把满足要求的 作为y(xn+1)的近似值yn+1，类似地可得出yn+2，yn+3，。01ny 01 ssny1 11snsnyy11sny2022-12-616011111111111122)2nnssnnsnnnnsnnyyhLyyhLyxfyxfhyy，（），（梯形公式（续）梯形公式（续）当f(x,y)关于y满足Lipschitz条件时，且步长h满足 时，上述迭代过程是收
12、敛的。这是因为：121hL2022-12-617实用中，h 取得较小时,为了简化计算，梯形公式第二式只迭代一次就结束，得到Euler预测校正方法(改进的Euler方法)：),(),(2),(011101nnnnnnnnnnyxfyxfhyyyxhfyy其中第一式称为预测算式，第二式称为校正算式。四、四、EulerEuler预测校正公式预测校正公式2022-12-618Euler预测校正公式（续）预测校正公式（续）若将Euler预-校方法中的第一式带入第二式，得),(,(),(211nnnnnnnnyxhfyxfyxfhyyRemark：这是一种显示的单步方法。有时为了计算方便，常将上式改写成：
13、),(),()(2121211hKyhxfKyxfKKKhyynnnnnn2022-12-619五、单步法的局部截断误差和阶五、单步法的局部截断误差和阶 设一般的单步法为：显式公式：），（hyxhyynnnn1隐式公式：），（hyyxhyynnnnn11设 为数值方法的精确值，y（）为微分方程的精确解。nxnynnnyxye）（nx定义1：为某一数值方法在 处的整体截断误差。Remark：整体截断误差不仅与 这步的计算有关，而且与前面所有点的计算的误差累计有关。为了简化误差分析，我们着重分析计算中的某一步。对一般的显式单步法，有如下定义：nx2022-12-620单步法的局部截断误差和阶单步法
14、的局部截断误差和阶（续）),(,()()(),()()(1111hxyxhxyxyhyxhyxyyxynnnnnnnnnn这就是上面定义中称Rn+1为“局部”的含义，我们一般用该式作为定义。这里应该注意，Rn+1和整体截断误差en+1是不同的。）（nnxyy 若设 ，即第n步及以前各步都没有误差，则由显示单步法计算一步所得之 与 之差为：1ny)(1nxy即在 的假设下，。111)(nnnyxyR）（nnxyy 定义2：对单步法，在 的假设下，称 为在 处的局部截断误差。）（nnxyy),(,()()(11hxyxhxyxyRnnnnn1nx2022-12-621单步法的局部截断误差和阶单步法
15、的局部截断误差和阶（续）Remark：由前面的定义可知，若某个单步方法是一种p阶方法，则有Rn+1=O(hp+1)，即p阶方法的局部截断误差为h的p1阶。我们往往比较关心Rn+1按h展开式的第一项。定义4：若一个单步方法是一种p阶方法，其局部截断误差可以写成：)()(,(211ppnnnhOhxyxR则(xn,y(xn)hp+1称为方法的主局部截断误差，或局部截断误差的主项。定义3：若一个单步方法的局部截断误差为O(hp+1)，即则称该方法为p阶方法（其中p为正整数）。)(),(,()()(111pnnnnhOhxyxhxyxyR2022-12-622例1：求显示Euler公式的局部截断误差。
16、)()(2321hOxyhRnn单步法的局部截断误差和阶单步法的局部截断误差和阶（续）故显示Euler公式是一阶方法,局部截断误差为：主局部截断误差为：。)(22nxyh)()(!31)(21)()()()(,()()(),(,()()(2321111hOxyhxyhxyhxyxyxyxhfxyxyhxyxhxyxyRnnnnnnnnnnnnnn 2022-12-623例2：求Euler预测校正公式),(),(2),(011101nnnnnnnnnnyxfyxfhyyyxhfyy的局部截断误差。单步法的局部截断误差和阶单步法的局部截断误差和阶（续）)()(,2)(2)()(,()(,)(,(2
17、)()(),(,()()(11111nnnnnnnnnnnnnnnnnnnxyhxyhxfhxyhxyxyxyxhfxyxfxyxfhxyxyhxyxhxyxyR2022-12-624单步法的局部截断误差和阶单步法的局部截断误差和阶（续）)()(,(2)()(,(2)(2)(,(2)(,()()(,()(,(2)(2)(32222222221hOyxyxfxyhyxxyxfxyhxxyxfhyxyxfxyhxxyxfhxyxfhxyhxyxynnnnnnnnnnnnnnnnnn又由)(,()(xyxfxyyxyxfxyxxyxfxy)(,()()(,()(故)()()(2)()(33211hO
18、hOxyhxyhxyxyRnnnnn 故Euler预测校正方法为二阶方法。2022-12-625单步法的局部截断误差和阶单步法的局部截断误差和阶（续）例3：求隐式Euler公式的局部截断误差。)()()(,()()(),(),(,()()(11111111nnnnnnnnnnnnnxyhxyxyxyxhfxyxyhxyxyxhxyxyR)()(2321hOxyhRnn 故隐式Euler公式是一阶方法。将上式中 、均在 处做Taylor展开，整理得)(1nxy)(1nxynx2022-12-626类似地可以证明，梯形公式的局部截断误差为：单步法的局部截断误差和阶单步法的局部截断误差和阶（续）)(
19、)(12431hOxyhRnn 即梯形公式为二阶方法。2022-12-627一、一、Taylor 方法方法 假定初值问题的解y(x)以及函数f(x,y)是足够光滑的，有)()(!)(!2)()(1)(21pnppnnnnhOxyphxyhxhyxyxy）（)()(,(!)(,(!2)(,()(1)1(2pnnppnnnnnhOxyxfphxyxfhxyxhfxy8.3 Runge Kutta 公式 Taylor展开法与数值积分法是推导高阶方法的常用手段。本节以Taylor展开法为基础，介绍如何推导高阶单步法。)(21!2pnpnnnnyphyhyhyy 当h充分小时，略去余项 ，将y(k)(x
20、n)用 来代替，则有p阶Taylor方法：)(1phO)(kny)(kny其中 (k=1,2,p)根据求导法则，其计算公式为：2022-12-628),(nnnyxfy),(),(),(nnynnnnxnyxfyxfyxfy),(2 22nnyyxyyxyxxnyxfffffffffyTaylor 方法（续）方法（续）Remark1：显然，p阶Taylor方法的局部截断误差为 。当p=1时，Taylor方法就是Euler方法。当p2时，需要计算公式中的高阶导数。)()()!1(12)1(11pnppnhOxyhpRRemark2：显然，Taylor方法可以得到任意阶精度的方法。但在实际计算中，
21、Taylor方法往往相当困难，因为公式中的高阶导数会很复杂。故Taylor方法很少单独使用，但常用它来启发思路。2022-12-629二、二、RungeKutta方法方法Euler公式：）（），（21111hORyxfkhkyynnnnn改进的Euler预测校正格式：)()231121211hORhkyhxfkyxfkkkhyynnnnnnn），（），（基本思想：用不同点的函数值作线性组合，构造高阶单步的近似公式，把近似公式和解的Taylor展开式比较，使前面尽可能多的项完全相同。这种方法间接应用Taylor展开的思想，避免了高阶导数计算的困难。2022-12-630一般的Runge-Kutt
22、a方法的形式为rikhyhxfkyxfkkchyxhyxhyyijjijnininnriiinnnnnn,3,2,),(),(11111，），（其中，为常数。选取这些常数的原则是使其截断误差阶尽可能高。iijiC，RungeKutta方法（续）方法（续）2022-12-631下面以二级RungeKutta公式为例进行具体推导。），（），（12122122111),(),(KhyhxfKyxfKKcKchyxhyxhyynnnnnnnnnn对 要求适当选取系数 ，使当 时，上式的局部截断误差为 ，即成为二阶方法。21221,cc）（nnxyy）3(hORungeKutta方法（续）方法（续）20
23、22-12-632RungeKutta方法（续）方法（续）1112221231222212()()(,(),)()()(),()()11()()()()()()2!3!1(,)(,)()(,)()(2nnnnnnnnnnnnnnnnxyxxy xy xhxy xhy xy xh c y xc f xh y xhy xy xhy xh y xh yxy xc hy xc hf x yhfx yhy x fx yhfx22322121(,(),)1()(,)()(,)()2nnxyyyxy xyhy x fx yhy xfx yO h按照局部截断误差的定义，有2022-12-633RungeKut
24、ta方法（续）方法（续）21222221(,()32222221224221(,()1(1)()()(,)()(,)21()(,)()(,)3!1()(,)()2nnnnnxyxy xxxxyyyxy xcc hy xhy xcfx ycy x fx yhyxcfx ycy x fx ycy xfx yO h 21222221(,()32222222212214(,()11(1)()2211111623621()()6nnnnnxyxy xxxxyyyyxyxy xcc hy xhcfcffhcfcffcf fffffO h 2022-12-634RungeKutta方法（续）方法（续）要使上
25、式等于 ，只需满足)(3hO021021012122221cccc在上式中有四个未知量，三个方程，故可以得到无穷多组解，也就是可以得到无穷多个二级二阶Runge-Kutta公式由于是二级方法，从而有 ，此时上式可以改写为02c221221211ccc在上述方程中，选取不同的c2，即可获得不同的二级二阶Runge-Kutta方法。常用的二阶Runge-Kutta方法有：2022-12-6351,2121221cc),(),()(2121211hKyhxfKyxfKKKhyynnnnnnEuler预测校正格式)2,2(),(12121KhyhxfKyxfKhKyynnnnnn中（间）点公式21,1
26、,021221cc RungeKutta方法（续）方法（续）2022-12-63632,43,4121221cc)32,32(),()3(4121211hKyhxfKyxfKKKhyynnnnnn二阶Heun（休恩）方法此时可以论证：,02161222c0312212c此时 ，即二阶Heun方法是局部截断误差项数最少的方法且二级Runge-Kutta方法不可能达到三阶。34,(,()()()6nnn hyxyxy xhRffffO hRungeKutta方法（续）方法（续）2022-12-637RungeKutta方法（续）方法（续）用类似的方法可以研究其它级的Runge-Kutta方法。常用
27、的其它级的Runge-Kutta方法有：)32,32()31,3(),()3(423121311hKyhxfKhKyhxfKyxfKKKhyynnnnnnnn三阶Heun（休恩）公式)2,()21,2(),()4(62131213211hKhKyhxfKhKyhxfKyxfKKKKhyynnnnnnnn三阶Kutta公式 2022-12-638标准（经典）四阶Runge-Kutta公式 基尔（Gill）公式RungeKutta方法（续）方法（续）),()21,2()21,2(),()22(6342312143211hKyhxfKhKyhxfKhKyhxfKyxfKKKKKhyynnnnnnnn
28、nn)221(22,()221(212,2()21,2(),()22()22(632421312143211hKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKKhyynnnnnnnnnn2022-12-639定理（整体截断误差与局部截断误差的关系）：若初值问题的一个单步方法之局部截断误差为 且单步法中函数 关于y满足Lipschitz条件，则)1(),(11phORpn),(hyx)()(111pnnnhOyxye定义：对于初值问题，如果一个单步显式方法产生近似解对于任一固定的 ,均有 则称该单步法是收敛的。nhxxbxx00,)(lim0 xyynh 一、收敛性8.4 单步法的
29、收敛性、相容性与稳定性Remark：该定义可类似地应用于单步隐式方法以及后面的线性多步法。从定义可知，若格式收敛，整体截断误差en=y(xn)-yn必然趋于零。2022-12-640)(),(,()()(111pnnnnnhOhxyxhxyxyR证明：根据局部截断误差的定义，有记),(,()(1hxyxhxyynnnn111)(pnnchyxy则有常数c0，使得 整体截断误差与局部截断误差的关系（续）由于 且 关于 满足Lipschitz条件，故有),(1hyxhyynnnn),(hyxy11nnyy),(),(,()(hyxhxyxhyxynnnnnnnnyxyhL)()1(2022-12-
30、641111)(nnnyxye故npehLCh)1(11111)(nnnnyyyxy整体截断误差与局部截断误差的关系（续）从而有递推关系：)1()1(1111nppnehLchhLche0111)1(1)1(1)1(ehLhLhLchnnp121)1()1(1 npehLhLch01)1(ehLn)1()1()1(1 21nphLhLhLch2022-12-642得：hLnnehL)1(1)1(1故11)1()1(11hLnPhLnpnehLchLeche当 固定时，1nxxabxxhnn01)1(若 则 且,)(00yxy00epLabpnhcehLce1)(11 所以整体截断误差与局部截断
31、误差的关系（续）证毕hLnehLnhLhLhL)(!1)(!2111122022-12-643整体截断误差与局部截断误差的关系（续）Remark2：在该定理的条件下，欧拉方法是一阶方法，欧拉预测校正方法是二阶方法。当f(x,y)关于y满足利普希茨条件时，r级Runge-Kutta方法中的 关于y也满足利普希茨条件，所以定理中的条件得到满足，解的收敛性得到了保证。),(hyxRemark1：该定理说明，数值方法的整体截断误差比局部截断误差低一阶。收敛的方法至少是一阶方法。要判断一个单步法的收敛性，一是要知道其局部截断误差，二是要判断其增量函数 是否关于y满足Lipschitz条件。),(hyx2
32、022-12-644整体截断误差与局部截断误差的关系（续）事实上，对于Euler预测校正方法，其增量函数),(,),(21),(yxhfyhxfyxfhyxyyhLLyxfyxfhyyLyyLyxhfyhxfyxhfyhxfyxfyxfhyxhyx2),(),(2121),(,(),(,(21),(),(21),(),(2故有从而Euler预测校正方法的增量函数 关于y满足Lipschitz条件。),(hyx2022-12-645二.相容性相容性)()0),(,()()0),(,()0),(,(2)()(),(,()()(2211hOxyxxyhhhxyxxyxhhxyhxyhxyxhxyxy
33、Rnnnnnnnnnnnnnn 故方法阶p1的充要条件是 ，而 ，从而可以给出如下的定义：0)0),(,()(xyxxy)(,()(xyxfxy由整体截断误差与局部截断误差的关系定理知，当 p1时单步方法收敛，且当y(x)是初值问题的准确解，显式单步方法有p阶精度时，若 关于h连续，则其局部截断误差),(hyx2022-12-646)(,()0),(,(xyxfxyx相容性（续）相容性（续）定义定义：如果单步方法的增量函数?(x,y,h)满足相容的充要条件为)(,()0),(,(xyxfxyx且?(x,y,h)关于h满足Lipschitz条件，则称单步方法与初值问题相容相容。Remark：以上
34、讨论表明，若?(x,y,h)关于h满足Lipschitz条件，p阶方法当p1时与初值问题相容，反之相容的方法至少是一阶方法Euler方法，Euler预测校正格式，Runge-Kutta方法等都与原微分方程相容。2022-12-647三、稳定性稳定性定义定义1*用一个数值方法求解微分方程初值问题时，对给定的步长h0，若在计算 时引入误差 （也称扰动），由此引起计算后面的 时的误差按绝对值均不增加，则称这个数值方法是绝对稳定的绝对稳定的。nyn),2,1(kykn关于单步法收敛性的概念和收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的。但实际计算时通常会有舍入误差及其累积。数值求解微分方
35、程的过程是一个递推公式，必须考虑误差积累是否得到控制。0Re）（为保证微分方程本身的稳定性，这里假定 。设f(x,y)关于y满足Lipschitz条件，我们总是针对模型方程 进行讨论，其中 为复常数。yy2022-12-648定义定义1 设步长为h0的单步法用于求解 时，yy中由 引起的误差 满足 ，则称nm)(nmnm单步法对于所用步长h和复数 是绝对稳定的。若在计算 时有误差 ,但在计算后面的)(nmymnny稳定性（续）稳定性（续）Remark1：在上面的定义中，可以取小于或等于关系符。取小于号是为了和线性多步法相一致。Remark2：从上面的定义可知，单步法是否稳定，与模型方程中的复数
36、以及所用步长h有关。若对复平面上的某个区域G，当 时，单步法绝对稳定，则称G为单步法的绝对稳定区域，G与实轴的交集为绝对稳定区间。Gh2022-12-649Euler显式公式显式公式将Euler显式公式用于模型方程 ，则有yynnnnyhyhyy)1(1稳定性（续）稳定性（续）下面考察几个常用公式的稳定性。设 有误差 ，nynnnnyy参与运算的量为 ，误差为 ，则实际得到近似 的量为1ny1n由此引起的要求误差不增加，即 ，nn11ny111nnnyy，即 。与前面公式)(1(11nnnnyhynnh)1(1相减，有 。必须 。111hnn是保证绝对稳定性对11h步长h加的限制。当为实数时，
37、可以得到用h表示的绝对稳定的区间(-2,0)。2022-12-650隐式隐式Euler公式公式后退的Euler公式用到 上，有 ，yy11nnnyhyy故 11nnnh即 nnh111令1111hnn11h得绝对稳定区域为 。稳定性（续）稳定性（续）可见，若取为实数，则对于任意的h都是绝对稳定的。2022-12-651梯形公式梯形公式 梯形公式用到 上，有yy)(211nnnnyyhyy故)(211nnnnh 即 nnhh21211 设iyxh稳定性（续）稳定性（续）2022-12-6522122224444yxxyxx212222)2()2(yxyxiyxiyx222121hh当x=Re()
38、0,上式右端总小于1，故梯形法的绝对稳hh定区域为Re()0，即左半平面。由于 nn1稳定性（续）稳定性（续）2022-12-653四阶经典四阶经典R-K方法方法四阶经典R-K方法用到 上，有 yynyk1)21()21(212hyhkyknn因此)21(21)21(223hhyhkyknn)4121()(32234hhhyhkyknn)4121()4121(243322322hhhhh)22(643211kkkkhyynn)21(262hyhynn)(241)(61)(211432hhhhyn稳定性（续）稳定性（续）2022-12-654扰动满足)(241)(61)(211 4321hhhh
39、nn令 ,得四阶经典RK方法的绝对稳定区域为11nn1)(241)(61)(211432hhhh稳定性（续）稳定性（续）2022-12-655步的信息来预测 ，则可以期望获得较高的精度。kny线性多步法的基本思想基本思想：如果充分利用前面多8.5 线性多步法一、线性多步法的一般公式)(nxy记 的近似值为 ,并记 ，则k步线性多步方法一般形式为 nhxxynn0,),(nnnyxffjj,11,kk其中 为常数，不全为零。1,.,1,11101Nkknfhyykiiniinkiin2022-12-656 若 则为隐式方法，若 则为显式方法。0101上式称为线性k步法步法。线性多步法的一般理论（
40、续）或写成Remark：前面的RK方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计算量增加了许多。而线性多步法每步只需要计算一个函数值。1010111kiiniinkiinnnnfhygfhgy2022-12-657线性多步法的一般理论（续）对于隐式情形（）的公式，由于f(x,y)一般是非线性函数，故难以求解到yn+1的显示表达式，故常用迭代法求解：01其中 任意给出，s0,1,2,迭代到满足给定的精度要求为止。)0(1ny可以证明，当f(x,y)满足Lipschitz条件（或 ）时，只要 ，迭代关系式就是收敛的。Lyf-11hL),()(1111)1(11snnnnsnny
41、xfhgfhgy2022-12-658定义 处的局部截断误差为 1nx 定义：对于线性多步法线性多步法的一般理论（续）1,.,1,11101Nkknfhyykiiniinkiin111101()()(,()kknnin iin in iiiRy xy xhf xy x若)()(2)1(111pnpppnhOxyhcR称 为线性多步法的主局部截断误差，并称该线性多步法为p阶方法。)()1(11npppxyhC2022-12-659设 ，则有 1,1,0i),(iikxyynn线性多步法的一般理论（续）111011)()()(kiiniinkiinnxyhxyRxy将上式与线性多步法公式相减，有1
42、1111111111111111111)(,),()(,()()(nnnnnnnnnnnnnnnyxyyxfhRyxfhxyxfhRyhxyhRyxy其中n+1介于y(xn+1)与yn+1之间，从而得到111111)(,1nnnnnyxyyxfhR2022-12-660线性多步法的一般理论（续）若线性多步法为显式方法（），则01111)(nnnyxyR这与显式单步方法的含义是一致的。如果线性多步法为隐式方法（），且该方法为p阶方法，则当y（x）充分可微时01)()()()(,1)()(,1)(2)1(112)1(111112)1(11111111pnppppnpppnnpnpppnnnnhOx
43、yhchOxyhcyxfhhOxyhcyxfhyxy即Rn+1与y(xn+1)-yn+1的局部截断误差主项相同。2022-12-661线性多步法的一般理论（续）Remark：可以证明，显式线性多步法的整体截断误差比局部截断误差低一阶。2022-12-662二、基于数值积分的构造方法数值积分的构造方法将 方程两端从 积分得)(,(xyxfy 1knnxx到（*）1)(,)()(1nknxxknndxxyxfxyxy记 ，利用Lagrang插值多项式来近似代替（*）式中的被积函数 ，则可以导出不同的线性多步法公式。下面导出常用的四阶Adams外插、内插公式。)(,()(xyxfxF)(xF2022
44、-12-663其插值余项为基于数值积分的构造方法（续）数值积分的构造方法（续）在（*）式中取k=0，并选择等距节点 作为插值节点，作函数F（x）的三次插值多项式321,nnnnxxxx33300()()njn iijn injj ixxL xF xxx nnnnnnxxxxxxxxxxFxR3321)4(3),)()()()(!41)(将 代入（*）式，得33()()()F xL xR x11133()()()()nnnnxxnnxxy xy xL x dxR x dx2022-12-664基于数值积分的构造方法（续）数值积分的构造方法（续）略去上式右端第三项，得113()()(),3,4,n
45、nxnnxy xy xL x dx n对上式积分部分作变量代换 ，并注意到thxxnhxxxxxxnnnnnn32211则1113023123()()()(1)(2)(3)(2)(3)3!2()()(1)(3)(1)(2)23!55()59()37()9()24nnxnnxnnnnnnF xF xL x dxtttt ttF xF xt ttt tthdthF xF xF xF x2022-12-665线性四步Adams显式公式基于数值积分的构造方法（续）数值积分的构造方法（续）3213322111937595524),(9),(37),(59),(5524nnnnnnnnnnnnnnnfff
46、fhyyxfyxfyxfyxfhyy其局部截断误差为1)()()()(!41321)4(1nnxxnnnnndxxxxxxxxxFR2022-12-666基于数值积分的构造方法（续）数值积分的构造方法（续）利用第二积分中值定理，得)(720251)(720251)()()()(!41)5(5)4(5321)4(11yhFhdxxxxxxxxxFRnnxxnnnnn,1nnxx由于插值多项式 是在 上作出的，而积分区间是 ，故上述线性四步Adams显式公式称为Adams外插公式外插公式。)(3xLnnxx,31,nnxx2022-12-667基于数值积分的构造方法（续）数值积分的构造方法（续）若
47、在（*）式中取k=0，但选择等距节点 作为插值节点，作函数F（x）的三次插值多项式，仿照上面的做法，可以得到线性三步的Adams隐式公式及其局部截断误差211,nnnnxxxx2111519924nnnnnnffffhyy)(72019)5(51yhRn上式称为Adams内插公式内插公式。由于Adams内插公式是隐式公式，故用它计算时也需用迭代法。通常把Adams外插公式与内插公式联合起来交替使用，先由前者提供初始值，再由后者进行修正，即 2022-12-668基于数值积分的构造方法（续）数值积分的构造方法（续）,2,1,0,519,92493759552421)(11)1(1321)0(1s
48、fffyxfhyyffffhyynnnsnnnsnnnnnnn21)0(111321)0(1519,924937595524nnnnnnnnnnnnnfffyxfhyyffffhyy可以看出，无论四步Adams显式公式还是三步Adams隐式公式，均为四阶公式。类似于以上讨论，还可以采用数值积分方法得到其它p阶的线性多步法公式。可以证明，当 时，迭代式收敛若上式中的第二式只迭代一次，则可以得到Adams预测校正公预测校正公式式183hL2022-12-669三、基于基于Taylor 展开的构造方法展开的构造方法 Taylor展开法更具一般性。下面构造线性两步方法利用线性多步法局部截断误差的定义，
49、有110111101nnnnnnfffhyyy111101101111011()()(,()()()()()()()kknnin iin in iiinnnnnnRy xy xhf xy xy xy xy xhy xhy xhy x将上式在xn处作Taylor展开，并按h的升幂整理排列，得到2022-12-670基于基于Taylor 展开的构造方法（续）展开的构造方法（续）)(24124112011201)(6161241241)(21216161)(2121)(1)(1)5(5111)4(4111311121111011101nnnnnnnxyhxyhxyhxyhxyhxyR2022-12-
50、671基于基于Taylor 展开的构造方法（续）展开的构造方法（续）241616124161212161212111111111111101110令求解上述方程组，得出 ，所得到的公式的局部截断误差至少为 。10110,)(5hO2022-12-672基于基于Taylor 展开的构造方法（续）展开的构造方法（续）可以只要求前面的几个方程成立，例如要求前面的四个方程成立时，所得公式的局部截断误差至少为 。由于此时方程个数少于未知量个数，故此种情形下方程组有无穷多组解。此时方程的解可以写成)(4hO0100010112531,3234,12131,1111143nnnnnfffhyy)()(901
