第9章 常微分方程数值解法ppt课件

1、第第9章章 常微分方程数值解常微分方程数值解/* Numerical Methods for Ordinary Differential Equations */ 待求解的问题：一阶常微分方程的初值问题待求解的问题：一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy 解的存在唯一性（解的存在唯一性（“常微分方程常微分方程理论）：只要理论）：只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件，即存条件，即存在与在与 x, y 无关的常数无关的常数 L 使使对任意

2、定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，则上述IVP存在唯一解。存在唯一解。| ),(),(|2121yyLyxfyxf 解析解法：（常微分方程理论）解析解法：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。定是解析表达式，而是函数表，无法用解析解法。如何求解如何求解计算解函数计算解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值),., 1()(nixyyii 节点间距节点间距 为步长，

3、通常采用等距节点，为步长，通常采用等距节点，即取即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii数值解法数值解法: 求解所有的常微分方程求解所有的常微分方程步进式：根据已知的或已求出的节点上的函数值计算步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。函数值的递推公式即可。111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f

4、xy1(,) 0, 1,.nnnnyyh f xyn-Eulers Method1 欧拉方法欧拉方法 /* Eulers Method */1 Eulers MethodTaylor展开法展开法几何意义几何意义亦称为欧拉折线法亦称为欧拉折线法 /* Eulers polygonal arc method*/ 几何直观是帮助我们寻几何直观是帮助我们寻找解决一个问题的思路找解决一个问题的思路的好办法哦的好办法哦定义定义在假设在假设 yn = y(xn)，即第，即第 n 步计算是精确的前提下，步计算是精确的前提下，考虑公式或方法本身带来的误差考虑公式或方法本身带来的误差: Rn = y(xn+1)

5、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展开。展开。 欧拉法的局部截断误欧拉法的局部截断误差：差：223111232() ( )( )( )( ) ( ,) ( )( )hnnnnnnnnnhnRy xyy xhy xy xO hyhf x yy xO hRn+1 的主项的主项/* leading term */1 Eulers Method定义定义若某算法的局部截

7、断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 欧拉法具有欧拉法具有 1 阶阶精度。精度。在在xn点用一阶向前差点用一阶向前差商近似一阶导数商近似一阶导数1()()()nnny xy xyxh 在第在第6章讨论牛顿插值公式时章讨论牛顿插值公式时介绍了差商的概念和性质，介绍了差商的概念和性质，各阶差商可以近似各阶导数，具有不同的精度，各阶差商可以近似各阶导数，具有不同的精度，且可以用函数值来表示。且可以用函数值来表示。上一章中数值微分的方法之一上一章中数值微分的方法之一就是用差商近似导数就是用差商近似导数111()()() ()()(,)nnnnn

8、nnnnny xy xhy xy xyy xyyh f xyEulers method1 Eulers Method1 Eulers Method 欧拉公式的改进：欧拉公式的改进：隐式欧拉法或后退欧拉法隐式欧拉法或后退欧拉法 /* implicit Euler method or backward Euler method*/11()()()nnny xy xy xhxn+1点向后差商近似导点向后差商近似导数数111111()()() ()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy隐式或后退欧拉公式隐式或后退欧拉公式由于未知数由于未知数 yn+1 同时出现

9、在等式的两边，故称为隐式同时出现在等式的两边，故称为隐式 /* implicit */ 欧拉公式，而前者称为显式欧拉公式，而前者称为显式 /* explicit */ 欧欧拉公式。隐式公式不能直接求解，一般需要用拉公式。隐式公式不能直接求解，一般需要用Euler显式显式公式得到初值，然后用公式得到初值，然后用Euler隐式公式迭代求解。因此隐隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好后面分析）。式公式较显式公式计算复杂，但稳定性好后面分析）。01(1)( )111(,)(,)nnnnkknnnnyyh f xyyyh f xy(1)()1111111()(0 )1111(1)

10、11111()1(,)(,) 1, ()(,)kknnnnnnkknnnnknnnnnnknyyh fxyfxyhL yyhLyyhLyykyyh fxyy 在 迭 代 公 式 中 取 极 限 ， 有因 此的 极 限 就 是 隐 式 方 程 的 解收敛性收敛性1 Eulers Method 见上图，见上图， 显然，这种近似也有一定误差，显然，这种近似也有一定误差，如何估计这种误差如何估计这种误差y(xn+1) yn+1 ？方法同上，基于方法同上，基于Taylor展开估计局部截断误差。展开估计局部截断误差。但是注意，隐式公式中右边含有但是注意，隐式公式中右边含有f(xn+1 , yn +1 )

11、,由于由于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 隐式欧拉法的局部截隐式欧拉法的局部截断误差：断误差：11111111121111111321, ,2

12、 , 2nnnnynnnnnnnnnnnnynnnnnnnnf xyf xy xfxyy xyy xhf xy xyxyxhyxyxyhfxyy xy xhhyxh yxyxy xy x由微分中值定理，得在，之间；又而 2326nnnnhhhyxyxyx1 Eulers Method111111232311(), 231,23nnnynnnnnynnnnRy xyhfxy xyhhyxyxhhhfxRyxyx 从而即2211121,1,(1)1ynynynhfxhfxhfxxxx 111 Eulers Method 2322111231221 1,23 3,226 ( )2nynynnnnyn

13、nnnnhhRhfxhfxy xyxhhy xfxy xy xhRy xo h 隐式欧拉法的局部截断误隐式欧拉法的局部截断误差：差：111()nnnRy xy232()()hny xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。1 Eulers Method1(,) 0, 1,.nnnnyyh f xyn比较尤拉显式公式和隐式公式及其局部截断误差231112()()()hnnnnRy xyy xO h显式公式111(,)nnnnyyh f xy隐式公式231112()()()hnnnnRy xyy xO h1 Eulers Method 若将这两种方法进行算术平均，即可消除误

14、差若将这两种方法进行算术平均，即可消除误差的主要部分的主要部分/*leading term*/而获得更高的精度而获得更高的精度,称为梯形法称为梯形法1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的平均显、隐式两种算法的平均111 (,)(,)2nnnnnnhyyf xyf xy注：的确有局部截断误差注：的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。但阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公

15、式相似。法，其迭代收敛性与欧拉公式相似。3111()()nnnRy xyO h梯形法的迭代计算和收敛性梯形法的迭代计算和收敛性01(1)( )111(,)(,)(,)2nnnnkknnnnnnyyh f xyhyyf xyf xy(1)()1111111()(0 )1111(1)11111()1(,)(,)2222 1, ()2(,)(,)2kknnnnnnkknnnnknnnnnnnnknhyyfxyfxyhhLL yyyyhLhyykLhyyfxyfxyy 当 时 ，在 迭 代 公 式 中 取 极 限 ， 有因 此的 极 限 就 是 隐 式 方 程 的 解收敛性收敛性1 Eulers Me

16、thod梯形法的简化计算梯形法的简化计算 迭代计算量大，且难以预测迭代次数。为了控制计迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法：统，也称为改进尤拉法： 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出),(1nnnnyxfhyy Step 2: 再将再将 代入隐式梯形公式的右边作校正，得到

17、代入隐式梯形公式的右边作校正，得到1 ny),(),(2111 nnnnnnyxfyxfhyy11( ,),( ,)2nnnnnnnnhyyf x yf xyh f x y1 Eulers Method注：此法亦称为预测注：此法亦称为预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有。可以证明该算法具有 2 阶精度，同时可以看到它是阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。面将看到，它的稳定性高于显式欧拉法。1 Eulers M

18、ethod+1+11+1 (,)(,(,)21 (,)(,)2nnnnnnnnnnnnnnnhyyf xyf xh yh f xyyyhf xyyhf xy或其它形式其它形式1+1(,) (,)12pnnncnnpnpcyyh f xyyyh f xyyyy或几何解释几何解释xnxn+1ABPn+1=(A+B)/2尤拉法尤拉法后退尤拉法后退尤拉法梯形法梯形法1 Eulers Method 0)(,),(yaybaxyxfdxdy 00( ),xxy xyft y tdt令令x=x1,得得 1010(),xxy xyft y tdtAnother point of view对右端积分采用左矩形、

19、右矩形、梯形积分公式，即对右端积分采用左矩形、右矩形、梯形积分公式，即可得尤拉显式、隐式、梯形公式可得尤拉显式、隐式、梯形公式1 Eulers Method 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点公式也具有即中点公式也具有 2 阶精度，且是显式的。阶精度，且是显式的。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2个

20、初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为双步法过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步法，而前面的三种算法都是单步法 /* single-step method */。1 Eulers Method几何解释几何解释xn xn+1尤拉法尤拉法后退尤拉法后退尤拉法中点法中点法xn-1 0)(,),(yaybaxyxfdxdy 00( ),xxy xyft y tdt令令x=x2,得得 2020(),xxy xyft y tdtAnother point of view对右端积分采用中矩形公式即得中点公式对右端积

21、分采用中矩形公式即得中点公式1 Eulers Method 预测预测-校正校正-改进系统改进系统中点法具有二阶精度，且是显式的，与梯形公中点法具有二阶精度，且是显式的，与梯形公式精度相匹配，用中点公式作预测，梯形公式式精度相匹配，用中点公式作预测，梯形公式作校正，得到如下预测校正系统：作校正，得到如下预测校正系统：111112(,)(,)(,)2nnnnnnnnnnyyh f xyhyyf xyf xy预测：校正： 3(3)1113(3)111312nnnnnnnnnnhyyy xyyxhyyy xyyx 预测误差(设, 准确): 校正误差(设,准确): 111114nnnny xyy xy

22、校正误差约为预测误校正误差约为预测误差的差的1/41 Eulers Method11111111111454515nnnnnnnnnnny xyyy xyyyy xyyy 预测误差和校正误差预测误差和校正误差的事后误差估计式的事后误差估计式利用上两式可以估计预测值和校正值与准确值的误差，可以利用上两式可以估计预测值和校正值与准确值的误差，可以期望，利用这两个误差分别作预测值和校正值的补偿，有可期望，利用这两个误差分别作预测值和校正值的补偿，有可能提高精度。能提高精度。 设设pn,cn分别为第分别为第n步的预测值和校正值，即步的预测值和校正值，即,nnnnpy cy1111111145 15nn

23、nnnnnnmppcycpc改进：此时此时cn+1未知，未知，故用故用pn -cn代替代替1 Eulers Method 预测预测-校正校正-改进公改进公式式1111111111111111245,215,nnnnnnnnnnnnnnnnnnnnnpyhymppcmf xmhcyymycpcyf xy预测：改进：计算：校正：改进：计算：注：利用该算法计算注：利用该算法计算yn+1时，需要时，需要11111110nnnnnyypcyypcypc,和，因此启动算法之前必须给出开始值 和，可用其它单步法计算，一般取为 。1 Eulers Method公式公式局部截断误差局部截断误差精精度度显显隐隐稳

24、稳定定性性步数步数尤拉显尤拉显式公式式公式1 1阶阶显显差差单步单步尤拉隐尤拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形公梯形公式式2 2阶阶隐隐差差单步单步中点法中点法2 2阶阶显显好好二步二步3(3)3nhyx3(3)12nhyx2(2)2nhyx2(2)2nhyxsummary两个预测两个预测-校正系统校正系统1111111111111111245,215,nnnnnnnnnnnnnnnnnnnnnpyhymppcmf xmhcyymycpcyf xy预测：改进：计算：校正：改进：计算：尤拉两步法和梯形公式尤拉两步法和梯形公式构成的预测构成的预测-校正校正-改进系改进系统统111111

25、,2nnnnnnnnnnyyhyyf xyhyyyy预测：计算：校正：尤拉公式和梯形公式构成尤拉公式和梯形公式构成的预测的预测-校正系统校正系统例例 HW: p.201 #1-5证明中点法和梯形公式的精度为证明中点法和梯形公式的精度为2阶阶2 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */建立高精度的单步递推格式建立高精度的单步递推格式: :在改进尤拉法和尤在改进尤拉法和尤拉两步法预测拉两步法预测- -校正系统中校正系统中, ,预测公式都是单步预测公式都是单步法法, ,如果预测误差很小如果预测误差很小, ,则通过校正后得到的近则通过校正后得到的近似值误差会更小似值

26、误差会更小, ,因此需要研究高精度的单步法因此需要研究高精度的单步法. . 1. Taylor级数法级数法 0)(,),(yaybaxyxfdxdyIVP:设其解为设其解为y=y(x)23126nnnnnhhy xy xhyxyxyx 由由Taylor展开，有展开，有(1)2 Runge-Kutta Method(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1)( ,)jjjjyf x yfff dyffyffxy dxxyffyffxyffyffxy (2)2 Runge-Kutta Method要使公式具有要使公式具有p阶精度，则在阶精度，则在(1)式中截取前式中截取前p+

27、1项，项，用用(2)式计算各阶导数，即得下面式计算各阶导数，即得下面Taylor公式：公式：2( )12!ppnnnnnhhyyhyyyp 局部截断误差局部截断误差(3) 1(1)111, 1 !ppnnnnhy xyyxxp 2 Runge-Kutta Method2.Taylor公式公式(3)表面上看形式简单，但具体构造时表面上看形式简单，但具体构造时往往很困难，因为按往往很困难，因为按(2)式求导，这一过程可能很复式求导，这一过程可能很复杂。因此通常不直接用杂。因此通常不直接用Taylor公式，而借鉴其思想公式，而借鉴其思想提出其它公式。提出其它公式。1. 由此看出，一种方法具有由此看出

28、，一种方法具有p阶精度阶精度公式对不公式对不超过超过p次的多项式准确成立局部截断误差为次的多项式准确成立局部截断误差为0）。）。这一等价条件也可以用来判断一种方法的精度。这一等价条件也可以用来判断一种方法的精度。2 Runge-Kutta Method单步递推法的基本思想是从单步递推法的基本思想是从 ( xn , yn ) 点出发，以某一点出发，以某一斜率沿直线达到斜率沿直线达到 ( xn+1 , yn+1 ) 点。欧拉法及其各种变点。欧拉法及其各种变形所能达到的最高精度为形所能达到的最高精度为2阶。阶。 2. RungeKutta Method由微分中值定理，有由微分中值定理，有*1*1()

29、( )(), 01 ()( )nnnnnnny xy xy xhf xh y xhkhy xy xhk k*称为区间称为区间xn, xn+1上的平均斜率，只要知道平均上的平均斜率，只要知道平均斜率，就可计算斜率，就可计算y(xn+1).因此只要对平均斜率提供一因此只要对平均斜率提供一种近似算法，则由种近似算法，则由(4)式可导出一种相应的求解公式。式可导出一种相应的求解公式。(4) 2 Runge-Kutta Method *1*112*121. 0, ,Eulers formula2. 1, ,Eulers implicit formula3. Trapezoid formula2nnnnk

30、f x y xkkf xy xkkkk例例121*12(,)4. (,) 2 modified Eulers method nnnnkf xykf xh yhkkkk2 Runge-Kutta Method 由此看出，改进的尤拉公式用由此看出，改进的尤拉公式用xn与与xn+1两个节两个节点的斜率的算术平均作为平均斜率，点的斜率的算术平均作为平均斜率， xn+1点的斜点的斜率通过已知信息率通过已知信息yn来预测。来预测。 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1， K

31、2 的平均值吗？的平均值吗？步长一定是步长一定是h 吗？即吗？即第二个节点一定是第二个节点一定是xn+1吗？吗？2 Runge-Kutta Method2 Runge-Kutta Method首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶精度，即在阶精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()

32、()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii l ll l),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii l ll ll ll ll l 2阶阶RungeKutta Method2 Runge-Kutta MethodStep 3: 将将 yi+1 与与 y( xi+

33、1 ) 在在 xi 点的泰勒展开作比较点的泰勒展开作比较)()()()(322211hOxyphxyhyyiiii l ll ll l)()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 pl ll ll l这里有这里有 个未知个未知数，数， 个方程。个方程。32存在无穷多个解。所有满足上式的格式统称为存在无穷多个解。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式。塔格式。21, 121 l ll lp注意到，注意到， 就是改进的欧拉法；就是改进的欧拉法；p=1/2, 1=0, 2=1, 变形尤拉

34、公式。变形尤拉公式。Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？改进的改进的Euler 公式推广为二阶公式推广为二阶Runge-Kutta公式公式带来这样的启示：带来这样的启示：若在若在xn, xn+1上多预测几个点的斜率值，然后上多预测几个点的斜率值，然后将它们的算术平均作为平均斜率，则有可能构造将它们的算术平均作为平均斜率，则有可能构造出具有更高精度的计算公式。出具有更高精度的计算公式。-Runge-Kutta方法的基本思想。方法的基本思想。注：二阶注：二阶Runge-Kutta公式用多算一次函数值公式用多算一次函数值f 的的办法避开了二阶办法避开了二

35、阶Taylor级数法所要计算的级数法所要计算的f 的导数。的导数。在这种意义上，可以说在这种意义上，可以说Runge-Kutta方法实质上是方法实质上是Taylor级数法的变形。级数法的变形。2 Runge-Kutta Method其中其中i ( i = 1, , m )，i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定系均为待定系数，确定这些系数的步骤与前面相似。数，确定这些系数的步骤与前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmm

36、mimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyyb bb bb ba ab bb ba ab ba al ll ll l 高阶高阶RungeKutta Method Gill公式：公式：4阶经典龙格阶经典龙格-库塔公式的一种改进库塔公式的一种改进112346121222 1231222222423222222( ,)(,)(,1)(,1)hiiiihhiihiiiiyyKKKKKf xyKf xyKKf xyhKhKKf xh yhKhK2 Runge-Kutta Method 最常用为四级最常用为四级4阶经典龙格阶经典龙格-库塔法库塔法 /

37、* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 2 Runge-Kutta Method注：注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4

38、hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算于光滑性不太好的解，最好采用低阶算法而将步长法而将步长h 取小。取小。2 Runge-Kutta Method 变步长的变步长的RungeKutta Method1(1)1ppnnRchyxQ: 由局部截断误差可以看出，步长由局部截断误差可以看出，步长 h 越小，局部截断越小，局部截断误差越小；但步长减小，在一定求解范围区间内误差越小；但步长减小，在一定求解范围区间内要完成

39、的步数就增加了，步数增加会引起计算量增大，要完成的步数就增加了，步数增加会引起计算量增大，导致舍入误差积累。因此要选取适当的步长。导致舍入误差积累。因此要选取适当的步长。选择步长时要考虑两个问题：选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？如何根据所获得的精度处理步长？HW: p.201 #6-83 单步法的收敛性与稳定性单步法的收敛性与稳定性 /* Convergency and Stability */ 前面介绍了两大类微分方程数值解法：一类是前面介绍了两大类微分方程数值解法：一类是用差商近似导数得到的尤拉

40、系列公式，另一类是用差商近似导数得到的尤拉系列公式，另一类是基于平均斜率概念的基于平均斜率概念的RungeKutta公式。基本思公式。基本思想都是通过某种离散化手续，将微分方程转化为想都是通过某种离散化手续，将微分方程转化为差分方程代数方程来求解。差分方程代数方程来求解。 Q1. 这种转化是否合理？要看差分问题的解这种转化是否合理？要看差分问题的解yn当当h0时是否收敛到微分方程的解时是否收敛到微分方程的解y(xn),即是否成立即是否成立 yn y(xn)， h0. -收敛性问题收敛性问题 Q2. 实际计算时，由于舍入误差的影响，差分方实际计算时，由于舍入误差的影响，差分方程的解本身也有误差，

41、这种误差在计算过程中会程的解本身也有误差，这种误差在计算过程中会不会扩大不会扩大 -稳定性问题稳定性问题3 Convergency and Stability 收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是收敛的。，则称该算法是收敛的。 例：就初值问题例：就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0)0(yyyyl l解：该问题的精确解为解：该问题的精确解为 xeyxyl l0)( 欧拉公式为欧拉公

42、式为iiiiyhyhyy)1 (1l ll l 0)1 (yhyiil l 对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyyl ll ll ll l)1()1(/10/0 ehhhl ll l/10)1(lim)(0ixxyeyi l l 3 Convergency and Stability显式单步法的收敛性显式单步法的收敛性(1)3 Convergency and Stability而整体截断误差为而整体截断误差为1111111nnnnnnney xyy xyyy证明证明: 设设 表示当表示当yn =y(xn)时，时， 由公式由公式(1)求得求得的结果，

43、即的结果，即1,nnnnyy xhxy xh1ny则局部截断误差为则局部截断误差为111pnny xych(2)(2)-(1),得得 11,11nnnnnnnnnnnnnnnyyy xyhx y xhx y hy xyhL y xyhLy xyhLe3 Convergency and Stability从而从而 11111211112010100111111111111111pppnnnppnnnnpnnnpnnpehLechhLhLechchhLehLchchehLehLhLhLchhLhLechhLhLhLechL3 Convergency and Stability000 ()1 ,1,

44、1011 0 1 ()nnnhLTLnxnxTLTLpnTLpnpnnxxnhTbahLeexRxexxeeee echeLceheLy xyo h 当时，有当时， 而 ，3 Convergency and Stability Euler法的收敛性法的收敛性 ： (x,y)=f(x,y),故当故当f(x,y)满足满足Lipschitz条件时，条件时，尤拉法收敛；尤拉法收敛；3 Convergency and Stability3 Convergency and Stability 稳定性稳定性 /* Stability */例：考察初值问题例：考察初值问题 在区间在区间0, 0.5上的解。上的

45、解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉显式欧拉显式 节点 xixey30 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.2341

46、1046.14421063.0590107What 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定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的算中都逐

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

48、 稳定，就是指稳定，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的大。的大。000yy hl l h3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability3 Convergency and Stability例：隐式龙格例：隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmiib bb ba al ll l而显式而显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii其

49、中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。小结小结 尤拉两步公式与尤拉单步公式相比，使用两尤拉两步公式与尤拉单步公式相比，使用两个节点上的已知信息将精度提高一阶。可以设想，个节点上的已知信息将精度提高一阶。可以设想，计算计算y(xn+1)时，充分利用前面已经求出的节点时，充分利用前面已经求出的节点上的上的 y 及及 y 值的线性组合来近似值的线性组合来近似y(xn+1)，精度，精度会大大

50、提高。会大大提高。-线性多步法的基本思想。线性多步法的基本思想。 构造线性多步法有多种途径，这里介绍两种：构造线性多步法有多种途径，这里介绍两种： 基于数值积分的构造方法；基于数值积分的构造方法； 基于基于Taylor展开的构造方法。展开的构造方法。4 线性多步法线性多步法 /* Multistep Method */).(.110111101rnrnnnrnrnnnffffhyyyy b bb bb bb ba aa aa a线性多步法的通式可写为：线性多步法的通式可写为：),(jjjyxff 当当 10 时，为隐时，为隐式公式式公式; 1=0 则为则为显式公式。显式公式。 基于数值积分的构

51、造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxfy 1,nnx x11()()( , ( )nnxnnxy xy xf x y x dx只要近似地算出右边的积分只要近似地算出右边的积分 ，则可通，则可通过过 近似近似y(xn+1) 。而选用不同近似式。而选用不同近似式 I，可得到，可得到不同的计算公式。例如利用左矩形积分公式得到尤拉公式；不同的计算公式。例如利用左矩形积分公式得到尤拉公式；梯形积分公式得到梯形公式。梯形积分公式得到梯形公式。1( , ( )nnxxIf x y x dx1nnyyI一般地，利用插值原理所建立的一系列数值积分一般地，利用插值原理所建立的一系

52、列数值积分方法也可以导出解微分方程的一系列计算公式。方法也可以导出解微分方程的一系列计算公式。运用插值方法的关键在于选取合适的插值节点。运用插值方法的关键在于选取合适的插值节点。假设已构造出假设已构造出f(x,y(x)的插值多项式的插值多项式Pr(x),那么那么1111( , ( )( )( )nnnnnnxxrxxxnnrxf x y x dxP x dxyyP x dx4 Multistep Method 亚当姆斯显式公式亚当姆斯显式公式 /* Adams explicit formulae */利用利用r+1 个节点上的被积函数值个节点上的被积函数值构造构造 r 阶牛顿后插多项式阶牛顿后

53、插多项式, 有有 11,nnnnn rn rxfxfxf11100( , ( )()()nnnx xthxrnrnxf x y x dxN xth hdtR xth hdtNewton插值余项插值余项1110001000()11rjjnnrnnn jjrrjjjnn jnjn jjjtyyhN xth dtyhfdtjtyhdtfyhfja0()11, rjjrnnjjjntNxthfjxxtjh ，其中0表示 阶向前差分，/*Adams 显式公式显式公式 */外推技术外推技术 /* extrapolation */10 is independent of n and r, 1jjjtdtja

54、awhere table 5-6, p.181ja aj0111/225/1233/8实际计算时，将差分展开实际计算时，将差分展开 000000111jijnjn iijrrrrriijjnjjn ijn irin ijjiij iijffijjffffiiaaab ijrj=i10rnnrin iiyyhfb局部截断误差为：局部截断误差为： 1110(2)12(2)102(2)()()()1 !()rnnrnrnrrrrnrrrnRy xyhR xth dtyht dtB hyrB hyx注：注：Br 与与yn+1 计算公式中计算公式中 fn , , fnk 各项的系数各项的系数 均可均可查

55、表得到查表得到 。 见下表。见下表。rib10123r2123122324552112162459125243724912583720251fnfn1fn2fn3Br例：例：1110,1,(3)2nnnnnnnryyhfhryyff35()12nRh yx常用的是常用的是 r = 3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式)9375955(243211 iiiiiiffffhyy5(5)251()720nRh yx21()2nRh y xAdams显式公式用显式公式用 作为插作为插值节点，在求积区间值节点，在求积区间xn, xn+1上用插值函数上用插值函数Nr(x)近似近似f(x,y(x)

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

57、隐式公式 /* Adams implicit formulae */利用利用r+1 个节点上的被积函数值个节点上的被积函数值 fn+1 , fn , , fnr+1 构造构造r 阶牛顿前插多项式。与显式多项式完全类似地可得到一系阶牛顿前插多项式。与显式多项式完全类似地可得到一系列隐式公式。列隐式公式。0*1110*110, 1, 1rjjnnjnjjjrrinnrinirijij ityyhfdtjjyyhfiaabba 差分展开，得j0 11 -1/22 -1/123 -1/24*ja局部截断误差为：局部截断误差为：2(2)112(2)()()()rrrnnrnrrrnRy xyB hyB

58、hyx10123k21 21125249211282419121 245 241121 241 72019 fi+1fifi1fi2Br小于小于Br注：注： 与与yn+1 计算公式中计算公式中 fn+1 , , fn+1k 各项的系数各项的系数 均可查表得到均可查表得到 。 见下表。见下表。*ribrB常用的是常用的是 k = 3 的的4阶亚当姆斯隐式公式阶亚当姆斯隐式公式)5199(242111 iiiiiiffffhyy较同阶显式较同阶显式稳定稳定例：例：11110,1,()2nnnnnnnryyhfhryyff31()12nRh yx 21()2nRh y x 5(5)19()720nR

59、h yx 4 Multistep Method 亚当姆斯预测亚当姆斯预测-校正系统校正系统 /* Adams predictor-corrector system */Step 1: 用用Runge-Kutta 法计算前法计算前 k 个初值；个初值；Step 2: 用用Adams 显式计算预测值；显式计算预测值；Step 3: 用同阶用同阶Adams 隐式计算校正值。隐式计算校正值。注意：三步所用公注意：三步所用公式的精度必须相同。式的精度必须相同。通常用经典通常用经典Runge-Kutta 法配合法配合4阶阶Adams 公式。公式。4阶阶Adams隐式公式的截断误差为隐式公式的截断误差为5(

60、5)1119()()720nnny xyh yx 4阶阶Adams显式公式的截断误差为显式公式的截断误差为5(5)11251()()720nnny xyh yx1111()19()251nnnny xyy xy Predicted value pn+1Corrected value cn+11111251()()270nnnny xyyy 111119()()270nnnny xyyyModified value mn+1预测值与校正预测值与校正值的事后误差值的事后误差估计估计1111251()()270nnnny xyyy111119()()270nnnny xyyy112311111211