第5章_常微分方程数值解法

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 xy1

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

5、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()()() ()()(,)nnnnnnnn

8、nny 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)11

10、111()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、 2nnnnynnnnnnnnnnnnynnnnnnnnfxyfxy xfxyy xyy xhfxy 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 ( )2nynynnnnynnnnnh

13、hRhfxhfxy 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 Method梯

16、形法的简化计算梯形法的简化计算 迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法预测校正系统，也称为改进尤拉法： 改进欧拉法改进欧拉法 /* 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 Method+1

18、+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个初值个初值 y

20、0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。1 Eulers Method几何解释几何解释xnxn+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 15nnnnnnnnmpp

23、cycpc改进：此时此时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,2nnnnnnn

25、nnnyyhyyf 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+1项，项，用用(2

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

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

29、 ()( )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 formula2nnnnkf x y xkk

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

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

32、phxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 2阶阶RungeKutta Method2 Runge-Kutta MethodStep 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)(

33、)()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法；就是改进的欧拉法；p=1/2, 1=0, 2=1, 变形尤拉公式变形尤拉公式。Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？改进的改

34、进的Euler 公式推广为二阶公式推广为二阶Runge-Kutta公式公式带来这样的启示：带来这样的启示：若在若在xn, xn+1上多预测几个点的斜率值，然后将上多预测几个点的斜率值，然后将它们的算术平均作为平均斜率，则有可能构造出它们的算术平均作为平均斜率，则有可能构造出具有更高精度的计算公式。具有更高精度的计算公式。-Runge-Kutta方法的基本思想。方法的基本思想。注：二阶注：二阶Runge-Kutta公式用多算一次函数值公式用多算一次函数值f 的的办法避开了二阶办法避开了二阶Taylor级数法所要计算的级数法所要计算的f 的导数。的导数。在这种意义上，可以说在这种意义上，可以说Ru

35、nge-Kutta方法实质上是方法实质上是Taylor级数法的变形。级数法的变形。2 Runge-Kutta Method其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定系数，确均为待定系数，确定这些系数的步骤与前面相似。定这些系数的步骤与前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 高

36、阶高阶RungeKutta Method Gill公式：公式：4阶经典龙格阶经典龙格-库塔公式的一种改进库塔公式的一种改进112346121222 1231222222423222222( ,)(,)(,1)(,1)hiiiihhiihiiiiyyKKKKKf xyKf xyKKf xyhKhKKf xh yhKhK2 Runge-Kutta Method 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxf

37、KKKKKyyiihihihihiiihii 2 Runge-Kutta Method注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性

38、不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。2 Runge-Kutta Method 变步长的变步长的RungeKutta Method1(1)1ppnnRchyxQ: 由局部截断误差可以看出，步长由局部截断误差可以看出，步长 h 越小，局部截断越小，局部截断误差越小；但步长减小，在一定求解范围（区间）内误差越小；但步长减小，在一定求解范围（区间）内要完成的步数就增加了，步数增加会引起计算量增大，要完成的步数就增加了，步数增加会引起计算量增大，导致舍入误差积累。因此要选取适当的步长。导致舍入误差积累。因此要选取适当的步长。

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

40、分方程转化为想都是通过某种离散化手续，将微分方程转化为差分方程（代数方程）来求解差分方程（代数方程）来求解。 Q1. 这种转化是否合理？要看差分问题的解这种转化是否合理？要看差分问题的解yn当当h0时是否收敛到微分方程的解时是否收敛到微分方程的解y(xn),即是否成立即是否成立 yn y(xn)， h0. -收敛性问题收敛性问题 Q2. 实际计算时，由于舍入误差的影响，差分方实际计算时，由于舍入误差的影响，差分方程的解本身也有误差，这种误差在计算过程中会程的解本身也有误差，这种误差在计算过程中会不会扩大不会扩大 -稳定性问题稳定性问题3 Convergency and Stability 收敛

41、性收敛性 /* Convergency */定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 0)1 (yhyii 对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 eh

42、hh /10)1(lim)(0ixxyeyi 3 Convergency and Stability显式单步法的收敛性显式单步法的收敛性(1)3 Convergency and Stability而整体截断误差为而整体截断误差为1111111nnnnnnney xyy xyyy证明证明: 设设 表示当表示当yn =y(xn)时，时， 由公式由公式(1)求得的求得的结果，即结果，即1,nnnnyy xhxy xh1ny则局部截断误差为则局部截断误差为111pnny xych(2)(2)-(1),得得 11,11nnnnnnnnnnnnnnnyyy xyhx y xhx y hy xyhL y x

43、yhLy xyhLe3 Convergency and Stability从而从而 11111211112010100111111111111111pppnnnppnnnnpnnnpnnpehLechhLhLechchhLehLchchehLehLhLhLchhLhLechhLhLhLechL3 Convergency and Stability000 ()1 ,1,1011 0 1 ()nnnhLTLnxnxTLTLpnTLpnpnnxxnhTbahLeexRxexxeeee echeLceheLy xyo h 当时，有当时， 而 ，3 Convergency and Stability E

44、uler法的收敛性法的收敛性 ： (x,y)=f(x,y),故当故当f(x,y)满足满足Lipschitz条件时，尤拉条件时，尤拉法收敛；法收敛；3 Convergency and Stability3 Convergency and Stability 稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式

45、欧拉显式 节点节点 xixey30 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 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so

46、unstable that a small perturbation on the conditions will cause a crash on the conclusions!3 Convergency and Stability3 Convergency and Stability定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。一般分析某算法的稳定性时，为简单起见，只考虑一般分析某算法的稳定性时，为简单起见，只考

47、虑模型方程模型方程或或试验方程试验方程 /* test equation */yy 常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。000yy h h h3 Convergency and Stability3 Converg

48、ency and Stability3 Convergency and Stability3 Convergency and Stability例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。一般来

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

50、ylorTaylor展开的构造方法展开的构造方法。4 线性多步法线性多步法 /* Multistep Method */).(.110111101rnrnnnrnrnnnffffhyyyy 线性多步法的通式可写为：线性多步法的通式可写为：),(jjjyxff 当当 1 0 时，为时，为隐式公隐式公式式; 1=0 则为则为显式公式显式公式。 基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxfy 1,nnx x11()()( , ( )nnxnnxy xy xf x y x dx只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通，则可通过过 近似近似y

51、(xn+1) 。而。而选用不同近似式选用不同近似式 I，可得到不，可得到不同的计算公式同的计算公式。例如利用左矩形积分公式得到尤拉公式；。例如利用左矩形积分公式得到尤拉公式；梯形积分公式得到梯形公式。梯形积分公式得到梯形公式。1( , ( )nnxxIf x y x dx1nnyyI一般地，利用插值原理所建立的一系列数值积分一般地，利用插值原理所建立的一系列数值积分方法也可以导出解微分方程的一系列计算公式。方法也可以导出解微分方程的一系列计算公式。运用插值方法的关键在于选取合适的插值节点。运用插值方法的关键在于选取合适的插值节点。假设已构造出假设已构造出f(x,y(x)的插值多项式的插值多项式

52、Pr(x),则则1111( , ( )( )( )nnnnnnxxrxxxnnrxf x y x dxP x dxyyP x dx4 Multistep Method 亚当姆斯显式公式亚当姆斯显式公式 /* Adams explicit formulae */利用利用r+1 个节点上的被积函数值个节点上的被积函数值构造构造 r 阶牛顿阶牛顿后插后插多项式多项式, 有有 11,nnnnn rn rxfxfxf11100( , ( )()()nnnx xthxrnrnxf x y x dxN xth hdtR xth hdtNewton插值余项插值余项1110001000()11rjjnnrnnn

53、 jjrrjjjnn jnjn jjjtyyhN xth dtyhfdtjtyhdtfyhfj0()11, rjjrnnjjjntNxthfjxxtjh ，其中0表示 阶向前差分，/*Adams 显式公式显式公式 */外推技术外推技术 /* extrapolation */10 is independent of n and r, 1jjjtdtjwhere table 5-6, p.181j j0111/225/1233/8实际计算时，将差分展开实际计算时，将差分展开 000000111jijnjn iijrrrrriijjnjjn ijn irin ijjiij iijffijjffffi

54、i ijrj=i10rnnrin iiyyhf局部截断误差为：局部截断误差为： 1110(2)12(2)102(2)()()()1 !()rnnrnrnrrrrnrrrnRy xyhR xth dtyht dtB hyrB hyx注：注：Br 与与yn+1 计算公式中计算公式中 fn , , fn k 各项的各项的系数系数 均可查表均可查表得到得到 。 见下表。见下表。ri10123r2123122324552112162459125243724912583720251fnfn 1fn 2fn 3Br例：例：1110,1,(3)2nnnnnnnryyhfhryyff35()12nRh yx常用

55、的是常用的是 r = 3 的的4阶亚当姆斯显式公式阶亚当姆斯显式公式)9375955(243211 iiiiiiffffhyy5(5)251()720nRh yx21()2nRh y xAdams显式公式用显式公式用 作为插作为插值节点，在求积区间值节点，在求积区间xn, xn+1上用插值函数上用插值函数Nr(x)近近似似f(x,y(x)，而，而xn+1不在插值节点内，因此是一个不在插值节点内，因此是一个外推外推的过程。虽然公式是显式的，便于计算，但的过程。虽然公式是显式的，便于计算，但效果并不理想，比如稳定性较差等。效果并不理想，比如稳定性较差等。因此改用通过因此改用通过r+1个节点个节点的

56、插值多项式的插值多项式Nr(x)近似近似f(x,y(x)，由于，由于xn+1是其中是其中一个插值点，因此是一个插值点，因此是内插多项式内插多项式，但导出的公式，但导出的公式是是隐式的隐式的。 11,nnnnn rn rxfxfxf 1111,nnnnn rn rxfxfxf 4 Multistep Method 亚当姆斯隐式公式亚当姆斯隐式公式 /* Adams implicit formulae */利用利用r+1 个节点上的被积函数值个节点上的被积函数值 fn+1 , fn , , fn r+1 构造构造r 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项式

57、完全类似地可得到一系列隐式公式。隐式公式。0*1110*110, 1, 1rjjnnjnjjjrrinnrinirijij ityyhfdtjjyyhfi 差分展开，得j0 11 -1/22 -1/123 -1/24*j局部截断误差为：局部截断误差为：2(2)112(2)()()()rrrnnrnrrrnRy xyB hyB hyx10123k21 21125249211282419121 245 241121 241 72019 fi+1fifi 1fi 2Br小于小于Br注：注： 与与yn+1 计算公式中计算公式中 fn+1 , , fn+1 k 各项的各项的系数系数 均均可查表得到可查表

58、得到 。 见下表。见下表。*rirB常用的是常用的是 k = 3 的的4阶亚当姆斯隐式公式阶亚当姆斯隐式公式)5199(242111 iiiiiiffffhyy较同阶显式较同阶显式稳定稳定例：例：11110,1,()2nnnnnnnryyhfhryyff31()12nRh yx 21()2nRh y x 5(5)19()720nRh yx 4 Multistep Method 亚当姆斯预测亚当姆斯预测-校正系统校正系统 /* Adams predictor-corrector system */Step 1: 用用Runge-Kutta 法法计算前计算前 k 个初值；个初值；Step 2: 用

59、用Adams 显式显式计算计算预测预测值；值；Step 3: 用同阶用同阶Adams 隐式隐式计算计算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必须相同。式的精度必须相同。通常用通常用经典经典Runge-Kutta 法法配合配合4阶阶Adams 公式公式。4阶阶Adams隐式公式的截断误差为隐式公式的截断误差为5(5)1119()()720nnny xyh yx 4阶阶Adams显式公式的截断误差为显式公式的截断误差为5(5)11251()()720nnny xyh yx1111()19()251nnnny xyy xy Predicted value pn+1Corrected

60、 value cn+11111251()()270nnnny xyyy 111119()()270nnnny xyyyModified value mn+1预测值与校正预测值与校正值的事后误差值的事后误差估计估计1111251()()270nnnny xyyy111119()()270nnnny xyyy11231111121111(5559379)24251()270(9195)2419()270nnnnnnnnnnnnnnnnnnnnhpyffffmppchcyffffycpc预测：改进：校正：改进：Adams显隐预测显隐预测-校正校正-改进系统改进系统4 Multistep Method