第四章常微分方程数值解法

1、四常微分方程数值解法2常微分方程数值解法引言（常微分方程数值解法概述）引言（常微分方程数值解法概述）显式欧拉法、隐式欧拉法、二步欧拉法显式欧拉法、隐式欧拉法、二步欧拉法局部截断误差与精度局部截断误差与精度改进的欧拉方法改进的欧拉方法龙格龙格-库塔方法库塔方法收敛性与稳定性简述收敛性与稳定性简述一阶常微分方程组与高阶常微分方程一阶常微分方程组与高阶常微分方程3引言引言一阶常微分方程初值问题：一阶常微分方程初值问题：00( , )()yf x yy xy 微分方程微分方程初始条件初始条件定理：定理：若若 f (x, y) 在某闭区域在某闭区域 r ： 00|,|(0,0)xxayybab上连续，且

2、在上连续，且在 r 域内满足李普希兹域内满足李普希兹 (lipschitz) 条件，条件，即存在正数即存在正数 l，使得对于，使得对于 r 域内的任意两值域内的任意两值 y1, y2，下，下列不等式成立：列不等式成立： 1221|( ,)( ,)|f x yf x yl yy则上述初值问题的连续可微的解则上述初值问题的连续可微的解 y(x) 存在并且唯一。存在并且唯一。4引言（续）引言（续）实际生产与科研中，除少数简单情况能获得初值问题实际生产与科研中，除少数简单情况能获得初值问题的初等解（用初等函数表示的解）外，绝大多数情况的初等解（用初等函数表示的解）外，绝大多数情况下是求不出初等解的。下

3、是求不出初等解的。有些初值问题即便有初等解，也往往由于形式过于复有些初值问题即便有初等解，也往往由于形式过于复杂而不便处理。杂而不便处理。实用的方法是在计算机上进行数值求解：即不直接求实用的方法是在计算机上进行数值求解：即不直接求 y(x) 的显式解，而是在解所存在的区间上，求得一系的显式解，而是在解所存在的区间上，求得一系列点列点 xn (n 0, 1, 2, ) 上解的近似值。上解的近似值。5欧拉欧拉（euler）方法方法方法一方法一化导数为差商的方法化导数为差商的方法10()()()()()limnnnnnhy xhy xy xy xy xhh 由于在逐步求解的过程中，由于在逐步求解的过

4、程中，y(xn) 的准确值无法求解的准确值无法求解出来，因此用其近似值代替。出来，因此用其近似值代替。为避免混淆，以下学习简记：为避免混淆，以下学习简记：y(xn)：待求函数：待求函数 y(x) 在在 xn 处的精确函数值处的精确函数值yn ：待求函数：待求函数 y(x) 在在 xn 处的近似函数值处的近似函数值6代入初值问题表达式可得：代入初值问题表达式可得：根据根据 y0 可以一步步计算出函数可以一步步计算出函数 y y(x) 在在 x1, x2, x3 x4, 上的近似值上的近似值 y1, y2, y3, y4 , 常微分方程数值解是一组离散的函数值数据，它的常微分方程数值解是一组离散的

5、函数值数据，它的精确表达式很难求解得到，但可以进行插值计算后精确表达式很难求解得到，但可以进行插值计算后用插值函数逼近用插值函数逼近 y(x)10()()()()()limnnnnnhy xhy xy xy xy xhh 100(,)0,1,2,()nnnnyyhf xynyy x 1(,)nnnnyyf xyh 7欧拉方法（续）欧拉方法（续）方法二方法二数值积分法数值积分法1111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 同样以近似值同样以近似值 yn 代替精确值代替精确值 y(xn) 可得可

6、得：100(,)()nnnnyyhf xyyy x 将微分方程将微分方程 y f (x, y) 在区间在区间 xn, xn+1 上积分：上积分：8欧拉方法的几何意义欧拉方法的几何意义xy00 x1x2x3x4x5x6x0p1p2p3p4p5p6p( )yy x9隐式欧拉法隐式欧拉法在数值积分法推导中，积分的近似值取为积分区间宽在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：度与右端点处的函数值乘积，即：这样便得到了隐式欧拉法：这样便得到了隐式欧拉法：11100(,)()nnnnyyh f xyyy x 含有未知含有未知的函数值的函数值隐式欧拉法没有显式欧拉法方便隐式

7、欧拉法没有显式欧拉法方便11111111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 10二步欧拉法二步欧拉法在数值积分法推导中，积分区间宽度选为两步步长，在数值积分法推导中，积分区间宽度选为两步步长，即积分区间为：即积分区间为：xn 1, xn 1，则：，则：11111111d()() , ( )d() , ()2, ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xh f xy x 以以 y(x) 在在 xn 1, xn 上的近似值上的近似值代替精确值代替精确值

8、可得可得：11002(,)()nnnnyyhf xyyy x 需要前两步需要前两步的计算结果的计算结果中矩形公式中矩形公式11梯形公式欧拉法梯形公式欧拉法在数值积分法中，如果用梯形公式近似计算在数值积分法中，如果用梯形公式近似计算 f (x, y) 在区间在区间 xn, xn+1 上的积分，即：上的积分，即：用近似值代替精确值可得梯形公式欧拉法：用近似值代替精确值可得梯形公式欧拉法：111 (,)(,)2nnnnnnhyyf xyf xy上式右端出现了未知项，可见梯形法是隐式欧拉法上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与的一种；实际上，梯形公式欧

9、拉法是显式欧拉法与隐式欧拉法的隐式欧拉法的算术平均算术平均。 111111, (), () , ( )d()2, (), ()2nnxnnnnnnxnnnnf xy xf xy xf x y xxxxhf xy xf xy x 12例例用显式欧拉法、隐式欧拉法、梯形法求解初值问题：用显式欧拉法、隐式欧拉法、梯形法求解初值问题：1(0)1yyxy 取取 h 0.1，计算到，计算到 x 0.5，并与精确解进行比较，并与精确解进行比较解：由已知条件可得：解：由已知条件可得：h 0.1，x0 0， y0 1， f (x, y) y x 1显式欧拉法：显式欧拉法：1(1)(1)0.90.10.1nnnn

10、nnnnyyhyxh yhxhyx 13例：（续）例：（续）隐式欧拉法：隐式欧拉法：111(1)nnnnyyhyx化简得：化简得：11(1)(1)(1)nnnnnh yyh xyh xh11(1)(0.10.11) 1.1(1)nnnnnyyh xhyxh 111(1)(1)2nnnnnnhyyyxyx 梯形公式欧拉法：梯形公式欧拉法：111(2)1(2)(2)(2)nnnnhyhyhxx11(2)(2) (2)(1.90.20.21) 2.1nnnnnnyh yh xxhyx14计算结果：计算结果：xn显式法显式法 yn隐式法隐式法 yn梯形法梯形法 yn精确解精确解 y (xn)0.011

11、110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.106531( )exy xx 本题的精确解为：本题的精确解为：10.90.10.1nnnyyx 1(0.10.11) 1.1nnnyyx 1(1.90.20.21) 2.1nnnyyx 显式法显式法隐式法隐式法梯形法梯形法15局部截断误差局部截断误差为了简

12、化分析某常微分方程数值算法的误差，现假为了简化分析某常微分方程数值算法的误差，现假设设 yn y(xn)，即，即在前一步在前一步 yn 准确的前提下准确的前提下，估计：，估计：111()nnnty xy称上述误差称上述误差 tn 1 为该常微分方程数值算法的为该常微分方程数值算法的局部截局部截断误差断误差如果某个常微分方程数值算法的局部截断误差可表如果某个常微分方程数值算法的局部截断误差可表示为示为 o(h p 1)，则称该数值算法的精度是，则称该数值算法的精度是 p 阶阶欧拉法的精度为一阶；二步欧拉法的精度为二阶；欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形公式欧拉法的精度为二阶。梯形公

13、式欧拉法的精度为二阶。16泰勒展开法泰勒展开法如果初值问题中的如果初值问题中的 f (x, y) 充分可微，则可将充分可微，则可将 y(xn 1) 在点在点 xn 处展开：处展开：21()()()()()2!nnnnnhy xy xhy xhy xyx 如果只保留线性项，忽略如果只保留线性项，忽略 h2 及以上各项，则：及以上各项，则：11()()()nnnny xy xhy xy (), ()nnny xf xy x ny1(,)nnnnyyhf xy 显式欧拉公式显式欧拉公式17局部截断误差的分析局部截断误差的分析利用泰勒公式展开，比较各算法与展开式的前几项利用泰勒公式展开，比较各算法与展

14、开式的前几项将将 y(xn 1) 在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：211( )()()()()2!(,)nnnnnnyy xy xhy xhy xhxx 22111( )()()2!nnnyty xyho h 显式欧拉法的局部截断误差：显式欧拉法的局部截断误差：1(), ()()()nnnnnnyy xhf xy xy xhy x 1(,)nnnnyyhf xy 欧拉法欧拉法()nnyy x 令令：1 阶精度阶精度18补充：二元函数微分中值定理补充：二元函数微分中值定理1100(,)(,)f xyf xy 1001010010()()()()xyxxfxxxyyfyyy 0

15、1 0000(,)(,)f xh ykf xy00()()xyhfxhkfyk000000(,)(,)()()xyf xh ykf xyhfxhkfyk0000(,)(,)f xyhkf xh ykxy191111111()(,),(, )()nnnynnnnyf xyf xxxfxyy 111()()( )(, )nnnynnyy xh y xo hfxt 21()()()()nnny xy xhy xo h y(xn 1) 在在 xn 点处展开：点处展开：211111()(, )()nnnynnty xyhfxto h 111(,)nnnnyyhf xy隐式欧拉法：隐式欧拉法：1 阶精度阶

16、精度()ny x2111()1(, )nynto hhfx 11 (),nny xy 111()(, )nynny xfxt 11()( )(, )nynny xo hfxt 211()()(, )()nnynny xhy xhfxto h 20111()2, ()()2()nnnnnnyy xhf xy xy xhy x 分别将分别将 y(xn 1), y(xn 1) 在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：231231()()()()()2!()()()()()()2!nnnnnnnnyxy xy xhy xho hyxy xy xhy xho h 111()nnnty xy二

17、步欧拉法的局部截断误差：二步欧拉法的局部截断误差：112(,)nnnnyyhf xy二步欧拉法：二步欧拉法：2 阶精度阶精度1()ny x ()ny x311()()2()()nnny xy xhy xo h 211111111111(,),(, )()()()(, )nnnynnnnynnnf xyf xfxyy xy xfty xx 111 (,)(,)2nnnnnnhyyf xyf xy梯形公式欧拉法：梯形公式欧拉法：()ny x2111()2()()()(, )2nnnnynnhyy xy xhyxo hfxt y(xn 1) 在在 xn 点处展开：点处展开：231()()()()()

18、2!nnnnhy xy xhy xyxo h 3111(, )()2nynnhtfxto h 31211()1(, )nhynto hfx 2 阶精度阶精度11 (),nny xy 2311()()()(, )()22nnnynnhhy xhy xyxfxto h 211()()()(, )nnynny xhyxo hfxt 22各种欧拉法的比较各种欧拉法的比较方法方法精度精度评述评述显式欧拉法显式欧拉法1最简单，精度低最简单，精度低隐式欧拉法隐式欧拉法1不便计算，稳定性好不便计算，稳定性好二步欧拉法二步欧拉法2需要两步初值，且第需要两步初值，且第 2 个初值只能由其个初值只能由其它方法给出，

19、可能对后面的递推精度有它方法给出，可能对后面的递推精度有影响影响梯形公式梯形公式欧欧 拉拉 法法2精度有所提高，但为隐式，需要迭代求精度有所提高，但为隐式，需要迭代求解，计算量大解，计算量大23改进的欧拉法改进的欧拉法从上述例子可以看到，梯形法由于具有二阶精度，其从上述例子可以看到，梯形法由于具有二阶精度，其局部截断误差比显式欧拉法和隐式欧拉法小，但梯形局部截断误差比显式欧拉法和隐式欧拉法小，但梯形法实质上是一种隐式算法法实质上是一种隐式算法显式欧拉法是一个显式算法，虽然计算量较小，但是显式欧拉法是一个显式算法，虽然计算量较小，但是精度不高精度不高综合两种方法的长处，可以先用显式欧拉法求出综合

20、两种方法的长处，可以先用显式欧拉法求出 y(xn+1) 的一个粗略近似值，然后用它代入梯形法公式的一个粗略近似值，然后用它代入梯形法公式的右端，用梯形法计算的右端，用梯形法计算 y(xn+1) 的较为精确的近似值。的较为精确的近似值。1ny 预预报报值值1ny 校校正正值值24改进的欧拉法（续）改进的欧拉法（续）按照上述思想，可以建立如下预报按照上述思想，可以建立如下预报-校正系统：校正系统：1(,)nnnnyyh f xy 预预报报: :按以上两式求解常微分方程的算法称为改进的欧拉按以上两式求解常微分方程的算法称为改进的欧拉法，它还可以表示为：法，它还可以表示为： 11(,),(,)2nnn

21、nnnnnhyyf xyf xyh f xy嵌套形式嵌套形式11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 平均化形式平均化形式2 阶精度阶精度111 (,)(,)2nnnnnnhyyf xyf xy校校正正: :25用改进欧拉法求上例所述的初值问题并与欧拉法和用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。梯形法比较误差的大小。解解：采用改进欧拉法的嵌套形式：采用改进欧拉法的嵌套形式： 11(,),(,)2nnnnnnnnhyyf xyf xyh f xy10.1,0.5(0)1yyxhxy 计计算算到到 1(1)(,)12nnnnnn

22、nhyyxyh f xyx (1)(1)12nnnnnnnhyyxyhyxxh( 2)(2)122nnhhhhyxh 0.9050.0950.1nnyx26计算结果计算结果xn改进欧改进欧拉法拉法 yn精确解精确解 y (xn) 误误 差差改进欧拉法改进欧拉法欧拉法欧拉法梯形法梯形法0.1 1.005000 1.0048371.6 10-44.8 10-37.5 10-50.2 1.019205 1.0197312.9 10-48.7 10-31.4 10-40.3 1.041218 1.0408184.0 10-41.2 10-21.9 10-40.4 1.070802 1.0703204.

23、8 10-41.4 10-22.2 10-40.5 1.107076 1.1065315.5 10-41.6 10-22.5 10-4可见，改进欧拉法的误差数量级与梯形法大致相同，可见，改进欧拉法的误差数量级与梯形法大致相同，而比欧拉法小得多。而比欧拉法小得多。 10.1,0.5(0)1yyxhxy 计计算算到到27改进的欧拉法的意义改进的欧拉法的意义11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改进的欧拉法改进的欧拉法的平均化形式的平均化形式1k1nyhk 2k12112222nnnhhkkyykkyh 11()()()( )(,)nnnnny xy

24、xhy xhyxx y (xn 1) 在点在点 xn 处的一阶展开式为：处的一阶展开式为：28改进的欧拉法的几何意义改进的欧拉法的几何意义2nyhkpnxny0 xyrh1211(,)(,)nnnnkf xykf xyhk 1nx1nyhkq斜率: k21ny 1ny 1222nhhykks斜率: k1291121211()/2( ,)(,)nnnnnnyyh kkkf x ykf xyhk 龙格龙格-库塔库塔(runge-kutta)方法方法11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改进的欧拉法改进的欧拉法（2 2 阶精度）阶精度）11()()()

25、( )() , ( )(,)nnnnnny xy xhy xhyy xh fyxx y (xn 1) 在点在点 xn 处的一阶泰勒展开式为：处的一阶泰勒展开式为：*k1(,)nnnnyyhf xy 显式欧拉法显式欧拉法（1 1 阶精度）阶精度）111(,)nnnnyyhkkf xy 30龙格龙格-库塔方法（续）库塔方法（续）显式欧拉法用一个点的值显式欧拉法用一个点的值 k1 作为作为 k* 的近似值的近似值改进的欧拉公式用二个点的值改进的欧拉公式用二个点的值 k1 和和 k2 的平均值作为的平均值作为 k* 近似值；近似值；改进的欧拉法比显式欧拉法精度高；改进的欧拉法比显式欧拉法精度高；在在

26、xn, xn 1 内多预报几个点的内多预报几个点的 ki 值，并用其加权平均值，并用其加权平均值作为值作为 k* 的近似值从而构造出具有更高精度的计算公的近似值从而构造出具有更高精度的计算公式，这就是式，这就是龙格龙格- -库塔方法的基本思想。库塔方法的基本思想。31二阶龙格二阶龙格-库塔方法库塔方法以以 k1 和和 k2 的加权平均来近似取代的加权平均来近似取代 k*12221121211()(,)(,)nnnnnncyyhkkkf xykf xh yhcakb 12221,ccab为为待待定定系系数数为分析局部截断误差，令为分析局部截断误差，令 yn y(xn)，由泰勒公式得：，由泰勒公式

27、得：231()()()()()()2!nnnnnhy xy xhy xhy xyxo h (), ()(,)nnnnny xf xy xf xy 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyo h ()(,)(,)(,) (,)nnnxnnynnnnyxfxyfxyfxyf xyny 32补充：二元泰勒展开式补充：二元泰勒展开式00000020000(,)(,)(,)1(,)2!1(,)!nf xh ykf xyhkf xyxyhkf xyxyhkf xynxy 20000200(,)2(,)(,)xxxyyyh fxyhk

28、 fxyk fxy 000nniin inin i x xyyifc h kx y 0000(,)(,)xyhfxyk fxy 3322211(,)nnkf xa h yb hk用二元泰勒公式展开用二元泰勒公式展开222211(,)(,)(,)()nnxnnynnkf xya h fxyb hk fxyo h123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyo h 11223211(,) (,)(,)(,)()nnnnnnxnnynnyyh c f xycf xya h fxyb hk fxyo h 将将 k1, k2 代入代入

29、中可得：中可得：11122()nnyyh c kc k 2122223221() (,)(,)(,) (,)()nnnxnnynnnnyh ccf xyh c a fxyh c b fxyf xyo h122222213() (,)2(,)2(,) (,)2()nnnxnnynnnnyh ccf xyhc a fxyc b fxyf xyo h 1(,)nnkf xy 341122322221() (,)2(,)2(,) (,)()2nnnnxnnynnnnyyh ccf xyhc a fxyc b fxyf xyo h 123()()(,)(,)(,) (,)()2nnnnxnnynnnny

30、 xy xhf xyhfxyfxyf xyo h 二阶龙格二阶龙格-库塔方法（续）库塔方法（续）122222112121ccc ac b 311()()nny xyo h2 阶精度阶精度35四个未知变量，只有三四个未知变量，只有三个方程，有无穷多组解个方程，有无穷多组解每组解的构成的龙格每组解的构成的龙格-库库塔方法均为二阶塔方法均为二阶11122122211()(,)(,)nnnnnnyyh c kc kkf xykf xa h yb hk 122211 21ccab 取取二阶龙格二阶龙格- -库塔方法库塔方法即为改进的欧拉方法即为改进的欧拉方法12221011 2ccab 取取12121(

31、,),22nnnnnnyyhkkf xyhhkfxyk 变形的欧拉法变形的欧拉法中中 点点 方方 法法122222112121ccc ac b36三阶龙格三阶龙格-库塔方法库塔方法三阶龙格三阶龙格-库塔方法是用三个值库塔方法是用三个值 k1, k2, k3 的加权平均来的加权平均来近似取代近似取代 k*111232213313231213122()(,)(,)(,)nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhcccababkb hk 要使三阶龙格要使三阶龙格-库塔方法具有三阶精度，必须使其局部库塔方法具有三阶精度，必须使其局部截断误差为截断误差为 o(h4)将将 k

32、1, k2, k3 代入代入 yn 1 的表达式中，在的表达式中，在 (xn, , yn) 处用二元处用二元泰勒公式展开，与泰勒公式展开，与 y(xn 1) 在在 xn 处的泰勒展开式比较处的泰勒展开式比较37三阶龙格三阶龙格-库塔方法（续）库塔方法（续）类似二阶龙格类似二阶龙格-库塔方法的推导过程，库塔方法的推导过程，8 个待定系数个待定系数 c1, c2, c3, a2, a3, b21, b31, b32 应满足：应满足： 123221331322233222233232211 21 31 6cccababbc ac ac ac ac b a 8 个未知参数，个未知参数，6 个个方程，有

33、无穷多组解方程，有无穷多组解11231213121416661122()(,)(,)()2,nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhkhk 库塔公式库塔公式38四阶龙格四阶龙格-库塔方法库塔方法类似可以推出四阶龙格类似可以推出四阶龙格-库塔公式，常用的有：库塔公式，常用的有：1123122166661122141213122243()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkkf xh yhk 标准四阶龙格标准四阶龙格- -库塔公式库塔公式39四阶龙格四阶龙格-库塔方法（续）库塔方法（续）222211

34、666611222 12212222222112341213124232()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkhkkf xh yhkhk 吉尔（吉尔（gill）公式）公式4 阶以上阶以上龙格龙格-库塔公式的计算量太大，并且精度不库塔公式的计算量太大，并且精度不一定提高，有时反而会降低，因此实际应用中一般一定提高，有时反而会降低，因此实际应用中一般选用四阶龙格选用四阶龙格-库塔已足可满足精度要求。库塔已足可满足精度要求。40用经典四阶龙格用经典四阶龙格- -库塔库塔方法求解前例的初值问题，并方法求解前例的初值问题，并与改进与改

35、进 欧拉欧拉 法、梯形法在法、梯形法在 x5 0.5 处比较其误差大小处比较其误差大小解解：采用经典四阶龙格：采用经典四阶龙格- -库塔库塔公式：公式：10000(,)10kf xyyx 10.1,0.5(0)1yyxhxy 计计算算到到112220011122010(,)()()10.05kf xh yhkyhkxh 112230021122020(,)()()10.0475kf xh yhkyhkxh 4003030(,)()()10.09525kf xh yhkyhkxh 411221666610123412216666()10.1(00.050.04750.09525)1.004837

36、5yyhkkkk 于于是是：23451.01873091.040818421.070320291.1065309yyyy 同同理理可可计计算算：四阶四阶r-k方法方法的精度比二阶的精度比二阶方法高得多方法高得多10.1,0.5(0)1yyxhxy 计计算算到到0.55( )e()0.5e1.106530660 xy xxy x 精确解为：精确解为：r-k方法的误差：方法的误差：755()2.4 10yy x 改进欧拉法的误差：改进欧拉法的误差：45.5 10 梯形法的误差：梯形法的误差：42.5 10 42变步长的龙格变步长的龙格-库塔方法库塔方法( )5111()hnnnnty xyc h设

37、设 y (xn) 在在 xn 处的值处的值 yn y (xn)，当，当 xn+1 xn+ h 时时 y (xn 1) 的近似值为的近似值为 ，由于四阶，由于四阶 r-k 方法的精度方法的精度为为 4 阶，故局部截断误差为：阶，故局部截断误差为：( )1hny 用四阶用四阶r-k方法求解初值问题精度较高，但要从理论方法求解初值问题精度较高，但要从理论上给出误差上给出误差 | | y (xn) yn | | 的估计式则比较困难；那么的估计式则比较困难；那么应如何判断计算结果的精度以及如何选择合适的步应如何判断计算结果的精度以及如何选择合适的步长长 h？通常是通过不同步长在计算机上的计算结果进行近通

38、常是通过不同步长在计算机上的计算结果进行近似估计。似估计。43若以若以 h/2 为步长，从为步长，从 xn 出发，经过两步计算，得到出发，经过两步计算，得到y(xn+1) 的近似值的近似值(2)1hny 变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）以上每步的截断误差约为以上每步的截断误差约为 cn(h/2)5，于是两步的局部，于是两步的局部截断误差为：截断误差为：(2)5111()2(2)hnnnnty xyc h(2)511( )511()2(2)1()16hnnnhnnny xyc hy xyc h 于是：于是：整理得：整理得：(2)(2)( )11111()15hhhnnnny

39、 xyyy44变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）记：记： ，给定的精度要求为，给定的精度要求为 e e(2)( )11115hhnnyy (2)1hny e e，反复将步长折半计算，直至，反复将步长折半计算，直至 e e，取最终得，取最终得到的到的 作为作为 y(xn+1) 的近似值。的近似值。 e e，反复将步长加倍计算，直至，反复将步长加倍计算，直至 e e，再将步长，再将步长折半一次计算，最终得到符合精度要求的折半一次计算，最终得到符合精度要求的 y(xn+1) 的的近似值。近似值。45单步法的收敛性单步法的收敛性显式单步法可统一写成：显式单步法可统一写成：1(, )

40、nnnnyyhxy h 增量函数，仅依赖于函数增量函数，仅依赖于函数 f，且仅仅是且仅仅是 xn, yn, h 的函数的函数求求 y y(x)求求 y(xn)，xn x0 nh离散化离散化求求 yn y(xn)某种数值方法某种数值方法h 0 时，近似解时，近似解是否收敛到精确解是否收敛到精确解，它应当，它应当是一个固定节点，因是一个固定节点，因此此 h 0 时应同时附时应同时附带带 n 0nxxnh( , , )x y h 46单步法的收敛性（续）单步法的收敛性（续）对于对于 p 阶的常微分方程数值算法，当阶的常微分方程数值算法，当 h 0， n 时，是否时，是否 yn+1 y(xn+1)？p

41、 阶算法的局部截断误差为：阶算法的局部截断误差为：111()()pnny xyo h 显然：显然：110,lim()0nnhny xy局部截断误差的前提假设是：局部截断误差的前提假设是： ()nnyy x 局部截断误差局部截断误差 0 并不能保证算法收敛并不能保证算法收敛47单步法的收敛性（续）单步法的收敛性（续）定义定义：若求解某初值问题的单步数值法，对于固定的：若求解某初值问题的单步数值法，对于固定的 当当 h 0 且且 n 时，它的近似时，它的近似 解趋向于精确解解趋向于精确解 y(xn)，即：，即：0nxxnh0,()limnnhnyy x 则称该单步法是收敛的。则称该单步法是收敛的。

42、定义定义：称：称 y(xn) yn 为单步法的近似解为单步法的近似解 yn 的的整体截断整体截断 误差误差。单步法收敛单步法收敛h 0 且且 n 时，时，yn 的整体截断误差的整体截断误差 0 48单步法的收敛性（续）单步法的收敛性（续）收敛性定理收敛性定理若某单步法满足以上条件，则该方法是收敛的若某单步法满足以上条件，则该方法是收敛的则该单步法的整体截断误差为：则该单步法的整体截断误差为：若单步法若单步法 具有具有 p 阶精度，阶精度，且增量函数且增量函数 关于关于 y 满足：满足：1(, )nnnnyyhxy h ( , , )x y h ()()pnny xyo h( , , )( ,

43、, )x y hx y hl yy lipschitz 条件：条件： 初值初值 y0 是准确的是准确的49假设在前一步假设在前一步 yn 准确的前提下求得的近似值为：准确的前提下求得的近似值为：( , , )( , , )x y hx y hl yy 1(, )nnnnyyhxy h 1,()()nnnnyhxhy xy x 算法精度为算法精度为 p 阶，局部截断误差：阶，局部截断误差：111()pnny xych 111111() ()nnnnnny xyy xyyy1111()nnnny xyyy1 (), (), )(, )pnnnnnnchy xyhxy xhhxy h 1(), ()

44、, )(, )pnnnnnnchy xyhxy xhxy h 1()()pnnnnchy xyhl y xy 1(1)()pnnhly xych ne e1ne e 5011(1)pnnhlch eeee 11(1)pnnhlch eeee 11(1) (1)ppnhlhlchche e 2 221(1)(1)1pnhlhlche e 2 22113(1)(1)(1)1ppnhlhlchhlche e 3213(1)(1)(1)1pnhlhlhlche e 110(1)(1)(1)1nnphlhlhlche e 101(1)1(1)(1)1nnphlhlchhl e e 0(1)(1)1pnn

45、chhlhll e e510(1)(1)1pnnnchhlhll eeee2ee112xxxx (1)ennxx即：即： (1)eeenhlnhltlnhl 0(1)e1ptlnnchhll eeee若初值是准确的，则若初值是准确的，则 e e 0 0 ，从而整体截断误差为：，从而整体截断误差为： e1ptlnchl e e0nxxnhty e x 为单调增函数，当为单调增函数，当 时时当当 h 0 且且 n 时时，则，则0ne e()po h52单步法的稳定性单步法的稳定性在讨论单步法收敛性时一般认为数值方法本身的计在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并非如此：

46、算过程是准确的，实际上并非如此：初始值初始值 y0 有误差有误差 d d y0 y(x0)后续的每一步计算均有舍入误差后续的每一步计算均有舍入误差这些初始和舍入误差在计算过程的传播中是逐步衰这些初始和舍入误差在计算过程的传播中是逐步衰减的还是恶性增长就是数值方法的稳定性问题减的还是恶性增长就是数值方法的稳定性问题53定义定义：若一种数值方法在节点：若一种数值方法在节点 xn 处的数值解处的数值解 yn 的扰的扰动动 ，而在以后各节点，而在以后各节点 ym (m n) 上产生的上产生的扰动为扰动为 ，如果：，如果：md d0nd d 单步法的稳定性（续）单步法的稳定性（续）(1,2,)mnmnn

47、dddd定义定义：设在节点：设在节点 xn 处用数值算法得到的理想数值解为处用数值算法得到的理想数值解为 yn，而实际计算得到的近似解为，而实际计算得到的近似解为 ，称差值：，称差值：ny nnnyyd d 为第为第 n 步的数值解的步的数值解的扰动扰动。则称该数值方法是则称该数值方法是稳定稳定的。的。54单步法的稳定性（续）单步法的稳定性（续）欧拉法：欧拉法：1(,)nnnnyyhf xy 由于函数由于函数 f (x, y) 的多样性，数值稳定性的分析相当复的多样性，数值稳定性的分析相当复杂，通常只研究模型方程杂，通常只研究模型方程(0)yy 考察模型方程：考察模型方程：(0)yy 1(1)

48、nnyhy 即：即：假设在节点值假设在节点值 yn 上有扰动上有扰动 d dn，在节点值，在节点值 yn 1 上有上有扰动扰动 d dn 1，且，且 d dn 1 仅由仅由 d dn 引起（即：计算过程中引起（即：计算过程中不再引起新的误差）不再引起新的误差）55111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d 欧拉法稳定欧拉法稳定1nndddd 11h 即：即：111h 欧拉法稳定的条件：欧拉法稳定的条件：20h 0 1(1)nnyhy 针对模型方程：针对模型方程：的显式欧拉法：的显式欧拉法：1(,)nnnnyyhf xy 化简得：化简得：20h 56隐式

49、欧拉法：隐式欧拉法：111(,)nnnnyyhf xy111111nnnnnnyyyyhhhd dd d 考察模型方程：考察模型方程：(0)yy 即：即：11nnnyyh y 化简为：化简为：11nnyyh 假设假设 yn 上有扰动上有扰动 ，则，则 yn 1 的扰动为：的扰动为：nnnyyd d 隐式欧拉法稳定隐式欧拉法稳定1nndddd 111h 0 0h ，上式均成立，所以：，上式均成立，所以：隐式欧拉法稳定是恒稳定的隐式欧拉法稳定是恒稳定的57一阶常微分方程组一阶常微分方程组11122212121012020( ,)( ,).( ,)(),(),()mmmmmmmyfx yyyyfx

50、yyyyfx yyyy xsyxsyxs 111222( ,)( )( )( ,)( )( ,)( )( ,)ttttmmmfx yy xsyxsfx yy xf x ysyxsfx y 0( )( ,)()tyxf x yy xs 5811121012212202( ,)()( ,)()yfx yyy xsyfx yyyxs 1(,)tnnnnyyhf x y 111111222212(,)(,)nnnnnnnnnnyyfxyyhyyfxyy 111222( )( ,)( )( ,)( )( ,)ttty xfx ysy xf x ysyxsfx y11(,)tnnnnyyf x y111111111112222112(,)(,)nnnnnnnnnnyyfxyyhyyfxyy 111111111121112222122112(,)(,