常微分方程数值解优秀PPT.ppt

常微分方程数值解1你现在浏览的是第一页，共151页2你现在浏览的是第二页，共151页9.1 引引 言言 考虑一阶常微分方程的初值问题（1.1）（1.2）如果存在实数 ，使得（1.3）则称 关于 满足利普希茨利普希茨(Lipschitz）条件条件，称为 的利普利普希茨常数希茨常数（简称Lips.常数）.3你现在浏览的是第三页，共151页 定理定理1 1 设 在区域 上连续，关于 满足利普希茨条件，则对任意 ，常微分方程(1.1),(1.2)式当 时存在唯一的连续可微解 .关于解对扰动的敏感性，有以下结论.定理定理2 2 设 在区域 （如定理1所定义）上连续，且关于 满足利普希茨条件，设初值问题的解为 ，则4你现在浏览的是第四页，共151页两者的区别:1.问题 5你现在浏览的是第五页，共151页6你现在浏览的是第六页，共151页7你现在浏览的是第七页，共151页左矩形右矩形梯形公式单步法:对初值问题,计算yn+1时只用到前一点的值yn,即yn+1=f(yn)k步法:计算yn+1时需要用到前k点的值yn,yn-1,yn-k+1,即yn+1=f(yn,yn-1,yn-k+1)对方程 离散化，建立求数值解的递推公式.描述这类算法，只要给出用已知信息 计算 的递推公式.8你现在浏览的是第八页，共151页9.2 简单的数值方法简单的数值方法 求解一阶微分方程初值问题:9你现在浏览的是第九页，共151页几何意义:它是用一条自点(x0,y0)出发的折线段去逼近积分曲线y=y(x)如下图9-1P28010你现在浏览的是第十页，共151页 例例1 1 求解初值问题（2.2）解解 欧拉公式的具体形式为 取步长 ，计算结果见表9-1.初值问题（2.2）的解为 ，按这个解析式子算出的准确值 同近似值 一起列在表9-1中，两者相比较可以看出欧拉方法的精度很差.11你现在浏览的是第十一页，共151页 还可以通过几何直观来考察欧拉方法的精度.假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图9-2）.12你现在浏览的是第十二页，共151页图9-2 从图形上看，这样定出的顶点 显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.误差分析：误差分析：为了分析计算公式的精度,通常可用泰勒展开将 在 处展开，则有 13你现在浏览的是第十三页，共151页在 的前提下，称为此方法的局部截断误差.于是可得欧拉法（2.1）的误差（2.3）（2.1）估算=精确14你现在浏览的是第十四页，共151页（2.5）称为后退的欧拉法后退的欧拉法(隐式欧拉公式隐式欧拉公式).欧拉公式是关于 的一个直接的计算公式，这类公式称作是显式的显式的；后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的隐式的.15你现在浏览的是第十五页，共151页 隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化.设用欧拉公式 给出迭代初值 ，用它代入（2.5）式的右端，使之转化为显式，直接计算得 然后再用 代入（2.5）式，又有 16你现在浏览的是第十六页，共151页如此反复进行，得（2.6）由于 对 满足利普希茨条件(1.3).由(2.6)减(2.5)得 由此可知，只要 迭代法(2.6)就收敛到解 .17你现在浏览的是第十七页，共151页 9.2.2 梯形方法梯形方法 （2.7）称为梯形方法梯形方法.梯形方法是隐式隐式单步法，可用迭代法求解.18你现在浏览的是第十八页，共151页 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得（2.8）同后退的欧拉方法一样，仍用欧拉方法欧拉方法提供迭代初值初值，则梯形法的迭代公式为（2.7）19你现在浏览的是第十九页，共151页如果选取 充分小，使得 则当 时有 ，这说明迭代过程(2.8)是收敛的.于是有 式中 为 关于 的利普希茨常数.20你现在浏览的是第二十页，共151页 9.2.3 改进欧拉公式改进欧拉公式 梯形方法虽然提高了精度，但其算法复杂.在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 ,而迭代又要反复进行若干次，计算量很大，而且往往难以预测.称之为预测值预测值，21你现在浏览的是第二十一页，共151页 这样建立的预测-校正系统通常称为改进的欧拉公式：改进的欧拉公式：预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正校正值值.预测校正（2.9）也可以表为下列平均化形式（2.7）（2.8）22你现在浏览的是第二十二页，共151页 例例2 2 用改进的欧拉方法求解初值问题（2.2）.解解 这里 改进的欧拉公式为（2.2）23你现在浏览的是第二十三页，共151页仍取 ，计算结果见表9-2.同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.24你现在浏览的是第二十四页，共151页 9.2.4 单步法的局部截断误差与阶单步法的局部截断误差与阶 初值问题(1.1)，(1.2)的单步法可用一般形式表示为（2.10）其中多元函数 与 有关，当 含有 时，方法是隐式的，若不含 则为显式方法，（2.11）称为增量函数，所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.（1.1）（1.2）（2.1）（2.3）25你现在浏览的是第二十五页，共151页 对一般显式单步法则可如下定义.定义定义1 1 设 是初值问题(1.1)，(1.2)的准确解，称（2.12）为显式单步法（2.11）的局部截断误差局部截断误差.之所以称为局部的，是假设在 前各步没有误差.当 时，计算一步，则有（1.1）（1.2）（2.11）26你现在浏览的是第二十六页，共151页在前一步精确的情况下用公式（2.11）计算产生的公式误差.根据定义，欧拉法的局部截断误差 即为（2.3）的结果.这里 称为局部截断误差主项.局部截断误差可理解为用方法（2.11）计算一步的误差，即显然（2.11）（2.3）27你现在浏览的是第二十七页，共151页 定义定义2 2 设 是初值问题(1.1)，(1.2)的准确解，若存在最大整数 使显式单步法(2.11)的局部截断误差满足（2.13）则称方法（2.11）具有 阶精度阶精度.若将（2.13）展开式写成 则 称为局部截断误差主项局部截断误差主项.以上定义对隐式单步法（2.10）也是适用的.（1.1）（1.2）（2.11）（2.10）28你现在浏览的是第二十八页，共151页 对后退欧拉法（2.5）其局部截断误差为 这里 ，是1 1阶方法，局部截断误差主项为 .（2.5）29你现在浏览的是第二十九页，共151页 对梯形法（2.7）有 所以梯形方法是二阶的，其局部误差主项为（2.7）30你现在浏览的是第三十页，共151页9.3 龙格龙格-库塔方法库塔方法31你现在浏览的是第三十一页，共151页一、Taylor展开法取等式右边前p+1项32你现在浏览的是第三十二页，共151页例例 取h=0.1,用三阶Taylor展开法求解33你现在浏览的是第三十三页，共151页从计算高阶导数的公式知道,方法的截断误差提高一阶,需要增加的计算量很大.下面我们用区间上若干点的导数f,而不是高阶导数,将它们作线性组合得到平均斜率,将其与解的Taylor展开相比较,使前面若干项吻合，从而得到提高阶的方法34你现在浏览的是第三十四页，共151页2 龙格龙格-库塔法库塔法(Runge-Kutta(Runge-Kutta法法)35你现在浏览的是第三十五页，共151页（2）龙格龙格-库塔法的一般形式库塔法的一般形式36你现在浏览的是第三十六页，共151页37你现在浏览的是第三十七页，共151页38你现在浏览的是第三十八页，共151页将以上结果代入局部截断误差公式则有 要使公式（3.6）具有 阶，必须使（3.6）39你现在浏览的是第三十九页，共151页即 非线性方程组（3.9）的解是不唯一的.令 ，则得 这样得到的公式称为二阶R-K方法，如取 ，则这就是改进欧拉法（3.1）.（3.9）40你现在浏览的是第四十页，共151页若取 ，则 得计算公式.称为中点公式中点公式，相当于数值积分的中矩形公式.（3.10）也可表示为（3.10）41你现在浏览的是第四十一页，共151页 9.3.3 三阶与四阶显式三阶与四阶显式R-K方法方法 要得到三阶显式R-K方法，必须 .（3.11）其中 及 均为待定参数.此时（3.4）,（3.5）的公式表示为 公式（3.11）的局部截断误差为（3.4）（3.5）42你现在浏览的是第四十二页，共151页只要将 按二元函数泰勒展开，使 ，可得待定参数满足方程（3.12）43你现在浏览的是第四十三页，共151页这是8个未知数6个方程的方程组，解也不是唯一的.所以这是一簇公式.满足条件（3.12）的公式（3.11）统称为三阶R-K公式.一个常见的公式为 此公式称为库塔库塔三阶方法.44你现在浏览的是第四十四页，共151页 继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为 .四阶龙格-库塔方法的每一步需要计算四次函数值 ，（3.13）45你现在浏览的是第四十五页，共151页 谢 谢！46你现在浏览的是第四十六页，共151页47你现在浏览的是第四十七页，共151页 9.3.4 变步长的龙格变步长的龙格-库塔方法库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了.步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累.因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题.在选择步长时，需要考虑两个问题：1 怎样衡量和检验计算结果的精度？48你现在浏览的是第四十八页，共151页 2 如何依据所获得的精度处理步长？考察经典的四阶龙格-库塔公式（3.13）从节点 出发，先以 为步长求出一个近似值 ，49你现在浏览的是第四十九页，共151页（3.14）然后将步长折半，即取 为步长从 跨两步到 ，再求得一个近似值 ，每跨一步的截断误差是 ，因此有（3.15）比较（3.14）式和（3.15）式我们看到，步长折半后，由于公式的局部截断误差为 ，故有 误差大约减少到 ，50你现在浏览的是第五十页，共151页由此易得下列事后估计式 这样，可以通过检查步长，折半前后两次计算结果的偏差 即有来判定所选的步长是否合适.具体地说，将区分以下两种情况处理：51你现在浏览的是第五十一页，共151页 1.对于给定的精度 ，如果 ，反复将步长折半进行计算，直至 为止.这时取最终得到的 作为结果；2.如果 ，反复将步长加倍，直到 为止，这种通过加倍或折半处理步长的方法称为变步长方法变步长方法.这时再将步长折半一次，就得到所要的结果.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.52你现在浏览的是第五十二页，共151页9.4 单步法的收敛性与稳定性单步法的收敛性与稳定性 9.4.1 收敛性与相容性收敛性与相容性 数值解法的基本思想是通过某种离散化手段将微分方程转化为差分方程，如单步法(2.11)，即（4.1）它在 处的解为 ，而初值问题（1.1）,（1.2）在 处的精确解为 ，记 称为整体截断误差.（1.1）（1.2）53你现在浏览的是第五十三页，共151页 收敛性就是讨论当 固定且 时的问题.定义定义3 3 若一种数值方法对于固定的 ,当 时有 ，其中 是(1.1),(1.2)的准确解，则称该方法是收敛收敛的.显然数值方法收敛是指 .对单步法（4.1）有下述收敛性定理：（1.1）（1.2）（4.1）54你现在浏览的是第五十四页，共151页 定理定理3 3 假设单步法（4.1）具有 阶精度，且增量函数 关于 满足利普希茨条件（4.2）又设初值 是准确的，即 ，则其整体截断误差整体截断误差（4.3）证明证明 设以 表示取 用公式（4.1）求得的结果，即（4.4）则 为局部截断误差，（4.1）（4.1）55你现在浏览的是第五十五页，共151页由于所给方法具有 阶精度，按定义2，存在定数 ，使又由式（4.4）与（4.1），得 利用假设条件（4.2），有 从而有（4.2）（4.1）（4.4）56你现在浏览的是第五十六页，共151页即对整体截断误差 成立下列递推关系式（4.5）反复递推，可得（4.6）再注意到当 时 最终得下列估计式（4.7）57你现在浏览的是第五十七页，共151页由此可以断定，如果初值是准确的，即 ，则（4.3）式成立.依据这一定理，判断单步法（4.1）的收敛性，归结为验证增量函数 能否满足利普希茨条件（4.2）.对于欧拉方法，由于其增量函数 就是 ，故当 关于 满足利普希茨条件时它是收敛的.再考察改进的欧拉方法，其增量函数给出，这时有（4.3）（4.2）（4.1）58你现在浏览的是第五十八页，共151页假设 关于 满足利普希茨条件，记利普希茨常数为 ,设 为定数），上式表明 关于 的利普希茨常数则由上式推得 因此改进的欧拉方法也是收敛的.59你现在浏览的是第五十九页，共151页 类似地，也可验证其他龙格-库塔方法的收敛性.定理3表明 时单步法收敛，并且当 是初值问题（1.1），（1.2）的解，（4.1）具有 阶精度时，有展开式 所以 的充要条件是 ，（4.1）（1.1）（1.2）60你现在浏览的是第六十页，共151页而 ，于是可给出如下定义：定义定义4 4 若单步法（4.1）的增量函数 满足 则称单步法（4.1）与初值问题（1.1），（1.2）相容相容.相容性是指数值方法逼近微分方程(1.1)，即微分方程(1.1)离散化得到的数值方法当 时可得到61你现在浏览的是第六十一页，共151页 定理定理4 4 阶方法(4.1)与初值问题(1.1)，(1.2)相容的充分必要条件是 由定理3可知单步法（4.1）收敛的充分必要条件是(4.1)是相容的.以上讨论表明 阶方法（4.1）当 时与(1.1),(1.2)相容，反之相容方法至少是1阶的.62你现在浏览的是第六十二页，共151页 9.4.2 绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域 定义定义5 5 若一种数值方法在节点值 上大小为 的扰动，于以后各节点值 上产生的偏差均不超过 ，则称该方法是稳定稳定的.以欧拉法为例考察计算稳定性.例例4 4 考察初值问题 其准确解 是一个按指数曲线衰减得很快的函数，如图9-3所示.63你现在浏览的是第六十三页，共151页若取 ，则欧拉公式的具体形式为 计算结果列于表9-4的第2列.可以看到，欧拉方法的解 （图9-3中用号标出）在准确值 的上下波动，计算过程明显地不稳定.图9-3 用欧拉法解方程 得 64你现在浏览的是第六十四页，共151页再考察后退的欧拉方法，取 时计算公式为 计算结果列于表9-4的第3列（图9-3中标以号），这时计算过程是稳定的.但若取 则计算过程稳定.65你现在浏览的是第六十五页，共151页 这表明稳定性不但与方法有关，也与步长 的大小有关，当然也与方程中的 有关.为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，（4.8）其中 为复数.例如在 的邻域，可展开为 模型方程为 对一般方程可以通过局部线性化化为这种形式.66你现在浏览的是第六十六页，共151页略去高阶项，再做变换即可得到 的形式.对于 个方程的方程组，也可线性化为 ，这里 为 的雅可比矩阵 .若 有 个特征值 ，则 还可能是复数，为保证微分方程本身的稳定性，还应假定 .先研究欧拉方法的稳定性.模型方程 的欧拉公式为 所以，为了使模型方程结果能推广到方程组，方程中 应为复数.67你现在浏览的是第六十七页，共151页（4.9）设在节点值 上有一扰动值 ，它的传播使节点值 产生大小为 的扰动值，假设用 按欧拉公式得出 的计算过程不再有新的误差，则扰动值满足 可见扰动值满足原来的差分方程（4.9）.如果差分方程的解是不增长的，即有 则它就是稳定的.68你现在浏览的是第六十八页，共151页即图9-4 显然，为要保证差分方程（4.9）的解不增长，只要选取 充分小，（4.10）在 的复平面上，这是以 为圆心，1为半径的单位圆内部（图9-4）.这个圆域称为欧拉法的绝对稳定域，一般情形可由下面定义.定义定义6 6 单步法（4.1）用于解模型方程（4.8），若得到的解 ，满足 ，则称方法（4.1）是绝对稳定绝对稳定的.（4.9）使（4.1）（4.8）69你现在浏览的是第六十九页，共151页 在 的平面上，使 的变量围成的区域，称为绝对稳定域绝对稳定域，对欧拉法，给出，绝对稳定区间为 .它与实轴的交称为绝对稳定区间绝对稳定区间.其绝对稳定域由 在例5中 ，即 为绝对稳定区间.当取 时 例4中取 故它是不稳定的，它是稳定的.70你现在浏览的是第七十页，共151页故 绝对稳定域由 得到.绝对稳定区间为 ，即 .类似可得三阶及四阶的R-K方法的 分别为 用二阶R-K方法解模型方程 可得到 71你现在浏览的是第七十一页，共151页由 可得到相应的绝对稳定域.当 为实数时则得绝对稳定区间.分别为 三阶显式R-K方法：即 四阶显式R-K方法：即 图9-5给出了R-K方法 到 的绝对稳定域.从以上讨论可知显式的R-K方法的绝对稳定域均为有限域，都对步长 有限制.如果 不在所给的绝对稳定区间内，方法就不稳定.图9-572你现在浏览的是第七十二页，共151页 例例5 5 分别取 ，及 用经典的四阶R-K方法（3.13）计算.解解本例分别为 及 .前者在绝对稳定区间内，后者则不在.经典的四阶R-K方法为（3.13）73你现在浏览的是第七十三页，共151页 以上结果看到，如果步长 不满足绝对稳定条件，误差增长很快.用四阶R-K方法计算其误差见下表：这里74你现在浏览的是第七十四页，共151页 对隐式单步法，可以同样讨论方法的绝对稳定性，例如对后退欧拉法，用它解模型方程可得故 由 可得绝对稳定域为 .它是以 为圆心，1为半径的单位圆外部，故绝对稳定区间为 .75你现在浏览的是第七十五页，共151页 当 时，则 ，即对任何步长均为稳定的.故 对 有 ，故绝对稳定域为 的左半平面，对隐式梯形法，解模型方程 得 76你现在浏览的是第七十六页，共151页绝对稳定区间为 ，即 时梯形法均是稳定的.定义定义7 7 如果数值方法的绝对稳定域包含了 那么称此方法是A-稳定的.隐式欧拉法与梯形方法的绝对稳定域均为 在具体计算中步长 的选择只需考虑计算精度及迭代收敛性要求而不必考虑稳定性，具有这种特点的方法特别重要.由定义知A-稳定方法对步长 没有限制.77你现在浏览的是第七十七页，共151页9.5 线性多步法线性多步法 在逐步推进的求解过程中，计算 之前事实上已经求出了一系列的近似值 ，如果充分利用前面多步的信息来预测 ，则可以期望会获得较高的精度.这就是构造所谓线性多步法的基本思想.本节主要介绍基于泰勒展开的构造方法.构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程 两端积分后利用插值求积公式得到.78你现在浏览的是第七十八页，共151页 9.5.1 线性多步法的一般公式线性多步法的一般公式 一般的线性多步法公式可表示为（5.1）其中 为 的近似，计算时需先给出前面 个近似值 ,再由（5.1）逐次求出 .如果计算 时，除了使用 的值，还用到的值，则称此方法为线性多步法线性多步法.为常数.及 不全为零，则称为线性线性 步法步法.若79你现在浏览的是第七十九页，共151页这时 可直接由（5.1）算出；如果 ，称（5.1）为显式显式 步法步法，如果 ，则（5.1）称为隐式隐式 步法步法，求解时与梯形法（2.7）相同，要用迭代法方可算出 .（5.1）（2.7）80你现在浏览的是第八十页，共151页 定义定义8 8 设 是初值问题(1.1)，(1.2)的准确解，线性多步法（5.1）在 上的局部截断误差为（5.2）（5.1）中系数 及 可根据方法的局部截断误差及阶确定，其定义为：（5.1）（1.1）（1.2）若 ，则称方法（5.1）是 阶的阶的，则称方法（5.1）与方程（1.1）是相容的相容的.81你现在浏览的是第八十一页，共151页 由定义8，对 在 处做泰勒展开，由于 代入（5.2）得（5.3）（5.2）82你现在浏览的是第八十二页，共151页其中（5.4）若在公式（5.1）中选择系数 及 ，使它满足 由定义可知此时所构造的多步法是 阶的.（5.1）83你现在浏览的是第八十三页，共151页（5.5）称右端第一项为局部截断误差主项局部截断误差主项，称为误差常数误差常数.根据相容性定义，,即 ，故方法（5.1）与微分方程（1.1）相容的充分必要条件是（5.6）成立.且（5.6）由（5.4）得（5.1）84你现在浏览的是第八十四页，共151页 当 时，若 ，则由（5.6）可求得 此时公式（5.1）为 即为欧拉法.从（5.4）可求得 ，故方法为1阶精度，这和第2节给出的定义及结果是一致的.且局部截断误差为（5.6）85你现在浏览的是第八十五页，共151页 对 ,若 ，方法为隐式公式.为了确定系数 ，可由 解得于是得到公式 即为梯形法.由（5.4）可求得 ,故 ，所以梯形法是二阶方法，其局部截断误差主项是 .这与第9.2节中的讨论也是一致的.86你现在浏览的是第八十六页，共151页 对 的多步法公式都可利用（5.4）确定系数 ，并由（5.5）给出局部截断误差.（5.5）87你现在浏览的是第八十七页，共151页 9.5.2 阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式 考虑形如（5.7）的 步法，称为阿当姆斯阿当姆斯（Adams)方法方法.为显式方法，为隐式方法，通常称为阿当姆斯显式与隐式公式，也称Adams-Bashforth公式与Adams-Monlton公式.这类公式可直接由对方程 两端从 到 积分求得.88你现在浏览的是第八十八页，共151页 也可以利用(5.4)由 推出，对比 与 可知此时系数 .显然 成立，下面只需确定系数 ，可令 ，求得 .若 ，则令 来求得 .89你现在浏览的是第八十九页，共151页 以 为例，由 ，根据 若 ，则由前三个方程解得得到 的阿当姆斯显式公式是（5.8）90你现在浏览的是第九十页，共151页由（5.4）求得 ，所以（5.8）是三阶方法，局部截断误差是 若 ，则可解得 于是得 的阿当姆斯隐式公式为（5.9）它是四阶方法，局部截断误差是（5.8）91你现在浏览的是第九十一页，共151页（5.10）类似的方法可求得阿当姆斯显式方法和隐式方法的公式，表9-5及表9-6分别列出了 时的阿当姆斯显式公式与阿当姆斯隐式公式，其中 为步数，为方法的阶，为误差常数.92你现在浏览的是第九十二页，共151页93你现在浏览的是第九十三页，共151页94你现在浏览的是第九十四页，共151页 例例6 6 用四阶阿当姆斯显式和隐式方法解初值问题 取步长 .解解本题 .从四阶阿当姆斯隐式公式得到 从四阶阿当姆斯显式公式得到 95你现在浏览的是第九十五页，共151页由此可直接解出 而不用迭代，得到 计算结果见表9-8，其中显式方法中的初值 及隐式方法中的 均用准确解 得到，对一般方程，可用四阶R-K方法计算初始近似.(表略，见教材.)96你现在浏览的是第九十六页，共151页 从以上例子看到同阶的阿当姆斯方法，隐式方法要比显式方法误差小.这可以从两种方法的局部截断误差主项的系数大小得到解释.这里 分别为 及 .97你现在浏览的是第九十七页，共151页 9.5.3 米尔尼方法与辛普森方法米尔尼方法与辛普森方法 考虑另一个 的显式公式 其中 为待定常数.由（5.4）可知 ，再令 得到 可根据使公式的阶尽可能高这一条件来确定其数值.98你现在浏览的是第九十八页，共151页解此方程组得 于是得到四步显式公式（5.11）称为米尔尼米尔尼(Milne）方法方法.由于 ，故方法为4阶，其局部截断误差为（5.12）米尔尼方法也可以通过方程 两端积分 99你现在浏览的是第九十九页，共151页得到.右端积分通过辛普森求积公式就有（5.13）称为辛普森方法辛普森方法.（5.14）它是隐式二步四阶方法，其局部截断误差为 若将方程 从 到 积分，可得 100你现在浏览的是第一百页，共151页 9.5.4 汉明方法汉明方法 辛普森公式是二步方法中阶数最高的，但它的稳定性差，为了改善稳定性，考察另一类三步法公式 其中系数 及 为常数.若希望导出的公式是四阶的，则系数中至少有一个自由参数.若取 ，则可得到辛普森公式.若取 ，仍利用泰勒展开，由（5.4），令101你现在浏览的是第一百零一页，共151页解此方程组得 于是有（5.15）则可得到 102你现在浏览的是第一百零二页，共151页称为汉明汉明(Hamming)方法方法.由于 ，故方法是四阶的，且局部截断误差（5.16）103你现在浏览的是第一百零三页，共151页 9.5.5 预测预测-校正方法校正方法 对于隐式的线性多步法，计算时要进行迭代，计算量较大.为了避免进行迭代，通常采用显式公式给出 的一个初始近似，记为 ，称为预测预测(predictor).接着计算 的值(evaluation)，再用隐式公式计算 ，称为校正校正（corrector）.在（2.13）中用欧拉法做预测，再用梯形法校正，得到改进欧拉法，它就是一个二阶预测-校正方法.（2.13）104你现在浏览的是第一百零四页，共151页 例如用四阶的阿当姆斯显式方法做预测，再用四阶阿当姆斯隐式公式做校正，得到以下格式：一般情况下，预测公式与校正公式都取同阶的显式方法与隐式方法相匹配.预测P：求值E：校正C：求值E：此公式称为阿当姆斯四阶预测阿当姆斯四阶预测-校正格式校正格式（PECE).105你现在浏览的是第一百零五页，共151页 依据四阶阿当姆斯公式的截断误差，对于PECE的预测步P有 对校正步C有 两式相减得 于是有下列事后误差估计 106你现在浏览的是第一百零六页，共151页容易看出 比 更好.（5.17）但在 的表达式中 是未知的，因此计算时用上一步代替，从而构造一种修正预测修正预测-校正格式校正格式(PMECME):107你现在浏览的是第一百零七页，共151页P：M：E：C：M：E：在PMECME格式中已将（5.17）的 及 分别改为 及 .108你现在浏览的是第一百零八页，共151页 利用米尔尼公式（5.11）和汉明公式（5.15）相匹配，并利用截断误差（5.12），（5.16）改进计算结果，可类似地建立四阶修正米尔尼修正米尔尼-汉明预测汉明预测-校正格式校正格式PMECME）：P：M：E：C：M：E：（5.15）（5.11）（5.16）（5.12）109你现在浏览的是第一百零九页，共151页 9.5.6 构造多步法公式的注记和例构造多步法公式的注记和例 构造多步法公式有基于数值积分和泰勒展开两种途径，数值积分方法只能用于可转化为等价的积分方程的微分方程，它是有局限性的.而用泰勒展开则可构造任意多步法公式，其做法是根据多步法公式的形式，直接在 处做泰勒展开.确定多步法（5.1）的系数 及 时不必套用系数公式（5.4），因为多步法公式不一定如（5.1）的形式，而且套用公式容易记错.（5.1）110你现在浏览的是第一百一十页，共151页 例例7 7 解初值问题 用显式二步法 其中 试确定参数 使方法阶数尽可能高,并求局部截断误差.解解 根据局部截断误差定义，用泰勒展开确定参数满足的方程.111你现在浏览的是第一百一十一页，共151页由于 112你现在浏览的是第一百一十二页，共151页为求参数 使方法阶数尽量高，可令 即得方程组 113你现在浏览的是第一百一十三页，共151页解得 ，此时公式为三阶，即为所求局部截断误差.而所得二步法为 而且 114你现在浏览的是第一百一十四页，共151页 例例8 8 证明存在 的一个值，使线性多步法 是四阶的.证明证明只要证明局部截断误差 ，则方法 仍用泰勒展开，由于 为四阶.115你现在浏览的是第一百一十五页，共151页116你现在浏览的是第一百一十六页，共151页当 时，故方法是四阶的.117你现在浏览的是第一百一十七页，共151页9.6 线性多步法的收敛性与稳定性线性多步法的收敛性与稳定性 9.6.1 相容性及收敛性相容性及收敛性 线性多步法(5.1)式的相容性：在定义8中给出的局部截断误差(5.2)中 ，若 称 步法(5.1)与微分方程(1.1)相容，它等价于 对多步法(5.1)可引入多项式（6.1）（6.2）118你现在浏览的是第一百一十八页，共151页和（6.3）分别称为线性多步法(5.1)的第一特征多项式和第二特征多项式.可以看出，如果(5.1)式给定，则 和 也完全确定.反之也成立.119你现在浏览的是第一百一十九页，共151页 定理定理5 5 线性多步法(5.1)式与微分方程(1.1)相容的充分必要条件是 关于多步法(5.1)的收敛性：由于用多步法求数值解需要 个初值，而微分方程(1.1)只给出一个初值 ，因此还要给出 个初值才能用多步法(5.1)进行求解，即（6.4）（6.5）其中 由微分方程初值给定，可由相应单步法给出.120你现在浏览的是第一百二十页，共151页 设由(6.5)式在 处得到的数值解为 ，这里 为固定点，于是有下面定义.定义定义9 9 设初值问题(1.1),(1.2)有精确解 .如果初始条件 满足条件 的线性 步法(6.5)在 处的解 有 则称线性 步法是收敛的.定理定理6 6 设线性多步法(6.5)是收敛的，则它是相容的.此定理的逆定理是不成立的.121你现在浏览的是第一百二十一页，共151页 例例9 9 用线性二步法 解初值问题 解解 此初值问题精确解 ，而由(6.6)式知 故有 ，故方法(6.6)是相容的，但方法(6.6)的解并不收敛，在方法(6.6)中取初值（6.6）（6.7）此时方法(6.6)为二阶差分方程122你现在浏览的是第一百二十二页，共151页其特征方程为解得其根为 及 .于是可求得(6.8)的解为故方法不收敛.多步法是否收敛与 的根有关.（6.8）123你现在浏览的是第一百二十三页，共151页 定义定义1010 如果多步法（5.1）式的第一特征多项式 的根都在单位圆内或圆上，且在单位圆上的根为单根，则称线性多步法（5.1）满足根条件.定理定理7 7 线性多步法是相容的，则线性多步法（6.5）收敛的充分必要条件是线性多步法（5.1）满足根条件.124你现在浏览的是第一百二十四页，共151页 9.6.2 稳定性与绝对稳定性稳定性与绝对稳定性 稳定性主要研究初始条件扰动与差分方程右端项扰动对数值解的影响，假设多步法（6.5）有扰动 ，则经过扰动后的解为 ，它满足方程（6.9）125你现在浏览的是第一百二十五页，共151页 定义定义1111 对初值问题(1.1)，(1.2)，由方法（6.5）得到的差分方程解 ，由于有扰动 ，使得方程（6.9）的解为 ，若存在常数 及 ，使对所有 ，当 时，有则称多步法（5.1）是稳定的或称为零稳定的.研究零稳定性就是研究 时差分方程（6.5）解 的稳定性.它表明当初始扰动或右端扰动不大时，解的误差也不大.126你现在浏览的是第一百二十六页，共151页 对多步法（5.1），当 时对应差分方程的特征方程为 .故有 定理定理8 8 线性多步法（5.1）是稳定的充分必要条件是它满足根条件.关于绝对稳定性只要将多步法（5.1）用于解模型方程（4.8），得到线性差分方程（6.10）127你现在浏览的是第一百二十七页，共151页利用线性多步法的第一、第二特征多项式 ，令此式称为线性多步法的稳定多项式，它是关于 的 次多项式.如果它的所有零点 满足 ，则(6.10)式的解 当 时，有 .（6.11）128你现在浏览的是第一百二十八页，共151页 定义定义1212 对于给定的 ，如果稳定多项式(6.11)的零点满足 ，则称线性多步法(5.1)关于此 值是绝对稳定的.若在 的复平面的某个区域 中所有 值线性多步法(5.1)都是绝对稳定的，而在区域 外，方法是不稳定的，则称 为多步法(5.1)的绝对稳定域绝对稳定域.与实轴的交集称为线性多步法(5.1)的绝对稳定区间绝对稳定区间.为实数时，可以只讨论绝对稳定区间.由于线性多步法的绝对稳定域较为复杂，通常采用根轨迹法.129你现在浏览的是第一百二十九页，共151页 阿当姆斯显式方法和隐式方法的绝对稳定域图形分别为图9-6和图9-7，绝对稳定区间见表9-9.图9-6图9-7130你现在浏览的是第一百三十页，共151页 例例1010 讨论辛普森方法的稳定性 解解 辛普森方法的第第二特征多项式为 的根分别为-1及1，它满足根条件，故方法是零稳定的.但它的稳定多项式为求绝对稳定域 的边界轨迹 .若 ，则可令 ，在 平面域 的边界轨迹 为131你现在浏览的是第一百三十一页，共151页可看出 在虚轴上，且对全部 从而可知 为虚轴上从 到 的线段，故辛普森公式的绝对稳定域为空集.即步长 ，此方法都不是绝对稳定的，故它不能用于求解.132你现在浏览的是第一百三十二页，共151页9.7 一阶方程组与刚性方程一阶方程组与刚性方程 9.7.1 一阶方程组一阶方程组 前面研究了单个方程 的数值解法，只要把 和 理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形.考察一阶方程组 的初值问题，初始条件给为 133你现在浏览的是第一百三十三页，共151页则上述方程组的初值问题可表示为（7.1）求解这一初值问题的四阶龙格-库塔公式为 式中 若采用向量的记号，记134你现在浏览的是第一百三十四页，共151页 考察两个方程的特殊情形：135你现在浏览的是第一百三十五页，共151页这时四阶龙格-库塔公式具有形式 其中（7.2）（7.3）136你现在浏览的是第一百三十六页，共151页 这是一步法，利用节点 上的值 ，由（7.3）式顺序计算 ，然后代入（7.2）式即可求得节点 上的 .（7.3）137你现在浏览的是第一百三十七页，共151页 9.7.2 化高阶方程为一阶方程组化高阶方程为一阶方程组 高阶微分方程（或方程组）的初值问题，原则上总可以归结为一阶方程组来求解.例如，考察下列 阶微分方程（7.4）初始条件为（7.5）只要引进新的变量 138你现在浏览的是第一百三十八页，共151页即可将 阶方程（7.4）化为如下的一阶方程组：（7.6）初始条件（7.5）则相应地化为（7.7）初值问题（7.4）,（7.5）和（7.6）,（7.7）是彼此等价的.（7.4）（7.5）139你现在浏览的是第一百三十九页，共151页 特别地，对于下列二阶方程的初值问题：引进新的变量 ，即可化为下列一阶方程组的初值问题：140你现在浏览的是第一百四十页，共151页针对这个问题应用四阶龙格-库塔公式（7.2），有 由（7.3）式 141你现在浏览的是第一百四十一页，共151页如果消去 ，则上述格式可表示为 这里142你现在浏览的是第一百四十二页，共151页 9.7.3 刚性方程组刚性方程组 在求解方程组（7.1）时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性刚性(stiff)问题问题.考察以下例子.给定系统（7.8）143你现在浏览的是第一百四十三页，共151页的特征值为 方程的准确解为 当 时，称为稳态解，中均含有快变分量 及慢变分量 .它可用解析方法求出准确解，方程右端系数矩阵 144你现在浏览的是第一百四十四页，共151页 对应于 的快速衰减的分量在 秒时已衰减到 .称 为时间常数时间常数.当 时快变分量即可被忽略，而对应于 的慢变分量，它的时间常数 ，它要计算到 时，才能衰减到 ，也就是说解 必须计算到 才能达到稳态解.它表明方程的解分量变化速度相差很大，是一个刚性方程组.145你现在浏览的是第一百四十五页，共151页 若