计算方法龙格库塔方法PPT课件

1、2022-7-61)(21! 2PiPiiiiyPhyhyhyy 其中ffffffffff ffyf ffyxfyfyyyxyyxyxxxyxiyxxiiii22)(2)(),(, P阶泰勒方法若取前三项，可得到截断误差为O(h3)的公式 )()(2)()()(321hOxyhxyhxyxyiiii )(),(),(),(2),(32hOyxfyxfyxfhyxhfyiiyiiiixiii类似地，若取前P+1项作为y(xi+1)的近似值，便得到第1页/共32页2022-7-62显然p=1时, y i+1=y i+hf(xi,y i)它即为我们熟悉的Euler方法。当p2时,要利用泰勒方法就需要

2、计算f(x,y)的高阶微商。这个计算量是很大的，尤其当f(x,y)较复杂时，其高阶导数会很复杂。因此，利用泰勒公式构造高阶公式是不实用的。但是泰勒级数展开法的基本思想是许多数值方法的基础。R-K方法不是直接使用Taylor级数,而是利用它的思想第2页/共32页2022-7-639.4.1 9.4.1 龙格龙格- -库塔库塔( (R-K) )法的基本思想法的基本思想Euler公式可改写成 ),(1iiiiyxhfKKyy则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同,即局部截断误差局部截断误差为O(h2)。Runge-Kutta 方法是一种高精度的单步法方法是一种高精度的

3、单步法, ,简称简称R-K法法第3页/共32页2022-7-64同理，改进Euler公式可改写成 ),(),(2121121211KyhxhfKyxhfKKKyyiiiiii 上述两组公式在形式上共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进欧拉法需计算两次f(x,y)的值，为二阶方法。局部截断误差局部截断误差为O(h3)第4页/共32页2022-7-65 于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公

4、式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在xi,xi+1这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格龙格库塔库塔（Runge-Kutta）法的基本思想法的基本思想。 第5页/共32页),(),(),(11,1112122122111ppppipipiiiippiiKbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyy一般龙格库塔方法的形式为2022-7-66其中ai,bij,ci为待

5、定参数，要求上式yi+1在点(xi,yi)处作Tailor展开，通过相同项的系数确定参数。称为P阶龙格库塔方法。第6页/共32页7Runge-Kutta方法的推导思想0)(),(yaybxayxfy对于常微分方程的初值问题的解y=y(x)，在区间xi, xi+1上使用微分中值定理，有)()()(1iiiyhxyxy即)()()(11iiiiixxyxyxy),(1iiixx其中2022-7-6第7页/共32页8上的平均斜率在区间可以认为是,)(1iixxxyyKKxyxyii)()(1引入记号)(,)(iiiyhfyhKKxxxyii上平均斜率的近似值间在区出只要使用适当的方法求,)(1就可得

6、到相应的Runge-Kutta方法ix1ixxy)(xyy Kyyii1K2022-7-6第8页/共32页9ix1ixxy)(xyy 如下图Kxxxyxxyiii上的平均斜率在处的斜率作为在如果以,)()(1即则上式化为),(1iiiiyxhfyy)(ixyhK)(,iixyxhf),(iiyxhf即Euler方法Euler方法也称为一阶一阶Runge-Kutta方法方法KK2022-7-6第9页/共32页9.4.2 9.4.2 二阶龙格二阶龙格库塔法库塔法 在xi, xi+1上取两点xi和xi+a2= xi +a2h,以该两点处的斜率值K1和K2的加权平均(或称为线性组合)来求取平均斜率k*

7、的近似值K，即 2211KcKcK式中:K1 1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy(xi) K2 2为xi +a2h点处的切线斜率值,比照改进的欧拉法,将xi+a2视为xi+1，即可得 ),(12122KbyhaxhfKii2022-7-610确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式第10页/共32页2022-7-611因此 Kxyxyii)()(1)()(2211KcKcxyi将y(xi+1)在x=xi处进行Taylor展开： )()(! 2)()()(321hOxyhxyhxyxyiiii )(! 2),()(32hOf ffhyxhfxyy

8、xiii将 在x=xi处进行Taylor展开： ),()(121222KbyhaxhfhaxyhKiii第11页/共32页2022-7-612)(),(212122hOfKbf hayxfhKyxiiK1 =hf(xi, yi)(),(3212hOfhfbf hayxfhyxii)(22111KcKcyyii)(),(32122hOfhfbf hayxfhcyxii),()(1iiiyxhfcxy),()()(21iiiyxhfccxy)(32221222hOf fhcbfhcayx第12页/共32页2022-7-61321,21,12212221cbcacc这里有 4 个未知数，3 个方程。

9、存在无穷多个解无穷多个解。所有满足上式的格式统称为2阶阶龙格龙格 - - 库塔格式库塔格式。令 11)(iiyxy对应项的系数相等，得到 第13页/共32页2022-7-614注意到，注意到， 就是二阶就是二阶龙格龙格 - - 库塔库塔公式，也就是公式，也就是改进的欧拉法。改进的欧拉法。 21, 121212ccba),(),(21121211KyhxhfKyxhfKKKyyiiiiiiix1ixxy)(xyyK1K2K 因此，凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格库塔格式。因此改进的欧拉格式是众多的二阶龙格库塔法中的一种特殊格式。 第14页/共32页若取若取 ，就是另一

10、种形式的二，就是另一种形式的二阶阶龙格龙格 - - 库塔公式库塔公式。 1, 0,2121212ccba2022-7-615)21,21(),(12121KyhxhfKyxhfKKyyiiiiii此计算公式称为变形的二阶龙格库塔法。式中 为区间 的中点。也称中点公式也称中点公式。 hxi211,iixxQ:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？第15页/共32页2022-7-616 二级R-K方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2, b21，可使每步计算两次函数值的二阶R-K方法达到二阶精度。能否

11、在计算函数值次数不变的情况下,通过选择不同的参数值,使得二阶R-K方法的精度再提高呢? 答案是否定的！无论四个参数怎样选择,都不能使公式的局部截断误差提高到三阶。 这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶。若要获得更高阶得数值方法若要获得更高阶得数值方法, ,就必须增加计算函就必须增加计算函数值的次数。数值的次数。第16页/共32页),(),(),(232131331212213322111KbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyyiiiiiiii9.4.3 三阶龙格三阶龙格库塔法库塔法2022-7-617 为进一步提高精度，在区间xi, xi+1上除两点

12、xi和xi+a2= xi +a2h,以外，再增加一点xi+a3= xi +a3h ，用这三点处的斜率值K1、K2和K3的加权平均得出平均斜率K*的近似值K，这时计算格式具有形式: ix3aix2aixxy)(xyy K1K2K3K第17页/共32页2022-7-618 同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进行Taylor展开，使局部截断误差达到O(h4)，使对应项的系数相等，得到系数方程组：61,6131)(31)(,3121)(,21133221232232323122213231332212323222332312213322321cbbcbacbbcbbbcacbac

13、acacbbcbcacaccc第18页/共32页参数的选择不唯一，从而构成一类不同的三阶R-K公式，下面给出一种常用的三阶R-K公式，形似simpson公式：)2,()21,2(),()4(612131213211KKyhxhfKKyhxhfKyxhfKKKKyyiiiiiiii2022-7-619第19页/共32页2022-7-6209.4.4 四阶四阶( (经典经典) )龙格龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，只需在区间xi,xi+1上用四个点处的斜率加权平均作为平均斜率K*的近似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部截断误差是O(h5)。 推导过程与

14、前面类似，由于过程复杂，这里从略，只介绍最常用的一种四阶经典龙格四阶经典龙格库塔公式库塔公式。 第20页/共32页2022-7-621 K1=hf (xi, yi) K2=hf (xi+a2h, yi+b21K1) K3=hf (xi+a3h, yi+b31K1+b32K2) K4=hf (xi+a4h, yi+b41K1+b42K2+b43K3) 其中c1、c2、c3、c4、a2、a3、a4、b21、b31、b32、b41、b42、b43均为待定系数。这里K1、K2、K3、K4为四个不同点上的函数值，分别设其为 设yi+1=yi+c1K1+c2K2+c3K3+c4K4第21页/共32页202

15、2-7-622 类似于前面的讨论，把K2、K3、K4分别在xi点展成h的幂级数，代入线性组合式中，将得到的公式与y(xi+1)在xi点上的泰勒展开式比较，使其两式右端直到h4的系数相等，经过较复杂的解方程过程便可得到关于ci，ai，bij的一组特解 a2=a3=b21=b32=1/2 b31=b41=b42=0 a4=b43=1 c1=c4=1/6 c2=c3=1/3 第22页/共32页23 四阶(经典)Runge-Kutta方法),()21,2()21,2(),()22(61342312143211KyhxhfKKyhxhfKKyhxhfKyxhfKKKKKyyiiiiiiiiii2022-

16、7-6第23页/共32页24例1. 使用高阶R-K方法计算初值问题1)0(5 . 002yxyy. 1 . 0h取解:(1) 使用三阶R-K方法时0i1 . 0201 hyK1103. 0)21(2102KyhK1256. 0)2(22103KKyhK1111. 1)4(6132101KKKyy2022-7-6第24页/共32页25其余结果如下:(2) 如果使用四阶R-K方法时0i1 . 0201 hyK1103. 0)21(2102KyhK i xi k1 k2 k3 yi 1.0000 0.1000 0.1000 0.1103 0.1256 1.1111 2.0000 0.2000 0.1

17、235 0.1376 0.1595 1.2499 3.0000 0.3000 0.1562 0.1764 0.2092 1.4284 4.0000 0.4000 0.2040 0.2342 0.2866 1.6664 5.0000 0.5000 0.2777 0.3259 0.4163 1.99932022-7-6第25页/共32页261113. 0)21(2203KyhK1235. 0)(2304KyhK)22(61 . 0432101KKKKyy1111. 1其余结果如下: i xi k1 k2 k3 k4 yi 1.0000 0.1000 0.1000 0.1103 0.1113 0.1

18、235 1.1111 2.0000 0.2000 0.1235 0.1376 0.1392 0.1563 1.2500 3.0000 0.3000 0.1562 0.1764 0.1791 0.2042 1.4286 4.0000 0.4000 0.2040 0.2342 0.2389 0.2781 1.6667 5.0000 0.5000 0.2777 0.3259 0.3348 0.4006 2.00002022-7-6第26页/共32页2022-7-627由上节分析常微分方程数值解法稳定性问题的方法，可得到各阶Runge-Kutta公式的稳定性条件：二阶二阶121122hh与欧拉预估校正公式一致三阶三阶1)(! 31)(! 21132hhh四阶四阶1)(! 41)(! 31)(! 211432hhhh9.4.5 龙格库塔方法的稳定性条件龙格库塔方法的稳定性条件第27页/共32页2022-7-628 龙格库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改