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

《常微分方程初值问题数值解法.pptx》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法.pptx（162页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、会计学1常微分方程初值问题数值解法常微分方程初值问题数值解法2 定理定理1 1 设 在区域 上连续，关于 满足利普希茨条件，则对任意 ，常微分方程(1.1),(1.2)式当 时存在唯一的连续可微解 .关于解对扰动的敏感性，有以下结论.定理定理2 2 设 在区域 （如定理1所定义）上连续，且关于 满足利普希茨条件，设初值问题的解为 ，则第1页/共161页3 这个定理表明解对初值的敏感性，它与右端函数 有关，当 的Lips.常数 比较小时，解对初值和右端函数相对不敏感，可视为好条件.若 较大则可认为坏条件，即为病态问题.如果右端函数可导，由中值定理有 若假定 在域 内有界，设 ，则第2页/共161

2、页4 求解常微分方程的解析方法只能用来求解一些特殊类型的方程，实际中归结出来的微分方程主要靠数值解法.它表明 满足利普希茨条件，且 的大小反映了右端函数 关于 变化的快慢，刻画了初值问题(1.1)和(1.2)式是否是好条件.第3页/共161页5 所谓数值解法数值解法，就是寻求解 在一系列离散节点 上的近似值 .相邻两个节点的间距 称为步长步长.如不特别说明，总是假定 为定数，这时节点为 .初值问题(1.1),(1.2)的数值解法的基本特点是采取“步进式”.即求解过程顺着节点排列的次序一步一步地向前推进.第4页/共161页6 描述这类算法，只要给出用已知信息 计算 的递推公式.一类是计算 时只用

3、到前一点的值 ，称为单步法单步法.另一类是用到 前面 点的值 ，称为 步法步法.其次，要研究公式的局部截断误差和阶，数值解 与精确解 的误差估计及收敛性，还有递推公式的计算稳定性等问题.首先对方程 离散化，建立求数值解的递推公式.第5页/共161页9.29.2 简单的数值方法简单的数值方法简单的数值方法简单的数值方法 7 9.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法 积分曲线上一点 的切线斜率等于函数 的值.如果按函数 在 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.在 平面上，微分方程 的解 称作它的积分曲线积分曲线.第6页/共161页8 基于上述

4、几何解释，从初始点 出发，先依方向场在该点的方向推进到 上一点 ，然后再从 依方向场的方向推进到 上一点 ，循此前进做出一条折线 （图9-1).图9-1第7页/共161页9 一般地，设已做出该折线的顶点 ，过 依方向场的方向再推进到 ，显然两个顶点 的坐标有关系 即（2.1）这就是著名的欧拉欧拉（Euler）方法方法.欧拉法实际上是对常微分方程中的导数用均差近似，即直接得到的.第8页/共161页10若初值 已知，则依公式（2.1）可逐步算出第9页/共161页11 例例1 1 求解初值问题（2.2）解解 欧拉公式的具体形式为 取步长 ，计算结果见表9-1.初值问题（2.2）的解为 ，按这个解析式

5、子算出的准确值 同近似值 一起列在表9-1中，两者相比较可以看出欧拉方法的精度很差.第10页/共161页12 还可以通过几何直观来考察欧拉方法的精度.假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图9-2）.第11页/共161页13图9-2 从图形上看，这样定出的顶点 显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.为了分析计算公式的精度，通常可用泰勒展开将 在 处展开，则有 第12页/共161页14在 的前提下，称为此方法的局部截断误差.于是可得欧拉法（2.1）的误差（2.3）（2.4）如果对方程 从 到 积分，得（2.1）第13页/共161页15右

6、端积分用左矩形公式 近似.再以 代替 如果在（2.4）中右端积分用右矩形公式 （2.5）称为后退的欧拉法后退的欧拉法.代替 也得到（2.1），局部截断误差也是（2.3）.近似，则得另一个公式 它也可以通过利用均差近似导数 ，即直接得到.第14页/共161页16 欧拉公式是关于 的一个直接的计算公式，这类公式称作是显式的显式的；后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的隐式的.第15页/共161页17 隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化.设用欧拉公式 给出迭代初值 ，用它代入（2.5）式的右端，使之转化为显式，直接计算得 然后再用 代入（2

7、.5）式，又有 第16页/共161页18如此反复进行，得（2.6）由于 对 满足利普希茨条件(1.3).由(2.6)减(2.5)得 由此可知，只要 迭代法(2.6)就收敛到解 .从积分公式可以看到后退欧拉方法的公式误差与欧拉法是相似的.第17页/共161页 9.2.2 9.2.2 梯形方法梯形方法梯形方法梯形方法 19 若用梯形求积公式近似等式（2.4）右端的积分并分别用 代替 则可得到比欧拉法精度高的计算公式（2.7）称为梯形方法梯形方法.梯形方法是隐式单步法，可用迭代法求解.（2.4）第18页/共161页20 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得（2.8）同后退的欧

8、拉方法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为（2.7）第19页/共161页21如果选取 充分小，使得 则当 时有 ，这说明迭代过程(2.8)是收敛的.于是有 式中 为 关于 的利普希茨常数.第20页/共161页 9.2.3 9.2.3 改进欧拉公式改进欧拉公式改进欧拉公式改进欧拉公式 22 梯形方法虽然提高了精度，但其算法复杂.在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 ,而迭代又要反复进行若干次，计算量很大，而且往往难以预测.称之为预测值

9、预测值，第21页/共161页23 这样建立的预测-校正系统通常称为改进的欧拉公式：改进的欧拉公式：预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正值校正值.预测校正（2.9）也可以表为下列平均化形式（2.7）（2.8）第22页/共161页24 例例2 2 用改进的欧拉方法求解初值问题（2.2）.解解 这里 改进的欧拉公式为（2.2）第23页/共161页25仍取 ，计算结果见表9-2.同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.第24页/共161页 9.2.49.2.4 单步法的局部截断误差与阶单步法的局部截断误差与阶单步法的局

10、部截断误差与阶单步法的局部截断误差与阶 26 初值问题(1.1)，(1.2)的单步法可用一般形式表示为（2.10）其中多元函数 与 有关，当 含有 时，方法是隐式的，若不含 则为显式方法，（2.11）称为增量函数，所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.（1.1）（1.2）（2.1）（2.3）第25页/共161页27 对一般显式单步法则可如下定义.定义定义1 1 设 是初值问题(1.1)，(1.2)的准确解，称（2.12）为显式单步法（2.11）的局部截断误差局部截断误差.之所以称为局部的，是假设在 前各步没有误差.当 时，计算一步，则有（1.1）

11、（1.2）（2.11）第26页/共161页28在前一步精确的情况下用公式（2.11）计算产生的公式误差.根据定义，欧拉法的局部截断误差 即为（2.3）的结果.这里 称为局部截断误差主项.局部截断误差可理解为用方法（2.11）计算一步的误差，即显然（2.11）（2.3）第27页/共161页29 定义定义2 2 设 是初值问题(1.1)，(1.2)的准确解，若存在最大整数 使显式单步法(2.11)的局部截断误差满足（2.13）则称方法（2.11）具有 阶精度阶精度.若将（2.13）展开式写成 则 称为局部截断误差主项局部截断误差主项.以上定义对隐式单步法（2.10）也是适用的.（1.1）（1.2）

12、（2.11）（2.10）第28页/共161页30 对后退欧拉法（2.5）其局部截断误差为 这里 ，是1阶方法，局部截断误差主项为 .（2.5）第29页/共161页31 对梯形法（2.7）有 所以梯形方法是二阶的，其局部误差主项为（2.7）第30页/共161页9.39.3 龙格龙格龙格龙格-库塔方法库塔方法库塔方法库塔方法32第31页/共161页 9.3.1 9.3.1 显式龙格显式龙格显式龙格显式龙格-库塔法的一般形式库塔法的一般形式库塔法的一般形式库塔法的一般形式33 上节给出了显式单步法的表达式其局部截断误差为对欧拉法 ，即方法为 阶，（3.1）若用改进欧拉法，它可表为 第32页/共161

13、页34此时增量函数（3.2）与欧拉法的 相比，增加了计算一个右函数 的值，可望 .若要使得到的公式阶数 更大，就必须包含更多的 值.（3.3）从方程 等价的积分形式（2.4），即 第33页/共161页35若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，必然要增加求积节点.为此可将（3.3）的右端用求积公式表示为点数 越多，精度越高，上式右端相当于增量函数 ，为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为（3.4）其中（3.3）第34页/共161页36（3.5）这里 均为常数.(3.4)和(3.5)称为 级显式龙格龙格-库塔库塔（Runge-Kutta）法法，简称R-K方

14、法.当 时，就是欧拉法，此时方法的阶为 .当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，（3.4）第35页/共161页37 下面将证明阶 .要使公式(3.4),(3.5)具有更高的阶 ,就要增加点数 .下面就 推导R-K方法.（3.4）（3.5）第36页/共161页 9.3.29.3.2 二阶二阶二阶二阶显式显式显式显式R-KR-K方法方法方法方法 38 对 的R-K方法，计算公式如下（3.6）这里 均为待定常数.希望适当选取这些系数，使公式阶数 尽量高.根据局部截断误差定义，（3.6）的局部截断误差为（3.7）第37页/共161页39这里 .为得到 的阶 ，要将上式各项在 处做泰勒

15、展开，由于 是二元函数，故要用到二元泰勒展开，各项展开式为其中（3.8）第38页/共161页40将以上结果代入局部截断误差公式则有 要使公式（3.6）具有 阶，必须使（3.6）第39页/共161页41即 非线性方程组（3.9）的解是不唯一的.令 ，则得 这样得到的公式称为二阶R-K方法，如取 ，则这就是改进欧拉法（3.1）.（3.9）第40页/共161页42若取 ，则 得计算公式.称为中点公式中点公式，相当于数值积分的中矩形公式.（3.10）也可表示为（3.10）的R-K公式(3.6)的局部误差不可能提高到 .（3.6）第41页/共161页43 把 多展开一项，从（3.8）的 看到展开式中 的

16、项是不能通过选择参数消掉的.实际上要使 的项为零，需增加3个方程，要确定4个参数 ，这是不可能的.故 的显式R-K方法的阶只能是 ，而不能得到三阶公式.第42页/共161页 9.3.39.3.3 三阶与四阶显式三阶与四阶显式三阶与四阶显式三阶与四阶显式R-KR-K方法方法方法方法 44 要得到三阶显式R-K方法，必须 .（3.11）其中 及 均为待定参数.此时（3.4）,（3.5）的公式表示为 公式（3.11）的局部截断误差为（3.4）（3.5）第43页/共161页45只要将 按二元函数泰勒展开，使 ，可得待定参数满足方程（3.12）第44页/共161页46这是8个未知数6个方程的方程组，解也

17、不是唯一的.所以这是一簇公式.满足条件（3.12）的公式（3.11）统称为三阶R-K公式.一个常见的公式为 此公式称为库塔库塔三阶方法.第45页/共161页47 继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为 .四阶龙格-库塔方法的每一步需要计算四次函数值 ，（3.13）第46页/共161页48 例例3 3 设取步长 ，从 到 用四阶龙格-库塔方法求解初值问题 解解 这里 ，经典的四阶龙格-库塔公式为第47页/共161页49 表9-3列出了计算结果，同时了出了相应的精确解.比较例3和例2的计算结果，显然以龙格-库塔方法的精

18、度为高.第48页/共161页50 但由于这里放大了步长 ，所以表9-3和表9-2所耗费的计算量几乎相同.龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质.反之，如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法.四阶龙格-库塔方法每一步要4次计算函数 ,其计算量比每一步只要2次计算函数 的二阶龙格-库塔方法-改进的欧拉方法大一倍.第49页/共161页 9.3.49.3.4 变步长的龙格变步长的龙格变步长的龙格变步长的龙格-库塔方法库塔方法库塔方法库塔方法 51 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求

19、解范围内所要完成的步数就增加了.步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累.因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题.在选择步长时，需要考虑两个问题：1 怎样衡量和检验计算结果的精度？第50页/共161页52 2 如何依据所获得的精度处理步长？考察经典的四阶龙格-库塔公式（3.13）从节点 出发，先以 为步长求出一个近似值 ，第51页/共161页53（3.14）然后将步长折半，即取 为步长从 跨两步到 ，再求得一个近似值 ，每跨一步的截断误差是 ，因此有（3.15）比较（3.14）式和（3.15）式我们看到，步长折半后，由于公式的局部截断误差为 ，故

20、有 误差大约减少到 ，第52页/共161页54由此易得下列事后估计式 这样，可以通过检查步长，折半前后两次计算结果的偏差 即有来判定所选的步长是否合适.具体地说，将区分以下两种情况处理：第53页/共161页55 1.对于给定的精度 ，如果 ，反复将步长折半进行计算，直至 为止.这时取最终得到的 作为结果；2.如果 ，反复将步长加倍，直到 为止，这种通过加倍或折半处理步长的方法称为变步长方法变步长方法.这时再将步长折半一次，就得到所要的结果.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.第54页/共161页9.49.4 单步法的收敛性与稳定性单步法的收敛性与稳定性单步法的

21、收敛性与稳定性单步法的收敛性与稳定性 56 9.4.1 收敛性与相容性收敛性与相容性 数值解法的基本思想是通过某种离散化手段将微分方程转化为差分方程，如单步法(2.11)，即（4.1）它在 处的解为 ，而初值问题（1.1）,（1.2）在 处的精确解为 ，记 称为整体截断误差.（1.1）（1.2）第55页/共161页57 收敛性就是讨论当 固定且 时的问题.定义定义3 3 若一种数值方法对于固定的 ,当 时有 ，其中 是(1.1),(1.2)的准确解，则称该方法是收敛收敛的.显然数值方法收敛是指 .对单步法（4.1）有下述收敛性定理：（1.1）（1.2）（4.1）第56页/共161页58 定理定

22、理3 3 假设单步法（4.1）具有 阶精度，且增量函数 关于 满足利普希茨条件（4.2）又设初值 是准确的，即 ，则其整体截断误差整体截断误差（4.3）证明证明 设以 表示取 用公式（4.1）求得的结果，即（4.4）则 为局部截断误差，（4.1）（4.1）第57页/共161页59由于所给方法具有 阶精度，按定义2，存在定数 ，使又由式（4.4）与（4.1），得 利用假设条件（4.2），有 从而有（4.2）（4.1）（4.4）第58页/共161页60即对整体截断误差 成立下列递推关系式（4.5）反复递推，可得（4.6）再注意到当 时 最终得下列估计式（4.7）第59页/共161页61由此可以断定

23、，如果初值是准确的，即 ，则（4.3）式成立.依据这一定理，判断单步法（4.1）的收敛性，归结为验证增量函数 能否满足利普希茨条件（4.2）.对于欧拉方法，由于其增量函数 就是 ，故当 关于 满足利普希茨条件时它是收敛的.再考察改进的欧拉方法，其增量函数给出，这时有（4.3）（4.2）（4.1）第60页/共161页62假设 关于 满足利普希茨条件，记利普希茨常数为 ,设 为定数），上式表明 关于 的利普希茨常数则由上式推得 因此改进的欧拉方法也是收敛的.第61页/共161页63 类似地，也可验证其他龙格-库塔方法的收敛性.定理3表明 时单步法收敛，并且当 是初值问题（1.1），（1.2）的解，

24、（4.1）具有 阶精度时，有展开式 所以 的充要条件是 ，（4.1）（1.1）（1.2）第62页/共161页64而 ，于是可给出如下定义：定义定义4 4 若单步法（4.1）的增量函数 满足 则称单步法（4.1）与初值问题（1.1），（1.2）相容相容.相容性是指数值方法逼近微分方程(1.1)，即微分方程(1.1)离散化得到的数值方法当 时可得到第63页/共161页65 定理定理4 4 阶方法(4.1)与初值问题(1.1)，(1.2)相容的充分必要条件是 由定理3可知单步法（4.1）收敛的充分必要条件是(4.1)是相容的.以上讨论表明 阶方法（4.1）当 时与(1.1),(1.2)相容，反之相容

25、方法至少是1阶的.第64页/共161页 9.4.2 9.4.2 绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域 66 定义定义5 5 若一种数值方法在节点值 上大小为 的扰动，于以后各节点值 上产生的偏差均不超过 ，则称该方法是稳定稳定的.以欧拉法为例考察计算稳定性.例例4 4 考察初值问题 其准确解 是一个按指数曲线衰减得很快的函数，如图9-3所示.第65页/共161页67若取 ，则欧拉公式的具体形式为 计算结果列于表9-4的第2列.可以看到，欧拉方法的解 （图9-3中用号标出）在准确值 的上下波动，计算过程明显地不稳定.图9-3 用欧拉法解方程 得

26、 第66页/共161页68再考察后退的欧拉方法，取 时计算公式为 计算结果列于表9-4的第3列（图9-3中标以号），这时计算过程是稳定的.但若取 则计算过程稳定.第67页/共161页69 这表明稳定性不但与方法有关，也与步长 的大小有关，当然也与方程中的 有关.为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，（4.8）其中 为复数.例如在 的邻域，可展开为 模型方程为 对一般方程可以通过局部线性化化为这种形式.第68页/共161页70略去高阶项，再做变换即可得到 的形式.对于 个方程的方程组，也可线性化为 ，这里 为 的雅可比矩阵 .若 有 个特征值 ，则 还可能是复数，为

27、保证微分方程本身的稳定性，还应假定 .先研究欧拉方法的稳定性.模型方程 的欧拉公式为 所以，为了使模型方程结果能推广到方程组，方程中 应为复数.第69页/共161页71（4.9）设在节点值 上有一扰动值 ，它的传播使节点值 产生大小为 的扰动值，假设用 按欧拉公式得出 的计算过程不再有新的误差，则扰动值满足 可见扰动值满足原来的差分方程（4.9）.如果差分方程的解是不增长的，即有 则它就是稳定的.第70页/共161页72即图9-4 显然，为要保证差分方程（4.9）的解不增长，只要选取 充分小，（4.10）在 的复平面上，这是以 为圆心，1为半径的单位圆内部（图9-4）.这个圆域称为欧拉法的绝对

28、稳定域，一般情形可由下面定义.定义定义6 6 单步法（4.1）用于解模型方程（4.8），若得到的解 ，满足 ，则称方法（4.1）是绝对稳定绝对稳定的.（4.9）使（4.1）（4.8）第71页/共161页73 在 的平面上，使 的变量围成的区域，称为绝对稳定域绝对稳定域，对欧拉法，给出，绝对稳定区间为 .它与实轴的交称为绝对稳定区间绝对稳定区间.其绝对稳定域由 在例5中 ，即 为绝对稳定区间.当取 时 例4中取 故它是不稳定的，它是稳定的.第72页/共161页74故 绝对稳定域由 得到.绝对稳定区间为 ，即 .类似可得三阶及四阶的R-K方法的 分别为 用二阶R-K方法解模型方程 可得到 第73页

29、/共161页75由 可得到相应的绝对稳定域.当 为实数时则得绝对稳定区间.分别为 三阶显式R-K方法：即 四阶显式R-K方法：即 图9-5给出了R-K方法 到 的绝对稳定域.从以上讨论可知显式的R-K方法的绝对稳定域均为有限域，都对步长 有限制.如果 不在所给的绝对稳定区间内，方法就不稳定.图9-5第74页/共161页76 例例5 5 分别取 ，及 用经典的四阶R-K方法（3.13）计算.解解本例分别为 及 .前者在绝对稳定区间内，后者则不在.经典的四阶R-K方法为（3.13）第75页/共161页77 以上结果看到，如果步长 不满足绝对稳定条件，误差增长很快.用四阶R-K方法计算其误差见下表：

30、这里第76页/共161页78 对隐式单步法，可以同样讨论方法的绝对稳定性，例如对后退欧拉法，用它解模型方程可得故 由 可得绝对稳定域为 .它是以 为圆心，1为半径的单位圆外部，故绝对稳定区间为 .第77页/共161页79 当 时，则 ，即对任何步长均为稳定的.故 对 有 ，故绝对稳定域为 的左半平面，对隐式梯形法，解模型方程 得 第78页/共161页80绝对稳定区间为 ，即 时梯形法均是稳定的.定义定义7 7 如果数值方法的绝对稳定域包含了 那么称此方法是A-稳定的.隐式欧拉法与梯形方法的绝对稳定域均为 在具体计算中步长 的选择只需考虑计算精度及迭代收敛性要求而不必考虑稳定性，具有这种特点的方

31、法特别重要.由定义知A-稳定方法对步长 没有限制.第79页/共161页9.59.5 线性多步法线性多步法线性多步法线性多步法 81 在逐步推进的求解过程中，计算 之前事实上已经求出了一系列的近似值 ，如果充分利用前面多步的信息来预测 ，则可以期望会获得较高的精度.这就是构造所谓线性多步法的基本思想.本节主要介绍基于泰勒展开的构造方法.构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程 两端积分后利用插值求积公式得到.第80页/共161页 9.5.19.5.1 线性多步法的一般公式线性多步法的一般公式线性多步法的一般公式线性多步法的一般公式82 一般的线性多步法公式可表示

32、为（5.1）其中 为 的近似，计算时需先给出前面 个近似值 ,再由（5.1）逐次求出 .如果计算 时，除了使用 的值，还用到的值，则称此方法为线性多步法线性多步法.为常数.及 不全为零，则称为线性线性 步法步法.若第81页/共161页83这时 可直接由（5.1）算出；如果 ，称（5.1）为显式显式 步法步法，如果 ，则（5.1）称为隐式隐式 步法步法，求解时与梯形法（2.7）相同，要用迭代法方可算出 .（5.1）（2.7）第82页/共161页84 定义定义8 8 设 是初值问题(1.1)，(1.2)的准确解，线性多步法（5.1）在 上的局部截断误差为（5.2）（5.1）中系数 及 可根据方法的

33、局部截断误差及阶确定，其定义为：（5.1）（1.1）（1.2）若 ，则称方法（5.1）是 阶的阶的，则称方法（5.1）与方程（1.1）是相容的相容的.第83页/共161页85 由定义8，对 在 处做泰勒展开，由于 代入（5.2）得（5.3）（5.2）第84页/共161页86其中（5.4）若在公式（5.1）中选择系数 及 ，使它满足 由定义可知此时所构造的多步法是 阶的.（5.1）第85页/共161页87（5.5）称右端第一项为局部截断误差主项局部截断误差主项，称为误差常数误差常数.根据相容性定义，,即 ，故方法（5.1）与微分方程（1.1）相容的充分必要条件是（5.6）成立.且（5.6）由（5

34、.4）得（5.1）第86页/共161页88 当 时，若 ，则由（5.6）可求得 此时公式（5.1）为 即为欧拉法.从（5.4）可求得 ，故方法为1阶精度，这和第2节给出的定义及结果是一致的.且局部截断误差为（5.6）第87页/共161页89 对 ,若 ，方法为隐式公式.为了确定系数 ，可由 解得于是得到公式 即为梯形法.由（5.4）可求得 ,故 ，所以梯形法是二阶方法，其局部截断误差主项是 .这与第9.2节中的讨论也是一致的.第88页/共161页90 对 的多步法公式都可利用（5.4）确定系数 ，并由（5.5）给出局部截断误差.（5.5）第89页/共161页 9.5.29.5.2 阿当姆斯显式

35、与隐式公式阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式 91 考虑形如（5.7）的 步法，称为阿当姆斯阿当姆斯（Adams)方法方法.为显式方法，为隐式方法，通常称为阿当姆斯显式与隐式公式，也称Adams-Bashforth 公式与Adams-Monlton公式.这类公式可直接由对方程 两端从 到 积分求得.第90页/共161页92 也可以利用(5.4)由 推出，对比 与 可知此时系数 .显然 成立，下面只需确定系数 ，可令 ，求得 .若 ，则令 来求得 .第91页/共161页93 以 为例，由 ，根据 若 ，则由前三个方程解得得到 的阿当姆斯显式公式是（5.8）第92页

36、/共161页94由（5.4）求得 ，所以（5.8）是三阶方法，局部截断误差是 若 ，则可解得 于是得 的阿当姆斯隐式公式为（5.9）它是四阶方法，局部截断误差是（5.8）第93页/共161页95（5.10）类似的方法可求得阿当姆斯显式方法和隐式方法的公式，表9-5及表9-6分别列出了 时的阿当姆斯显式公式与阿当姆斯隐式公式，其中 为步数，为方法的阶，为误差常数.第94页/共161页96第95页/共161页97第96页/共161页98 例例6 6 用四阶阿当姆斯显式和隐式方法解初值问题 取步长 .解解本题 .从四阶阿当姆斯隐式公式得到 从四阶阿当姆斯显式公式得到 第97页/共161页99由此可直

37、接解出 而不用迭代，得到 计算结果见表9-8，其中显式方法中的初值 及隐式方法中的 均用准确解 得到，对一般方程，可用四阶R-K方法计算初始近似.(表略，见教材.)第98页/共161页100 从以上例子看到同阶的阿当姆斯方法，隐式方法要比显式方法误差小.这可以从两种方法的局部截断误差主项的系数大小得到解释.这里 分别为 及 .第99页/共161页 9.5.39.5.3 米尔尼方法与辛普森方法米尔尼方法与辛普森方法米尔尼方法与辛普森方法米尔尼方法与辛普森方法 101 考虑另一个 的显式公式 其中 为待定常数.由（5.4）可知 ，再令 得到 可根据使公式的阶尽可能高这一条件来确定其数值.第100页

38、/共161页102解此方程组得 于是得到四步显式公式（5.11）称为米尔尼米尔尼(Milne）方法方法.由于 ，故方法为4阶，其局部截断误差为（5.12）米尔尼方法也可以通过方程 两端积分 第101页/共161页103得到.右端积分通过辛普森求积公式就有（5.13）称为辛普森方法辛普森方法.（5.14）它是隐式二步四阶方法，其局部截断误差为 若将方程 从 到 积分，可得 第102页/共161页 9.5.49.5.4 汉明方法汉明方法汉明方法汉明方法 104 辛普森公式是二步方法中阶数最高的，但它的稳定性差，为了改善稳定性，考察另一类三步法公式 其中系数 及 为常数.若希望导出的公式是四阶的，则

39、系数中至少有一个自由参数.若取 ，则可得到辛普森公式.若取 ，仍利用泰勒展开，由（5.4），令第103页/共161页105解此方程组得 于是有（5.15）则可得到 第104页/共161页106称为汉明汉明(Hamming)方法方法.由于 ，故方法是四阶的，且局部截断误差（5.16）第105页/共161页 9.5.59.5.5 预测预测预测预测-校正方法校正方法校正方法校正方法 107 对于隐式的线性多步法，计算时要进行迭代，计算量较大.为了避免进行迭代，通常采用显式公式给出 的一个初始近似，记为 ，称为预测预测(predictor).接着计算 的值(evaluation)，再用隐式公式计算 ，

40、称为校正校正（corrector）.在（2.13）中用欧拉法做预测，再用梯形法校正，得到改进欧拉法，它就是一个二阶预测-校正方法.（2.13）第106页/共161页108 例如用四阶的阿当姆斯显式方法做预测，再用四阶阿当姆斯隐式公式做校正，得到以下格式：一般情况下，预测公式与校正公式都取同阶的显式方法与隐式方法相匹配.预测P：求值E：校正C：求值E：此公式称为阿当姆斯四阶预测阿当姆斯四阶预测-校正格式校正格式（PECE).第107页/共161页109 依据四阶阿当姆斯公式的截断误差，对于PECE的预测步P有 对校正步C有 两式相减得 于是有下列事后误差估计 第108页/共161页110容易看出

41、 比 更好.（5.17）但在 的表达式中 是未知的，因此计算时用上一步代替，从而构造一种修正预测修正预测-校正格式校正格式(PMECME):第109页/共161页111P：M：E：C：M：E：在PMECME格式中已将（5.17）的 及 分别改为 及 .第110页/共161页112 利用米尔尼公式（5.11）和汉明公式（5.15）相匹配，并利用截断误差（5.12），（5.16）改进计算结果，可类似地建立四阶修正米尔尼修正米尔尼-汉明预测汉明预测-校正格式校正格式PMECME）：P：M：E：C：M：E：（5.15）（5.11）（5.16）（5.12）第111页/共161页 9.5.69.5.6 构

42、造多步法公式的注记和例构造多步法公式的注记和例构造多步法公式的注记和例构造多步法公式的注记和例 113 构造多步法公式有基于数值积分和泰勒展开两种途径，数值积分方法只能用于可转化为等价的积分方程的微分方程，它是有局限性的.而用泰勒展开则可构造任意多步法公式，其做法是根据多步法公式的形式，直接在 处做泰勒展开.确定多步法（5.1）的系数 及 时不必套用系数公式（5.4），因为多步法公式不一定如（5.1）的形式，而且套用公式容易记错.（5.1）第112页/共161页114 例例7 7 解初值问题 用显式二步法 其中 试确定参数 使方法阶数尽可能高,并求局部截断误差.解解 根据局部截断误差定义，用泰

43、勒展开确定参数满足的方程.第113页/共161页115由于 第114页/共161页116为求参数 使方法阶数尽量高，可令 即得方程组 第115页/共161页117解得 ，此时公式为三阶，即为所求局部截断误差.而所得二步法为 而且 第116页/共161页118 例例8 8 证明存在 的一个值，使线性多步法 是四阶的.证明证明只要证明局部截断误差 ，则方法 仍用泰勒展开，由于 为四阶.第117页/共161页119第118页/共161页120当 时，故方法是四阶的.第119页/共161页1219.6 线性多步法的收敛性与稳定性线性多步法的收敛性与稳定性 9.6.1 相容性及收敛性相容性及收敛性 线性

44、多步法(5.1)式的相容性：在定义8中给出的局部截断误差(5.2)中 ，若 称 步法(5.1)与微分方程(1.1)相容，它等价于 对多步法(5.1)可引入多项式（6.1）（6.2）第120页/共161页122和（6.3）分别称为线性多步法(5.1)的第一特征多项式和第二特征多项式.可以看出，如果(5.1)式给定，则 和 也完全确定.反之也成立.第121页/共161页123 定理定理5 5 线性多步法(5.1)式与微分方程(1.1)相容的充分必要条件是 关于多步法(5.1)的收敛性：由于用多步法求数值解需要 个初值，而微分方程(1.1)只给出一个初值 ，因此还要给出 个初值才能用多步法(5.1)

45、进行求解，即（6.4）（6.5）其中 由微分方程初值给定，可由相应单步法给出.第122页/共161页124 设由(6.5)式在 处得到的数值解为 ，这里 为固定点，于是有下面定义.定义定义9 9 设初值问题(1.1),(1.2)有精确解 .如果初始条件 满足条件 的线性 步法(6.5)在 处的解 有 则称线性 步法是收敛的.定理定理6 6 设线性多步法(6.5)是收敛的，则它是相容的.此定理的逆定理是不成立的.第123页/共161页125 例例9 9 用线性二步法 解初值问题 解解 此初值问题精确解 ，而由(6.6)式知 故有 ，故方法(6.6)是相容的，但方法(6.6)的解并不收敛，在方法(

46、6.6)中取初值（6.6）（6.7）此时方法(6.6)为二阶差分方程第124页/共161页126其特征方程为解得其根为 及 .于是可求得(6.8)的解为故方法不收敛.多步法是否收敛与 的根有关.（6.8）第125页/共161页127 定义定义1010 如果多步法（5.1）式的第一特征多项式 的根都在单位圆内或圆上，且在单位圆上的根为单根，则称线性多步法（5.1）满足根条件.定理定理7 7 线性多步法是相容的，则线性多步法（6.5）收敛的充分必要条件是线性多步法（5.1）满足根条件.第126页/共161页128 9.6.2 稳定性与绝对稳定性稳定性与绝对稳定性 稳定性主要研究初始条件扰动与差分方

47、程右端项扰动对数值解的影响，假设多步法（6.5）有扰动 ，则经过扰动后的解为 ，它满足方程（6.9）第127页/共161页129 定义定义1111 对初值问题(1.1)，(1.2)，由方法（6.5）得到的差分方程解 ，由于有扰动 ，使得方程（6.9）的解为 ，若存在常数 及 ，使对所有 ，当 时，有则称多步法（5.1）是稳定的或称为零稳定的.研究零稳定性就是研究 时差分方程（6.5）解 的稳定性.它表明当初始扰动或右端扰动不大时，解的误差也不大.第128页/共161页130 对多步法（5.1），当 时对应差分方程的特征方程为 .故有 定理定理8 8 线性多步法（5.1）是稳定的充分必要条件是它

48、满足根条件.关于绝对稳定性只要将多步法（5.1）用于解模型方程（4.8），得到线性差分方程（6.10）第129页/共161页131利用线性多步法的第一、第二特征多项式 ，令此式称为线性多步法的稳定多项式，它是关于 的 次多项式.如果它的所有零点 满足 ，则(6.10)式的解 当 时，有 .（6.11）第130页/共161页132 定义定义1212 对于给定的 ，如果稳定多项式(6.11)的零点满足 ，则称线性多步法(5.1)关于此 值是绝对稳定的.若在 的复平面的某个区域 中所有 值线性多步法(5.1)都是绝对稳定的，而在区域 外，方法是不稳定的，则称 为多步法(5.1)的绝对稳定域绝对稳定域

49、.与实轴的交集称为线性多步法(5.1)的绝对稳定区间绝对稳定区间.为实数时，可以只讨论绝对稳定区间.由于线性多步法的绝对稳定域较为复杂，通常采用根轨迹法.第131页/共161页133 阿当姆斯显式方法和隐式方法的绝对稳定域图形分别为图9-6和图9-7，绝对稳定区间见表9-9.图9-6图9-7第132页/共161页134 例例1010 讨论辛普森方法的稳定性 解解 辛普森方法的第第二特征多项式为 的根分别为-1及1，它满足根条件，故方法是零稳定的.但它的稳定多项式为求绝对稳定域 的边界轨迹 .若 ，则可令 ，在 平面域 的边界轨迹 为第133页/共161页135可看出 在虚轴上，且对全部 从而可

50、知 为虚轴上从 到 的线段，故辛普森公式的绝对稳定域为空集.即步长 ，此方法都不是绝对稳定的，故它不能用于求解.第134页/共161页9.79.7 一阶方程组与刚性方程一阶方程组与刚性方程一阶方程组与刚性方程一阶方程组与刚性方程 136 9.7.1 一阶方程组一阶方程组 前面研究了单个方程 的数值解法，只要把 和 理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形.考察一阶方程组 的初值问题，初始条件给为 第135页/共161页137则上述方程组的初值问题可表示为（7.1）求解这一初值问题的四阶龙格-库塔公式为 式中 若采用向量的记号，记第136页/共161页138 考察两个方程