1、第第5章章 常微分方程数值解法常微分方程数值解法5.1.2 局部误差和方法的阶局部误差和方法的阶5.1.1 Euler 方法及其有关的方法方法及其有关的方法5.1 Euler 方法方法5.1.1 Euler 方法及其有关的方法方法及其有关的方法考虑一阶常微分方程初值的问题：考虑一阶常微分方程初值的问题：00)(),(yxyyxfy设设f（x,y）是连续函数，对）是连续函数，对y满足满足Lipschitz条件，这样初值问题的解是存条件，这样初值问题的解是存在唯一的，而且连续依赖于初始条件。在唯一的，而且连续依赖于初始条件。为了求得离散点上的函数值，将微分方程的连续问题进行离散化。为了求得离散点上
2、的函数值，将微分方程的连续问题进行离散化。一般是引入点列一般是引入点列 ，这里为步长，这里为步长，经常考虑定长的情形，即经常考虑定长的情形，即 。记记 为初始问题的问题准确解为初始问题的问题准确解 在在 处的值，用均差近似代处的值，用均差近似代替的导数得替的导数得 nxnnnnhnhxx。称称,.2,1,1 ，100 nnhxxhhnn)(xy)(nxynx。，)()()()()()(11nnnnnnnnxyxfhxyhxyxyxfhxyhxy令令 为为 的近似值，将上面两个近似写成等式，整理后得的近似值，将上面两个近似写成等式，整理后得ny)(nxy。，10)(10)(1111nyxhfyy
3、nyxhfyynnnnnnnn（5.1.2）（5.1.3）从从 处的初值处的初值 开始，按（开始，按（5.1.2）可逐步计算以后各点上的值。称）可逐步计算以后各点上的值。称（5.1.2）式为）式为显式显式Euler。由于（。由于（5.1.3）式的右端隐含有待求函数值）式的右端隐含有待求函数值 ，不能逐步显式计算，称（不能逐步显式计算，称（5.1.3）式为）式为隐式隐式Euler公式公式或或后退后退Euler公式公式。如果。如果将（将（5.1.2）和（）和（5.1.3）两式作算术平均，就得）两式作算术平均，就得梯形公式。梯形公式。0 x0y1ny梯形公式也是隐式公式。以上公式都是由梯形公式也是隐
4、式公式。以上公式都是由 去计算去计算 ，故称它们为单步法。，故称它们为单步法。例例5.1 取取h=0.1，用，用Euler方法、隐式方法、隐式Euler方法和梯形方法解方法和梯形方法解ny1ny。，1)0(1yyxy 解解 本题有本题有 如果用如果用Euler方法，由（方法，由（5.1.2）并代入并代入h=0.1得得 。，11)(0yyxyxf。1.09.01.01nnnyxy同理，用隐式同理，用隐式Euler方法有方法有。)1.01.0(1.1111nnnyxy。，1)()(2111onyxfyxfhyynnnnnn（5.1.4）用梯形公式有用梯形公式有。)105.095.01.0(05.1
5、11nnnyxy三种方法及准确解三种方法及准确解 的数值结果如表的数值结果如表5-1所示。从表中看所示。从表中看 到，在到，在 处，处，Euler方法和隐式方法和隐式Euler方法的误差方法的误差 分分别是别是 和和 ，而梯形方法的误差却是，而梯形方法的误差却是 。xexxy)(5.0nxnnyxy)(2104.12106.14105.2 在例在例5.1中，由于中，由于f（x，y）对）对y是线性的，所以对隐式公式也可以方便地计是线性的，所以对隐式公式也可以方便地计算算 。但是，当。但是，当f（x，y）是）是y的非线性函数时，如的非线性函数时，如 ，其隐式，其隐式Euler公式为公式为 。显然，
6、它是。显然，它是 的非线性方程，可以选择的非线性方程，可以选择非线性方程求根的迭代求解非线性方程求根的迭代求解 。以梯形公式为例，可用显式。以梯形公式为例，可用显式Euler公式提供公式提供迭代初值迭代初值 ，用公式，用公式1ny35yxy)5(3111nnnnyxhyy1ny)0(1ny1ny表表5-1 Euler方法方法 隐式隐式Euler方法方法 梯形法梯形法 准确解准确解nx 0 1 1 1 10.1 1.000000 1.009091 1.004762 1.0048370.2 1.010000 1.026446 1.018549 1.018731 0.3 1.029000 1.051
7、315 1.040633 1.040818 0.4 1.056100 1.083013 1.070096 1.070320 0.5 1.090490 1.120921 1.106278 1.106531，102)()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn反复迭式，直到反复迭式，直到，)(1)1(1knknyy其中，步长其中，步长h成为迭代参数，它需要满足一定的条件，才能收敛。若将成为迭代参数，它需要满足一定的条件，才能收敛。若将（5.1.4）式减去该迭代公式，得）式减去该迭代公式，得)(1111)1(112knnnnknnyxfyxfhyy，假设假设f（
8、x，y）关于）关于y满足满足Lipschiz条件，则有条件，则有，)(11)1(112knnknnyyhLyy这里，这里，L是是Lipschiz常数。当常数。当hL/21即即h2/L时，迭代序列时，迭代序列 收敛收敛 。)(1kny 1 ny 对于隐式公式，通常采用估计对于隐式公式，通常采用估计-校正技术，即先用显式公式计算，得到校正技术，即先用显式公式计算，得到预估值，然后以预估值作为隐式公式的迭代初值，用隐式公式迭代一次得到预估值，然后以预估值作为隐式公式的迭代初值，用隐式公式迭代一次得到校正值，称为校正值，称为预估预估-校正技术校正技术。例如，用显式。例如，用显式Euler公式作预估，用
9、梯形公式公式作预估，用梯形公式作校正，即作校正，即。，1021111nyxfyxfhyyyxhfyynnnnnnnnnn称该公式为称该公式为改进的改进的Euler公式公式。它显然等价于显式公式为。它显然等价于显式公式为nnnnnnnnyxhfyxfyxfhyy,211，（5.1.6）也可以表示为下列平均化的形式也可以表示为下列平均化的形式。，qpnpnnqnnnpyyyyxhfyyyxhfyy2111例例5.2 取取h=0.1，用改进的，用改进的Euler方法解方法解。，102yyxyy解解 按（按（5.1.5），改进的），改进的Euler方法解方法解。，10)2()2(2)2(11111ny
10、xyyxyhyyyxyhyynnnnnnnnnnnnn由由 得计算结果如表得计算结果如表5-2。该初值问题的准确解为。该初值问题的准确解为 。1.010hy，xxy21表表 5-2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6153 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6165 nxny)(nxy5.1.2 局部误差和方法的阶局部误差和方法的阶 初值问题（初值问题（5.1.1）的单步法可以写成如下统一形式）
11、的单步法可以写成如下统一形式，)(111hyyxxhyynnnnnn（5.1.7）f其中其中 与与 有关。若有关。若 中不含中不含 ，则方法是显式的，否则是隐式的，所以一，则方法是显式的，否则是隐式的，所以一般显式单步法表示为般显式单步法表示为1ny。，hyxhyynnnn1（5.1.8）例如，例如，Euler方法中，有方法中，有 yxfhyx，对于不同的方法，计算值对于不同的方法，计算值 与准确解与准确解 的误差各不相同。所以有必要讨论方的误差各不相同。所以有必要讨论方法的截断误差。我们称法的截断误差。我们称 为某一方法在为某一方法在 点的点的整体截断误差整体截断误差。显。显然，然，不单与不
12、单与 这步的计算有关，它与以前各步的计算也有关，所以误差被称这步的计算有关，它与以前各步的计算也有关，所以误差被称为整体的。分析和估计整体截断误差为整体的。分析和估计整体截断误差 是复杂的。为此，我们假设是复杂的。为此，我们假设 处的处的 没有误差，即没有误差，即 ，考虑从，考虑从 到到 这一步的误差，这就是如下的局部这一步的误差，这就是如下的局部误差的概念。误差的概念。ny nxynnnyxyenxnenxnenxnynnxyy nx1nx hxyxyxxhxyxyTnnnnnnn，1111定义定义5.1 设设 是初值问题（是初值问题（5.1.1）的准确解，则称）的准确解，则称 xy为单步法
13、（为单步法（5.1.7）的）的局部截断误差局部截断误差。定义定义5.2 如果给定方法的如果给定方法的 局部截断误差局部截断误差 ，其中，其中 为整数，则称该方法是为整数，则称该方法是p阶阶的，或具有的，或具有p阶精度阶精度。若一个。若一个p阶单步法的局部阶单步法的局部截断误差为截断误差为)(11pnhOT1p，211ppnnnhohxyxgT则称其第一个非零项则称其第一个非零项 为该方法的为该方法的局部截断误差的主项局部截断误差的主项。对于对于Euler方法，有方法，有Taylor展开有展开有1pnnhxyxg nnnnnxyxhfxyxyT，11)(1nnnxyhxyxy 2432)(62h
14、ohoxyhxyhnn 对于隐式对于隐式Euler方法，其局部截断误差为方法，其局部截断误差为所以所以Euler方法是一种一阶方法，其局部截断误差的主项为方法是一种一阶方法，其局部截断误差的主项为 。nxyh 221111nnnnnxyxhfxyxyT，)(11nnnxyhxyxy2322hOhOxyhn 梯形方法也是一种隐式单步法，类似可得其局部截断误差梯形方法也是一种隐式单步法，类似可得其局部截断误差所以隐式所以隐式Euler方法也是一种一阶方法，该方法的局部截断误差的主项方法也是一种一阶方法，该方法的局部截断误差的主项为为 ，仅与显式，仅与显式Euler方法的局部截断误差的主项反一个符号
15、。方法的局部截断误差的主项反一个符号。nxyh 22 11112nnnnnnnxyxfxyxfhxyxyT，。34312hOhOxyhn 可见，梯形方法是二阶精度的。可见，梯形方法是二阶精度的。5.2 Runge-Kutta 方法方法5.2.2 几类显式几类显式Runge-Kutta方法方法5.2.1 Runge-Kutta方法的基本思想方法的基本思想5.2.1 Runge-Kutta方法的基本思想方法的基本思想 显式显式Euler方法是最简单的单步法，它是一阶的，它可以看作方法是最简单的单步法，它是一阶的，它可以看作Talylor展开后展开后取前两项。因此，得到高阶方法的一个直接想法是用取前
16、两项。因此，得到高阶方法的一个直接想法是用Talylor展开，如果能计算展开，如果能计算 的高阶导数，则可写出的高阶导数，则可写出p阶方法的计算方法阶方法的计算方法 xy，pnpnnnnyphyhyhyy!221 其中其中 是是 的近似值，的近似值，若将若将 分别记成分别记成 则对于二阶和三阶导数可表示为则对于二阶和三阶导数可表示为 jny njxy。，pj210，yfxfyxf，yxfff。，fffffffffyfffyyyyyxxyxxyx222 这个方法并不实用，因为一般情况下，求这个方法并不实用，因为一般情况下，求 的导数相当麻烦。从计算高的导数相当麻烦。从计算高阶导数的公式知道，方法
17、的截断误差提高一阶，需要增加的计算量很大。但阶导数的公式知道，方法的截断误差提高一阶，需要增加的计算量很大。但是由此启发我们用区间上若干个点的导数是由此启发我们用区间上若干个点的导数 ，而不是高阶导数，将它们作线，而不是高阶导数，将它们作线性组合得到平均斜率，将其与解的性组合得到平均斜率，将其与解的Taylor展开相比较，使前面若干项吻合，从展开相比较，使前面若干项吻合，从而得到具有一定阶的方法。这就是而得到具有一定阶的方法。这就是Runge-Kutta方法的基本思想，其一般形式方法的基本思想，其一般形式为为yxf，f，LiKahcyhcxfKyxfKKhyyijjijinininniLiin
18、n32,11111(5.2.1)其中，其中，与与 的区别在于：用微分方程准确解的区别在于：用微分方程准确解 代替代替 中的中的 就得到就得到 。参数参数 和和 待定，确定它们的原则和方法是：将（待定，确定它们的原则和方法是：将（5.2.2）式中）式中的的 在在 处作处作Taylor展开，将展开，将 在在 处作二元处作二元Taylor展开，将展开展开，将展开式按式按H的幂次整理后，令的幂次整理后，令 中中h的低次幂的系数为零，使的低次幂的系数为零，使 首项中首项中h的幂次尽的幂次尽量高，比如使量高，比如使 ，则称（，则称（5.2.1）式为）式为L级级p阶阶Runge-Kutta方法方法（简称（简
19、称R-K法法）。）。*iKiK nxyiKny*iKiic，ija 1 nxynx nnxyx，1 nT1 nT 11 pnhOT*iK其中，其中，。它的局部截断误差是。它的局部截断误差是111111 ijijLiiiac，*111iLiinnnKhxyxyT ，（5.2.2）它与显式它与显式R-K公式的区别在于：显式公式中，对系数公式的区别在于：显式公式中，对系数 求和的上限是求和的上限是 ，从，从而而 构成的矩阵是一个严格下三角阵。而在隐式公式中，对系数构成的矩阵是一个严格下三角阵。而在隐式公式中，对系数 求和的上求和的上限是限是L，从而，从而 构成的矩阵是方阵，需要用迭代法求出近似斜率构
20、成的矩阵是方阵，需要用迭代法求出近似斜率 推导隐式公式的思路和方法与显式公式法类似。推导隐式公式的思路和方法与显式公式法类似。ija1 iijaijaija LiKi，21 类似于显式类似于显式R-K公式（公式（5.2.1），稍加改变，就得到隐式），稍加改变，就得到隐式R-K公式公式，iLiinnKhyy11。，LiKahcyhcxfKLjjijinini2115.2.2 几类显式几类显式Runge-Kutta方法方法 对于对于L=2，则，则。，1222122111hKcyhcxfKyxfKKKhyynnnnnn其局部截断误差是其局部截断误差是 *22*1111KKhxyxyTnnn（5.2.
21、3）将将 中的各项作中的各项作Taylor展开，并利用展开，并利用 则有则有1nT ，fffyxyxfxyyxnnn，nnnnnnnnnnnxyhcxyhcxfKxyxyxfKhOxyhxyhxyhxyxy 22*2*1432162，32222 222hOfffffhcxhycxyyyxyxxnn将它们代入（将它们代入（5.2.3）式，整理后得）式，整理后得 nnnxyhcxyhT 222211211 4222232261hOfffffcxyhyyxyxxn 选取选取 和和 ，使方法的阶尽可能高，就是使，使方法的阶尽可能高，就是使 h 和和 的系数为零，因为的系数为零，因为 的系数一般不为零。
22、于是得到方程组的系数一般不为零。于是得到方程组21，2c2h3h。，2112221c显然，该方程组有无穷多组解，从而得到一族显然，该方程组有无穷多组解，从而得到一族二级二阶二级二阶R-K方法方法。2c若以若以 为自由参数，取为自由参数，取 得得中点公式中点公式212c，nnnnnnyxfhyhxhfyy221 （5.2.4）取取c=2/3得得Heun公式公式nnnnnnnnyxhfyhxfyxfhyy，3232341，（5.2.5）取取c=1得改进的得改进的Euler公式（公式（5.1.6）。）。对于对于L=3的情形，要计算三个斜率的近似值：的情形，要计算三个斜率的近似值：。，23213133
23、312221KaKahcyhcxfKhKcyhcxfKyxfKnnnnnn 类似于二阶方法的推导，可以得三阶的方法，所得系数应满足的方程组是类似于二阶方法的推导，可以得三阶的方法，所得系数应满足的方程组是。，161312111323132323233222332221321aaacccccca该方程组的解也是该方程组的解也是 不唯一的。常见的一种不唯一的。常见的一种三级三阶方法三级三阶方法是是 。，2131221321122246hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnn 对于对于L=4的情形，可进行类似推导。最常用的四级四阶方法是如下的情形，可进行类似推导。最常用的四
24、级四阶方法是如下经典经典R-K方法方法 。，3423121432112222226hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn （5.2.6）为了分析经典为了分析经典R-K公式的计算量和计算精度，将四阶经典公式的计算量和计算精度，将四阶经典R-K公式（公式（5.2.6）与一阶显式与一阶显式Euler公式（公式（5.1.2）及二阶改进的）及二阶改进的Euler公式相比较。一般说来，公公式相比较。一般说来，公式的级数越大，计算右端项式的级数越大，计算右端项 f 的次数越多，计算量越大。在同样步长的情况的次数越多，计算量越大。在同样步长的情况下，下，Euler
25、方法每步只计算一个函数值，而经典方法要计算方法每步只计算一个函数值，而经典方法要计算4个函数值。四阶个函数值。四阶R-K法的计算量差不多是改进的法的计算量差不多是改进的Euler公式的公式的2倍，是显式倍，是显式Euler公式的公式的4倍。下面倍。下面的例子中的例子中Euler方法用步长方法用步长 ，二阶改进的，二阶改进的Euler法用步长法用步长 ，而四阶经典公，而四阶经典公式用步长式用步长 。这样，从。这样，从 到到 三种方法都计算了三种方法都计算了4个函数制，计算个函数制，计算量大体相当。量大体相当。12h14h1hnx14hxn例例5.3 考虑初值问题考虑初值问题 。，001 yyy其
26、解析解为其解析解为 。分别用。分别用h=0.025的显式的显式Euler方法，方法，h=0.05改进改进Euler法和法和h=0.1的经典的经典R-K方法计算到方法计算到x=0.5。三种方法在。三种方法在x方方向每前进向每前进0.1都要计算都要计算4个右端函数值，计算量相当。计算结果列于个右端函数值，计算量相当。计算结果列于表表5-3。从计算结果看，在工作量大致相同的情况下，还是经典方法。从计算结果看，在工作量大致相同的情况下，还是经典方法比其他两种方法的结果好得多。在比其他两种方法的结果好得多。在x=0.5处，三种方法的误差分别处，三种方法的误差分别是是 。经典。经典R-K法对多数好条件问法
27、对多数好条件问题（题（，参见下节单步法的稳定性），能获得好的效果。，参见下节单步法的稳定性），能获得好的效果。xexy 1743108.2103.1108.3 ，0yf表表 5-3Euler法法 改进改进Euler法法 经典经典R-K法法 准确解准确解h=0.025 h=0.05 h=0.1 nx 0.1 0.096312 0.095123 0.09516250 0.09516258 0.2 0.183348 0.181193 0.18126910 0.18126925 0.3 0.262001 0.259085 0.25918158 0.25918178 0.4 0.333079 0.329
28、563 0.32967971 0.32967995 0.5 0.397312 0.393337 0.39346906 0.39346934 nxy 在微分方程数值解法的实际计算中，有个如何选择步长的问题。在微分方程数值解法的实际计算中，有个如何选择步长的问题。因为单从每一步看，步长越小，截断误差越小。但随着步长的缩小，因为单从每一步看，步长越小，截断误差越小。但随着步长的缩小，在一定求解范围内所要完成的步数就增加了。步数的增加不但引起计在一定求解范围内所要完成的步数就增加了。步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累。算量的增大，而且可能导致舍入误差的严重积累。在选择步长时
29、，我们需要衡量和检验计算结果的精度，并依据所获得的精度在选择步长时，我们需要衡量和检验计算结果的精度，并依据所获得的精度处理步长。下面以经典处理步长。下面以经典R-K方法为例进行说明。方法为例进行说明。从节点从节点 出发，先以出发，先以h为步长求出一个近似值为步长求出一个近似值 ，由于公式的局部截断，由于公式的局部截断误差为误差为 ，故有，故有nx hny1 5hOnx然后将步长折半，既然后将步长折半，既 h/2 为步长，从为步长，从 跨两步到跨两步到 ，再求得一个近似，再求得一个近似值值 ，每跨一步的截断误差约为，每跨一步的截断误差约为 ，因此有，因此有1nx 21hnh 52hc 511c
30、hyxyhnn521122hcyxyhnn比较上述二式，有比较上述二式，有。16111211hnnhnnyxyyxy由此易得下列事后估计式由此易得下列事后估计式。hnhnhnnyyyxy121211151这样，我们可以通过检查步长折半前后两次计算结果的偏差这样，我们可以通过检查步长折半前后两次计算结果的偏差 hnhnyy121来判定所选的步长是否合适。来判定所选的步长是否合适。具体地说，对于给定的精度具体地说，对于给定的精度 ，将按两种情况处理。如果，将按两种情况处理。如果 ，我，我们反复将步长折半进行计算，直到们反复将步长折半进行计算，直到 为止，这时取最终得到的为止，这时取最终得到的 作作
31、为结果。如果为结果。如果 ，我们反复将步长加倍，直到，我们反复将步长加倍，直到 为止，这时再将为止，这时再将前一次步长折半的结果作为所要的结果。这种通过加倍或折半处理步长的方法前一次步长折半的结果作为所要的结果。这种通过加倍或折半处理步长的方法称作称作变步长方法变步长方法。虽然为了选择步长，每一步的计算量有所增加，但总体考虑。虽然为了选择步长，每一步的计算量有所增加，但总体考虑是值得的。是值得的。21hny5.55.5 线性多步法线性多步法5.5.2 基于基于Taylor展开的方法展开的方法5.5.1 基于数值积分的方法基于数值积分的方法5.55.5 线性多步法线性多步法 常微分方程初值问题（
32、常微分方程初值问题（5.1.1）的数值解法中，除了）的数值解法中，除了Runge-Kutta型公式等单步法之外，还有另一种类型的解法，即某一步的型公式等单步法之外，还有另一种类型的解法，即某一步的公式不仅与前一步解的值有关，而且与前若干步解的值有关，利公式不仅与前一步解的值有关，而且与前若干步解的值有关，利用前面多步的信息预测下一步的值，这就是多步法的基本思想，用前面多步的信息预测下一步的值，这就是多步法的基本思想，可以期望获得较高的精度。构造多步法有多种途径，下面先讨论可以期望获得较高的精度。构造多步法有多种途径，下面先讨论基于数值积分的方法。基于数值积分的方法。5.5.1 基于数值积分的方
33、法基于数值积分的方法将（将（5.1.1）中的方程在区间）中的方程在区间 上积分，可以得到上积分，可以得到1nnxx，。，dxxyxfxyxynnxxnn 11（5.5.1）如推导如推导Newton-Cotes求积公式一样，用等距节点的插值多项式来替代被积函求积公式一样，用等距节点的插值多项式来替代被积函数，再对插值多项式积分，这样就得到一系列求积公式。数，再对插值多项式积分，这样就得到一系列求积公式。例如，用梯形方法计算积分项例如，用梯形方法计算积分项 ，1121 nnnnxxxyxfxyxfhdxxyxfnn代入（代入（5.5.1）式有）式有 。，1112 nnnnnnxyxfxyxfhxy
34、xy据此即可导出公式（据此即可导出公式（5.1.4）。）。一般地，设由一般地，设由 个数据点个数据点 构造插值多项式构造插值多项式 ，这里，这里，。运用插值运用插值公式有公式有 1 r rnrnnnnnfxfxfx ，11 xPr khxxyxffkkkk 0，。，rjkkknjnknjjrjjnrxxxxxlxlfxP00将（将（5.5.1）离散化即得下列计算公式）离散化即得下列计算公式，jnrjrjnnfhyy 01（5.5.2）其中其中 。，rjdtjkktdxxlhrjkkxxjrjnn1011001 由此可得（由此可得（5.5.2）中的系数，其具体数值见表）中的系数，其具体数值见表5
35、-6。公式（。公式（5.5.2）是一个）是一个r+1步的显式公式，称为步的显式公式，称为Adams显式公式显式公式。r=0时，即为时，即为Euler公式。公式。表表 5-6j 0 1 2 3 4 1 3 -1 23 -16 5 55 -59 37 -9 1901 -2774 2616 -1274 251 jjjjj4321072024122应用实例应用实例：考虑跳伞员的下落速度。考虑跳伞员的下落速度。自由落体运动可用牛顿第二定律描述：自由落体运动可用牛顿第二定律描述：F=ma。实验表明，空气阻力模型。实验表明，空气阻力模型为为 ，其中，其中 ，比例系数，比例系数 k 依赖于物体的大小、形状，空
36、气依赖于物体的大小、形状，空气的密度和粘度。跳伞员下落的速度可描述为下列模型：的密度和粘度。跳伞员下落的速度可描述为下列模型：pkvF 21 p，00 vmgvkdtdvmp负号表示下降。显然，当负号表示下降。显然，当 1 p 2 时，适合于数值方法求解。时，适合于数值方法求解。设设 k/m =1.5，g=32，先用中点法提供开始值，再用下列两步而阶方法，先用中点法提供开始值，再用下列两步而阶方法，21211iffhvviiii求其他需要计算的值。当求其他需要计算的值。当 p=1 时，取时，取 h=0.2 有有。，0497.219552.208292.206611.204371.201383.
37、207400.192088.195007.185564.172975.166187.143816.123920.94400.501514131211109876543210vvvvvvvvvvvvvvvv可见，三秒末跳伞员的末速度约有可见，三秒末跳伞员的末速度约有 21 。若将模型修改为若将模型修改为 p=1.1，取，取 h=0.2，则有计算结果：，则有计算结果：sec/ft。，1113.160912.160612.160165.169500.158508.157030.154830.151552.151466749411.138630.122565.118911.83216.50151413
38、1211109876543210vvvvvvvvvvvvvvvv可见三秒末跳伞员的末速度减慢了。计算结果如下图所示可见三秒末跳伞员的末速度减慢了。计算结果如下图所示00.511.522.53-25-20-15-10-50+表示表示 p=1时的解，时的解，*表示表示 p=1.1时的解时的解 在上述在上述Adams显式公式的推导中，选用了显式公式的推导中，选用了 作为插值作为插值节点。这样的插值多项式节点。这样的插值多项式 在求积区间在求积区间 上逼近上逼近 是一是一个外推结果。为了改善逼近效果，我们变外推为内推，即改用个外推结果。为了改善逼近效果，我们变外推为内推，即改用 为插值节点，用数据点为
39、插值节点，用数据点 构造插值构造插值多项式多项式 ，则有，则有 rnnnxxx ，1 xPr 1 nnxx，xyxf，11 rnnnxxx，1111 rnrnnnnnfxfxfx，xPr rjkkknjnknjjrjjnrxxxxxlxlfxP，。，011101于是我们有如下的计算公式于是我们有如下的计算公式，101jnrjrjnnfhyy（5.5.3）其中其中，rjdtjkktdxxlhrjkkxxjrjnn1010101 其具体数值见表其具体数值见表 5-7。公式（。公式（5.5.3）是隐式公式，称为）是隐式公式，称为 Adams隐式公式。隐式公式。r=0，1时分别为隐式时分别为隐式 Eu
40、ler公式和梯形公式。公式和梯形公式。表表5-7j 0 1 2 3 4 1 1 1 5 8 -1 9 19 -5 1 251 646 -264 106 -19 jjjjj4321072024122对于隐式公式（对于隐式公式（5.5.3），需要用迭代求解。确定），需要用迭代求解。确定 的迭代公式为的迭代公式为1ny ，101111011 sfyxfhyyrjjnrjsnnrnsn 迭代收敛条件为迭代收敛条件为 ，其中，其中 的的Lipschitz常数常数 利用插值多项式的余项，可以求出利用插值多项式的余项，可以求出 Adams方法的局部截断误差。当然方法的局部截断误差。当然也可以从得到的显式和隐
41、式也可以从得到的显式和隐式 Adama公式，有局部截断误差的定义来求出方公式，有局部截断误差的定义来求出方法的局部截断误差。表法的局部截断误差。表 5-8中列出了它们的局部截断误差的主项，有表中列出了它们的局部截断误差的主项，有表 5-8可以看出，可以看出，Adams隐式方法的局部截断误差要小。隐式方法的局部截断误差要小。10 Lhr yfL关关于于为为r 0 1 2 3表表 5-8Adams显式公式显式公式Adams隐式公式隐式公式 nxyh 221 nxyh 221 nxyh 3125 nxyh 3121 nxyh4483 nxyh44241 nxyh55720251 nxyh557201
42、9 5.5.2 基于基于Taylor展开的方法展开的方法 基于数值积分可以构造出一系列求解常微分方程的计算公式，下面介绍基于数值积分可以构造出一系列求解常微分方程的计算公式，下面介绍基于基于 Taylor 展开的待定系数法，它可灵活地构造出线性多步法。对固定的系展开的待定系数法，它可灵活地构造出线性多步法。对固定的系数，可以选取待定系数使线性多步法的阶尽可能高。还可以根据需要，确定数，可以选取待定系数使线性多步法的阶尽可能高。还可以根据需要，确定显式还是隐式。显式还是隐式。设构造如下具有设构造如下具有 p 阶精度的线性多步公式阶精度的线性多步公式 。rnrnnrnrnnnfffhyyyy 01
43、11101(8.4.4)当当 时，则（时，则（8.4.4）为显式多步式。当）为显式多步式。当 时，（时，（8.4.4）为隐式）为隐式多步式。它们的局部截断误差为多步式。它们的局部截断误差为01 01 ，rkrkknknkknknnxyxfhxyxyT0111 利用原微分方程，有利用原微分方程，有 。rkrkknkknknnxyhxyxyT0111 （5.5.5）现利用现利用Taylor展开定理，确定线性多步公式（展开定理，确定线性多步公式（5.5.4）中的待定参数）中的待定参数 ，使她达到使她达到 阶精度，即阶精度，即 。对（对（5.5.5）式的右端各项在）式的右端各项在 点处作点处作Tayl
44、or展开有展开有 kk ，p 11 pnhOTnx 。，11112110!1!1!pnppnjpjjknpnppnjpjjknhOxypkhxyjkhxyhOxypkhxyjkhxy将它们代入（将它们代入（5.5.5）式整理后得）式整理后得 )()()()1()(1!1)()(1!)1(21111111110 pnpkrkpkprkpnjkjrkkjrkpjnrkkhOxykpkphxykjkjhxy 1nT使使 的系数为零，得到关于的系数为零，得到关于 和和 的线性方程组的线性方程组 pnhhhxy，2kk 。，pjkjkkjrkkjrkrkk21111110 （5.5.6）而且得到线性多步
45、法的局部截断误差而且得到线性多步法的局部截断误差 1nT 。21111111!1 pnpkprkkprkphOxykpkph 下面我们构造几个著名的四阶线性多步公式，考虑下列形式的公式下面我们构造几个著名的四阶线性多步公式，考虑下列形式的公式 1ny ，3322110113322110 nnnnnnnnnfffffhyyyy （5.5.7）。6543153151511!5hOxykhTnkkkkn （5.5.8）由于由于 r=3，p=4，由（，由（5.5.6）得到）得到5个方程，而（个方程，而（5.5.7）中有）中有9个为知量，个为知量，因此，（因此，（5.5.7）中有）中有4个自由度。个自由
46、度。若取若取 ，由（，由（5.5.6）式得到其他）式得到其他5个待定参数的个待定参数的方程组，解之得方程组，解之得003211 ，。，249243724592455132100 代入（代入（5.5.7）和（）和（5.5.8）式，得到常用的）式，得到常用的四步四阶显式四步四阶显式Admas公式公式和它的余和它的余项：项：。，655132131720251937595524hOxyhTffffhyynnnnnnnn （5.5.9）（5.5.10）若若 取取 ，由（，由（5.5.6）式得到其他）式得到其他5个待定参数的个待定参数的方程组，解之得方程组，解之得002101 ，。，03834381321
47、03 由此构造成著名的由此构造成著名的四步四阶显式四步四阶显式Milne公式公式和它的余项和它的余项 。，655121314514)22(34hOxyhTfffhyynnnnnnn （5.5.11）（5.5.12）若取若取 由（由（8.4.6）式得到其他）式得到其他5个待定参数的方程个待定参数的方程组，从而得组，从而得三步四阶隐式三步四阶隐式Admas公式公式及余项：及余项：003321，。，6551211172019)5199(24hOxyhTffffhyynnnnnnnn（5.5.13）（5.5.14）若取若取 ，求解（，求解（5.5.6）得著名的三步四阶隐式）得著名的三步四阶隐式Hamm
48、ing公式及其余项：公式及其余项：003230，。，65511121401)2(83981hOxyhTfffhyyynnnnnnnn （5.5.15）（5.5.16）若取若取 ，求解（，求解（5.5.6）得到隐式）得到隐式Simpson公式及其公式及其余项：余项：003320，。，65511111901)4(3hOxyhTfffhyynnnnnnn 例例5.5 分别取分别取 h=0.2，2，用四阶显式，用四阶显式 Milne公式和四阶隐式公式和四阶隐式 Hamming公式公式求解例求解例5.4所给的初值问题。所给的初值问题。解解 我们用单步法提供多步法的初值。由我们用单步法提供多步法的初值。由
49、4阶经典阶经典R-K公式为公式为Milne公式提公式提供初值供初值 ，为，为 Hamming公式提供公式提供 。h=0.2和和h=2时的时的计算结果及准确解之间的误差分别列于表计算结果及准确解之间的误差分别列于表8-9和表和表8-10。从表从表5-9看出，两种多步法的计算精度都很高，看出，两种多步法的计算精度都很高，Hamming公式化比公式化比 Milne公公式更精确。这是因为式更精确。这是因为 Hamming公式的截断误差主项的系数比公式的截断误差主项的系数比 Milne公式小。从公式小。从表表5-10看到，当计算步长变大后，显式多步法看到，当计算步长变大后，显式多步法 Milne公式的计
50、算结果误差增大，公式的计算结果误差增大，不稳定，而隐式多步法不稳定，而隐式多步法 Hamming公式的计算结果仍然是稳定的，这说明隐式公公式的计算结果仍然是稳定的，这说明隐式公式的稳定性比同阶的显式公式好。式的稳定性比同阶的显式公式好。3210yyyy，210yyy，表表5-955565107.4104.1108.3100.5108.166666106.4108.4108.4106.4102.4 Milne方法方法 误差误差 Hammins方法方法 误差误差2.2 0.94294268 0.94291955 2.4 1.12283349 1.12283386 2.6 1.30643214 1.
