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

1、第7章 常微分方程初值问题的数值解法7.1 引 言科学研究和工程技术中的许多问题在数学上往往归结为微分方程的求解问题。为了确定微分方程的解，一般要加上定解条件，根据不同的情况，这些定解条件主要有初始条件(initial condition)和边界条件(boundary condition). 只含初始条件作为定解条件的微分方程求解问题称为初值问题(initial-value problem); 例如天文学中研究星体运动,空间技术中研究物体飞行等，都需要求解常微分方程初值问题(initial-value problem for ordinary differential equations).

2、只含边界条件作为定解条件的微分方程求解问题称为边值问题(boundary-value problem).除特殊情形外，微分方程一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，而且实际问题往往只要求在某一时刻解的函数值. 为了解决这个问题，有两种方法可以逼近原方程的解。第一种方法是：将原微分方程化简为可以准确求解的微分方程，然后使用化简后的方程的解近似原方程的解；第二种方法是：将求原微分方程的解析解转化为求原方程的数值解，这是实际中最常用的方法。本章将介绍求解常微分方程初值问题的常用的数值方法。第8章将介绍常微分方程的边值问题的常用的数值方法。为简明起见，本章主要介

3、绍形如 (7.1.1)的初值问题的数值解法. 在介绍这些方法之前，需要了解常微分方程的一些相关定义和结果.定义7.1.1 函数称为在集合上关于变量满足Lipschitz条件(Lipschitz condition)，简称Lip条件，如果存在常数，使得 (7.1.2)对所有都成立. 常数称为Lipschitz常数(Lipschitz constant)，简称Lip常数.定义7.1.2 如果对所有，都有， (7.1.3)其中，那么称集合为凸集(convex set).注 如果没有特别说明，本章中的集合为，显然它是凸集。定理7.1.1 设函数在上连续. 如果在上关于变量满足Lipschitz条件，那

4、么初值问题(7.1.1)对有唯一解.定理7.1.2 设函数定义在凸集上. 如果存在常数，使得 (7.1.4)对一切成立，那么在上关于变量满足Lipschitz条件，且为Lipschitz常数. 注 因为Lipschitz条件一般不易验证，所以往往用条件(7.1.4)代替Lipschitz条件(7.1.2)。由定理7.1.1立得：定理7.1.3 设函数在上连续. 如果存在常数，使得对一切成立，那么初值问题(7.1.1)对有唯一解.理论上，初始值是准确的，微分方程的计算过程也是准确的；而实际情况并非如此. 事实上，初始数据可能有误差，计算过程中的数字的舍入也会产生误差。对初值问题(7.1.1)，若

5、初值有误差，即实际得到的初值是，再考虑到计算过程中的数字的舍入也会产生误差，则得到的初值问题变为 (7.1.5)这种误差的扰动在计算过程中是否会增长很快，以至影响结果，这就是初值问题的适定性问题.定义7.1.3 设初值问题(7.1.1)和(7.1.5)的均有唯一解，分别是和，且连续。若对任意，存在常数，当时，对任意，恒有，则称初值问题(7.1.1)是适定的(well-posed)。定理7.1.6 设函数在上连续. 如果在上关于变量满足Lipschitz条件，那么初值问题(7.1.1)是适定的。注 初值问题只有是适定的才有意义，因此本章考察的所有初值问题都是适定的。求微分方程数值解的主要问题：(

6、1) 如何将原微分方程离散化, 并建立求其数值解的递推公式;(2) 如何求递推公式的局部截断误差, 数值解与精确解的误差估计；(3) 研究递推公式的稳定性与收敛性.7.2 Euler方法及改进的Euler方法7.2.1 Euler格式与梯形格式考虑一阶常微分方程的初值问题 (7.2.1)设，其中为等距节点，步长. 在上对两边积分，得 (7.2.2)7.2.1.1 Euler格式用左矩形求积公式计算式(7.2.2)右端积分项，得,代入式(7.2.2)右端，得. (7.2.3)用的近似值代入式(7.2.3)右端，记所得结果为，得 (7.2.4)并称式(7.2.4)为求解初值问题(7.2.1)的Eu

7、ler方法(Eulers method)或Euler格式(Euler scheme)，称为差分方程(difference equation).注 Euler方法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它的精度不高，很少被采用，但是它反映了微分方程的数值解法的基本思想和特征.若式(7.2.2)右边的积分由数值积分的右矩形公式近似，并用近似值替代，近似值替代，则可得到 (7.2.5)并称式(7.2.5)为后退的Euler方法(backward Eulers method) 或后退的Euler格式(backward Eulers scheme). 是差分方程。注 在平面上，微分方程的解

8、称为积分曲线，积分曲线上一点的切线斜率等于函数在点的值. 如果在中每一点都画一条以在点的值为斜欧拉(Leonhard Paul Euler 1707年4月15日 1783年9月18日) 是瑞士杰出的数学家，自然科学家. 他在数学分析，图论，力学，光学和天文学等很多领域都做出了卓越的贡献，被公认为有史以来最伟大的数学家之一.率并指向增加方向的有向线段(即：在上作出了一个由方程确定的方向场），那么从几何上看，微分方程的解就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向相一致.从初始点出发，过这点的积分曲线为，斜率为. 设在附近可用过点的切线近似表示，切线方程为. 当时，的近似值为

9、，并记为，这就是得到时计算的近似公式当时，的近似值为，并记为. 于是就得到当时计算的近似公式重复上面方法，一般可得当的计算的近似公式如果，则上面公式就是式(7.2.4). 将连起来，就得到一条折线，所以Euler方法又称为折线法(polygon method), 见图7.2.1.由公式(7.2.4)知，已知便可算出；已知，便可算出；如此继续下去，这种只用前一步的值便可计算出的递推公式称为单步法(one-step method).图7.2.17.2.1.2 梯形格式用梯形求积公式计算式(7.2.2)右端积分项，得,代入式(7.2.2)右端，得. (7.2.6)在式(7.2.6)中，将用近似替代，

10、所得结果记为，得 (7.2.7)并称(7.2.7)为求解初值问题(7.2.1)的梯形方法(trapezoidal method), 或梯形格式(trapezoidal scheme).是差分方程.7.2.2 改进的Euler方法因为梯形求积公式的精度比矩形求积公式的精度高，所以求解初值问题(7.2.1)的梯形方法比Euler方法精度高；但梯形方法是隐式方法，不便于计算。为了克服这个困难，先用Euler格式(7.2.4)求得一个初步的近似值，称之为预测值；然后用公式(7.2.6)作一次迭代得，即将校正一次，得如下方法： 定义7.2.1 对，称预测校正系统(predictor-corrector

11、system) (7.2.8)为改进的Euler格式(modified Euler scheme), 称这种方法为求解初值问题(7.2.1)的改进的Euler方法(modified Euler method).注 改进的Euler格式(7.2.8)还可以表示为下列平均化形式： (7.2.9)例7.2.1 取步长，分别用Euler方法及改进的Euler方法求解初值问题解 这个初值问题的准确解为. 根据题设知Euler方法的计算格式为改进的Euler方法的计算格式为或初始值均为，将计算结果列于表7.2.1.表7.2.1Euler方法改进的Euler方法准确值0100011.10.27182820.

12、34237780.345919921.20.68475560.85831450.866642531.31.27697831.59274961.607215141.42.09354772.59829822.620359651.53.18744513.93644413.967666361.64.62081785.67890715.720961571.76.46639647.90920927.963873581.88.809119710.724467110.793624791.911.747996514.237441814.3230815102.015.398235718.578882518.683

13、0971从表7.2.1可以看出, Euler方法的计算结果与准确值差距较大, 最多只有1位有效数字；而改进的Euler方法却至少有2位有效数字, 这表明改进的Euler方法的精度比Euler方法高。7.2.3 局部截断误差与阶无论是Euler方法，还是改进的Euler方法，得到都是近似解，因此首先要考虑的是这些数值解的误差. 为使讨论简单起见，不妨假定在计算时用到前面一步的值是准确值，即.定义7.2.2 在的前提下，称误差为局部截断误差(local truncation error). 若一种数值方法的局部截断误差， (7.2.10)则称这种方法的精度为阶的(accuracy of pth o

14、rder).定理7.2.2 Euler方法是具有一阶精度的数值方法; 即其局部截断误差为.证 设，由(7.2.1)知：. 将在处作Taylor展开. (7.2.11)由Euler方法(7.2.4)得. (7.2.12)将式(7.2.11)与式(7.2.12)相减得. (7.2.13)于是, 即. 由定义8.2.1知：Euler方法是具有一阶精度的数值方法. 证毕！仿定理7.2.2的证明，并利用二元Taylor展开式得：定理7.2.3 改进的Euler方法是具有2阶精度的数值方法; 即其局部截断误差为.注 从实际情况来看，在计算时用到前面一步的值不可能是准确值，而应该是的前面每一步都有误差，这样

15、得到的与准确值的误差称为整体截断误差. 作为误差估计，最根本的是估计整体截断误差. 而估计整体截断误差并不容易，好在整体截断误差与局部截断误差之间通常有一定的联系(见定理8.4.4):整体截断误差局部截断误差 (7.2.14)即：单步法的整体截断误差比局部截断误差低一阶. 因此要构造高精度的计算方法，必须设法提高该方法的局部截断误差中的. 一般地，一个方法的整体截断误差阶越高，该方法的精度也越高. 由式(7.2.14)知：Euler方法的整体截断误差为, 改进的Euler方法的整体截断误差为.7.3 Runge-Kutta方法7.3.1 Runge-Kutta方法的基本思想对初值问题 (7.3

16、.1)由7.2节知：Euler方法每一步只计算一次的值，其整体截断误差为，它的计算格式可以写成 (7.3.2)而改进的Euler方法，每迭代一步要计算两个函数值，其整体截断误差为，它的计算格式可改写成 (7.3.3)上述2种情况启发我们，能否通过在每步迭代中多计算几次函数的值，然后对这些值作线性组合，使相应的数值方法的整体截断误差的阶数升高，从而提高该方法的精度, 而库塔(Martin Wilhelm Kutta 1867年11月3日-1944年12月25日) 是德国数学家。这正是龙格-库塔(Runge-Kutta)方法的基本思想. 由7.2节可知：因为整体截断误差与局部截断误差的关系，所以要

17、构造高精度的计算方法，只要设法提高该方法的局部截断误差即可；因此下面利用局部截断误差讨论Runge-Kutta方法。7.3.2 二阶Runge-Kutta方法对初值问题(7.3.1)，以计算两次函数值为例，设一般计算格式为 (7.3.4)其中。适当选择参数的值，使得在的假设下，局部截断误差的阶尽可能高. 为此, 将在点作Taylor展开，因为 的展开式为考虑到，所以 (7.3.5)再将在作Taylor展开 (7.3.6)若要使局部截断误差，则通过比较式(8.3.5)和(8.3.6)知，参数必须满足, (7.3.7)这是3个未知量，2个方程式的方程组，因此解不惟一. 定义7.3.1 满足条件(7

18、.3.7)的数值方法(7.3.4)称为二阶Runge-Kutta方法 (Runge-Kutta method of order two) (简记为二阶R-K方法)。注 二阶Runge-Kutta方法的局部截断误差都是.(1) 若取, 则式(7.3.4)化为 (7.3.8)这就是改进的Euler方法(7.2.8).(2) 若取, 则式(8.3.4)化为 (7.3.9)并称之为中点方法(Midpoint method).(3) 若取, 则式(7.3.4)化为 (7.3.10)并称之为Heun方法(Heun method).例7.3.1 取步长，分别用改进的Euler方法、中点方法及Heun方法求解

19、初值问题解 根据题设知，改进的Euler方法的计算格式为中点方法的计算格式为Heun方法的计算格式为上述计算格式的初始值均为，将计算结果列于表7.3.1.表7.3.1改进的Euler方法中点方法Heun方法准确值01000011.10.34237780.34094440.34139210.345919921.20.85831450.85490960.85597430.866642531.31.59274961.58670781.58860011.607215141.42.59829822.58880772.59178592.620359651.53.93644413.92252503.9269

20、0203.967666361.65.67890715.65938695.66553815.720961571.77.90920927.88269477.89106767.963873581.810.724467110.689312410.700436710.793624791.914.237441814.191712614.206212214.3230815102.018.578882518.520314618.538920818.6830971从表7.3.1可以看出, 三种R-K方法均至少有2位有效数字。7.3.3 三阶及四阶Runge-Kutta方法为了提高方法的精度,考虑每步计算三次函数

21、值. 根据两次计算函数值的做法，很自然地取的形式为 (7.3.11)其中。适当选择参数, , , , , , ，使截断误差的阶尽可能高. 由于在推导公式时只考虑局部截断误差，故设. 类似前面二阶R-K公式的推导方法，将在处作Taylor展开，然后再将在处作Taylor展开，只要两个展开式的前四项相同，便有，而要两个展开式的前四项相同，参数必须满足 (7.3.12)这是7个未知数5个方程的方程组，解不是惟一的，可以得到很多公式. 定义7.3.2 满足条件(7.3.12)的数值方法(7.3.11)称为三阶Runge-Kutta方法 (Runge-Kutta method of order two)

22、 (简记为三阶R-K方法)。注 三阶Runge-Kutta方法的局部截断误差都是.若在(7.3.12)中取则得到一种常用的三阶Runge-Kutta公式 (Runger-Kutta formula of order three) (7.3.13)类似地，如果每步计算4次函数的值，那么可以导出局部截断误差为的四阶Runge-Kutta公式. 详细推导过程这里略去，本节只给出一种最常用的四阶R-K公式： (7.3.14)公式(7.3.14)称为四阶经典Runge-Kutta方法(classical Runge-Kutta method of order four ).三阶、四阶Runge-Kutt

23、a公式都是显式单步公式. 例7.3.2 取步长，分别用三阶Runge-Kutta方法(7.3.13)和四阶经典Runge-Kutta方法(7.3.14)求解初值问题解 根据题设知，三阶Runge-Kutta方法(8.3.13)的计算格式为四阶经典R-K方法(7.3.14)的计算格式为初始值均为，将计算结果列于表7.3.2.表7.3.2三阶R-K方法四阶经典R-K方法准确值0100011.10.34561110.34591030.345919921.20.8594490.86662170.866642531.31.60604471.60718131.607215141.42.61862802.6

24、2031132.620359651.53.96528053.96760193.967666361.65.71782345.72087935.720961571.77.95987977.96377187.963873581.810.788666110.793501810.793624791.914.317042914.322935714.3230815102.018.675856818.682926618.6830971从表7.3.2可以看出, 三阶R-K方法至少有3位有效数字；四阶经典R-K方法至少有5位有效数字。例7.3.3 试用Euler方法、改进的Euler方法及四阶经典Runge-Ku

25、tta方法在不同步长下计算初值问题在1.2, 1.4, 1.8, 2.0处的近似值，并比较它们的数值结果。解 对上述三种方法，每执行一步所需计算的次数分别为1, 2, 4. 为了公平起见，上述三种方法的步长之比应为. 因此，在用Euler方法、改进的Euler方法及四阶经典R-K方法计算1.2, 1.4, 1.8, 2.0处的近似值时，它们的步长应分别取为0.05, 0.1, 0.2，以使三种方法的计算量大致相等.Euler方法的计算格式为改进的Eluer方法的计算格式为四阶经典R-K方法的计算格式为上述计算格式的初始值均为，将计算结果列于表7.3.3.表7.3.3Euler方法(步长h=0.

26、05)改进的Euler方法(步长h=0.1)四阶经典R-K方法(步长h=0.2)准确解1.20.76969600.85831450.86637910.86664251.42.34022362.59829822.61974052.62035961.65.13737865.67890715.71989535.72096151.89.743489410.724467110.792017610.79362472.016.949013318.578882518.680852418.6830971从表7.3.3可以看出，在计算量大致相等的情况下，Euler方法计算的结果只有1位有效数字，改进的Euler方

27、法计算的结果有2位有效数字，而四阶经典R-K方法计算的结果却有4位有效数字，这与理论分析是一致的. 例7.2.1和例7.3.3的计算结果说明，在解决实际问题时，选择恰当的算法是非常必要的。需要指出的是Runge-Kutta方法是基于Taylor展开法，因而要求方程的解具有足够的光滑性. 如果解的光滑性差，使用四阶Runge-Kutta方法求得数值解的精度，可能不如改进的Euler方法精度高. 因此，在实际计算时，要根据具体问题的特性，选择合适的算法。7.3.4 变步长的Runge-Kutta方法上面导出的Runge-Kutta方法都是定步长的，单从每一步来看，步长h越小，局部截断误差也越少.

28、但随着步长的减小，在一定范围内要进行的步数就会增加，而步数增加不仅增加计算量，还有可能导致舍入误差的积累过大. 由于Runge-Kutta方法是单步法，每一步计算步长都是独立的，所以步长的选择具有较大的灵活性. 因此根据实际问题的具体情况合理选择每一步的步长是非常有意义的. 正如第7章建立变步长的Simpson求积公式那样，我们来建立变步长的Runge-Kutaa公式.以四阶经典Runge-Kutta公式为例进行说明. 从初始点出发，先选一个步长，利用公式(7.3.14)求出的近似值记为, 由于公式的局部截断误差为, 所以有, (7.3.15)其中为常数. 然后将步长进行折半, 即取步长为,

29、由基点出发, 到算一步, 由到再算一步, 将求得的近似值记为, 因此有 (7.3.16)由式(7.3.15)和(7.3.16)得一般地, 从出发, 按上述做法可得到或写成 (7.3.17)式(7.3.17)是事后估计式, 记. (7.3.18)对给定的步长h，若是预先指定的精度)，说明步长太大，应折半进行计算，直至为止，这时取作为近似值. 若，则将步长加倍，直到为止，这时再将步长折半一次，就把这次所得结果作为近似值. 变步长的Runge-Kutta方法的计算步骤如下：设误差上限为，误差最小下限为，步长最大值为，从出发进行计算, 步长为.(1) 用步长和Runge-Kutta公式计算，用步长计算

30、两步得，并计算；(2) 若，说明步长过大，应将折半，返回步骤(1)重新计算；(3) 若，说明步长过小，在下一步将放大，但不超过.这种通过加倍和减半处理步长的Runge-Kutta方法就称为变步长Runge-Kutta方法(Runge-Kutta varing step-size method). 从表面上看，为了选择步长，每一步的计算量增加了，但从总体考虑还是合算的.7.4 单步法的相容性、收敛性与稳定性常微分方程的初值问题 (7.4.1)的数值解法的基本思想是通过某种离散化手段将微分方程转化为差分方程，很自然地，我们必须回答以下三个问题：1、这种转化是否合理？即，相容性问题。2、由数值解法得

31、到的数值解是否收敛到值问题(7.4.1)的解析解？即，收敛性问题。3、数值解法是否稳定？即，稳定性问题。本节仅介绍单步法的相容性、收敛性和稳定性，而多步法的的相容性、收敛性和稳定性的证明比较复杂，这里就不介绍了，有兴趣的读者请参阅文献1.下面均设，其中为等距节点，步长.7.5.1 相容性定义7.4.1 若解初值问题(7.4.1)数值方法的计算格式可以写成 (7.4.2)则称该方法为显式单步法(explicit one-step method). 其中称为增量函数(increment function).注 增量函数是依赖于微分方程中的的一个泛函. 例如：Euler方法的增量函数是；改进的Eul

32、er方法的增量函数是.定义7.4.2 设求解处初值问题(7.4.1)的单步法为(7.4.2)，若关于连续，且对任意固定的，都有 (7.4.3)则称数值方法(7.4.2)与初值问题(7.4.1)是相容；或数值方法(7.4.2)与初值问题(7.4.1)是相容的(consistent)。因为，所以有， (7.4.4)式(7.4.4)称为相容性条件(consistent condition)。注1 若关于连续，则式(7.4.3)与式(7.4.4)式等价的；因此，若数值方法(7.4.2)的增量函数满足相容性条件(7.4.4)，则称数值方法(7.4.2)与初值问题(7.4.1)相容。注2 相容性保证了，当

33、时，数值方法的差分方程收敛于原初值问题的微分方程，说明数值解法通过某种离散化手段将微分方程转化为差分方程求解原初值问题是合理的。注3 容易验证，Euler方法和改进的Euler方法都与初值问题(8.4.1)相容。7.5.2 收敛性定义7.4.3 对单步法(7.4.2)，若对任意固定的，都有，则称方法(7.4.2)是收敛的(convergent).定义7.4.4 若在计算(除外)时, 用到的不是，而是近似值；即: 每一步计算除局部截断误差外，还考虑前一步不准确而引起的误差，则称这种误差为整体截断误差(global truncation error). 整体截断误差包含各步的局部误差，也包括了以前

34、各步的局部误差在逐步计算中的积累. 因此作为误差估计，最根本的是估计整体截断误差. 但估计整体截断误差并不容易，好在整体截断误差与局部截断误差之间通常有一定的联系(见定理7.4.4), 因此要构造高精度的计算方法，只要设法提高该方法的局部截断误差中的阶数即可.定理7.4.1 若 (a) 单步法(7.4.2)具有阶精度；即: 它的局部截断误差为， (7.4.5)其中.(b) 增量函数在区域内关于满足Lipschitz条件；即：存在，对于，有. (7.4.6)(c) 初始值是准确的：,则(1) 当时，方法(7.4.2)是收敛的.(2) 方法(7.4.2)的整体截断误差为， (7.4.7)其中.证

35、由条件(a)知，存在常数，使得,且由条件(b)得因此即有递推关系. (7.4.8)利用递推关系(7.4.8)反复递推，得注意到，当(是某个固定的常数)时，于是有 (7.4.9)因此当时，由式(7.4.9)知，方法(7.4.2)是收敛的，且(7.4.7)成立. 证毕！注1 在定理7.4.1的条件下，单步法(7.4.2)的整体截断误差比局部截断误差低一阶；即：整体截断误差局部截断误差.注2 要判断某个单步法是否收敛，只要判断增量函数是否关于满足Lipschitz条件.推论1 Euler方法是收敛的.证明 因为Euler方法满足；且由7.1节的假定知：Euler方法的增量函数关于满足Lipschit

36、z条件，故由定理7.4.4知：Euler方法是收敛的. 证毕！推论2 改进的Euler方法也是收敛的.证明 因为改进的Euler方法的增量函数是，由7.1节的假定知：存在，对于，有；于是其中(是某个固定的常数), 所以改进的Euler方法的增量函数关于满足Lipschitz条件.又因为改进的Euler方法满足；所以由定理7.4.1知：改进的Euler方法也是收敛的. 证毕！注 相容性保证了当步长时，数值方法的差分方程收敛到初值问题的微分方程；收敛性保证了当步长时，数值方法的数值解收敛到初值问题的解析解。 下面的定理阐述了单步法的收敛性与相容性的关系：定理7.4.2 若单步法(7.4.2)的增量

37、函数在区域上关于满足Lipschitz条件；即：存在，对于，有，则当精度时，方法(7.4.2)收敛的充分必要条件是相容性条件(7.4.4)成立。证明参见文献13。7.4.3 稳定性在理论上，讨论方法的收敛性时，总是假定初始值是准确的，数值方法本身的计算也是准确的，而实际情况并非如此. 事实上，初始数据可能有误差，计算过程中的数字的舍入也会产生误差. 这种误差的扰动在计算过程中是否会增长很快，以至影响结果，这就是数值方法的稳定性问题. 定义7.4.5 若一个数值方法在节点处的值有扰动，在以后的各点处的值产生的偏差均不超过，则称这种方法是稳定的(stable).在选择数值方法时，我们总是希望选择稳

38、定的数值方法.定理7.4.3 对初值问题(7.4.1)，若增量函数是连续的，在集合上关于满足Lipschitz条件；即：存在，对于，有， (7.4.10)其中是常数, 则方法(7.4.2)是稳定的.当考察某数值方法的稳定性时，为简明起见，通常将此方法用于简单的试验方程以检验该方法的稳定性. 试验(test equation)方程为 (7.4.11)其中是常数，可以是实数，也可以是复数. 定义7.4.6 将单步法(7.4.2)用于试验方程(7.4.11)，设得到的迭代格式为，其中步长是固定的. 若增量函数满足，则称此方法是绝对稳定的(absolute stable)。记，则在平面上，使的变量围成

39、的区域成为绝对稳定域(region of absolute stable)，它与实轴的交集称为绝对稳定区间(interval of absolute stable).下面分别将Euler方法，后退的Euler方法，改进的Euler方法和梯形方法分别应用于试验方程(7.4.11)，并记，确定上述方法的绝对稳定区域.(1) 将Euler方法应用于试验方程(7.4.9)，得, (7.4.12)其中. 若有误差而变为，相应地有误差而变为，则. (7.4.13)将式(7.4.12)与式(7.4.13)相减, 并记，得因此. (7.4.14)由此可见，若要求误差不增长，则比值必须满足，即 . (7.4.1

40、5)因为可以是复数，所以在的复平面上，式(7.4.15)表示Euler方法的绝对稳定区域是半径为1，圆心为的圆域内部，不包括圆周(如图7.4.1). 图7.4.1绝对稳定区间是(此时显然要求)，即. 这说明：当步长满足(其中)时， Euler方法是稳定的.因为稳定性对步长有限制，因此称这样的稳定为条件稳定(conditional stable).(2) 将后退的Euler方法应用于试验方程(7.4.11)，并记, 得 即 , (7.4.16)其中. 若有误差而变为，相应地有误差而变为，则. (7.4.17)将式(7.4.16)与式(7.4.17)相减, 并记，得,由此得. (7.4.18)由此

41、可见，若要求误差不增长，则比值必须满足 或 , (7.4.19)即. (7.4.20)因此在复平面上，式(7.4.20)表示后退的Euler方法的绝对稳定区域是半径为1，圆心为的圆域的外部，不包括圆周(如图7.4.2).图7.4.2当时，在复平面上，式(7.4.20)表示后退的Euler方法的绝对稳定区域是左半平面，绝对稳定区间是，即. 这说明后退的Euler方法对任意步长都是稳定的，此时称后退的Euler方法是无条件稳定的(unconditional stable).(3) 将改进的Euler方法(一种二阶Runge-Kutta方法)应用于试验方程(7.4.11)，得, (7.4.21)其中

42、. 仿上述过程，得.若要求误差不增长，则比值必须满足, 即 或 . (7.4.22)于是改进的Euler方法的绝对稳定区域由确定，由得改进的Euler方法的绝对稳定区间(此时显然要求)，即. 这说明：当步长满足(其中)时，改进的Euler方法是稳定的，且是条件稳定的.(4) 将梯形方法应用于试验方程(7.4.11)，得, (7.4.23)其中. 仿上述过程，得.若要求误差不增长，则比值必须满足 即 . (7.4.24)当时，在复平面上，式(7.4.24)表示隐式梯形方法的绝对稳定区域是左半平面，绝对稳定区间是，即. 这说明隐式梯形方法对任意步长都是稳定的，即隐式梯形方法也是无条件稳定的。从前面

43、的讨论可以看出，隐式方法的稳定区域(如后退Euler方法和隐式梯形方法)分别比显式方法(如Euler方法和改进的Euler方法)的稳定区域大，而且在条件下后退Euler方法和隐式梯形方法都是无条件稳定的，而Euler方法和改进的Euler方法却是条件稳定的。一般地，显式Runge-Kutta方法皆具有有限的绝对稳定区域，即对步长都有限制，因此它们都是条件稳定的。而隐式方法在某些情况下是无条件稳定的。 7.5 线性多步法7.3节介绍的Runge-Kutta法是单步法，即：计算时，只用到前面一步的信息. 为了提高精度，三阶及三阶以上的Runge-Kutta法用的都是内部的点处的函数值信息，这样做的

44、不足是：计算量增加；而且这些新增加的函数值信息不保留到以后步骤的运算，在算法上似乎有点浪费. 这就启发我们：如果能充分利用已知信息, , 来计算，那么不但有可能提高精度，而且大大减少了计算量；这就是构造所谓线性多步法的基本思想. 又因为误差随着的增大而趋于增大，所以当求时，使用这些以前更精确的数据信息的方法也是合理的.定义7.5.1 求解初值问题 (7.5.1)的k步线性多步法 (linear k-step multistep method) 的一般公式是 (7.5.2)其中和是与无关的常数, ，.当时，式(7.5.2)中不含，称方法(7.5.2)为显式的 (explicit 或 open)；

45、当时，式(7.5.2)中含，称方法(7.5.2)为隐式的 (implicit 或 closed).为清晰简明起见，本节主要介绍最常用的线性多步法四阶Adams方法；本节最后还简介了其他一些较常用的线性多步法。7.5.1 Adams方法定义7.5.2 对, 称 (7.5.3)为四阶显式Adams方法 (explicit Adams method of order four).对, 称 (7.5.4)为四阶隐式Adams方法 (implicit Adams method of order four). 推导式(7.5.3)和式(7.5.4)：将区间等分,其中为等距节点，步长. 阿达姆斯(John

46、Couch Adams 1819年6月5日 1892年1月21日) 英国数学家和天文学家. 他最著名的成就是仅用数学理论和方法准确地预测了海王星的存在和位置。在每个小区间上，对两边积分，得. (7.5.5)对运用插值型求积公式，便可得到相应的线性多步法公式.将式(7.5.5)右边积分的被积函数用插值多项式来近似代替. 这里，我们取三次插值多项式. 插值节点除,外，通常还要在内再取两点作为插值节点. 但区间内的点处的函数值是未知的, 而却是已知的，因此可取插值节点, 则三次插值多项式为插值余项为 (7.5.6)其中. 于是有 (7.5.7)将式(7.5.5)右端的积分中的用代替，令，然后积分，并

47、用的近似值替代，得即得式(7.5.3).若插值节点取为作三次插值多项式其中.插值余项为, (7.5.8)其中. 于是有 (7.5.9)将式(7.5.5)右端的积分中的用代替，令，然后积分，再分别将, , , 用近似值, , , 表示，得即得式(7.5.4).定理7.5.1 若，则四阶显式Adams方法(7.5.3)的局部截断误差为. (7.5.10)四阶隐式Adams方法(7.5.4)的局部截断误差为 (7.5.11)证明 由余项表示式(7.5.6)得，四阶显式Adams方法的局部截断误差为 , (7.5.12)其中.因为，所以在上连续, 亦在上连续. 由介值定理知：在上存在最大值(记为)和最

48、小值(记为)，于是有. 记,则有令，上式化为即 或 .由介值定理知：存在，使得,即.代入式(7.5.12)，得 (7.5.13)注 这里不能用积分第二中值定理证明，因为未必是连续的.由余项表示式(7.5.8)得，四阶隐式Adams方法的局部截断误差为由介值定理，仿式(7.5.13)的证明，得证毕！7.4.2 初始值的计算四阶显式Adam方法 (7.5.3)计算时需要知道, , , , 而微分方程(7.5.1)只提供了一个初始值, 还有3个初值需要用其他方法来计算.四阶隐式Adam方法(7.5.4)计算时需要知道, , , 而微分方程(7.5.1)只提供了一个初始值，还有2个初值需要用其他方法来

49、计算. 因此，线性多步法不是能自行启动的方法，它需要其他方法提供初始值. 因为微分方程的初始值在定解时起着重要作用，所以补充的初始值对精度要求较高. 对于四阶Adams公式（显式的或隐式的），其局部截断误差为，因此用于确定补充初始值的其他方法的局部截断误差至少应不低于，否则由于初始值不准确，会影响到最后的结果. 对四阶Adams公式，最常用的补充初始值的方法是用四阶Runge-Kutta方法，计算四阶Adams方法中所缺少的初值. 还可以用Taylor展开法，如将在处作Taylor展开取它的前若干项，使余项满足精度要求.7.5.3 Admas预测-校正方法四阶显式Adams方法，在补充缺少的初

50、值后就可计算了. 每步只算一个函数值，局部截断误差可达到，这是单步法难以达到的. 而四阶隐式Adams方法，每步需要解一个非线性方程. 与四阶显式Adams方法相比，在计算上增加了许多困难. 但是同阶的隐式Adams方法的精度及稳定性比同阶的显式Adams方法好得多.在实际应用中，为了保留隐式Adams方法的优点，一般将显式公式与隐式公式联合使用，用显式公式提供预测值，用隐式公式加以校正; 这样，既不用求解非线性方程，又能达到较高的精度. 这种方法称为预测-校正方法. 定义7.5.4 对，称 (7.5.14)为四阶Adams预测-校正迭代系统 (forth-order Adams predictor-corrector iterative system).式(