数值分析第9章

1、1第第9章章 常微分方程初值问题数值解常微分方程初值问题数值解法法9.1 引引 言言2本章研究的问题本章研究的问题:3n9.2 欧欧 拉拉 方方 法法n9.2.1 欧拉公式欧拉公式1 欧拉公式欧拉公式4图9.1 欧拉折线法5(2)(3)672 欧拉公式的截断误差欧拉公式的截断误差(4)83单步法的局部截断误差与阶单步法的局部截断误差与阶局部截断误差可以理解为计算一步的误差局部截断误差可以理解为计算一步的误差.9局部截断误差可以理解为计算一步的误差局部截断误差可以理解为计算一步的误差.10则称该方法具有则称该方法具有P阶精度阶精度.)(nxyp定义定义2 设设是初值问题的准确解，若存在最大整数是

2、初值问题的准确解，若存在最大整数使显式单步法的局部截断误差满足使显式单步法的局部截断误差满足1111),()()(pnnnnnnnhOhyxhxyhxyyxyT若把上式展开写成若把上式展开写成)()(,(211ppnnnhOhxyxT1)(,(pnnhxyx则则 称为局部截断误差主项称为局部截断误差主项 。 114 后退的欧拉方法后退的欧拉方法(5)(6)(6)式称为式称为后退的欧拉方法后退的欧拉方法,它是它是隐式的隐式的,欧拉公式欧拉公式(2)是是显式的显式的,1213(7)(6)1415后退的欧拉方法的局部截断误差后退的欧拉方法的局部截断误差:165 梯形方法梯形方法(8)式称为梯形方法式

3、称为梯形方法.(8)17梯形方法的局部截断误差梯形方法的局部截断误差:189.2.4 改进欧拉法及局部截断误差改进欧拉法及局部截断误差nnnnyxhfyy,1111,2nnnnnnyxfyxfhyy预测步预测步校正步校正步cpnpnncnnnpyyyyxhfyyyxhfyy21,11或者写成或者写成1.改进的欧拉公式改进的欧拉公式:192.改进的欧拉方法的局部截断误差改进的欧拉方法的局部截断误差20),(1kkkkyxhfyy),(),(2111kkkkkkyxfyxfhyy考虑改进考虑改进Euler法法如果将其改成如果将其改成)(2211KKhyykk),(1kkyxfK ),(12hKyh

4、xfKkk)(00 xyy -(1)9.3 Runge-Kutta法法21改进改进Euler法是由梯形公式和法是由梯形公式和Euler公式复合而成公式复合而成梯形公式具有梯形公式具有2阶精度阶精度(1)式为一种式为一种二阶二阶Runge-Kutta法法同样可以证明同样可以证明,改进改进Euler法也具有法也具有2阶精度阶精度22Runge-Kutta方法的推导方法的推导23Runge-Kutta方法的一般形式：方法的一般形式：riKhyhxFKyxfKKchyyijjijnininnriiinn, 3 , 2),(),(11111确定了阶数之后，再通过确定了阶数之后，再通过Taylor展开、比

5、较两边系数的方展开、比较两边系数的方法，确定各待定系数：法，确定各待定系数：,iiijc 24二阶显式二阶显式Runge-Kutta方法方法2526272829例例.1)0(,:ODE2yydxdy求求解解初初值值问问题题.11xy易易知知其其精精确确解解为为：21 . 01 . 0. 1,21Euler,21. 1222121221nnnnnyyyyycca积积分分公公式式：法法：改改进进的的321 . 0321 . 0.23,31,32,31. 2222121221nnnnnyyyyycca积积分分公公式式：分分别别用用以以下下两两种种系系数数：步步长长都都取取为为1 . 0h30结果及比

6、较结果及比较31三阶显式三阶显式Runge-Kutta方法方法nnnYnntnnYnnnYYnnntYnnttn)F,Y(tF),Y(tF),Y(tF)F,Y(tF)F,Y(tF),Y(tFY 22在推导二阶显式方法的过程中，注意到局部截断误差表达式中在推导二阶显式方法的过程中，注意到局部截断误差表达式中h3项包项包含了以下表达式：含了以下表达式：因此若要在局部截断误差中消去因此若要在局部截断误差中消去h3项，必须增加包含了以上各项的多项，必须增加包含了以上各项的多个方程，同时我们注意到个方程，同时我们注意到r=2时，只有时，只有 等四个待定系数，少等四个待定系数，少于方程的数目，所以这样的系

7、数不存在。故：于方程的数目，所以这样的系数不存在。故： r=2时时Runge-Kutta方方法只能是二阶的法只能是二阶的。要得到三阶的方法，则必须有。要得到三阶的方法，则必须有r=3。),(),(),()(232131331212213322111hKhKYhtFKhKYhtFKYtFKKcKcKchYYnnnnnnnn其其局局部部截截断断误误差差为为：).()()(33221111KcKcKchtYtYdnnn1,2121,c 32三阶显式三阶显式Runge-Kutta方法方法程程如如下下：可可得得待待定定系系数数满满足足的的方方展展开开，使使得得作作以以及及将将),(Taylor)(,41

8、132hOdtYKKnn61312113223233222332232313212321cccccccc)2,()2,2(),()4(62131213211hKhKYhtFKKhYhtFKYtFKKKKhYYnnnnnnnn三三阶阶公公式式，如如下下：常常见见的的公公式式称称为为公公式式。特特别别地地，一一个个三三阶阶，它它们们统统称称为为因因此此可可以以得得到到众众多多公公式式个个未未知知数数，解解不不唯唯一一。个个方方程程要要决决定定KuttaKuttaRunge8633四阶显式四阶显式Runge-Kutta方法方法个个：。下下面面列列出出最最常常见见的的一一的的局局部部截截断断误误差差满

9、满足足公公式式，它它们们导导出出各各种种四四阶阶的的类类似似前前面面的的推推导导，可可以以)(KuttaRunge51hOdn),()2,21()2,21(),(226342312143211hKYhtFKKhYhtFKKhYhtFKYtFKKKKKhYYnnnnnnnnnn34.1)0(,:ODE2yxydxdy求求解解初初值值问问题题xexxy222易易知知其其精精确确解解为为：方方法法求求解解：分分别别用用二二阶阶、四四阶阶步步长长都都取取为为KR1 . 0hx四阶二阶真解四阶误差二阶误差0.01.0000001.0000001.0000000.00000.0000000.11.1048

10、291.1024501.1048291.60E-72.38E-30.21.2185971.2115071.2185973.40E-77.09E-30.31.3401411.3257661.3401415.48E-71.44E-20.41.4681751.4436711.4681757.69E-72.45E-20.51.6012781.5635061.6012799.95E-73.78E-20.61.7378801.6833741.7378811.20E-65.45E-20.71.8762461.8011791.8762471.42E-67.51E-20.82.0144571.9146032.0

11、144591.68E-69.99E-20.92.1503952.0210862.1503971.96E-61.29E-11.02.2817162.1178002.2817182.32E-61.64E-1例例35结果及比较结果及比较36结果及比较结果及比较.1045.1505.0).10(1060.1101040.2510.01514时时误误差差为为而而二二阶阶公公式式，相相对对误误差差仅仅为为仍仍然然是是相相当当精精确确的的结结果果时时误误差差为为时时误误差差为为对对四四阶阶公公式式，xhxxh37关于关于Runge-Kutta方法方法RungeKuttaRungeKutta类似前面的推导,可

12、以导出更高阶的公式.关于方法,有以下几点需要特别指出：。解解曲曲线线比比较较光光滑滑的的情情形形别别适适用用于于展展开开的的方方法法，因因此此它它特特方方法法的的推推导导基基于于TaylorKuttaRunge.1次次右右函函数数。、分分别别须须计计算算阶阶数数相相同同，即即它它们们每每步步数数的的次次数数和和方方法法，每每一一步步计计算算右右函函二二阶阶、三三阶阶、四四阶阶的的432Kutta-Runge. 27)9(, 6)8(, 6)7(, 5)6(, 4)5()(KuttaRunge)4(.3NNNNNvNvN阶阶数数，则则有有：次次右右函函数数可可达达到到的的最最高高表表示示计计算算

13、若若用用比比阶阶数数大大。次次的的次次数数方方法法每每步步须须计计算算右右函函数数阶阶的的的的波波动动。，局局部部误误差差会会有有比比较较大大如如果果采采用用固固定定步步长长计计算算的的步步长长等等等等因因素素相相关关。微微分分方方程程的的性性质质、采采用用方方法法具具体体的的系系数数、待待解解阶阶数数、比比较较复复杂杂，它它和和方方法法的的法法的的局局部部截截断断误误差差估估计计KuttaRunge.438提高提高Runge-Kutta方法的精度的方法方法的精度的方法提高积分方法的精度,我们最熟悉的(不一定是最好的)措施是1( )212()2212Euler( )( )211( )( )24

14、nnnhhyyhyyxy xc hc hhyxy xc hc h我们用一个例子予以说明如下法的近似解：将步长减半为时，有外外推推法法Richardson0h将上二式作适当线性组合,可使的一次项为 ：)()(2)()()2(xyxyxyhh。计计算算量量同同时时增增加加了了一一倍倍但但我我们们注注意意到到右右函函数数的的得得精精度度提提高高一一阶阶。利利用用此此式式计计算算时时，可可使使提高精度最简单的方提高精度最简单的方法是缩短步长，但要法是缩短步长，但要以牺牲计算速度和积以牺牲计算速度和积累舍入误差为代价。累舍入误差为代价。39变步长的变步长的Runge-Kutta方法方法作为妥协，如果能在

15、计算过程中实时控制步长的大小，就可以在获得较高作为妥协，如果能在计算过程中实时控制步长的大小，就可以在获得较高的计算速度的同时，保证较高的精度。的计算速度的同时，保证较高的精度。( )( )( )111111()()()()2222111pRungeKutta()2()222hhhpnnnphhhhpnnnpdY tYchhhhdY tYcc一般地,设有 阶的公式,其局部截断误差为将步长减半为时,有()()2()()1211()()()()()()2221111111(1)211;.1221hhhhpnnphhhhhhnnnnnnppcccYYchdYYdYY假定,记为 ,则可以估计误差如下：

16、因此可以从两次计算当中估计出每一步的截断误差,有了这个误差估计之后,通过与控制误差限比较,就可以控制步注意这个方法增加了长.计算量.40Runge-Kutta-Fehlberg方法方法Fehlberg设计了一个更加精巧的嵌套方法如下：设计了一个更加精巧的嵌套方法如下：pp+1:在采用一个 阶方法的同时,计算一个阶的结果,并由此给出误差估计11p()pnnYY thch阶的方法：#21p1()pnnYY thch阶的方法：p因此 阶方法的局部截断误差可以近似为：.1#11nnnYYd制制步步长长。可可以以用用这这个个估估计计式式来来控控RungeKutta,pp+1iiijc我们记得方法中的待定

17、系数不是唯一确定的.因此就有可能利用这个特点,选择适当的系数使得 阶和阶公式中尽可能多的系数相同,从而达到减少计算量的目的.pp+1ic实际上,可以找到仅有权系数不同的 阶和阶算法,如此一来,右函数的计算就可以大大减少.41Runge-Kutta-Fehlberg方法方法Fehlberg给出的四阶、五阶公式给出的四阶、五阶公式RKF4(5)如下：如下：552040114104185925653544227821509514104845513368082164391564302856141042197217972962179720021791932131212825665625651408329

18、32383004141135162162500#54321iiiiiiiicc61#1111RKF4(5);( ,);(,). (2,6)inniinnininijjijYYc KKhF t YKhF th YKi实际使用时,采用五阶公式计算下一步点的值：6#11()niiiidcc K计算中并不计算四阶的结果,只计算误差估计控制步长：42Runge-Kutta-Fehlberg方法方法七阶、八阶七阶、八阶RKF7(8)8404101041126433825141002193822891025449616434100410017771840410041641341320534160000205

19、3008404141181644582454100213382301102544961643410041002383128091216176019543111359761082300108913128093906791074570465300232359900139222561000300316135954125276510812500108256510534514100201210162516250125125081024161012136191027227208404100#121110987654321iiiiiiiiiiiiiiicc4384041010411264338251410021938228910254496164341004100177718404100416413413205341600002053008404141181644582