常微分方程数值解法[高教课堂]

1、 第九章第九章 常微分方程的数值解法常微分方程的数值解法 1、引言引言 2、初值问题的数值解法、初值问题的数值解法-单步法单步法 3、龙格、龙格-库塔方法库塔方法 4、收敛性与稳定性收敛性与稳定性 5、初值问题的数值解法、初值问题的数值解法多步法多步法 6、方程组和刚性方程、方程组和刚性方程 7、习题和、习题和总结总结 主主 要要 内内 容容 1教育教学 主主 要要 内内 容容 研究的问题研究的问题 数值解法的意义数值解法的意义 1、 引引 言言 2教育教学 现实世界中大多数事物 内部联系非常复杂 找出其状态和状态变化规律之间的相互联系， 也即一个或一些函数与他们的导数之间的关系 此种关系的数

2、学表达就为 1.1.什么是微分方程什么是微分方程 ? ? 3教育教学 如何建立数学模型已在建模课程中得到讨论， 各类微分方程本身和他们的解所具有的特性 已在常微分方程及数学物理方程中得以解释， 4教育教学 寻找解析解的过程称为求解微分方程组。寻找解析解的过程称为求解微分方程组。 一个或一组具有所要求阶连续导数的解析函数，将 它代入微分方程（组），恰使其所有条件都得到满 足的解称为解析解解析解( (或古典解或古典解),),称为真解或解。称为真解或解。 5教育教学 4.4.什么是微分方程的数值解什么是微分方程的数值解? ? 虽然求解微分方程有许多解析方法,但解析方法 只能够求解一些特殊类型的方程,

3、从实际意义 上来讲 我们更关心的是某些 特定的自变量在某一个 定义范围内的一系列离散点一系列离散点上的近似值. 寻找数值解的过程称为数值求解微分方程。寻找数值解的过程称为数值求解微分方程。 数值解数值解 6教育教学 在大量的实际方程中出现的函数起码的连续性都 无法保证，更何况要求阶的导数 求解数值解 很多微分方程很多微分方程 根本求不到根本求不到 问题的解析解！问题的解析解！ 重要手段。 7教育教学 常微分方程的数值解法常用来求近似解 根据提供的算法 通过计算机通过计算机 便捷地实现便捷地实现 5.常微分方程数值解法的特点常微分方程数值解法的特点 8教育教学 00 ( , )(1.1) ( )

4、 (1.2) yf x yaxb y xy 本章主要讨论一阶常微方程的初值问题 6.基本知识基本知识 其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。初值条件。 9教育教学 则称 f (x,y) 对y 满足李普希兹条件，L 称为 Lipschitz常数. 常微分方程的理论指出：常微分方程的理论指出： 当 f (x,y) 定义在区域 G=(axb,y) 若存在正的常数 L 使： 1212 |( ,)( ,)| (1.3)f x yf x yL yy 12 , ,xa byy使得对任意的及都成立 就可保证方程解的存在唯一性就可保证方程解的存在唯一性 (Lipschitz)条件条

5、件 10教育教学 我们以下的讨论，都在满足上述条件下进行. 若 f (x,y) 在区域 G连续，关于y 一阶常微分方程的初值问题的解存在且唯一问题的解存在且唯一. 一阶常微分方程组常表述为： 111 1 ( ,) ( ) ( ,) m mmm yf x yy axb yfx yy 101 0 () () mm yx yx 11教育教学 写成向量形式为 (1)() 00000 ( , ) (),(,) mT yf x yaxb y xyxxx 高阶常微分方程定解问题如二阶定解问题： (,) () () yfxyyaxb ya y a 12教育教学 这些解法都可以写成向量形式 13教育教学 简单的

6、数值方法与基本概念简单的数值方法与基本概念 1. 简单简单欧拉法（欧拉法（Euler） 2后退的欧拉法后退的欧拉法 3梯形法梯形法 4改进改进EulerEuler法法 2、初值问题的数值解法、初值问题的数值解法单步法单步法 14教育教学 1. 简单的欧拉简单的欧拉(Euler)方法方法 考虑模型： 00 ( , ) (1.1) () (1.2) yf x yaxb y xy 在精度要求不高时 通过欧拉方法的讨论 欧拉方法欧拉方法 15教育教学 2. 欧拉方法的导出欧拉方法的导出 把区间a,b 分为n个小区间 步长为 要计算出解函数 y(x) 在一系列节点 () ii yy x )/()( ii

7、i xaihhhban,一般取即等距节点 处的近似值 01 n axxxb 1 (-) iii hxx N N等分等分 16教育教学 00 ( ,)(1.1) () (1.2) yfx yaxb y xy 对微分方程(1.1)两端从 1nn xx 到 进行积分 11 ( ,( ) nn nn xx xx y dxfx y xdx 1 1 ()()( , ( ) n n x nn x y xy xf x y x dx 17教育教学 右端积分用右端积分用 左矩形数值左矩形数值 求积公式求积公式 2 () ( )()()() 2 b a g g x dxba g aba ( )( , ( )g xf

8、 x y x令 1 ) 1 (, () (, () nnnn xx nn f xy x nn yyf xy x h 得 18教育教学 x0 x1 1 ()()()(,) nnnnnn y xy xhy xyh f xy )1,., 0(),( 1 niyxfhyy iiii 亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 每步计算只用到1n y 或用向前差 商近似导数 1 ()() () nn n y xy x yx h 1000 2111 1 (,) (,) (,) nnnn yyhf xy yyhf xy yyhf xy 依上述公式逐次计

9、算可得： n y 例题 19教育教学 3.欧拉公式有明显的几何意义 依此类推得到一折线 故也称Euler为单步法。n y 公式右端只含有已知项 所以又称为显格式的单步法。 0000 111 22 (,) ( )(,) ( ) (,) (,) xyy xxy y xxxxy xy 过 点 的 曲 线 是 解 在 作 的 切 线 与 直 线 交 于 再 作 切 线 交 于 20教育教学 也称欧拉折线法. 就是用这条折线近似地代替曲线 ( )y x ( )yy x x y n x 1n x n p 1n p 1n p x 21教育教学 从上述几何意义上得知，由Euler法所得的 折线明显偏离了积分曲

10、线，可见此方法 非常粗糙。非常粗糙。 4.欧拉法的局部截断误差： 在假设第 i 步计算是精确的前提下，考虑 截断误差 称为局部截断误差 /* local truncation error */。 1 () iii Ryxy 22教育教学 定义定义若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。 Ri 的的主项主项 /* leading term */ 欧拉法的局部截断误差： 2 3 2 ()() h i yxO h 欧拉法具有 1 阶精度。 ),()()()()()( 3 211 2 iiii h iiiii yxhfyhOxyxyhxyyxyR 定义定义在假设 yi = y(

11、xi)，即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。 23教育教学 如果单步差分公式的局部截断误差为O(hp+1), 则称该公式为p p阶方法阶方法.这里p为非负整数. 显然,阶数越高,方法的精度越高. Taylor展开式,一元函数的Taylor展开式为: 32 1 ! 3 )( ! 2 )( )()()()(h xy h xy hxyxyhxyxy nn nnnn 2 3 () ( )( )( )() ( ,) 112 h Ry xyy xhy xy xO hyhf x y

12、iiiiiiiii 2 3 2 ()() h i yxO h 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。 Ri 的的主项主项 /* leading term */ 24教育教学 5. 欧拉公式的改进欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */ 向后差商近似导数向后差商近似导数 h xyxy xy )()( )( 01 1 x0 x1 )(,()( 1101 xyxfhyxy )1,., 0(),( 111 niyxfhyy iiii 由于未知数由于未知数 yi+1 同时出现在等

13、式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故 称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。 一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差： 11) ( iii yxyR)()( 3 2 2 hOxy i h 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 25教育教学 6.梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均

14、平均 )1,., 0(),(),( 2 111 niyxfyxf h yy iiiiii 注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。 但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到 迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。 )()( 3 11 hOyxyR iii 中点欧拉公式中点欧拉公式 /* midpoint formula */ 中心差商近似导数中心差商近似导数 h xyxy xy 2 )()( )( 02 1 x0 x2x

15、1 )(,(2)()( 1102 xyxfhxyxy 1,., 1),(2 11 niyxfhyy iiii 假设假设 ，则可以导出，则可以导出 即中点公式具有即中点公式具有 2 阶精度。阶精度。 )(),( 11iiii xyyxyy )()( 3 11 hOyxyR iii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推 过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /* single-step method */。 26教育教学 方方 法法 显式欧拉显式欧拉 隐式

16、欧拉隐式欧拉 梯形公式梯形公式 中点公式中点公式 简单简单精度低精度低 稳定性最好稳定性最好精度低精度低, 计算量大计算量大 精度提高精度提高计算量大计算量大 精度提高精度提高, 显式显式 多一个初值多一个初值, 可能影响精度可能影响精度 27教育教学 改进欧拉法改进欧拉法 /* modified Eulers method */ Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),( 1iiii yxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到 1 i y ),(),( 2 111 iiiiii yxfyxf

17、 h yy 注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。 可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单 步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将 看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。 )1,., 0(),(,),( 2 11 niyxfhyxfyxf h yy iiiiiiii 28教育教学 3 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */ 考察改进的欧

18、拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为： ),( ),( 2 1 2 1 12 1 211 hKyhxfK yxfK KKhyy ii ii ii 斜率斜率 一定取一定取K1 K2 的的平均值平均值吗？吗？ 步长一定是一个步长一定是一个h 吗？吗？ 单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜 率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所 能达到的最高精度为能达到的最高精度为2阶阶。 建立高精度的单步递推格式。建立高精度的单步递推格式。 29教育教学 首先

19、希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算法格式有2阶阶 精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )( ii xyy )()( 3 11 hOyxyR iii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开 )(),(),(),( ),( 2 1 12 hOyxfphKyxphfyxf phKyphxfK iiyiixii ii )()()( 2 hOxyphxy ii 将改进欧拉法推广为：将改进欧拉法推广为： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii

20、 ii ii ),(),(),( ),(),( ),()( yxfyxfyxf dx dy yxfyxf yxf dx d xy yx yx Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()( )()()()( 32 221 2 211 hOxyphxyhy hOxyphxyxyhyy iii iiiii 30教育教学 Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较 )()()()( 32 2211 hOxyphxyhyy iiii )()( 2 )()()( 3 2 1 hOxy h xyhxyxy iii

21、i 要求要求 ，则必须有：，则必须有：)()( 3 11 hOyxyR iii 2 1 ,1 221 p 这里有这里有 个未知个未知 数，数， 个方程。个方程。 3 2 存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库 塔格式塔格式。 2 1 , 1 21 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？ 31教育教学 其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1

22、, , i 1 ) 均为待定均为待定 系数，确定这些系数的系数，确定这些系数的 步骤与前面相似。步骤与前面相似。 ).,( . ),( ),( ),( . 112211 23213133 12122 1 22111 mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ： ),( ),( ),( ),( )22( 34 2223 1222 1 432161 hKyhxfK KyxfK K

23、yxfK yxfK KKKKyy ii h i h i h i h i ii h ii 32教育教学 注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的 值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精 度阶数的关系：度阶数的关系： 753 可达到的最高精度可达到的最高精度 642 每步须算每步须算Ki 的个数的个数 )( 2 hO)( 3 hO)( 4 hO)( 5 hO)( 6 hO)( 4 hO)( 2n hO 8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的

24、导出基于泰勒展开，故精度主要受 解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好 采用采用低阶算法低阶算法而将步长而将步长h 取小取小。 深入研究龙格深入研究龙格- 库塔法请看库塔法请看此处此处！ 33教育教学 4 收敛性与稳定性收敛性与稳定性 /* Convergency and Stability */ 1.收敛性收敛性 /* Convergency */ 定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛收

25、敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。 0 )0(yy yy 解：解：该问题的精确解为该问题的精确解为 x eyxy 0 )( 欧拉公式为欧拉公式为 iiii yhyhyy)1 ( 1 0 )1 (yhy i i 对任意固定的对任意固定的 x = xi = i h ，有，有 ii xhhx i hyhyy )1()1( /1 0 / 0 eh h h /1 0 )1(lim )( 0i x xyey i 34教育教学 2.稳定性稳定性 /* Stability */ 例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解

26、。 分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0( )(30)( y xyxy 0.0 0.1 0.2 0.3 0.4 0.5 精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xi x ey 30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.0000 2.5000 10 1 6.2500 10 2 1.5625 10 2 3.9063 10 3 9.7656 10 4 1.0000 2.5000 6.2500 1.5626 101 3.906

27、3 101 9.7656 101 1.0000 4.9787 10 2 2.4788 10 3 1.2341 10 4 6.1442 10 6 3.0590 10 7 What is wrong ?! 35教育教学 定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计 算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。 一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */ yy 常数，可以常数，可以 是复数是复数

28、当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值 产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该 算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。 我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。 000 yy h h h 36教育教学 例：例：考察显式欧拉法考察显式欧拉法 0 1 1 )1(yhyhyy i iii 000 yy 0 1 1 )1(yhy i i 0 1 111 )1( i ii

29、i hyy 由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减， 必须满足：必须满足： hh 1|1| h 0-1-2 Re Img 例：例：考察隐式欧拉法考察隐式欧拉法 11 iii yhyy ii y h y 1 1 10 1 1 1 1 i i h 可见绝对稳定区域为：可见绝对稳定区域为：1|1| h 210Re Img 注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式一般来说，隐式欧拉法的绝对稳定性比同阶的显式 法的好。法的好。 37教育教学 例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1( ).,( . 11 111 mj hKhKyhxfK

30、KKhyy mmjjijij mmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定 区域为区域为 ) 2 , 2 ( 11 11 K h y h xfK hKyy ii ii 其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为 0 Re Img k=1 k=2 k=3 k=4 -1-2-3 - - - 1 2 3 Re Img 无条件稳定无条件稳定 38教育教学 例例1 用欧拉方法用欧拉方法，隐式欧拉方法隐式欧拉方法和和欧拉中点公式欧拉中点公式计算计算 01 (0 )1 yyx y 的近似解，取步长的近似解，取步长h=0.1，并与精确值比较，并与精确值比较 解解:欧拉方法的算

31、式为：欧拉方法的算式为： 1 0 (1)0,1, 9 . 1 ii yyhi y 欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为： 1 0 0,1, 9 . 1 1 i i y yi h y 欧拉中点公式是两步法，第一步欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式，以后用公式 11 0 21, 2, 9 . 1 iii yyh yi y 39教育教学 本题精确解为本题精确解为y=e-x,计算结果见表计算结果见表9-1 例例2 用用欧拉公式欧拉公式和和梯形公式梯形公式的预估校正法计算：的预估校正法计算： 的数值解，取的数值解，取

32、h=0.1，梯形公式只迭代一次，并与精确值比较，梯形公式只迭代一次，并与精确值比较. 方程的解析解为方程的解析解为: 2 01 (0 )1 x yyx y y xy21 解解: 本例中欧拉公式为：本例中欧拉公式为： 1 0 20.2 0.11.10,1,9. 1 ii iiii ii xx yyyyi yy y 40教育教学 梯形公式只校正一次的格式为梯形公式只校正一次的格式为 (0) 1 (0) 1 11 (0) 1 0 0.2 1.1 22 0.05 1 i ii i ii iiii ii x yy y xx yyyy yy y 结果列入表结果列入表9.2 41教育教学 1. Runge-

33、Kutta 1. Runge-Kutta 法的一般形式法的一般形式 2. 22. 2阶阶Runge-Kutta Runge-Kutta 方法方法 3. 3. 经典经典Runge-Kutta Runge-Kutta 方法方法 4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法 5. 5. 隐式隐式R-KR-K方法方法 6. 6. 变步长方法变步长方法 42教育教学 1.Runge-1.Runge-Kutta Kutta 法的一般形式法的一般形式 Runge-Kutta 法是用区间上若干个点上的导数 的线性组合得到平均斜率，以提高方法的阶。 其一般形式为 ： 1 1 L nni

34、i i yyhk 1 2221 333311322 1 1 ,1,2, (,) (,) (,() nn nn nn i ininii jj j kfxc h yc ha kiL kf x y kf xc h yc hk kf xc h yc h a ka k 43教育教学 式(9.11) 称L级级p阶阶Runge-Kutta方法方法(简称R-K法法)。 当L1就是欧拉法，当L2时为改进的欧拉法。 1 11 1 , 1 , 1 Li iii j ij ca 。 11 1 ()() L nnnii i Ty xy xhk 其中 它的局部截断误差是 (9.11) 2. 22. 2级级2 2阶阶Run

35、ge-KuttaRunge-Kutta方法方法 令令 L=2L=2，则，则 11 122 () nn yyhkk 44教育教学 1 2221 (,) (,) nn nn kfxy kfxc h yc hk 111 122 ()()() nnn Ty xy xhkk 其局部截断误差是其局部截断误差是： 这是3个未知量的两个方程，其中有一个自由参数， 方程组有无穷多解，从而得到一族2级2阶R-KR-K方法 。 12 22 1 1 2 c 45教育教学 3. 3. 经典经典Runge-KuttaRunge-Kutta方法方法 我们可以构造出一族3级3阶，一族4级4阶和一族5级4 阶等R-K方法。最常

36、用的4级4阶是如下的经典经典R-KR-K方法方法： 11234 1 21 32 43 (22) 6 (,) (,) 22 (,) 22 (,) nn nn nn nn nn h yykkkk kfxy hh kfxyk hh kfxyk kfxh yhk 46教育教学 4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法 R-K-FehhlbergR-K-Fehhlberg方法方法是在R-K方法的基础上引进误差和步长 控制的办法。即利用5阶R-K法 112456 1666562856192 13512825564305055 ii yyKKKKK 估计4阶R-K的局部误差，其

37、中 54311 5 1 4104 2197 2565 1048 216 25 KKKKyy ii ) 32 9 32 3 , 8 3 ( ) 4 1 , 4 ( ),( 213 12 1 KKyhxhfK Ky h xhfK yxhfK ii ii ii 47教育教学 ) 40 11 4104 1859 2565 5344 2 27 8 , 2 ( ) 4104 854 513 3680 8 216 439 ,( ) 2197 7296 2197 7200 2197 1932 , 13 12 ( 543216 43215 3214 KKKKKy h xhfK KKKKyhxhfK KKKyhx

38、hfK ii ii ii 注：注：R-K-Fehhlberg比4阶R-K方法具有更大的优 越性，他是计算稳定，高精度的方法，他的显著 优点是，每一步仅需计算f的6个值；若用4阶R- K-L与5阶R-K-L在一起使用，则每步需要计算f的 10个值。推荐使用！ 48教育教学 5. 5. 隐式隐式R-KR-K方法方法 类似于显式R-K公式，稍加改变，就得到隐式隐式R-K方法方法。 1 1 L nnii i yyhk 1 , L ininii jj j kfxc h yc ha k 它与显式R-K公式的区别在于：显式公式中对系数求 和的上限是i-1，从而构成的矩阵是一个严格下三角阵。 而在隐式公式中对

39、系数求和的上限是L，从而构成的 矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式 的思路和方法与显式R-K法类似。通常，同级的隐式 公式获得比显式公式更高的阶。 49教育教学 通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式 R-KR-K法有法有: 2 , 2 1 1 1 nn nnn yy hxhfyy ),(),( 2 111 nnnnnn yxfyxf h yy 112 112 212 () 2 331323 , 6412 333231 , 6124 nn nn nn h yykk kfxh yk hk h kfxh yk hk

40、 h 1级级2阶中点公式阶中点公式 : 2级级2阶梯形公式：阶梯形公式： 2级级4阶阶R-K公式：公式： 50教育教学 6.6.变步长方法变步长方法 在单步法中每一积分步步长实际上是相互独立的，步长的 选择具有了灵活性。根据合理地选择每一积分步的步长， 既保证精度的要求，又可以减少计算量，从而减小舍入误 差。其方便的控制手段是基于误差的事后估计式。 对于给定的精度 ,如果 ,反复将步长折半进行计算， 直至 为止,这时取最终得到的作为结果； 如果 为止，这时再 将步长折半一次，就得到所要的结果。 这种通过加倍或折半处理步长的计算方法称为变步长方法。 注注：推荐使用精度好计算量低的变步长方法。 用

41、四阶显式R-K方法做变步长方法是实践中较好的方法！ 51教育教学 例 分别用改进的欧拉格式和四阶龙 格库塔格式解初值问题（取步长 h=0.2）： 2 (0)1 x yy y y (01)x 52教育教学 节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.183216 0.4 1.348312 1.341667 1.341641 0.6 1.493704 1.483281 1.483240 0.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051 表 （注：已指出过准确解12yx

42、 ） 53教育教学 单步法的相容性、收敛性和稳定性单步法的相容性、收敛性和稳定性 1.1.单步法的相容性单步法的相容性 2.2.单步法的收敛性单步法的收敛性 3.3.单步法的稳定性单步法的稳定性 4.4.相容性和收敛性的关系相容性和收敛性的关系 5.5.相容性和方法阶的关系相容性和方法阶的关系 6.6.稳定性和收敛性稳定性和收敛性 7.7.绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 54教育教学 定义一定义一：对于（：对于（9.1.1）常微分方程初值问题）常微分方程初值问题 单步法的形式可以变表示为单步法的形式可以变表示为 （9.2.19） 其中其中 h为步长为步长 若对求解区间中任一固定的

43、若对求解区间中任一固定的x，当，当 时皆成立时皆成立 1 00 ( ,; ) nnnn yyh x y h yy 0h 00 ()() limlim( ,;) hh y xhy x x y h h 则称由（9.2.19）确定的单步法与微分方程初值问题是相容相容的 注注意到上式左边极限为由（由（9.1.19.1.1）知它应等于从而由相容性 定义得 称相容性条件相容性条件。 ( , ;0)( , )x yf x y 55教育教学 单步法的收敛性单步法的收敛性 定义二：定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19） 产生的数值解，对于每一个固定的 , ， 当 即 。若成立

44、， 则称该方法是收敛的。 n y n x 0h 1nn xx 1 () nn yy x 1nn xxh 56教育教学 定义三定义三： 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （mn）处的值产生的偏差均不超 过 ，则称该方法是的。 单步法的稳定性有以下定理 n y m y 定理二定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是的。 57教育教学 相容性和收敛性的关系相容性和收敛性的关系 若单步法的增量函数对变量y满足Lipschtiz 条件， 即存在与 h , x 无关的常数 L， 对区域D= 任意两点 （x，y1）,(x,y2)成立，则单步法收

45、敛的充分必要 条件是相容性条件成立。（读者自证） ,axby 58教育教学 相容性和方法阶的关系相容性和方法阶的关系 若单步法是若单步法是p p阶方法则成立阶方法则成立 若单步法满足相容性条件，得若单步法满足相容性条件，得 所以所以 =0=0 也就是说单步法的阶数一定要是正数。由于我们考虑也就是说单步法的阶数一定要是正数。由于我们考虑 的单步法皆为正整数，的单步法皆为正整数，p p至少为一。因此我们考虑的至少为一。因此我们考虑的 单步法都满足相容性条件。单步法都满足相容性条件。 1 ()( )( , ; )() p y x hy xhx y hO h 1 0 1 ( )( , )() lim

46、p h y xf x yO h h 0 () lim p h O h 59教育教学 稳定性和收敛性的关系稳定性和收敛性的关系 若单步法的增量函数满足定理二的条件即单若单步法的增量函数满足定理二的条件即单 步法是稳定的则单步法收敛的充分必要条件步法是稳定的则单步法收敛的充分必要条件 是是 相容性条件成立。相容性条件成立。 ( , ;0)( , )x yf x y 60教育教学 绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 稳定性问题是一个比较复杂的问题。为了简化讨论一稳定性问题是一个比较复杂的问题。为了简化讨论一 般仅对试验方程般仅对试验方程 进行考察。这里假定进行考察。这里假定 Re0,Re0

47、h0和的允许范围，称为该方法的和的允许范围，称为该方法的 绝对稳定域绝对稳定域。 yy n n k y 61教育教学 绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域2 2 将Euler方法应用到试验方程得 误差方程是 要求误差不增长则 nnnn yhhyyy)1 ( 1 nnn h 1 1|1| 1 h n n 62教育教学 则Euler 方法的绝对稳定域是以 为半径的圆的内部。 为了保证稳定性步长有所控制。 假如当 时h应满足 ，当 时， h 应满 足 等等。 同样我们可以将试验方程用到其它各种单步法当中，求出其 绝对稳定域。在实际应用中必须在这个范围内，否则误差传 播相当大，即不稳定。 1h

48、 5 5 2 0 h 100 50 1 100 2 0 h 63教育教学 1.Adams1.Adams方法方法 2.2.米尔尼方法、汉明方法及辛普森方法米尔尼方法、汉明方法及辛普森方法 3.3.预测校正方法预测校正方法 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性 5 初值问题的数值解法初值问题的数值解法 多步法多步法 64教育教学 考虑型如考虑型如 的的k步法，步法， 称为称为阿当姆斯阿当姆斯（Adams）方法）方法 1 0 k nknkini i yyhf 11233 0,0,0 0 1012 112 112 112 1 1 2241 33121 44321 1.A

49、dams1.Adams方法方法 拿其中一种为例，推导其公式。对上面所说的法一般形拿其中一种为例，推导其公式。对上面所说的法一般形 式式若取若取 ，有，有方程组方程组 同样解之，得到同样解之，得到3步步4阶隐式阶隐式Adams公式和它的余项。公式和它的余项。 65教育教学 1112 (9195) 24 nnnnnn h yyffff 5(5) 121 19 (),(,) 720 nnnnn Th yxx 其它请读者自证。我们仅将结论列于下表。其它请读者自证。我们仅将结论列于下表。 Adams显式公式显式公式 1p c 1nnn yyhf 1 2 211 ( 3) 2 nnnn h yyff 5

50、12 3221 (23165) 12 nnnnn h yyfff 3 8 43321 (5559379) 24 nnnnnn h yyffff 251 720 kp公式 11 2 2 33 44 66教育教学 Adams隐式公式隐式公式 1p c 11 () 2 nnnn h yyff 1 12 2121 (58) 12 nnnnn h yyfff 1 24 32321 (9195) 24 nnnnnn h yyffff 19 720 434321 (25164626410619 ) 720 nnnnnnn h yyfffff 3 160 k p公式 12 23 34 45 注：注：在这些方法

51、中，在这些方法中，4阶的阶的Adams预测校正方法具有方法预测校正方法具有方法 简洁、使用方便、易排序、高精度等优点。尤其当函数简洁、使用方便、易排序、高精度等优点。尤其当函数 f比较复杂时更能显出它的优越性。比较复杂时更能显出它的优越性。 67教育教学 同 理 得 到同 理 得 到 5 个 待 定 参 数 方 程 组 。 解 之个 待 定 参 数 方 程 组 。 解 之 得得 ， ， ， ， 。 构成著名的构成著名的Miline 4Miline 4步步4 4阶显式公式和它的余阶显式公式和它的余 项。项。 1 3 3 8 0 3 4 1 ， 3 8 2 0 3 13nn yy 12 4 (22

52、) 3 nnn h fff 5(5) 131 14 ( ),(,) 45 nnnnn Th yxx 2.2.米尔尼方法、汉明方法及米尔尼方法、汉明方法及 辛普森方法辛普森方法 68教育教学 同理得到同理得到5 5个参数方程组。求解后就构成著名的个参数方程组。求解后就构成著名的 3 3步步4 4阶隐式阶隐式HammingHamming公式和它的余项。公式和它的余项。 若取，也得到若取，也得到5个参数方程组。求解后就构成个参数方程组。求解后就构成 Simpson隐式公式和其余项。隐式公式和其余项。 12 1 (9) 8 nnn yyy 11 3 (2) 8 nnn h fff 5(5) 121 1

53、 (),(,) 40 nnnnn Th yxx 米尔尼方法、汉明方法及米尔尼方法、汉明方法及 辛普森方法辛普森方法 1111 (4) 3 nnnnn h yyfff 5(5) 111 1 (),(,) 90 nnnnn Th yxx 69教育教学 不论单步法或多步法，隐式公式比显式公式稳定 性好，但在实际使用隐式公式时，都会遇到两个 问题：一个是隐式公式如何能方便地进行计算； 另一个是实际计算步长取多大。如隐式梯形公式， 每往前推进一步，不必进行多次迭代，而是采用 一阶显式Euler公式预测，二阶隐式梯形公式校正 一次，构成显式改进Euler公式，能达到与梯形公 式同阶的精度，即二阶精度。 3

54、. 预测校正方法预测校正方法 70教育教学 对于线性多步公式，一般地，不难验证，如果 预测公式是阶或阶精度，校正公式为阶精度， 则用预测公式提供初值，校正公式迭代一次的 效果也能达到阶精度，再迭代下去，效果就不 明显了。预测-校正技术即保证了计算精度，又 使隐式计算显式化，克服了隐式公式要反复迭 代的困难。至于实际计算步长的选取，我们对 预测-校正公式使用外推原理，得到误差估计式 ，用来调整计算步长，使达到控制误差要求。 对于线性多步方法常用的预测-校正公式有Miline- Hamming方法和4阶阶Adams显隐式显隐式预测-校正公式 71教育教学 将Miline Miline 公式公式和H

55、ammingHamming公式公式结合，构成预测-校正公式 Miline公式公式和Hamming公式公式的局部截断误差分别为 1312 4 (22) 3 p nnnnn h yyfff 12111 13 (9)( (,) 2) 88 p nnnnnnn h yyyf xyff 5(5)6 111 14 ()()() 45 MM nnnn Ty xyh yxO h 5(5)6 111 1 ()( )() 40 HH nnnn Ty xyh yxO h 修正修正Miline-Hamming公式公式 72教育教学 利用外推原理，即加上后消去局部截断误差的主项。使 说明经过外推后的算法其精度提高一阶。

56、忽略误差项，上 式可改写为 11 6 1 114 4045 ()() 114 4045 MH nn n yy y xO h 111 114114 () ()0 40454045 MH nnn y xyy 111 91129112 () ()0 360360360360 MH nnn y xyy 73教育教学 由于 和是 在计算过程中获得的数据，称为 Miline Miline 公式公式和HammingHamming公式公式的事后误差估计式。我们可 以用它们来调节计算步长的大小，即选择一个合适的的 步长，使 ， 其中的是要求达到的计算精度。 111 121 ()91120 MH nnn y xy

57、y 1111 112 ()() 121 MHM nnnn y xyyy 1111 9 ()() 121 HHM nnnn y xyyy 1 M n y 1 H n y 11 9 () 121 HM nn yy 74教育教学 又可得到 MilineMiline公式公式和HammingHamming公式公式的修正公 式，它们分别是 从而构成如下的修正从而构成如下的修正HammingHamming预测预测- -校正公式，校正公式， 简称修正简称修正HammingHamming公式：公式： 1111 112 () 121 MmMHM nnnn yyyy 1111 9 () 121 HmHHM nnnn

58、 yyyy 75教育教学 预测： 修正： 校正： 修正： 在应用这套公式时，先由同阶单步法提供初 值 ， ， ， 。计算 时可取。 31 n p n yy)22( 3 4 21 nnn fff h )( 121 112 11 p n c n p n pm n yyyy )2),( 8 3 )9 ( 8 1 11121 nn pm nnnn c n ffyxf h yyy )( 121 9 1111 p n c n c nn yyyy 0 y 1 y 2 y 3 y4 y pc yy 33 76教育教学 将将4步步4阶显式阶显式Adams公式作为预测公式和公式作为预测公式和3步步4阶阶 式式Ad

59、ams公式作为校正公式，构成公式作为校正公式，构成Adams预测预测- 校正公式。校正公式。 它们的局部截断误差分别是它们的局部截断误差分别是 1123 (5559379) 24 p nnnnnn h yyffff 1112 (9195) 24 c nnnnnn h yyffff 5(5)6 111 251 ()()() 720 pp nnnn Ty xyh yxO h 5(5)6 111 19 ()( )() 720 cc nnnn Ty xyh yxO h Adams预测预测-校正公式校正公式 77教育教学 利用外推原理，将上两式作线性组合消去局部截 断误差的主项。使计算精度至少提高一阶，

60、同 时得到两个修正公式，将它们和上两式构成了 如下修正修正Adams公式公式： 预测： 修正： 校正: 1123 (5559379) 24 p nnnnnn h yyffff 11 251 () 270 pmpcp nnnn yyyy 11112 (9 (,) 195) 24 cpm nnnnnnn h yyf xyfff 78教育教学 修正： 同理，在计算时，调节计算步长 使 , 由同阶单步法提供初值 , ， ， 。计算 时，可 取 。 理论分析得出它们的绝对稳定区域，修正Hamming法： ; 4阶Adams预测预测-校正法校正法： 其中 )( 270 19 1111 p n c n c