常微分方程数值解法(IV).ppt

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

1、绪论 在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究微分方程的数值解法是有重要意义的。所谓微分方程数值解法，就是研究利用计算机求解微分方程的近似解的数值方法及相关理论。微分方程数值解法是“信息与计算科学”专业的专业基础课程之一，与数值代数、数值逼近和计算几何统称三大核心课程。第一部分第一部分 常微分方程数值解常微分方程数值解/*Numerical Methods for Ordinary Differe

2、ntial Equations*/待求解的问题待求解的问题：一阶一阶常微分方程的常微分方程的初值问题初值问题/*Initial-Value Problem*/:解的存在唯一性解的存在唯一性（“常微分方程常微分方程”理论）：只要理论）：只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存，即存在与在与 x,y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a,b 上的上的 y1(x)和和 y2(x)都成立，则上述都成立，则上述IVP存在唯一解存在唯一解。解析解法解析解法：（常微分方程理论）：（常微分方程理论）只能求解极少一类

3、常微分方程；实际中给定的问题不一只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。定是解析表达式，而是函数表，无法用解析解法。如何求解如何求解计算解函数计算解函数 y(x)在一系列节点在一系列节点 a=x0 x1 xn=b 处的近似值处的近似值节点间距节点间距 为步长，通常采用为步长，通常采用等距节点等距节点，即取，即取 hi=h(常数常数)。数值解法数值解法:求解所有的常微分方程求解所有的常微分方程步进式步进式：根据已知的或已求出的节点上的函数值计算：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需当前节点上的函数

4、值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。函数值的递推公式即可。-Eulers Method1 欧拉方法欧拉方法/*Eulers Method*/1 Eulers MethodTaylor展开法展开法几何意义几何意义亦称为亦称为欧拉折线法欧拉折线法/*Eulers polygonal arc method*/几何直观是帮助我们寻几何直观是帮助我们寻找解决一个问题的思路找解决一个问题的思路的好办法哦的好办法哦在在假假设设 yn=y(xn)，即即第第 n 步步计计算算是是精精确确的的前前提提下下，

5、考考虑虑公公式式或或方方法法本本身身带带来来的的误误差差:Rn=y(xn+1)yn+1,称称为为局局部部截断误差截断误差/*local truncation error*/。说明 显然，这种近似有一定误差，显然，这种近似有一定误差，而且步长越大，误差越大，而且步长越大，误差越大，如何估计这种误差如何估计这种误差y(xn+1)yn+1？1 Eulers Method截断误差截断误差:实际上，实际上，y(xn)yn，yn 也有误差，它对也有误差，它对yn+1的误的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此

6、称为方法误差或截断误差。方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析局部截断误差的分析：由于假设：由于假设yn=y(xn)，即，即yn准确，因此准确，因此分析局部截断误差时将分析局部截断误差时将y(xn+1)和和 yn+1都用点都用点xn上的信息来表上的信息来表示，工具：示，工具：Taylor展开。展开。欧拉法的局部截断误差：欧拉法的局部截断误差：Rn+1 的的主项主项/*leading term*/1 Eulers Method若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度。阶精度。欧拉法具有欧拉法具有 1 阶精

7、度。阶精度。在在xn点用一阶向前差点用一阶向前差商近似一阶导数商近似一阶导数Eulers method1 Eulers Method1 Eulers Method 欧拉公式的改进欧拉公式的改进：隐式欧拉法或后退欧拉法隐式欧拉法或后退欧拉法/*implicit Euler method or backward Euler method*/xn+1点向后差商近似导数点向后差商近似导数隐式或后退欧拉公式隐式或后退欧拉公式由于未知数由于未知数 yn+1 同时出现在等式的两边，故称为同时出现在等式的两边，故称为隐式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显式显式/*expli

8、cit*/欧拉公欧拉公式。隐式公式不能直接求解，一般需要用式。隐式公式不能直接求解，一般需要用Euler显式公式显式公式得到初值，然后用得到初值，然后用Euler隐式公式迭代求解。因此隐式公隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。式较显式公式计算复杂，但稳定性好（后面分析）。收敛性收敛性1 Eulers Method 见上图，见上图，显然，这种近似也有一定误差，显然，这种近似也有一定误差，如何估计这种误差如何估计这种误差y(xn+1)yn+1？方法同上，基于方法同上，基于Taylor展开估计局部截断误差。展开估计局部截断误差。但是注意，隐式公式中右边含有但是注

9、意，隐式公式中右边含有f(xn+1,yn+1),由于由于yn+1不准确，所以不能直接用不准确，所以不能直接用y(xn+1)代替代替f(xn+1,yn+1)设已知曲线上一点设已知曲线上一点 Pn(xn,yn),过过该点作弦线，斜率为该点作弦线，斜率为(xn+1,yn+1)点的方向场点的方向场f(x,y)方向方向,若步长若步长h充分充分小，可用弦线和垂线小，可用弦线和垂线x=xn+1的交点的交点近似曲线与垂线的交点。近似曲线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)1 Eulers Method 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：1 Eulers Method1

10、Eulers Method 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。1 Eulers Method比较尤拉显式公式和隐式公式及其局部截断误差显式公式隐式公式1 Eulers Method 若将这两种方法进行算术平均，即可消除误差若将这两种方法进行算术平均，即可消除误差的主要部分的主要部分/*leading term*/而获得更高的精度而获得更高的精度,称为梯形法称为梯形法1 Eulers Method 梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法的显、隐式两种算法的平均平均注：注：的确有局部截断

11、误差的确有局部截断误差 ，即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。但阶精度，比欧拉方法有了进步。但注意到该公式是注意到该公式是隐式隐式公式，计算时不得不用到迭代法，公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。其迭代收敛性与欧拉公式相似。梯形法的迭代计算和收敛性梯形法的迭代计算和收敛性收敛性收敛性1 Eulers Method梯形法的简化计算梯形法的简化计算 迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预次就转入下一点的

12、计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法测校正系统，也称为改进尤拉法：改进欧拉法改进欧拉法/*modified Eulers method*/Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1nnnnyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+ny),(),(2111+=nnnnnnyxfyxfhyy1 Eulers Method注：注：此法亦称为此法亦称为预测预测-校正法校正法/*predictor-corrector method*/。可以证明该算法具有可

13、以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。1 Eulers Method其它形式其它形式几何解释几何解释xnxn+1ABPn+1=(A+B)/2尤拉法尤拉法后退尤拉法后退尤拉法梯形法梯形法1 Eulers Method令令x=x1,得得Another point of view对右端积分采用左矩形、右矩形、梯形积分公式，即对右端积分采用左矩形、右矩形、梯形积分公式，即可得尤拉显式、隐式、梯形公式可得尤

14、拉显式、隐式、梯形公式1 Eulers Method 中点欧拉公式中点欧拉公式/*midpoint formula*/中心差商近似导数中心差商近似导数x0 x2x1假设假设 ，则可以导出，则可以导出即中点公式也具有即中点公式也具有 2 阶精度，且是显式的。阶精度，且是显式的。需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法/*double-step method*/，而前面的三种算法都是，而前面的三种算法都是单步法单步法/*single-step method*/。1 Eulers Method几何解释几何解释xnxn+1尤拉法尤拉

15、法后退尤拉法后退尤拉法中点法中点法xn-1令令x=x2,得得Another point of view对右端积分采用中矩形公式即得中点公式对右端积分采用中矩形公式即得中点公式1 Eulers Method 预测预测-校正校正-改进系统改进系统中点法具有二阶精度，且是显式的，与梯形公式精度相匹配，用中点公式作中点法具有二阶精度，且是显式的，与梯形公式精度相匹配，用中点公式作预测，梯形公式作校正，得到如下预测校正系统预测，梯形公式作校正，得到如下预测校正系统：校正误差约为预测误校正误差约为预测误差的差的1/41 Eulers Method预测误差和校正误差预测误差和校正误差的事后误差估计式的事后误

16、差估计式利用上两式可以估计预测值和校正值与准确值的误差，可以利用上两式可以估计预测值和校正值与准确值的误差，可以期望，利用这两个误差分别作预测值和校正值的补偿，有可期望，利用这两个误差分别作预测值和校正值的补偿，有可能提高精度。能提高精度。设设pn,cn分别为第分别为第n步的预测值和校正值，即步的预测值和校正值，即此时此时cn+1未知，未知，故用故用pn-cn代替代替1 Eulers Method 预测预测-校正校正-改进公式改进公式注：利用该算法计算注：利用该算法计算yn+1时，需要时，需要1 Eulers Method公式公式局部截断误差局部截断误差精精度度显显隐隐稳稳定定性性步数步数尤拉

17、显尤拉显式公式式公式1 1阶阶显显差差单步单步尤拉隐尤拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形公梯形公式式2 2阶阶隐隐差差单步单步中点法中点法2 2阶阶显显好好二步二步summary两个预测两个预测-校正系校正系统统尤拉两步法和梯形公式构尤拉两步法和梯形公式构成的预测成的预测-校正校正-改进系统改进系统尤拉公式和梯形公式构成尤拉公式和梯形公式构成的预测的预测-校正系统校正系统例例 HW:p.201#1-5证明中点法和梯形公式的精度为证明中点法和梯形公式的精度为2阶阶2 龙格龙格-库塔法库塔法/*Runge-Kutta Method*/建立高精度的单步递推格式建立高精度的单步递推格式:在

18、改进尤拉法和尤拉在改进尤拉法和尤拉两步法预测两步法预测-校正系统中校正系统中,预测公式都是单步法预测公式都是单步法,如果预如果预测误差很小测误差很小,则通过校正后得到的近似值误差会更小则通过校正后得到的近似值误差会更小,因因此需要研究高精度的单步法此需要研究高精度的单步法.1.Taylor级数法级数法IVP:设其解为设其解为y=y(x)由由Taylor展开，有展开，有(1)2 Runge-Kutta Method(2)2 Runge-Kutta Method要使公式具有要使公式具有p阶精度，则在阶精度，则在(1)式中截取前式中截取前p+1项，项，用用(2)式计算各阶导数，即得下面式计算各阶导数

19、，即得下面Taylor公式：公式：局部截断误差局部截断误差(3)2 Runge-Kutta Method2.Taylor公式公式(3)表面上看形式简单，但具体构造时表面上看形式简单，但具体构造时往往很困难，因为按往往很困难，因为按(2)式求导，这一过程可能很式求导，这一过程可能很复杂。因此通常不直接用复杂。因此通常不直接用Taylor公式，而借鉴其思公式，而借鉴其思想提出其它公式。想提出其它公式。1.由此看出，一种方法具有由此看出，一种方法具有p阶精度阶精度公式对不超公式对不超过过p次的多项式准确成立（局部截断误差为次的多项式准确成立（局部截断误差为0）。）。这一等价条件也可以用来判断一种方法

20、的精度。这一等价条件也可以用来判断一种方法的精度。2 Runge-Kutta Method单步递推法的单步递推法的基本思想基本思想是从是从(xn,yn)点出发，以点出发，以某一某一斜率斜率沿直线达到沿直线达到(xn+1,yn+1)点。欧拉法及其各种变形点。欧拉法及其各种变形所能达到的最高精度为所能达到的最高精度为2阶阶。2.RungeKutta Method由微分中值定理，有由微分中值定理，有k*称为区间称为区间xn,xn+1上的上的平均斜率平均斜率，只要知道平均，只要知道平均斜率，就可计算斜率，就可计算y(xn+1).因此只要对平均斜率提供一因此只要对平均斜率提供一种近似算法，则由种近似算法

21、，则由(4)式可导出一种相应的求解公式。式可导出一种相应的求解公式。(4)2 Runge-Kutta Method例例2 Runge-Kutta Method 由此看出，改进的尤拉公式用由此看出，改进的尤拉公式用xn与与xn+1两个节点两个节点的斜率的算术平均作为平均斜率，的斜率的算术平均作为平均斜率，xn+1点的斜率通点的斜率通过已知信息过已知信息yn来预测。来预测。考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1，K2 的的平均值平均值吗？吗？步长一定是步长一定是h 吗？即吗？即第二个节点一定是第二个节点一定是xn+1吗？吗？2 Runge-Ku

22、tta Method2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入第代入第1式，得到式，得到 2阶阶RungeKutta Method2 Runge-Kutta MethodStep 3:将将 yi+1 与与 y(xi+

23、1)在在 xi 点的点的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式塔格式。注意到，注意到，就是改进的欧拉法；就是改进的欧拉法；p=1/2,1=0,2=1,变形尤拉公式变形尤拉公式。Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？改进的改进的Euler 公式推广为二阶公式推广为二阶Runge-Kutta公式公式带来这样的启示：带来这样的启示：若在若在xn,xn+1上多预测几个点的斜率值，

24、然后将上多预测几个点的斜率值，然后将它们的算术平均作为平均斜率，则有可能构造出它们的算术平均作为平均斜率，则有可能构造出具有更高精度的计算公式。具有更高精度的计算公式。-Runge-Kutta方法的基本思想。方法的基本思想。注：二阶注：二阶Runge-Kutta公式用多算一次函数值公式用多算一次函数值f 的的办法避开了二阶办法避开了二阶Taylor级数法所要计算的级数法所要计算的f 的导数。的导数。在这种意义上，可以说在这种意义上，可以说Runge-Kutta方法实质上是方法实质上是Taylor级数法的变形。级数法的变形。2 Runge-Kutta Method其中其中 i (i=1,m)，i

25、 (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待定系数，均为待定系数，确定这些系数的步骤与前面相似。确定这些系数的步骤与前面相似。2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 高阶高阶RungeKutta Method Gill公式：公式：4阶经典龙格阶经典龙格-库塔公式的一种改进库塔公式的一种改进2 Runge-Kutta Method 最常用为四级最常用为四级4阶阶经典龙格

26、经典龙格-库塔法库塔法/*Classical Runge-Kutta Method*/：2 Runge-Kutta Method注：注：龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。

27、对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。2 Runge-Kutta Method 变步长的变步长的RungeKutta MethodQ:由局部截断误差可以看出，步长由局部截断误差可以看出，步长 h 越小，局部截断越小，局部截断误差越小；但步长减小，在一定求解范围（区间）内误差越小；但步长减小，在一定求解范围（区间）内要完成的步数就增加了，步数增加会引起计算量增大，要完成的步数就增加了，步数增加会引起计算量增大，导致舍入误差积累。因此要选取适当的步长。导致舍入误差积累。因此要选取适当的步长。选择步长时要考虑两个问题：选择步长时要考虑两个问题：1.如何衡量和

28、检验计算结果的精度？如何衡量和检验计算结果的精度？2.如何根据所获得的精度处理步长？如何根据所获得的精度处理步长？HW:p.201#6-83 单步法的单步法的收敛性与稳定性收敛性与稳定性 /*Convergency and Stability*/前面介绍了两大类微分方程数值解法：一类是前面介绍了两大类微分方程数值解法：一类是用差商近似导数得到的尤拉系列公式，另一类是用差商近似导数得到的尤拉系列公式，另一类是基于平均斜率概念的基于平均斜率概念的RungeKutta公公式。式。基本思基本思想都是通过某种离散化手续，将微分方程转化为想都是通过某种离散化手续，将微分方程转化为差分方程（代数方程）来求解

29、差分方程（代数方程）来求解。Q1.这种转化是否合理？要看差分问题的解这种转化是否合理？要看差分问题的解yn当当h0时是否收敛到微分方程的解时是否收敛到微分方程的解y(xn),即是否成立即是否成立 yn y(xn)，h0.-收敛性问题收敛性问题 Q2.实际计算时，由于舍入误差的影响，差分方实际计算时，由于舍入误差的影响，差分方程的解本身也有误差，这种误差在计算过程中会程的解本身也有误差，这种误差在计算过程中会不会扩大不会扩大 -稳定性问题稳定性问题3 Convergency and Stability 收敛性收敛性/*Convergency*/若若某某算算法法对对于于任任意意固固定定的的 x=x

30、i=x0+i h，当当 h0(同时同时 i )时有时有 yi y(xi)，则称该算法是，则称该算法是收敛收敛的。的。例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 3 Convergency and Stability显式单步法的收敛性显式单步法的收敛性(1)3 Convergency and Stability而整体截断误差为而整体截断误差为证明证明:设设 表示当表示当yn=y(xn)时，时，由公式由公式(1)求得的求得的结果，即结果，即则局部截

31、断误差为则局部截断误差为(2)(2)-(1),得得3 Convergency and Stability从而从而3 Convergency and Stability3 Convergency and Stability Euler法的收敛性法的收敛性：(x,y)=f(x,y),故当故当f(x,y)满足满足Lipschitz条件时，条件时，尤拉法收敛；尤拉法收敛；3 Convergency and Stability3 Convergency and Stability 稳定性稳定性/*Stability*/例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、

32、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10

33、 46.1442 10 63.0590 10 7What is wrong?!An Engineer complains:Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!3 Convergency and Stability3 Convergency and Stability若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是

34、绝绝对对稳稳定定的的/*absolutely stable*/。一般分析某算法的稳定性时，为简单起见，只考虑一般分析某算法的稳定性时，为简单起见，只考虑模型方程模型方程或或试验方程试验方程/*test equation*/常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该算法，则若此误差以后逐步衰减，就称该算法相对于相对于 绝对稳定绝对稳定，的全体构成的全体构成绝对稳定区域绝对稳定区域。我们。我们称称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝

35、对稳定区域比的绝对稳定区域比 B 的的大大。h h=h3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。一般来说，隐式欧拉法的绝

36、对稳定性比同阶的显式法的好。小结小结 尤拉两步公式与尤拉单步公式相比，使用两尤拉两步公式与尤拉单步公式相比，使用两个节点上的已知信息将精度提高一阶。可以设想，个节点上的已知信息将精度提高一阶。可以设想，计算计算y(xn+1)时，充分利用前面已经求出的节点上时，充分利用前面已经求出的节点上的的 y 及及 y 值的值的线性组合线性组合来近似来近似y(xn+1)，精度会大，精度会大大提高。大提高。-线性多步法的基本思想线性多步法的基本思想。构造线性多步法有多种途径，这里介绍两种：构造线性多步法有多种途径，这里介绍两种：基于数值积分的构造方法基于数值积分的构造方法；基于基于TaylorTaylor展开

37、的构造方法展开的构造方法。4 线性多步法线性多步法 /*Multistep Method*/).(.110111101rnrnnnrnrnnnffffhyyyy +=线性多步法的通式可写为：线性多步法的通式可写为：当当 1 0 时，为时，为隐式公式隐式公式;1=0 则为则为显式公式显式公式。基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通过，则可通过 近似近似y(xn+1)。而。而选用不同近似式选用不同近似式 I，可得到不同的计算公式，可得到不同的计算公式。例如利用左矩形积分公式得到尤拉公式；梯形积分公式得例

38、如利用左矩形积分公式得到尤拉公式；梯形积分公式得到梯形公式。到梯形公式。一般地，利用插值原理所建立的一系列数值积分一般地，利用插值原理所建立的一系列数值积分方法也可以导出解微分方程的一系列计算公式。方法也可以导出解微分方程的一系列计算公式。运用插值方法的关键在于选取合适的插值节点。运用插值方法的关键在于选取合适的插值节点。假设已构造出假设已构造出f(x,y(x)的插值多项式的插值多项式Pr(x),则则4 Multistep Method 亚当姆斯显式公式亚当姆斯显式公式/*Adams explicit formulae*/利用利用r+1 个节点上的被积函数值个节点上的被积函数值构造构造 r 阶

39、牛顿阶牛顿后插后插多项式多项式,有有Newton插值余项插值余项/*Adams 显式公式显式公式*/外推技术外推技术/*extrapolation*/table 5-6,p.181jj0111/225/1233/8实际计算时，将差分展开实际计算时，将差分展开 ijrj=i局部截断误差为：局部截断误差为：注：注：Br 与与yn+1 计算公式中计算公式中 fn,fn k 各项的各项的系数系数 均可查表均可查表得到得到。见下表。见下表。10123rfnfn 1fn 2fn 3Br例：例：常用的是常用的是 r=3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式Adams显式公式用显式公式用 作为插值节点，

40、作为插值节点，在求积区间在求积区间xn,xn+1上用插值函数上用插值函数Nr(x)近似近似f(x,y(x)，而，而xn+1不在插值节点内，因此是一个不在插值节点内，因此是一个外推外推的过程。虽然公式是显式的，便于计算，但效果的过程。虽然公式是显式的，便于计算，但效果并不理想，比如稳定性较差等。并不理想，比如稳定性较差等。因此改用通过因此改用通过r+1个节点个节点的插值多项式的插值多项式Nr(x)近似近似f(x,y(x)，由于，由于xn+1是其中是其中一个插值点，因此是一个插值点，因此是内插多项式内插多项式，但导出的公式，但导出的公式是是隐式的隐式的。4 Multistep Method 亚当姆

41、斯隐式公式亚当姆斯隐式公式/*Adams implicit formulae*/利用利用r+1 个节点上的被积函数值个节点上的被积函数值 fn+1,fn,fn r+1 构造构造r 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项式完全类似地可得到一系列隐式公式。隐式公式。j0 11-1/22-1/123-1/24局部截断误差为：局部截断误差为：10123kfi+1fifi 1fi 2Br小于小于Br注：注：与与yn+1 计算公式中计算公式中 fn+1,fn+1 k 各项的各项的系数系数 均均可查表得到可查表得到。见下表。见下表。常用的是常用的是 k=3 的的4阶

42、亚当姆斯隐式公式阶亚当姆斯隐式公式较同阶显式较同阶显式稳定稳定例：例：4 Multistep Method 亚当姆斯预测亚当姆斯预测-校正系统校正系统 /*Adams predictor-corrector system*/Step 1:用用Runge-Kutta 法法计算前计算前 k 个初值；个初值；Step 2:用用Adams 显式显式计算计算预测预测值；值；Step 3:用同阶用同阶Adams 隐式隐式计算计算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必须相同。式的精度必须相同。通常用通常用经典经典Runge-Kutta 法法配合配合4阶阶Adams 公式公式。4阶阶Ada

43、ms隐式公式的截断误差为隐式公式的截断误差为4阶阶Adams显式公式的截断误差为显式公式的截断误差为Predicted value pn+1Corrected value cn+1Modified value mn+1预测值与校正预测值与校正值的事后误差值的事后误差估计估计Adams显隐预测显隐预测-校正校正-改进系统改进系统4 Multistep Method Adams 4th-Order predictor-corrector AlgorithmTo approximate the the solution of the initial-value problemAt(N+1)equal

44、ly spaced numbers in the interval a,b.Input:endpoints a,b;integer N;initial value y0.Output:approximation y at the(N+1)values of x.Step 1 Set h=(b a)/N;x0=a;y0=y0;Output(x0,y0);Step 2 For i=1,2,3 Compute yi using classical Runge-Kutta method;Output(xi,yi);Step 3 For i=4,N do steps 4-10Step 5 ;/*pred

45、ict*/Step 6 ;/*modify*/Step 7 ;/*correct*/Step 8 ;/*modify the final value*/Step 9 Output(xi+1,yi+1);Step 10 For j=0,1,2,3 Set xi =xi+1;yi =yi+1;/*Prepare for next iteration*/Step 11 STOP.应为应为(ci+1 pi+1),但因但因ci+1 尚未尚未算出，只好用算出，只好用(ci pi)取代之。取代之。4 Multistep Method 基于泰勒展开的构造法基于泰勒展开的构造法).(.110111101rnrn

46、nnrnrnnnffffhyyyy +=将通式中的右端各项将通式中的右端各项 yn 1,yn k;fn+1,fn 1,fn k(yn+1,yn k)分别在分别在 xn 点作点作泰勒展开泰勒展开，与精确解，与精确解 y(xn+1)在在 xn 点的泰勒展开作点的泰勒展开作比较比较。通过令。通过令同类项系数相等同类项系数相等，得到足以确定待，得到足以确定待定系数定系数 0,r;1,0,r 的等式，则可构的等式，则可构造出线性多步法的公式。造出线性多步法的公式。当当 1 0 时，为时，为隐式公式隐式公式;1=0 则为则为显式公式显式公式。局部截断误差为局部截断误差为例：例：设设)(3322110221

47、101 +=nnnnnnnnyyyyhyyyy 确定式中待定系数确定式中待定系数 0,1,2,0,1,2,3,使得公式具有使得公式具有4阶精度。阶精度。4 Multistep Method解：解：/*y(xn)=yn*/个未知数个未知数个方程个方程754 Multistep Method 令令 1=2=0Adams 显式显式公式公式 以以 y n+1 取代取代 y n 3，并取，并取 1=2=0Adams 隐式隐式公式公式 以以 yn 3 取代取代 y n 3，则可导出另一组，则可导出另一组4 阶显式算法，其阶显式算法，其中包含了著名的中包含了著名的米尔尼米尔尼/*Milne*/公式公式(4步

48、步4阶显式公式阶显式公式)其局部截断误差为其局部截断误差为注：注：上式也可通过上式也可通过数值积分数值积分导出，即将导出，即将 在区间在区间 上积分，得到上积分，得到 再过再过 做做 f 的的2次插值多项式即可。次插值多项式即可。取取 1=1,2=0得到得到辛甫生辛甫生/*Simpson*/公式公式与与Milne 公式匹配使用公式匹配使用,构成预测构成预测-校正系统校正系统F辛甫生辛甫生/*Simpson*/公公式式（2步步4阶）阶）在区间在区间xn 1,xn+1上积分，并上积分，并用用Simpson数值积分数值积分公式来近似公式来近似积分项，亦可得此积分项，亦可得此Simpson公式。公式。

49、其局部截断误差为其局部截断误差为Milne-Simpson预测校正系统的局部截断误差预测校正系统的局部截断误差(14/450.31,-1/90 -0.01)都比同价的都比同价的Adams公式公式(251/720 0.34,-19/720 -0.02)小，计算量又比小，计算量又比4阶阶Runge-Kutta公式少。缺点是公式少。缺点是,校正公式是弱稳定的，校正公式是弱稳定的，即对某些问题来说，其稳定性不好。即对某些问题来说，其稳定性不好。Hamming对对Simpson公式的稳定性作了改进，得到公式的稳定性作了改进，得到Hamming公式。公式。4 Multistep Method Milne-

50、Simpson 系统的缺点是系统的缺点是稳定性差稳定性差，为改善稳定性，为改善稳定性，考虑另一种隐式校正公式：考虑另一种隐式校正公式：要求公式具有要求公式具有4 阶精度。通过泰勒展开，可得到阶精度。通过泰勒展开，可得到 个等式，个等式，从中解出从中解出 个未知数，则有个未知数，则有 个自由度。个自由度。561取取 1=1 得得Simpson 公式公式哈明哈明/*Hamming*/用用 1 的不同数值进行试验，发现当的不同数值进行试验，发现当 1=0 时，公式的稳定性较好，即：时，公式的稳定性较好，即：其局部截断误差为其局部截断误差为注：注：哈明公式不能用数值积分方法推导出来。哈明公式不能用数值