第三章 常微分方程的差分法

1、 第三章第三章 常微分方程的差分法常微分方程的差分法 1、引言引言 2、初值问题的数值解法、初值问题的数值解法-单步法单步法 3、龙格、龙格-库塔方法库塔方法 4、收敛性与稳定性收敛性与稳定性 5、方程组和高阶方程的情形方程组和高阶方程的情形 6、边值、边值问题问题主主 要要 内内 容容主主 要要 内内 容容研究的问题研究的问题数值解法的意义数值解法的意义1、 引引 言言现实世界中大多数事物内部联系非常复杂找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系此种关系的数学表达就为1.1.什么是微分方程什么是微分方程 ? ?如何建立数学模型已在建模课程中得到讨论，各类

2、微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，寻找解析解的过程称为求解微分方程组。寻找解析解的过程称为求解微分方程组。 一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解解析解( (或古典解或古典解),),称为真解或解。称为真解或解。4.4.什么是微分方程的数值解什么是微分方程的数值解? ?虽然求解微分方程有许多解析方法,但解析方法只能够求解一些特殊类型的方程,从实际意义上来讲我们更关心的是某些 特定的自变量在某一个定义范围内的一系列离散点一系列离散点上的近似值.寻找数值解的过程称为数值求解微分方程。寻找数值解的

3、过程称为数值求解微分方程。数值解数值解在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数求解数值解很多微分方程很多微分方程根本求不到根本求不到问题的解析解！问题的解析解！重要手段。常微分方程的数值解法常用来求近似解根据提供的算法通过计算机通过计算机便捷地实现便捷地实现5.常微分方程数值解法的特点常微分方程数值解法的特点00( , )(1.1)( ) (1.2)yf x yaxby xy本章主要讨论一阶常微方程的初值问题6.基本知识基本知识其中f (x,y)是已知函数，(1.2)是定解条件也称为初值条件。初值条件。则称 f (x,y) 对y 满足李普希兹条件，L 称为Lips

4、chitz常数.常微分方程的理论指出：常微分方程的理论指出：当 f (x,y) 定义在区域 G=(axb,y)若存在正的常数 L 使：1212|( ,)( ,)| (1.3)f x yf x yL yy12 , ,xa byy使得对任意的及都成立就可保证方程解的存在唯一性就可保证方程解的存在唯一性(Lipschitz)条件条件我们以下的讨论，都在满足上述条件下进行.若 f (x,y) 在区域 G连续，关于y一阶常微分方程的初值问题的解存在且唯一问题的解存在且唯一.一阶常微分方程组常表述为：1111( ,) ( ) ( ,)mmmmyf x yyaxbyfx yy 1010() ()mmyxyx

5、写成向量形式为(1)()00000( , )(),(,)mTyf x yaxby xyxxx高阶常微分方程定解问题如二阶定解问题：(,) ()()yfxyyaxbyay a这些解法都可以写成向量形式简单的数值方法与基本概念简单的数值方法与基本概念 1. 简单简单欧拉法（欧拉法（Euler） 2后退的欧拉法后退的欧拉法 3梯形法梯形法 4改进改进EulerEuler法法 2、初值问题的数值解法、初值问题的数值解法单步法单步法1. 简单的欧拉简单的欧拉(Euler)方法方法考虑模型：00( , ) (1.1)() (1.2)yf x yaxby xy 在精度要求不高时通过欧拉方法的讨论欧拉方法欧拉

6、方法2. 欧拉方法的导出欧拉方法的导出把区间a,b分为n个小区间步长为要计算出解函数 y(x) 在一系列节点()iiyy x)/()( iiixaihhhban,一般取即等距节点处的近似值01 naxxxb1(-)iiihxxN N等分等分00( ,)(1.1)() (1.2)yfx yaxby xy 对微分方程(1.1)两端从1nnxx到进行积分11( ,( )nnnnxxxxy dxfx y xdx11()()( , ( )nnxnnxy xy xf x y x dx右端积分用右端积分用左矩形数值左矩形数值求积公式求积公式2()( )()()()2bagg x dxba g aba( )(

7、 , ( )g xf x y x令1)1(, ()(, () nnnnxxnnf xy xnnyyf xy xh得x0 x11()()()(,)nnnnnny xy xhy xyh f xy)1,., 0(),(1 niyxfhyyiiii亦称为亦称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 每步计算只用到1ny或用向前差商近似导数1()()()nnny xy xyxh100021111(,)(,) (,)nnnnyyhf xyyyhf xyyyhf xy依上述公式逐次计算可得：ny例题3.欧拉公式有明显的几何意义依此类推得到一折线故也称Euler

8、为单步法。ny公式右端只含有已知项所以又称为显格式的单步法。000011122(,) ( )(,)( ) (,)(,) xyy xxyy xxxxyxy 过 点 的 曲 线 是 解 在 作 的 切 线 与 直 线 交 于 再 作 切 线 交 于 也称欧拉折线法.就是用这条折线近似地代替曲线( )y x( )yy xxynx1nxnp1np1np x从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。非常粗糙。4.欧拉法的局部截断误差：在假设第 i 步计算是精确的前提下，考虑截断误差称为局部截断误差/* local truncation error */。1()

9、iiiRyxy定义定义若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。Ri 的的主项主项/* leading term */ 欧拉法的局部截断误差：232()()hiyxO h欧拉法具有 1 阶精度。),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR 定义定义在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p p阶方法阶方法.这里p为非负

10、整数.显然,阶数越高,方法的精度越高.Taylor展开式,一元函数的Taylor展开式为: 321! 3)(! 2)()()()()(hxyhxyhxyxyhxyxynnnnnn23() ( )( )( )() ( ,)112hRy xyy xhy xy xO hyhf x yiiiiiiiii232()()hiyxO h若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有，则称该算法有p 阶精度。阶精度。Ri 的的主项主项/* leading term */5. 欧拉公式的改进欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method *

11、/向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxy

12、ih 即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。6.梯形公式梯形公式 /* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii 中点欧拉公式中点欧拉公

13、式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2个初值个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是，而前面的三种算法都是单步法单步法 /

14、* single-step method */。方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(

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

16、，可以将其改写为：考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度的单步递推格式。首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格

17、式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx St

18、ep 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开作比较展开作比较)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii21,1221 p 这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有

19、满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。21, 121 p注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定均为待定系数，确定这些系数的系数，确定这些系数的步骤与前面相似。步骤与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhK

20、hKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的

21、最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！4 收敛性与稳定性收敛性与稳定性 /* Convergency and Stability */ 1.收敛性收敛性 /* Convergency */定义定义 若某算法对于任意固

22、定的若某算法对于任意固定的 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 ehhh /10)1(lim)(0ixxyeyi 2.稳定性稳定性 /* S

23、tability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 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.5

24、0006.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 ?!定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */yy 常数，可以常数，可以是复数

25、是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。000yy h h h例：例：考察显式欧拉法考察显式欧拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可见，要保证初始误差由此可见，要保证初始

26、误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh 1|1| h0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法11 iiiyhyy iiyhy 11101111 iih可见绝对稳定区域为：可见绝对稳定区域为：1|1| h210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。法的好。例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111Khyhxf

27、KhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定例例1 用欧拉方法用欧拉方法，隐式欧拉方法隐式欧拉方法和和欧拉中点公式欧拉中点公式计算计算01(0 )1yyxy 的近似解，取步长的近似解，取步长h=0.1，并与精确值比较，并与精确值比较解解:欧拉方法的算式为：欧拉方法的算式为：10(1)0,1, 9 .1iiyyhiy欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：100,1, 9 .11iiyyihy欧拉中点公式是两步法，第一步欧拉中

28、点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式，以后用公式11021, 2, 9 .1iiiyyh yiy本题精确解为本题精确解为y=e-x,计算结果略。计算结果略。例例2 用用欧拉公式欧拉公式和和梯形公式梯形公式的预估校正法计算：的预估校正法计算：的数值解，取的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较，梯形公式只迭代一次，并与精确值比较.方程的解析解为方程的解析解为:201(0 )1xyyxyy xy21 解解: 本例中欧拉公式为：本例中欧拉公式为：1020.20.11.10,1,9.1iiiiiiiixxyyyyiyyy梯形公式只校正一次的格式为梯形公式只校正一

29、次的格式为(0)1(0)111(0)100.21.1220.051iiiiiiiiiiiixyyyxxyyyyyyy结果结果略。结果结果略。1. Runge-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. 变步长方法变步长方法 1.Runge-1.Runge-Kutta Kutta 法的一般形式法

30、的一般形式Runge-Kutta 法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为 ：11Lnniiiyyhk1222133331132211,1,2,(,)(,)(,()nnnnnniininii jjjkfxc h yc ha kiLkf x ykf xc h yc hkkf xc h yc h a ka k式(9.11) 称L级级p阶阶Runge-Kutta方法方法(简称R-K法法)。当L1就是欧拉法，当L2时为改进的欧拉法。 1111 , 1 , 1 Liiii jijca。111()()LnnniiiTy xy xhk其中它的局部截断误差是(9.11)

31、2. 22. 2级级2 2阶阶Runge-KuttaRunge-Kutta方法方法令令 L=2L=2，则，则 11 122()nnyyhkk12221(,)(,)nnnnkfxykfxc h yc hk111 122()()()nnnTy xy xhkk其局部截断误差是其局部截断误差是： 这是3个未知量的两个方程，其中有一个自由参数，方程组有无穷多解，从而得到一族2级2阶R-KR-K方法 。1222112c3. 3. 经典经典Runge-KuttaRunge-Kutta方法方法我们可以构造出一族3级3阶，一族4级4阶和一族5级4阶等R-K方法。最常用的4级4阶是如下的经典经典R-KR-K方法方

32、法：112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyykkkkkfxyhhkfxykhhkfxykkfxh yhk4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法R-K-FehhlbergR-K-Fehhlberg方法方法是在R-K方法的基础上引进误差和步长控制的办法。即利用5阶R-K法 112456166656285619213512825564305055iiyyKKKKK估计4阶R-K的局部误差，其中 5431151410421972565104821625KKKKyyii)329323,83()41,4(),(213121

33、KKyhxhfKKyhxhfKyxhfKiiiiii)401141041859256553442278,2()410485451336808216439,()219772962197720021971932,1312(543216432153214KKKKKyhxhfKKKKKyhxhfKKKKyhxhfKiiiiii注：注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！ 5. 5. 隐式隐式R-KR-K方法方法类似于显式R-K

34、公式，稍加改变，就得到隐式隐式R-K方法方法。11Lnniiiyyhk1,Lininii jjjkfxc h yc ha k它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-KR-K法有法有:2,2111nnnnnyyhxhfyy),(),(2111nnnnnnyxf

35、yxfhyy112112212()2331323,6412333231,6124nnnnnnhyykkkfxh yk hk hkfxh yk hk h1级级2阶中点公式阶中点公式 :2级级2阶梯形公式：阶梯形公式： 2级级4阶阶R-K公式：公式： 6.6.变步长方法变步长方法在单步法中每一积分步步长实际上是相互独立的，步长的选择具有了灵活性。根据合理地选择每一积分步的步长，既保证精度的要求，又可以减少计算量，从而减小舍入误差。其方便的控制手段是基于误差的事后估计式。对于给定的精度 ,如果 ,反复将步长折半进行计算，直至 为止,这时取最终得到的作为结果；如果 为止，这时再将步长折半一次，就得到所

36、要的结果。这种通过加倍或折半处理步长的计算方法称为变步长方法。 注注：推荐使用精度好计算量低的变步长方法。用四阶显式R-K方法做变步长方法是实践中较好的方法！ 例 分别用改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：2(0)1xyyyy (01)x节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410.6 1.493704 1.483281 1.4832400.8 1.627861 1.612514 1.612452 1 1.754205 1.7321

37、42 1.732051表（注：已指出过准确解12yx）单步法的相容性、收敛性和稳定性单步法的相容性、收敛性和稳定性 1.1.单步法的相容性单步法的相容性 2.2.单步法的收敛性单步法的收敛性 3.3.单步法的稳定性单步法的稳定性 4.4.相容性和收敛性的关系相容性和收敛性的关系 5.5.相容性和方法阶的关系相容性和方法阶的关系 6.6.稳定性和收敛性稳定性和收敛性 7.7.绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 定义一定义一：对于常微分方程初值问题：对于常微分方程初值问题单步法的形式可以变表示为单步法的形式可以变表示为 其中其中 h为步长为步长若对求解区间中任一固定的若对求解区间中任一固

38、定的x，当，当 时皆成立时皆成立 100( ,; )nnnnyyh x y hyy0h00()()limlim( ,;)hhy xhy xx y hh则称由上式确定的单步法与微分方程初值问题是相容相容的注注意到上式左边极限为由原方程由原方程知它应等于从而由相容性定义得称相容性条件相容性条件。 ( , ;0)( , )x yf x y单步法的收敛性单步法的收敛性定义二：定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19）产生的数值解，对于每一个固定的 , ， 当 即 。若成立 ， 则称该方法是收敛的。nynx0h1nnxx1()nnyy x1nnxxh定义三定义三： 若一个数

39、值方法在基点 处的值有 的扰动，在 此后各基点 （mn）处的值产生的偏差均不超 过 ，则称该方法是的。 单步法的稳定性有以下定理nymy定理二定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是的。相容性和收敛性的关系相容性和收敛性的关系 若单步法的增量函数对变量y满足Lipschtiz 条件，即存在与 h , x 无关的常数 L，对区域D= 任意两点（x，y1）,(x,y2)成立，则单步法收敛的充分必要条件是相容性条件成立。（读者自证）,axby 相容性和方法阶的关系相容性和方法阶的关系 若单步法若单步法是是p p阶方法则成立阶方法则成立 若单步法满足相容性条件

40、，得若单步法满足相容性条件，得 所以所以 =0=0也就是说单步法的阶数一定要是正数。由于我们考虑也就是说单步法的阶数一定要是正数。由于我们考虑的单步法皆为正整数，的单步法皆为正整数，p p至少为一。因此我们考虑的至少为一。因此我们考虑的单步法都满足相容性条件。单步法都满足相容性条件。 1()( )( , ; )()py x hy xhx y hO h101( )( , )()limphy xf x yO hh0()limphO h 稳定性和收敛性的关系稳定性和收敛性的关系若单步法的增量函数满足定理二的条件即单若单步法的增量函数满足定理二的条件即单步法是稳定的则单步法收敛的充分必要条件步法是稳定

41、的则单步法收敛的充分必要条件是是相容性条件成立相容性条件成立。 ( , ;0)( , )x yf x y绝对稳定性和绝对稳定绝对稳定性和绝对稳定域域 稳定性问题是一个比较复杂的问题。为了简化讨论一稳定性问题是一个比较复杂的问题。为了简化讨论一般仅对试验方程般仅对试验方程 进行考察。这里假定进行考察。这里假定Re0,Re0h0和的允许范围，称为该方法的和的允许范围，称为该方法的 绝对稳定域绝对稳定域。 yynn ky绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域2 2将Euler方法应用到试验方程得误差方程是要求误差不增长则nnnnyhhyyy)1 (1nnnh11|1|1hnn则Euler 方法

42、的绝对稳定域是以 为半径的圆的内部。为了保证稳定性步长有所控制。假如当 时h应满足 ，当 时， h 应满足 等等。同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。 1h5520 h10050110020 h 1.Adams1.Adams方法方法 2.2.米尔尼方法、汉明方法及辛普森方法米尔尼方法、汉明方法及辛普森方法 3.3.预测校正方法预测校正方法 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性 5 初值问题的数值解法初值问题的数值解法多步法多步法 考虑型如考虑型如 的的k步法，步法， 称为称为

43、阿当姆斯阿当姆斯（Adams）方法）方法10knknkiniiyyhf112330,0,00101211211211211224133121443211.Adams1.Adams方法方法拿其中一种为例，推导其公式。对上面所说的法一般形拿其中一种为例，推导其公式。对上面所说的法一般形式式若取若取 ，有，有方程组方程组同样解之，得到同样解之，得到3步步4阶隐式阶隐式Adams公式和它的余项。公式和它的余项。1112(9195)24nnnnnnhyyffff5(5)12119(),(,)720nnnnnTh yxx 其它请读者自证。我们仅将结论列于下表。其它请读者自证。我们仅将结论列于下表。 Ada

44、ms显式公式显式公式1pc1nnnyyhf12211( 3)2nnnnhyyff5123221(23165)12nnnnnhyyfff3843321(5559379)24nnnnnnhyyffff251720kp公式11223344Adams隐式公式隐式公式1pc11()2nnnnhyyff1122121(58)12nnnnnhyyfff12432321(9195)24nnnnnnhyyffff19720434321(25164626410619 )720nnnnnnnhyyfffff3160kp公式12233445注：注：在这些方法中，在这些方法中，4阶的阶的Adams预测校正方法具有方法预

45、测校正方法具有方法简洁、使用方便、易排序、高精度等优点。尤其当函数简洁、使用方便、易排序、高精度等优点。尤其当函数f比较复杂时更能显出它的优越性。比较复杂时更能显出它的优越性。 同 理 得 到同 理 得 到 5 个 待 定 参 数 方 程 组 。 解 之个 待 定 参 数 方 程 组 。 解 之得得 ， ， ， ， 。 构成著名的构成著名的Miline 4Miline 4步步4 4阶显式公式和它的余阶显式公式和它的余项。项。 13380341， 3820313nnyy124(22)3nnnhfff5(5)13114( ),(,)45nnnnnTh yxx2.2.米尔尼方法、汉明方法及米尔尼方法

46、、汉明方法及辛普森方法辛普森方法 同理得到同理得到5 5个参数方程组。求解后就构成著名的个参数方程组。求解后就构成著名的3 3步步4 4阶隐式阶隐式HammingHamming公式和它的余项。公式和它的余项。 若取，也得到若取，也得到5个参数方程组。求解后就构成个参数方程组。求解后就构成Simpson隐式公式和其余项。隐式公式和其余项。 121(9)8nnnyyy113(2)8nnnhfff5(5)1211(),(,)40nnnnnTh yxx米尔尼方法、汉明方法及米尔尼方法、汉明方法及辛普森方法辛普森方法1111(4)3nnnnnhyyfff5(5)1111(),(,)90nnnnnTh y

47、xx 不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。 3. 预测校正方法预测校正方法 对于线性多步公式，一般地，不难验证，如果 预测公式是阶或阶精度，校正公式为阶精度， 则用预测公式提供初值，校正公式迭代一次的 效果也能达到阶精度，再迭代下去，效果就不 明显了。预测-校正技术即保证了计算精度，又 使隐式计算显式

48、化，克服了隐式公式要反复迭 代的困难。至于实际计算步长的选取，我们对 预测-校正公式使用外推原理，得到误差估计式 ，用来调整计算步长，使达到控制误差要求。对于线性多步方法常用的预测-校正公式有Miline-Hamming方法和4阶阶Adams显隐式显隐式预测-校正公式 将Miline Miline 公式公式和HammingHamming公式公式结合，构成预测-校正公式 Miline公式公式和Hamming公式公式的局部截断误差分别为13124(22)3pnnnnnhyyfff1211113(9)( (,) 2)88pnnnnnnnhyyyf xyff5(5)611114()()()45MMnn

49、nnTy xyh yxO h5(5)61111()( )()40HHnnnnTy xyh yxO h修正修正Miline-Hamming公式公式利用外推原理，即加上后消去局部截断误差的主项。使 说明经过外推后的算法其精度提高一阶。忽略误差项，上式可改写为 11611144045()()1144045MHnnnyyy xO h111114114() ()040454045MHnnny xyy11191129112() ()0360360360360MHnnny xyy 由于 和是 在计算过程中获得的数据，称为Miline Miline 公式公式和HammingHamming公式公式的事后误差估计

50、式。我们可以用它们来调节计算步长的大小，即选择一个合适的的步长，使 ，其中的是要求达到的计算精度。 111121 ()91120MHnnny xyy1111112()()121MHMnnnny xyyy11119()()121HHMnnnny xyyy1Mny1Hny119()121HMnnyy 又可得到 MilineMiline公式公式和HammingHamming公式公式的修正公式，它们分别是 从而构成如下的修正从而构成如下的修正HammingHamming预测预测- -校正公式，校正公式，简称修正简称修正HammingHamming公式：公式： 1111112()121MmMHMnnnn

51、yyyy11119()121HmHHMnnnnyyyy预测：修正：校正： 修正： 在应用这套公式时，先由同阶单步法提供初值 ， ， ， 。计算 时可取。31npnyy)22(3421nnnfffh)(12111211pncnpnpmnyyyy)2),(83)9 (8111121nnpmnnnncnffyxfhyyy)(12191111pncncnnyyyy0y1y2y3y4ypcyy33 将将4步步4阶显式阶显式Adams公式作为预测公式和公式作为预测公式和3步步4阶阶式式Adams公式作为校正公式，构成公式作为校正公式，构成Adams预测预测-校正公式。校正公式。 它们的局部截断误差分别是它

52、们的局部截断误差分别是 1123(5559379)24pnnnnnnhyyffff1112(9195)24cnnnnnnhyyffff5(5)6111251()()()720ppnnnnTy xyh yxO h5(5)611119()( )()720ccnnnnTy xyh yxO hAdams预测预测-校正公式校正公式利用外推原理，将上两式作线性组合消去局部截断误差的主项。使计算精度至少提高一阶，同时得到两个修正公式，将它们和上两式构成了如下修正修正Adams公式公式：预测： 修正： 校正: 1123(5559379)24pnnnnnnhyyffff11251()270pmpcpnnnnyy

53、yy11112(9 (,) 195)24cpmnnnnnnnhyyf xyfff修正： 同理，在计算时，调节计算步长 使 , 由同阶单步法提供初值 , ， ， 。计算 时，可取 。 理论分析得出它们的绝对稳定区域，修正Hamming法： ; 4阶Adams预测预测-校正法校正法： 其中 )(270191111pncncnnyyyyh)(2701911pncnyy0y1y2y3y4ypcyy33069. 0h075.0h)(,(xyxyf我们也可以在教学演示上看出修正的预测校正格式的误差非常的小。 多步法的相容性条件多步法的相容性条件 k步法的一般形式为 其中 由对单步法的讨论可知,当 时,方法

54、阶数至少为1。即对 y1,y=x 应精确成立。令 y 1，所以令y=x得 101111011 ()nnnrn rnnnrn ryyyyhffff0| , 000k0 ,y100 ,kk1y4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性所以所以我们称为我们称为线性多步法的相容性条件。线性多步法的相容性条件。 1010101110() (1) ()(1)nkkkkkkkkkkxh kkhkk 即 k步法需要k个出发值，而初值问题只提供了一个初值，在使用k步法时尚需要其它方法作补充k-1个出发值，今假定它们 是 ，当 这k-1个值都 应收敛于共同的极限y（a）即 在讨论多步法收

55、敛性时我们总假定（9.3.12）成立 定义四定义四： 1,.2, 1),(klhyll0h 1,.2 , 1,)(limlim00klayhylhlhnknkknknknkknkfffhyayaya011011. 的解 收敛于初值问题的解 y(x)。这里 xa+nh. 定义五定义五：称多项式 为k步法的特征多项式。如果特征多项式的零点皆不大于1，且等 于1的零点是单重的，称根条件成立。称多项式 为第二特征多项式。 显然根条件可以表示 定理定理三三：k步法收敛的充要条件为： （1）相容性条件成立。 （2）特征多项式的零点皆不大于1，且等于1的零点是单重的。 （称2为）特征根条件。 ny,baxn

56、x011.)(kkkk011.)(kkkk) 1 () 1 (, 0) 1 (多步法的稳定性 关于线性多步法成立以下定理：关于线性多步法成立以下定理： 定理四：定理四：若函数若函数f(x,y)f(x,y)对变量对变量y满足满足Lipschtiz 条件在条件在与与h,x无关的常数无关的常数L，对区域对区域D= D= 任意两点任意两点（x x，y1y1）,(x,y2),(x,y2)成立成立 k步法的相容性、收敛性、稳定性有以下关系步法的相容性、收敛性、稳定性有以下关系 对于常微分方程右端函数对于常微分方程右端函数f(x,y)，在相容性条件成立情，在相容性条件成立情况下，况下，k步法的收敛性和稳定性

57、是等价的。步法的收敛性和稳定性是等价的。 事实上在相容性条件成立时，收敛性和稳定性的充要条事实上在相容性条件成立时，收敛性和稳定性的充要条件都是特征根条件成立。件都是特征根条件成立。ybxaD,| 12| ); 2,(); 1,(|yyLhyxhyx多步法的绝对稳定性和绝对稳定域多步法的绝对稳定性和绝对稳定域 将线性多步法的公式应用到试验方程 进行考察。这里假定Re 0,即试验方程本身是稳定的。 记 得 是常系数差分方程，其特征方程为 记它的 k个特征根为 并设它们是互异的。显然根与 有关，不妨记为 注意到当 时 这时由特征方程得 由线性多步法的相容性条件得 是一个根。不妨设, 差分方程的解为

58、 其中系数由线性多步法的出发值确定。yyhu inkiikiinkyy000)()(kii,.2, 1, hu kii,.2 , 1),(0h00)(10, 0)(1nkknnndddy.2211另一方面，y(0)=1的试验方程的精确解为 ， 设多步法截断误差为 ，由此可得 ，我们称 为主根，其它根都为增根。 定义五定义五：线性多步法的绝对稳定区域对给定的 ，如果特征方程的特征根 皆按模小于1，则线性多步法关于u是绝对稳定的。使得 成立的 构成绝对稳定区域。注：从误差角度来看绝对稳定区域的方法是一个理想的方法。这样，绝对稳定区域越大，方法适用性越广，因而越优越。 xexy)()(1phO)(1

59、1phhOe1ui1|ihu 6 6 方程组和刚性方程方程组和刚性方程 1 一阶方程组 2 化高阶方程为一阶方程组 3 刚性方程组1. 一阶方程组 前面我们研究了单个方程y=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形。 考虑一阶方程组的出值问题 若采用向量的记号，记miyayyyyxfyiimii,.,2 , 1,)(),.,(021TmTmTmyxfyxfyxfyxfyyyyyyy),(),.,(),(),().,(),.,(y210,2010021 则常微分方程组的出值问题可以表示为 前几节我们主要讨论了常微分方程的出值问题的数值解法。只要将

60、 y 和f 改写为向量，那么前面提供的各种计算公式即可适用于一阶常微分方程组的出值问题。 Runge-Kutta方法 对于方程组 四级四阶显示Runge- Kutta公式 0)(),(yyayyxf0)(),(yyayyxf1123412132431(22)6(,) (,)21 (,)22 (,)nnnnnnnnnnyykkkkkhf xyhkhf xykhkhf xykkhf xh yk若写成分量形式就是 i=1,2,m 为了帮助理解这一公式的计算过程，我们再考虑两个方程的特殊情形： ),()21,2()21,2(),()22(61342312143211kyhxhfkkyhxhfkkyhx