第7章-2 收敛性稳定性R-K方法

1、常微分方程的数值解法的常微分方程的数值解法的 收敛性、稳定性收敛性、稳定性第第7 7章章-2-2 以上我们讨论了求解问题（以上我们讨论了求解问题（7-1）, ,（7-2）的单步）的单步法法和多步法。和多步法。应关注三个问题：应关注三个问题：、数值方法的局部截断误差和阶、数值方法的局部截断误差和阶二、在离散点二、在离散点tn处的数值解处的数值解un是否收敛到精确解是否收敛到精确解u( (tn) ) 三、数值方法的稳定性三、数值方法的稳定性具体说，具体说，对于上述两类方法求近似解（数值解）还对于上述两类方法求近似解（数值解）还误差估计误差估计、收敛性和稳定性收敛性和稳定性。 对于第一个问题前面我们

2、已经讨论过，而关于数值对于第一个问题前面我们已经讨论过，而关于数值方法收敛性问题我们在这里不详细讨论，只给出一些基方法收敛性问题我们在这里不详细讨论，只给出一些基本结论性的结果，即：本结论性的结果，即：对单步法，当方法的阶对单步法，当方法的阶p11时，有整体误差时，有整体误差)()(pnnnhOutuE故有故有 0lim0nhE，因此方法是收敛的。，因此方法是收敛的。 对于多步法，若方法是对于多步法，若方法是k 步步p 阶法，那么（阶法，那么（7-24）是）是 一个一个k阶差分方程，引入多步法（阶差分方程，引入多步法（7-24）的）的第一特征第一特征多多项项式式和和第二特征多项式第二特征多项式

3、： 0( ),kjjj kjjj0)(定义定义7.17.1 若（若（7-24）的第一特征多项式）的第一特征多项式( () )的所有的所有 根在单位圆内或圆上（根在单位圆内或圆上（11），且位于单位圆周上），且位于单位圆周上的根都是单根，称多步法（的根都是单根，称多步法（7-24）满足）满足根条件根条件。第二特征多项式第二特征多项式第一特征多项式第一特征多项式定理定理7.27.2 若线性多步法若线性多步法( (7-24) )的阶的阶p11，且满足，且满足根条件，则方法是收敛的。根条件，则方法是收敛的。 对于常用的数值方法都是满足收敛性条件的对于常用的数值方法都是满足收敛性条件的。 下面我们着重讨

4、论第三个问题，即下面我们着重讨论第三个问题，即数值方法的稳数值方法的稳是有误差的，且这些误差将在计算中传递下去。是有误差的，且这些误差将在计算中传递下去。定性定性问题。问题。110,kuuu误差积累无限增长，则会歪曲真解，这样的算法是不误差积累无限增长，则会歪曲真解，这样的算法是不如果如果能用的。能用的。用多步法计算时，各种因素如初值用多步法计算时，各种因素如初值精确解为精确解为 2 2( ) (1)u tt 。考虑二步考虑二步三阶显式法三阶显式法： 21145(42)nnnnnuuuhff例如例如 初值问题初值问题124,utu 02;t (0)1u取步长取步长h=0.1=0.1，初值初值u

5、0=1=1，附加值：附加值： 2 21(1) (0.1)uhh。数值结果表数值结果表在开始几步数值解与精确解符合，但在再往后算，数值在开始几步数值解与精确解符合，但在再往后算，数值解的误差则急剧增长，完全歪曲了真解解的误差则急剧增长，完全歪曲了真解. .通常人们都是通过模型方程来讨论方法的通常人们都是通过模型方程来讨论方法的数值稳定性数值稳定性。uu( (7-32) )而一般形式的一阶微分方程总能化成（而一般形式的一阶微分方程总能化成（7-32）的形式。）的形式。因为实际计算时，。因为实际计算时，h是固定的。是固定的。 当某一步当某一步un有舍入误差时，有舍入误差时， 若以后的计算中不会逐步扩

6、大，称这种稳定性为若以后的计算中不会逐步扩大，称这种稳定性为绝对稳定性绝对稳定性。此后，若不做特殊说明，都是指绝对稳定性此后，若不做特殊说明，都是指绝对稳定性 。模型方程为：模型方程为：本书中数值方法的稳定性也是如此。前提是求解好条件问题，本书中数值方法的稳定性也是如此。前提是求解好条件问题，其中其中Re( () )0 0。另外，我们也不考虑另外，我们也不考虑h00时方法的渐近稳定性时方法的渐近稳定性 例如例如，对最简单的，对最简单的Euler法法2, 1, 0,1nhfuunnn( (7-33) ) 用其求解模型方程（用其求解模型方程（7-32）得到）得到1nnnuuh u 当当un有舍入误

7、差时，其近似解为有舍入误差时，其近似解为 nu，从而有，从而有 nnuhu)1(1取取 nnnuu，得到误差传播方程，得到误差传播方程 ,)1 (1nnh(1),0,1, 2nh un记记 hh，只要，只要 11h都不会恶性发展，此时方法绝对稳定。都不会恶性发展，此时方法绝对稳定。，则显式，则显式Euler方法的解和误差方法的解和误差从从 11,h可得可得 20h 。即即20 h 时，时，(-1,0)(-1,0)为圆心，为圆心，1 1为半径的单位圆。为半径的单位圆。 又由于实数又由于实数0 0，（7-33）绝对稳定，绝对稳定，若若为复数，在为复数，在 hh的复平面上，则的复平面上，则 11h表

8、示为以表示为以20110, 1绝绝对对稳稳定定区区域域绝对稳定区间绝对稳定区间 定义定义7.27.2 一个数值方法用于求解模型问题（一个数值方法用于求解模型问题（7-32），若在），若在 平面中的某一区域平面中的某一区域D中方法都是绝对稳定的，而在区域中方法都是绝对稳定的，而在区域D外，方法外，方法是不稳定的是不稳定的，则称则称D是方法的是方法的绝对稳定区域；绝对稳定区域；绝对稳定区间绝对稳定区间。它与实轴的交称为它与实轴的交称为例如，显式例如，显式Euler方法的方法的绝对稳定区域、区间。如图绝对稳定区域、区间。如图现在考察多步法现在考察多步法( (7-24) )，将它用于解模型方程，将它用

9、于解模型方程( (7-32) )得到得到k阶线性差分方程阶线性差分方程0kjnjju(7-34) 若取若取 hh，则记（，则记（7-34）的特征方程为）的特征方程为0)()(h(7-35) 其中其中 kjjj0)(kjjj0)(0kjnjjhu由由k阶线性差分方程的性质我们可以得到如下结论，阶线性差分方程的性质我们可以得到如下结论，区域：区域： ()1,1,2,jhhjkD例如，对于例如，对于k=1=1时，考虑隐式方法中最简单的后退时，考虑隐式方法中最简单的后退Euler法法111(,)0, 1,nnnnuuh f tun方程（方程（7-35）的根都在单位圆内（）的根都在单位圆内（1) 1)

10、，则线性多步法，则线性多步法（7-4）关于）关于 hh绝对稳定绝对稳定, ,其绝对稳定域是复平面其绝对稳定域是复平面 h上的上的其特征方程为：其特征方程为：( )( )h (1)10h 若特征若特征得得 11,1h当当 11 h时，时， 11 ,故故 11h就是就是隐式隐式Euler法的绝对稳定区域。法的绝对稳定区域。当当0 0为实数时，绝对稳定区间为为实数时，绝对稳定区间为 (-,0)(-,0)。h平面上以（平面上以（1.01.0）为圆心的单位圆外区域。）为圆心的单位圆外区域。它是它是11h11h0Re2,0h0,当当Re 0 0时，它位于时，它位于 平面上平面上y轴左侧区域轴左侧区域。h又

11、如，梯形法又如，梯形法, 1 , 0)(2111nffhuunnnn其特征方程为：其特征方程为：( )( )h 其根其根1( )h当当Re0 0时，时，121,12hh故梯形公式故梯形公式h平面的左半平面。绝对稳定区间为平面的左半平面。绝对稳定区间为(-(-，0 0) )。的绝对的绝对稳定域是稳定域是 这样检验绝对稳定性归结为检验特征方程这样检验绝对稳定性归结为检验特征方程( (7-35) )的根是否在单位的根是否在单位圆内（圆内（1)1)。对此有对此有很多判别法，如很多判别法，如Schur准则、轨迹法。准则、轨迹法。12,12hh11022hhk=1=14 4的隐式的隐式Adams类方法的绝

12、对稳定区间（类方法的绝对稳定区间（0 0为实数）。为实数）。实系数二次方程实系数二次方程2 2- -b - -c=0=0的根在单位圆内的充要条件为的根在单位圆内的充要条件为: : 21cb这里我们给出一种简单的、常用的判别法：这里我们给出一种简单的、常用的判别法：例例 证明求解一阶常微分方程初值问题：证明求解一阶常微分方程初值问题：),( utfu 0)0(uu的差分格式的差分格式)85 (121212nnnnnfffhuu收敛收敛并求其并求其局部截断误差主项局部截断误差主项、绝对稳定区间绝对稳定区间。解解：由差分格式可知，：由差分格式可知，2( ), ()(1)0 , 则则其特征值其特征值满

13、足根条件满足根条件。令令得得1=0,=0,2=1=1。2581( ),121212 故此为故此为隐式二步三阶隐式二步三阶法，其局部截断误差主项为：法，其局部截断误差主项为：4(4)1()24nh ut。01 10C 1C 2C 3C 4C 1852 1121212 013( 14)22 032118512262 1212 776604311851224!3! 12120注意，注意，0120,1,1,从而从而012185,121212由定理由定理7.27.2 可知，此方法收敛。可知，此方法收敛。( )( )h 12 812 5hh而而 1212 5hh自然成立。自然成立

14、。hhhh512412512812得得 12 4h 12 812 4hh即有即有 332,hh 可得其可得其绝对稳定区间绝对稳定区间： 60h 。又其特征方程为又其特征方程为21280125125hhhh而使得而使得 1 1的充要条件为：的充要条件为：124125hh12125hh现在再由现在再由128 h124 h2581110121212hhh进一步进一步而而 自然成立。自然成立。 显式显式Runge-KuttaRunge-Kutta法法第第7 7章章-3-37.1.47.1.4 四阶显式四阶显式Runge-Kutta法法 通过观察我们发现显式通过观察我们发现显式Euler法和隐法和隐Eu

15、ler法各用到了法各用到了u( (t) )在在 t,t+h 上的一个一阶导数值，它们都是一阶方法。上的一个一阶导数值，它们都是一阶方法。改进的改进的Euler法用到了法用到了u( (t) )在在 t,t+h 上的两个一阶导数值，它上的两个一阶导数值，它 梯形法和梯形法和们都是二阶方法。们都是二阶方法。 我们要研究的我们要研究的Runge-Kutta方法是一种高阶单步法，它使用方法是一种高阶单步法，它使用u( (t) )在在 t,t+h 上的斜率上的斜率f 在一些点的值非线性表示在一些点的值非线性表示 使得其局部截断误差的阶和使得其局部截断误差的阶和Taylor展开法相等。展开法相等。( , (

16、 ), )t u t h Euler是最简单的单步法。单步法不需要附加初值，所需的是最简单的单步法。单步法不需要附加初值，所需的存储量小，改变步长灵活，但线性单步法的阶最高为存储量小，改变步长灵活，但线性单步法的阶最高为2 2，Taylor展开法，用在同一点展开法，用在同一点( (tn,un) )的高阶导数表示的高阶导数表示 ，这不便，这不便于计算。于计算。( , ( ), )t u t h先引进若干记号，首先先引进若干记号，首先 t, t+h 取上的取上的m个点个点: :123,mtttttth令令, 2,1mihathattiii 11,2,iijijbaimmiiicc11, 02131

17、32121mmmmbbbbbb Runge-Kutta矩阵矩阵B为严格下三角矩阵：为严格下三角矩阵：满足满足显式显式 Runge-Kutta Runge-Kutta 公式公式假设三组系数已给定，则求解（假设三组系数已给定，则求解（7-1），（），（7-2）的一般）的一般1( , ),nnnnuuht u h （7-12）其中其中1( , )mnniiituhc k( (7-13) )( (7-14) )显式显式Runge-Kutta法的计算过程如下：法的计算过程如下：1( ,),nnkf t u2221 1(,) ,nnkf tha uhb k3331 1322(,(),nnkf tha uh

18、 b kb k31323,bba11(,),mmnmnmjjjkf thauhb k 212,ba11mmjmjba0,1,n 现在推导一些常用的计算方案，特别地，给出现在推导一些常用的计算方案，特别地，给出 m=3 =3 显式显式首先将首先将u( (t+h) )在在t处展开到处展开到h的三次幂的三次幂, ,即：即：3( )41()( )( )()!lllhu thu tutO hl( (7-15) )其中其中( (7-16) ) Runge-Kutta法的推导。法的推导。( )( , ( ), ),u tht u t h2311( , , )()()26ut u hfhfhfffO htuf

19、fff22tttuuuffffff其次，由二元函数其次，由二元函数f(t,u(t)在在( (t, ,u) )点处的点处的Taylor展开式可得：展开式可得： ,)(,(1ftutfk),(1222khauhatfk21()tufha fk f)(2132222hOfahfhaf)()21(323322233hOfaffbahfhafku于是，将于是，将k1 1, ,k2 2, ,k3 3代入（代入（7-13）中，即）中，即31( , , )iiit u hc k1c f22332 32312ucfa f ha b f fa fh1 12233c kc kc k222321112()2tttuu

20、uh afk fk fO h222312cfha fa f h3()O h1232 23 3cccfh a ca cf31( , , )iiit u hck22232 32 3223312()()2uha b c ffa ca cfO h( (7-17) ) 由由( (7-16) )已得已得(3)211( , , )( )( )( )26t u hu tu th uth其中其中 ,tuffff22tttuuufffff f。2311()26ufh fhfffO h合并合并f( (t, ,u, ,h) )展开式中的各阶展开式中的各阶hl( (l= =0,1,2 2) )的系数，得的系数，得比较比

21、较 ( , , )t u h和和 ( , , )t u h的同次幂系数，的同次幂系数，231111()2233ufh fhfffO h可得可得（一）（一）m=1 =1 此时此时 c2= =c3= =0, , f( (t, ,u, ,h)=)=c1 f , , 比较比较h的零次幂，知的零次幂，知( , , ),t u hf方法（方法（7-21）为一级一阶）为一级一阶Runge-Kutta法，实际上为法，实际上为Euler法。法。 （二）（二）m=2,=2,此时此时 c3=0，则则223122 2221( , , )()2t u hccfha c fh a c fO h与与 ( , , )t u

22、h比较比较1, ,h的系数，则的系数，则 2112221cacc它有无穷多组解，从而有无穷多个二级二阶方法。它有无穷多组解，从而有无穷多个二级二阶方法。2311()26ufh fhfffO h( , , )t u h,21, 1, 0221acc（1 1）称为称为中点法中点法。此时此时12,nnuuhk1( ,),nnkf t u2111,22nnkfth uhk(7-18) 三个常见的方法是：三个常见的方法是：1221,1,2cca（2 2）称为称为改进的改进的Euler法法。此时此时112,2nnhuukk1( ,),nnkf t u21(,)nnkf th uhk(7-19) 12213

23、2,443cca（3 3）此时此时1123,4nnhuukk1( ,),nnkf t u2122,33nnkfth uhk（三）（三）m=3 3 比较（比较（7-16）和（）和（7-17），令），令 f,f,h, ,h2 2 的系数的系数相等，并注意相等，并注意 ff,的任意性，得的任意性，得 21, 13322321cacaccc.61,313322323222cbacaca四个方程不能完全确定六个系数，因此这是含两个参数的四个方程不能完全确定六个系数，因此这是含两个参数的三级三级三阶三阶方法类。方法类。常见方案有：常见方案有：12313,0,44ccc2332122,333aab。 Heu

24、n三阶方法三阶方法。 此时取此时取1133,4nnhuukk1( ,),nnkf t u2111,33nnkfth uhk(7-20) 3222,33nnkfth uhk（2 2）Kutta三阶方法三阶方法，123121,636ccc23321,1,22aab 。),4(613211kkkhuunn),(1nnutfk ),21,21(12hkuhtfknn).2,(213hkhkuhtfknn(7-21) (7-21) 此时此时（四）（四）m=4=4将（将（7-167-16）和（）和（7-177-17）展开到）展开到h3, ,比较比较)3 , 2 , 1 , 0( ihi的系数，则含的系数，

25、则含1313个待定系数的个待定系数的1111个方程，由此得到含两个参数的个方程，由此得到含两个参数的四级四阶四级四阶Runge-Kutta方法方法类，类，其中最常用的有以下两个方法：其中最常用的有以下两个方法： 经典四阶经典四阶Runge-Kutta方法方法：),22(643211kkkkhuunn),(1nnutfk 2111,22nnkf th uhk3211,22nnkf th uhk).,(34hkuhtfknn（7-22）),33(843211kkkkhuunn),(1nnutfk 2111,33nnkf th uhk31221,33nnkf th uhkhk).,(34hkuhtf

26、knnButcher表分别为：表分别为：10021021613131611212101001313181838381132310以上讨论的是以上讨论的是m级级Runge-Kutta法在法在m=1,2,3,4=1,2,3,4时，可分别时，可分别得到最高阶级一、二、三、四阶，但是，通常得到最高阶级一、二、三、四阶，但是，通常m级级Runge-Kutta Runge-Kutta 方法最高阶不一定是方法最高阶不一定是m阶。若设阶。若设p( (m) )是是m级级Runge-Kutta方法可方法可达到的最高阶，可证达到的最高阶，可证：(5)4,(6)5,(7)6,(8)6,(9)7ppppp 。1221,

27、1nnnnnt uuuht改进的改进的Euler法计算公式为：法计算公式为：经典经典Runge-Kutta法计算公式为：法计算公式为：例例1 分别用分别用EulerEuler法，改进的法，改进的EulerEuler法法(7-27)(7-27)和经典和经典Runge-Kutta法法(7-30)(7-30)求解初值问题求解初值问题: :221,1tuut (0)0u02,t 解解：Euler法计算公式为：法计算公式为： 1121(),2nnuuh kk1221,1nnnt ukt 1222()()1,1 ()nnnth uhkkth ),22(643211kkkkhuunn1221,1nnnt ukt ,)21(1)21)(21( 21212hthkuhtknnn232112221,1 ()nnnthuhkkth .)(1)(21234hthkuhtknnn三个方法计算结果比较表三个方法计算结果比较表作比较作比较 ，计算结果见下表：，计算结果见下表：取步长取步长h=0.5,=0.5,tn=0.5=0.5n, ,n