第五章常微分方程数值解法PPT课件

1、第五章常微分方程数值解法2常微分方程数值解法引言（常微分方程数值解法概述）引言（常微分方程数值解法概述）显式欧拉法、隐式欧拉法、二步显式欧拉法、隐式欧拉法、二步欧拉法欧拉法局部截断误差与精度局部截断误差与精度改进的欧拉方法改进的欧拉方法3常微分方程数值解法龙格龙格-库塔方法库塔方法收敛性与稳定性简述收敛性与稳定性简述一阶常微分方程组与高阶常微一阶常微分方程组与高阶常微分方程分方程4differential equati.()on1微分凡含有变量和未知函数的导数 或微分的方程定称为(方程)义tanyyyxx例：引言引言5ordinary differential equationpartial

2、differential equation 未知函数是一元函数的微分方程称为简称为微分方程;未知函数常微分方程偏微分方程()是多元函数的微分方程称为,().6微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数都是一次的,则称它是线性的,否则称为非线性的)(,nyyy 7 在高等数学中，对于常微分方程的在高等数学中，对于常微分方程的求解，给出了一些典型方程求求解，给出了一些典型方程求解析解解析解的的基本方法，基本方法，如如: :可可分离变量法、常系数分离变量法、常系数齐次线性方程的解法、常系数非齐次线齐次线性方程的解法、常系数非齐次线性方程的解法等性方程的解法

3、等。8但能求解的常微分方程仍然是有限的，但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析大多数的常微分方程是不可能给出解析解。解。 例如例如 22yxy 这个一阶微分方程就不能用初等函这个一阶微分方程就不能用初等函数及其积分来表达它的解。数及其积分来表达它的解。 9引言引言一阶常微分方程初值问题：一阶常微分方程初值问题：00( , )()yf x yy xy 微分方程微分方程初始条件初始条件10定理：定理：若若 f (x, y) 在某闭区域在某闭区域 R ： 00|,|(0,0)xxayybab上 连 续 ， 且 在上 连 续 ， 且 在 R 域 内 满 足 李 普 希 兹

4、域 内 满 足 李 普 希 兹 (Lipschitz) 条件，即存在正数条件，即存在正数 L，使得对，使得对于于 R 域内的任意两值域内的任意两值 y1, y2，下列不等式，下列不等式成立：成立： 1221|( ,)( ,)|f x yf x yL yy则上述初值问题的连续可微的解则上述初值问题的连续可微的解 y(x) 存在存在并且唯一。并且唯一。11 虽然求解常微分方程有各种各样虽然求解常微分方程有各种各样的的解析方法解析方法，但解析方法只能用来求，但解析方法只能用来求解一些特殊类型的方程，实际问题中解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠归结出来的微分方程主要靠数值解法数值

5、解法.12所谓所谓数值解法数值解法, 就是寻求解就是寻求解y(x)在一系列在一系列离散节点离散节点 121nnxxxx上的近似值上的近似值 y1, y2, , yn, yn+1,. 相邻两个相邻两个节点的间距节点的间距hn=xn+1- -xn称为称为步长步长. 今后如不今后如不特别说明，总是假定特别说明，总是假定 hi=h(i=1,2,)为为定数定数, 这时节点为这时节点为xn=x0+nh(i=0,1,2,) (等距节点等距节点).13解析解与数值解解析解与数值解解的类型解的类型解的形式解的形式解的区间解的区间变量变量解的性质解的性质解析解y(x)连续精确解数值解y0,y1,.yna=x0,

6、x1, xn=b离散近似解14一阶微分方程的初值问题：一阶微分方程的初值问题：)2 . 1 (.)() 1 . 1 (),(00yxyyxfy两边取定积分有两边取定积分有可见可见,一阶微分方程的初值问题转换为计算一阶微分方程的初值问题转换为计算定积分问题定积分问题,采用不同的定积分数值解法就采用不同的定积分数值解法就会得到不同的一阶微分方程的数值解法会得到不同的一阶微分方程的数值解法。, bax.)(,()()(11nnxxnndttytfxyxy15 初值问题的初值问题的数值解法数值解法有个有个基本特点基本特点，他们都采取他们都采取“步进式步进式”，即求解过程顺着，即求解过程顺着节点排列的次

7、序一步一步地向前推进节点排列的次序一步一步地向前推进. 描描述这类算法，只要给出用已知信息述这类算法，只要给出用已知信息yn,yn- -1,yn- -2,计算计算yn+1的递推公式的递推公式.数值方法的基本思想数值方法的基本思想16首先首先，要对微分方程离散化，建立求解，要对微分方程离散化，建立求解数值解的数值解的递推公式递推公式. p一类：是一类：是计算计算yn+1时只用到前一点的值时只用到前一点的值yn，称为，称为单步法单步法. p另一类：是另一类：是用到用到yn+1前面前面 k 点的值点的值yn, yn-1, yn-k+1，称为，称为k步法步法. 数值方法的基本思想数值方法的基本思想17

8、其次，要研究公式的其次，要研究公式的局部截断误差局部截断误差和和阶阶，数值解数值解yn与精确解与精确解y(xn)的的误差估计误差估计及及收敛性收敛性，还有递推公式的，还有递推公式的计算稳定性计算稳定性等等问题问题.数值方法的基本思想数值方法的基本思想18一阶微分方程的初值问题：一阶微分方程的初值问题：)2 . 1 (.)() 1 . 1 (),(00yxyyxfy 121nnxxxx求对应的：求对应的：欧拉方法欧拉方法19数值积分法数值积分法1111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 将微

9、分方程将微分方程 y f (x, y) 在区间在区间 xn, xn+1 上积分：上积分：欧拉方法欧拉方法20由于在逐步求解的过程中，由于在逐步求解的过程中，y(xn) 的准确的准确值无法求解出来，因此用其近似值代替。值无法求解出来，因此用其近似值代替。为避免混淆，以下学习简记：为避免混淆，以下学习简记：欧拉方法欧拉方法y(xn)：待求函数：待求函数 y(x) 在在 xn 处的精确函数处的精确函数值值yn ：待求函数：待求函数 y(x) 在在 xn 处的近似函数值处的近似函数值21以以近似值近似值 yn 代替精确值代替精确值 y(xn) 可得可得：100(,)0,1,2,()nnnnyyhf x

10、ynyy x 根据根据 y0 可以一步步计算出函数可以一步步计算出函数 y y(x) 在在 x1, x2, x3 x4, 上的近似值上的近似值 y1, y2, y3, y4 , 欧拉公式欧拉公式欧拉方法欧拉方法22常微分方程数值解是一组离散的函数值常微分方程数值解是一组离散的函数值数据，它的精确表达式很难求解得到，数据，它的精确表达式很难求解得到，但可以进行插值计算后用但可以进行插值计算后用插值函数插值函数逼近逼近 y(x)欧拉方法欧拉方法23欧拉方法的几何意义欧拉方法的几何意义xy00 x1x2x3x4x5x6x0P1P2P3P4P5P6P( )yy x24隐式欧拉法隐式欧拉法在数值积分法推

11、导中，积分的近似值在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数取为积分区间宽度与右端点处的函数值乘积，即：值乘积，即：11111111d()( ) , ( )d() , (), ()nnnnxxnnxxnnnnnny xy xy xf x y xxxx f xy xhf xy x 25这样便得到了隐式欧拉法：这样便得到了隐式欧拉法：11100(,)()nnnnyyh f xyyy x 含有未知含有未知的函数值的函数值隐式欧拉法没有显式欧拉法方隐式欧拉法没有显式欧拉法方便便隐式欧拉法隐式欧拉法26二步欧拉法二步欧拉法在数值积分法推导中，积分区间宽度选在数值积分法推导中，积分区

12、间宽度选为两步步长，即积分区间为：为两步步长，即积分区间为：xn 1, xn 1，则：则：11111111d()() , ( )d() , ()2, ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xh f xy x 中矩形公式中矩形公式11002(,)()nnnnyyhf xyyy x 27以以 y(x) 在在 xn 1, xn 上的近似值上的近似值代替代替精确精确值值可得可得：需要前两需要前两步的计算步的计算结果结果二步欧拉法二步欧拉法28梯形公式欧拉法梯形公式欧拉法在数值积分法中，如果用在数值积分法中，如果用梯形公式梯形公式近似计近似计算算 f (x,

13、y) 在区间在区间 xn, xn+1 上的积分，即：上的积分，即：用近似值代替精确值可得梯形公式欧拉法：用近似值代替精确值可得梯形公式欧拉法：111 (,)(,)2nnnnnnhyyf xyf xy 111111, (), () , ( )d()2, (), ()2nnxnnnnnnxnnnnf xy xf xy xf x y xxxxhf xy xf xy x 29上式右端出现了未知项，可见梯形法是上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的欧拉法是显式欧拉法与隐式欧拉法的算算术平均术平均。梯形公式欧拉

14、法梯形公式欧拉法111 (,)(,)2nnnnnnhyyf xyf xy30例例用显式欧拉法、隐式欧拉法、梯形法求解用显式欧拉法、隐式欧拉法、梯形法求解初值问题：初值问题：1(0)1yyxy 取取 h 0.1，计算到，计算到 x 0.5，并与，并与精确解进精确解进行行比较比较31解：由已知条件可得：解：由已知条件可得：h 0.1，x0 0， y0 1，f (x, y) y x 1显式欧拉法：显式欧拉法：1(1)(1)0.90.10.1nnnnnnnnyyhyxh yhxhyx 32例：（续）例：（续）隐式欧拉法：隐式欧拉法：111(1)nnnnyyhyx化简得：化简得：11(1)(1)(1)n

15、nnnnh yyh xyh xh11(1)(0.10.11) 1.1(1)nnnnnyyh xhyxh 33111(1)(1)2nnnnnnhyyyxyx 梯形公式欧拉法：梯形公式欧拉法：111(2)1(2)(2)(2)nnnnhyhyhxx11(2)(2) (2)(1.90.20.21) 2.1nnnnnnyh yh xxhyx34计算结果：计算结果：xn显式法显式法 yn隐式法隐式法 yn梯形法梯形法 yn精确解精确解 y (xn)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.3

16、1.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 显式法显式法隐式法隐式法梯形法梯形法35局部截断误差局部截断误差为了简化分析某常微分方程数值算法的为了简化分析某常微分方程数值算法的误差，现假设误差，现假设 yn y(xn)，即，即在前一步在前一步 yn 准确的前提下准

17、确的前提下，估计：，估计：111()nnnTy xy称上述误差称上述误差 Tn 1 为该常微分方程数值算为该常微分方程数值算法的法的局部截断误差局部截断误差36如果某个常微分方程数值算法的局部截如果某个常微分方程数值算法的局部截断误差可表示为断误差可表示为 O(h p 1)，则称该数值算，则称该数值算法的精度是法的精度是 p 阶阶欧拉法的精度为一阶；二步欧拉法的精欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形公式欧拉法的精度为二度为二阶；梯形公式欧拉法的精度为二阶。阶。局部截断误差局部截断误差37泰勒展开法泰勒展开法如果初值问题中的如果初值问题中的 f (x, y) 充分可微，则可将充分可微

18、，则可将 y(xn 1) 在点在点 xn 处展开：处展开：21()()()()()2!nnnnnhy xy xhy xhy xyx 如果只保留线性项，忽略如果只保留线性项，忽略 h2 及以上各项，则：及以上各项，则：11()()()nnnny xy xhy xy (), ()nnny xf xy x ny1(,)nnnnyyhf xy 显式欧拉公式显式欧拉公式38局部截断误差的分析局部截断误差的分析利用泰勒公式展开，比较各算法与展开利用泰勒公式展开，比较各算法与展开式的前几项式的前几项显式欧拉法的局部截断误差：显式欧拉法的局部截断误差：1(), ()()()nnnnnnyy xhf xy xy

19、 xhy x 1(,)nnnnyyhf xy 欧拉法欧拉法()nnyy x 令令：39将将 y(xn 1) 在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：211( )()()()()2!(,)nnnnnnyy xy xhy xhy xhxx 22111( )()()2!nnnyTy xyhO h 1 阶精度阶精度1(), ()()()nnnnnnyy xhf xy xy xhy x 局部截断误差局部截断误差显式欧拉法的局部截断误差显式欧拉法的局部截断误差40所以，所以， 隐式隐式Euler方法也具有方法也具有 1 阶精度。阶精度。 , ! 2)()()()()(2111hyxyhxyhx

20、yxynnnn )(nxy1nx将将 在在点点 处处一阶一阶Taylor展开展开 )()()(211hOxyhxynn)( )()()()()(2121111hOxyhhOxyhxyxyTnnnnn111(,)nnnnyyhf xy隐式欧拉法：隐式欧拉法：1 阶精度阶精度隐式欧拉法隐式欧拉法的局部截断误差：的局部截断误差：41111()2, ()()2()nnnnnnyy xhf xy xy xhy x 分别将分别将 y(xn 1), y(xn 1) 在在 xn 点处用泰勒公式展开：点处用泰勒公式展开：231231()()()()()2!()()()()()()2!nnnnnnnnyxy xy

21、 xhy xhO hyxy xy xhy xhO h 111()nnnTy xy二步欧拉法的局部截断误差：二步欧拉法的局部截断误差：112(,)nnnnyyhf xy二步欧拉法：二步欧拉法：2 阶精度阶精度1()ny x ()ny x311()()2()()nnny xy xhy xO h 42111 (,)(,)2nnnnnnhyyf xyf xy梯形公式欧拉法：梯形公式欧拉法：()ny x2 阶精度阶精度112312312111: ()()( )( )()( )( )2!3!( )( )()( )( )2!3!( )()( )( )2!()( )( )()2iiiiiiiiiiiiiiii

22、iiy xy xxTaylory xyy xy xhy xhhy xyy xy xhy xhhyy xy xhy xhhy xy xy xy xhy梯形公式和在点 处展开:则2323( )( )( )2!3!( )( )( )( )( )22!12iiiiiy xyxhhhyyy xy xhy xhh 43各种欧拉法的比较各种欧拉法的比较方法方法精度精度评述评述显式欧拉法显式欧拉法1最简单，精度低最简单，精度低隐式欧拉法隐式欧拉法1不便计算，稳定性好不便计算，稳定性好二步欧拉法二步欧拉法2需要两步初值，且第需要两步初值，且第 2 个初值只能由其个初值只能由其它方法给出，可能对后面的递推精度有它

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

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

25、的欧拉法欧拉法2 阶精度阶精度47改进改进的的欧拉法还欧拉法还可以表示为：可以表示为： 11(,),(,)2nnnnnnnnhyyf xyf xyh f xy嵌套形式嵌套形式11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 平均化形式平均化形式48用改进欧拉法求上例所述的初值问题并与欧拉法和用改进欧拉法求上例所述的初值问题并与欧拉法和梯形法比较误差的大小。梯形法比较误差的大小。解解：采用改进欧拉法的嵌套形式：采用改进欧拉法的嵌套形式： 11(,),(,)2nnnnnnnnhyyf xyf xyh f xy10.1,0.5(0)1yyxhxy 计计算算到到 1

26、(1)(,)12nnnnnnnhyyxyh f xyx (1)(1)12nnnnnnnhyyxyhyxxh( 2)(2)122nnhhhhyxh 0.9050.0950.1nnyx49计算结果计算结果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.070

27、802 1.0703204.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 计计算算到到502(01)(0)10.1,.():12EulerEulerxyyyxyhyx 用公式和改进的公式计算:的近似解,取并与精确值比较本方程的精确解为例25110(0)1(0)11102( , )2()12() ( ,)(,)21iiiiiiiiiiiiiii

28、ixf x yyyxyyh yyyxyyh yyhyyf x yf xyy解：故欧拉公式为梯形公式的预报校正公式为5253改进的欧拉法的意义改进的欧拉法的意义11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改 进 的 欧改 进 的 欧拉 法 的 平拉 法 的 平均 化 形 式均 化 形 式1k1nyhk 2k12112222nnnhhkkyykkyh 11()()()( )(,)nnnnny xy xhy xhyxx y (xn 1) 在点在点 xn 处的一阶展开式为：处的一阶展开式为：54改进的欧拉法的几何意义改进的欧拉法的几何意义2nyhkPnxny0

29、 xyRh1211(,)(,)nnnnkf xykf xyhk 1nx1nyhkQ斜率: k21ny 1ny 1222nhhykkS斜率: k1551121211()/2( ,)(,)nnnnnnyyh kkkf x ykf xyhk 龙格龙格-库塔库塔(Runge-Kutta)方法方法11(,)(,)() 2pnnncnnpnpcyyh f xyyyh f xyyyy 改进的欧拉法改进的欧拉法（2 2 阶精度）阶精度）11()()()( )() , ( )(,)nnnnnny xy xhy xhyy xh fyxx y (xn 1) 在点在点 xn 处的一阶泰勒展开式为：处的一阶泰勒展开式为

30、：*k1(,)nnnnyyhf xy 显式欧拉法显式欧拉法（1 1 阶精度）阶精度）111(,)nnnnyyhkkf xy 56龙格龙格-库塔方法（续）库塔方法（续）显式欧拉法用一个点的值显式欧拉法用一个点的值 k1 作为作为 k* 的近似值的近似值改进的欧拉公式用二个点的值改进的欧拉公式用二个点的值 k1 和和 k2 的平均的平均值作为值作为 k* 近似值；近似值；改进的欧拉法比显式欧拉法精度高；改进的欧拉法比显式欧拉法精度高；在在 xn, xn 1 内多预报几个点的内多预报几个点的 ki 值，并用其加值，并用其加权平均值作为权平均值作为 k* 的近似值从而构造出具有更高的近似值从而构造出具

31、有更高精度的计算公式，这就是精度的计算公式，这就是龙格龙格- -库塔方法的基库塔方法的基本思想。本思想。57二阶龙格二阶龙格-库塔方法库塔方法以以 k1 和和 k2 的加权平均来近似取代的加权平均来近似取代 k*12221121211()(,)(,)nnnnnncyyhkkkf xykf xh yhcakb 12221,ccab为为待待定定系系数数确定这些系数，使得该公式具有确定这些系数，使得该公式具有2阶精度阶精度58为分析局部截断误差，令为分析局部截断误差，令 yn y(xn)，由泰，由泰勒公式得：勒公式得：231()()()()()()2!nnnnnhy xy xhy xhy xyxO

32、h (), ()(,)nnnnny xf xy xf xy 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h ()(,)(,)(,) (,)nnnxnnynnnnyxfxyfxyfxyf xyny 二阶龙格二阶龙格-库塔方法库塔方法59补充：二元泰勒展开式补充：二元泰勒展开式00000020000(,)(,)(,)1(,)2!1(,)!nf xh ykf xyhkf xyxyhkf xyxyhkf xynxy 20000200(,)2(,)(,)xxxyyyh fxyhk fxyk fxy 000nniin inin i

33、x xyyifC h kx y 0000(,)(,)xyhfxyk fxy 6022211(,)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 代入代入 中中可得：可得：11122()nnyyh c kc

34、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 611122322221() (,)2(,)2(,) (,)()2nnnnxnnynnnnyyh ccf xyhc a fxyc b fxyf xyO h 123()()(,)(,)(,) (,)()2nnnnxnnynnnny xy xhf xyhfxyfxyf xyO h 二

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

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

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

38、16661122()(,)(,)()2,nnnnnnnnyyhkkkkf xykf xh yhkkf xh yhkhk 库塔公式库塔公式65四阶龙格四阶龙格-库塔方法库塔方法类似可以推出四阶龙格类似可以推出四阶龙格-库塔公式，常用的有：库塔公式，常用的有：1123122166661122141213122243()(,)(,)(,)(,)nnnnnnnnnnyyhkkkkkf xykf xh yhkkf xh yhkkf xh yhk 标准四阶龙格标准四阶龙格- -库塔公式库塔公式66四阶龙格四阶龙格-库塔方法（续）库塔方法（续）222211666611222 1221222222211234

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

40、 欧拉欧拉 法、梯形法在法、梯形法在 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 1221666610123412216666()10.1(00.050.04750.09525)1.0048375yy

41、hkkkk 于于是是：23451.01873091.040818421.070320291.1065309yyyy 同同理理可可计计算算：10.1,0.5(0)1yyxhxy 计计算算到到0.55( )e()0.5e1.106530660 xy xxy x 精确解为：精确解为：70四阶四阶R-K方法方法的精度比二阶的精度比二阶方法高得多方法高得多R-K方法的误差：方法的误差：755()2.4 10yy x 改进欧拉法的误差：改进欧拉法的误差：45.5 10 梯形法的误差：梯形法的误差：42.5 10 71变步长的龙格变步长的龙格-库塔方法库塔方法用四阶用四阶R-K方法求解初值问题精度较高，方法

42、求解初值问题精度较高，但要从理论上给出误差但要从理论上给出误差 | | y (xn) yn | | 的估的估计式则比较困难计式则比较困难；单从每一步看，步长越小，截断误差就单从每一步看，步长越小，截断误差就越小。但随着步长的缩小，计算量增加，越小。但随着步长的缩小，计算量增加，舍入误差可能积累。舍入误差可能积累。72变步长的龙格变步长的龙格-库塔方法库塔方法那么应如何判断计算结果的精度以及如何那么应如何判断计算结果的精度以及如何选择合适的步长选择合适的步长 h？通常是通过不同步长在计算机上的计算结通常是通过不同步长在计算机上的计算结果进行近似估计。果进行近似估计。73( )5111()hnnn

43、nTy xyc h设设 y (xn) 在在 xn 处的值处的值 yn y (xn)，当，当 xn+1 xn+ h 时时, y (xn 1) 的近似值为的近似值为 ，由于四，由于四阶阶 R-K 方法的精度为方法的精度为 4 阶，故局部截断阶，故局部截断误差为：误差为：( )1hny 变步长的龙格变步长的龙格-库塔方法库塔方法若以若以 h/2 为步长，从为步长，从 xn 出发，经过两步计出发，经过两步计算，算，得到得到y(xn+1) 的近似值的近似值(2)1hny 74变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）以上每步的截断误差约为以上每步的截断误差约为 cn(h/2)5，于是，于是两

44、步的局部截断误差为：两步的局部截断误差为：(2)5111()2(2)hnnnnTy xyc h(2)511( )511()2(2)1()16hnnnhnnny xyc hy xyc h 于是：于是：整理得：整理得：(2)(2)( )11111()15hhhnnnny xyyy75变步长的龙格变步长的龙格-库塔方法（续）库塔方法（续）记：记： ，给定的精度要求给定的精度要求为为 e e(2)( )11115hhnnyy e e，反复将步长折半计算，直至，反复将步长折半计算，直至 e e，取，取最终得到最终得到的的 作为作为 y(xn+1) 的近似值。的近似值。 (2)1hny e e，反复将步长

45、加倍计算，直至，反复将步长加倍计算，直至 e e，再，再将步长折半一次计算，最终得到符合精度要求将步长折半一次计算，最终得到符合精度要求的的 y(xn+1) 的近似值。的近似值。76单步法的收敛性单步法的收敛性显式单步法可统一写成：显式单步法可统一写成：1(, )nnnnyyhxy h 增量函数，仅依赖于函数增量函数，仅依赖于函数 f，且仅仅是且仅仅是 xn, yn, h 的函数的函数求求 y y(x)求求 y(xn)，xn x0 nh离散化离散化求求 yn y(xn)某种数值方法某种数值方法h 0 时，近似解时，近似解是否收敛到精确解是否收敛到精确解 ，它应当它应当是一个固定节点，因是一个固

46、定节点，因此此 h 0 时应同时附时应同时附带带 n 0nxxnh ( , , )x y h 77单步法的收敛性（续）单步法的收敛性（续）对于对于 p 阶的常微分方程数值算法，当阶的常微分方程数值算法，当 h 0， n 时，是否时，是否 yn+1 y(xn+1)？p 阶算法的局部截断误差为：阶算法的局部截断误差为：111()()pnny xyO h 显然：显然：110,lim()0nnhny xy78局部截断误差的前提假设局部截断误差的前提假设是：是： ()nnyy x 局部截断误差局部截断误差 0 并不能保证算法收敛并不能保证算法收敛单步法的收敛性（续）单步法的收敛性（续）79单步法的收敛性

47、（续）单步法的收敛性（续）定义定义：若求解某初值问题的单步数值法，：若求解某初值问题的单步数值法，对于固定对于固定的的 当当 h 0 且且 n 时，它的时，它的近似解近似解趋向趋向于精确解于精确解 y(xn)，即：，即：0nxxnh0,()limnnhnyy x 则称该单步法是收敛的。则称该单步法是收敛的。80定义定义：称：称 y(xn) yn 为单步法的近似解为单步法的近似解 yn 的的整体整体截断误差截断误差。单步法收敛单步法收敛h 0 且且 n 时，时，yn 的整体截的整体截断误差断误差 0 单步法的收敛性（续）单步法的收敛性（续）81例：例：就初值问题就初值问题 考察考察欧拉显式格式的

48、收敛性。欧拉显式格式的收敛性。0)0(yyyy 解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 82欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 对任意固定的对任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 83单步法的收敛性（续）单步法的收敛性（续）收敛性定理收敛性定理若某单步法满足以上条件，则该方法是收敛的若某单步法满足以上条件，则该方法是收敛的则该单步法的整体截断误差为：则该单步法的整体截断误差为：若若单步法单步法 具有具有 p 阶阶精度，且增量函数精度，且增量函数 关于

49、关于 y 满足：满足：1(, )nnnnyyhxy h ( , , )x y h ()()pnny xyO h( , , )( , , )x y hx y hL yy Lipschitz 条件：条件： 初值初值 y0 是准确的是准确的84假设在前一步假设在前一步 yn 准确的前提下求得的近似值为：准确的前提下求得的近似值为：( , , )( , , )x y hx y hL yy 1(, )nnnnyyhxy h 1,()()nnnnyhxhy xy x 算法精度为算法精度为 p 阶，局部截断误差：阶，局部截断误差：111()pnny xych 111111() ()nnnnnny xyy x

50、yyy1111()nnnny xyyy1 (), (), )(, )pnnnnnnchy xyhxy xhhxy h 1(), (), )(, )pnnnnnnchy xyhxy xhxy h 1()()pnnnnchy xyhL y xy 1(1)()pnnhLy xych ne e1ne e 8511(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

51、(1)(1)(1)1nnphLhLhLche e 101(1)1(1)(1)1nnphLhLchhL e e 0(1)(1)1pnnchhLhLL e e860(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 h87根据这一定理，判

52、断单步法根据这一定理，判断单步法 的收敛性，归结为验证增量函数的收敛性，归结为验证增量函数 能否满足能否满足Lipschitz条件：条件：1(, )nnnnyyhxy h ( , , )( , , )x y hx y hL yy Lipschitz 条件：条件： 88100(,)0,1,2,()nnnnyyhf xynyy x Euler 方法：方法：收敛条件：收敛条件：89单步法的稳定性单步法的稳定性在讨论单步法收敛性时一般认为数值方在讨论单步法收敛性时一般认为数值方法本身的计算过程是准确的，实际上并法本身的计算过程是准确的，实际上并非如此：非如此：初始值初始值 y0 有误差有误差 d d

53、y0 y(x0)后续的每一步计算均有舍入误差后续的每一步计算均有舍入误差90这些初始和舍入误差在计算过程的传播这些初始和舍入误差在计算过程的传播中是逐步衰减的还是恶性增长就是数值中是逐步衰减的还是恶性增长就是数值方法的稳定性问题方法的稳定性问题单步法的稳定性单步法的稳定性91例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解上的解。分别。分别用欧拉法、隐式欧拉法和改进的欧拉用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法隐式欧拉法隐式欧拉法欧拉法欧拉法 节点节点 x

54、ixey30 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 792单步法的稳定性（续）单步法的稳定性（续）定义定义：设在节点：设在节点 xn 处用数值算法得到的处用数值算法得到的理想数值解为理想数值解为 y

55、n，而实际计算得到，而实际计算得到的近似解为的近似解为 ，称差值：，称差值：ny nnnyyd d 为第为第 n 步的数值解的步的数值解的扰动扰动。93定义定义：若一种数值方法在节点：若一种数值方法在节点 xn 处的数处的数值解值解 yn 的扰动的扰动 ，而在以后各节，而在以后各节点点 ym (m n) 上产生的扰动为上产生的扰动为 ，如果：如果：md d0nd d (1,2,)mnmnndddd则称该数值方法是则称该数值方法是稳定稳定的。的。94单步法的稳定性（续）单步法的稳定性（续）欧拉法：欧拉法：1(,)nnnnyyhf x y 由于函数由于函数 f (x, y) 的多样性，数值稳定性的

56、多样性，数值稳定性的分析相当复杂，通常只研究模型方程的分析相当复杂，通常只研究模型方程(0)yy 考察模型方程：考察模型方程：(0)yy 1(1)nnyhy 即：即：95假设在节点值假设在节点值 yn 上有扰动上有扰动 d dn，在节点值，在节点值 yn 1 上有扰动上有扰动 d dn 1，且，且 d dn 1 仅由仅由 d dn 引起引起（即：计算过程中不再引起新的误差）（即：计算过程中不再引起新的误差）1(1)nnyhy 即：即：96111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d 欧拉法稳定欧拉法稳定1nndddd 11h 即：即：111h 欧拉法稳定

57、的条件：欧拉法稳定的条件：20h 0 1(1)nnyhy 针对模型方程：针对模型方程：的显式欧拉法：的显式欧拉法：1(,)nnnnyyhf xy 化简得：化简得：20h 97隐式欧拉法：隐式欧拉法：111(,)nnnnyyhf xy考察模型方程：考察模型方程：(0)yy 即：即：11nnnyyh y 化简为：化简为：11nnyyh 98111111nnnnnnyyyyhhhd dd d 假设假设 yn 上有扰动上有扰动 ，则，则 yn 1 的的扰动为：扰动为：nnnyyd d 隐式欧拉隐式欧拉法稳定法稳定1nndddd 111h 0 0h ，上式均成立，所以：，上式均成立，所以：隐式欧拉法稳定

58、是恒稳定的隐式欧拉法稳定是恒稳定的一些单步法的绝对稳定区间见下表一些单步法的绝对稳定区间见下表方方 法法绝对稳定区间绝对稳定区间Euler方法方法改进的改进的Euler方法方法三阶三阶R-K法法四阶四阶R-K法法隐式隐式Euler法法梯形法梯形法- -2 z 0- -2 z 0- -2.51 z 0- -2.785 z 0- - z 0- - z 0100 对二阶对二阶R- -K方法方法，解模型方程，解模型方程(4.1)可得到可得到,2)(121nnyhhy 故故.2)(1)(2 hhhE 绝对稳定域绝对稳定域由由|E(h)|1得到，于是可得得到，于是可得绝对稳定区间绝对稳定区间为为- -2h

59、0，即，即0h2/.12121(,),22nnnnnnyyhkkf xyhhkfxyk 101 类似可得三阶及四阶类似可得三阶及四阶R- -K方法方法的的E(h)分别为分别为.!3)(!2)(1)(32 hhhhE .!4)(!3)(!2)(1)(432 hhhhhE 由由|1+h|1可得到相应的可得到相应的绝对稳定域绝对稳定域. 当当为实数时，为实数时，则得则得绝对稳定区间绝对稳定区间，它们分别为，它们分别为 三阶显式三阶显式R- -K方法方法： 四阶显式四阶显式R- -K方法方法：./51. 20, 051. 2 hh即即./78. 20, 078. 2 hh即即102从以上讨论可知显式从

60、以上讨论可知显式R- -K方法方法的的绝对绝对稳定域稳定域均为有限域，都对步长均为有限域，都对步长h有限制有限制. 如果如果h不在所给的不在所给的绝对稳定区间绝对稳定区间内，方法内，方法就不稳定就不稳定.103例例 分别取分别取h=1，2，4，用经典四阶，用经典四阶R-K方法求方法求解解, 1)0(yyy讨论稳定性。讨论稳定性。解：解：421, 1，分分别别为为则则hz 由上表知道，四阶由上表知道，四阶R-K法的绝对稳定区间是法的绝对稳定区间是 -2.785 z 0，即：当即：当h 2.785时，四阶时，四阶R-K方法才稳定。方法才稳定。故，当故，当h =1，2时，四阶时，四阶R-K方法的计算