龙格-库塔方法基本原理3

在连续系统的仿真中，主要的计算工作是求解一阶微分方程y'=f(x,y)y（x0）=y0解析法只能用来求解一些特殊类型的方程，实际仿真问题中归结出来的微分方程主要靠数值解法。由于实际运算只能完成有限项或有限步运算，因此要将有些需用极限或无穷过程进行的运算有限化，对无穷过程进行截断，这样产生的误差成为截断误差。根据实际情况建立的数学模型往往难以求解。通常需要通过近似替代，将所求解的数学模型简化为易求解的数值计算问题后再进行求解。数学模型的理论解与数值计算问题的精确解之间的误差称为截断误差。这是计算方法本身带来的误差，所以也成为方法误差。2020/2/103得到高精度方法的一个直接想法是利用Taylor展开假设式y'=f(x,y)(a≤x≤b)中的f(x,y)充分光滑,将y(xi+1)在xi点作Taylor展开，若取右端不同的有限项作为y(xi+1)的近似值，就可得到计算y(xi+1)的各种不同截断误差的数值公式。例如：取前两项可得到)()()()(21hOxyhxyxyiii)(),()())(,()(22hOyxhfyhOxyxhfxyiiiiii3.2龙格－库塔方法§2020/2/104)(21!!2PiPiiiiyPhyhyhyy其中ffffffffffffyfffyxfyfyyyxyyxyxxxyxiyxxiiii22)(2)(),(,P阶泰勒方法若取前三项，可得到截断误差为O(h3)的公式)()(2)()()(321hOxyhxyhxyxyiiii)(),(),(),(2),(32hOyxfyxfyxfhyxhfyiiyiiiixiii类似地，若取前P+1项作为y(xi+1)的近似值，便得到2020/2/105显然p=1时,yi+1=yi+hf(xi,yi)它即为我们熟悉的Euler方法。当p≥2时,要利用泰勒方法就需要计算f(x,y)的高阶微商。这个计算量是很大的，尤其当f(x,y)较复杂时，其高阶导数会很复杂。因此，利用泰勒公式构造高阶公式是不实用的。但是泰勒级数展开法的基本思想是许多数值方法的基础。R-K方法不是直接使用Taylor级数,而是利用它的思想2020/2/106龙格-库塔(R-K)法的基本思想Euler公式可改写成),(1iiiiyxhfKKyy则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同,即局部截断误差为O(h2)。Runge-Kutta方法是一种高精度的单步法,简称R-K法2020/2/107同理，改进Euler公式可改写成),(),(2121121211KyhxhfKyxhfKKKyyiiiiii上述两组公式在形式上共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进欧拉法需计算两次f(x,y)的值，为二阶方法。局部截断误差为O(h3)2020/2/108于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在[xi,xi+1]这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格—库塔（Runge-Kutta）法的基本思想。),(),(),(11,1112122122111ppppipipiiiippiiKbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyy一般龙格－库塔方法的形式为2020/2/109其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作Tailor展开，通过相同项的系数确定参数。称为P阶龙格－库塔方法。10Runge-Kutta方法的推导思想0)(),(yaybxayxfy对于常微分方程的初值问题的解y=y(x)，在区间[xi,xi+1]上使用微分中值定理，有)()()(1iiiyhxyxy即))(()()(11iiiiixxyxyxy),(1iiixx其中2020/2/1011上的平均斜率在区间可以认为是],[)(1iixxxyyKKxyxyii)()(1引入记号)(,)(iiiyhfyhKKxxxyii上平均斜率的近似值间在区出只要使用适当的方法求],[)(1就可得到相应的Runge-Kutta方法ix1ixxy)(xyyKyyii1K2020/2/1012ix1ixxy)(xyy如下图Kxxxyxxyiii上的平均斜率在处的斜率作为在如果以],[)()(1即则上式化为),(1iiiiyxhfyy)(ixyhK)](,[iixyxhf),(iiyxhf即Euler方法Euler方法也称为一阶Runge-Kutta方法KK2020/2/10二阶龙格—库塔法在[xi,xi+1]上取两点xi和xi+a2=xi+a2h,以该两点处的斜率值K1和K2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即2211KcKcK式中:K1为xi点处的切线斜率值K1=hf(xi,yi)=hy'(xi)K2为xi+a2h点处的切线斜率值,比照改进的欧拉法,将xi+a2视为xi+1，即可得),(12122KbyhaxhfKii2020/2/1013确定系数c1、c2、a2、b21，可得到有2阶精度的算法格式2020/2/1014因此Kxyxyii)()(1)()(2211KcKcxyi将y(xi+1)在x=xi处进行Taylor展开：)()(!2)()()(321hOxyhxyhxyxyiiii)(!2),()(32hOfffhyxhfxyyxiii将在x=xi处进行Taylor展开：),()(121222KbyhaxhfhaxyhKiii2020/2/1015)(),(212122hOfKbfhayxfhKyxiiK1=hf(xi,yi))(),(3212hOfhfbfhayxfhyxii)(22111KcKcyyii)(),(32122hOfhfbfhayxfhcyxii),()(1iiiyxhfcxy),()()(21iiiyxhfccxy)(32221222hOffhcbfhcayx2020/2/101621,21,12212221cbcacc这里有4个未知数，3个方程。存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。令11)(iiyxy对应项的系数相等，得到2020/2/1017注意到，就是二阶龙格-库塔公式，也就是改进的欧拉法。21,121212ccba),(),(21121211KyhxhfKyxhfKKKyyiiiiiiix1ixxy)(xyyK1K2K因此，凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。若取，就是另一种形式的二阶龙格-库塔公式。1,0,2121212ccba2020/2/1018)21,21(),(12121KyhxhfKyxhfKKyyiiiiii此计算公式称为变形的二阶龙格—库塔法。式中为区间的中点。也称中点公式。hxi211,iixxQ:为获得更高的精度，应该如何进一步推广？2020/2/1019二级R-K方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2,b21，可使每步计算两次函数值的二阶R-K方法达到二阶精度。能否在计算函数值次数不变的情况下,通过选择不同的参数值,使得二阶R-K方法的精度再提高呢?答案是否定的！无论四个参数怎样选择,都不能使公式的局部截断误差提高到三阶。这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶。若要获得更高阶得数值方法,就必须增加计算函数值的次数。),(),(),(232131331212213322111KbKbyhaxhfKKbyhaxhfKyxhfKKcKcKcyyiiiiiiii三阶龙格—库塔法2020/2/1020为进一步提高精度，在区间[xi,xi+1]上除两点xi和xi+a2=xi+a2h,以外，再增加一点xi+a3=xi+a3h，用这三点处的斜率值K1、K2和K3的加权平均得出平均斜率K*的近似值K，这时计算格式具有形式:ix3aix2aixxy)(xyyK1K2K3K2020/2/1021同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进行Taylor展开，使局部截断误差达到O(h4)，使对应项的系数相等，得到系数方程组：61,6131)(31)(,3121)(,21133221232232323122213231332212323222332312213322321cbbcbacbbcbbbcacbacacacbbcbcacaccc参数的选择不唯一，从而构成一类不同的三阶R-K公式，下面给出一种常用的三阶R-K公式，形似simpson公式：)2,()21,2(),()4(612131213211KKyhxhfKKyhxhfKyxhfKKKKyyiiiiiiii2020/2/10222020/2/1023四阶(经典)龙格—库塔法如果需要再提高精度，用类似上述的处理方法，只需在区间[xi,xi+1]上用四个点处的斜率加权平均作为平均斜率K*的近似值，构成一系列四阶龙格—库塔公式。具有四阶精度，即局部截断误差是O(h5)。推导过程与前面类似，由于过程复杂，这里从略，只介绍最常用的一种四阶经典龙格—库塔公式。2020/2/1024K1=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+c4K42020/2/1025类似于前面的讨论，把K2、K3、K4分别在xi点展成h的幂级数，代入线性组合式中，将得到的公式与y(xi+1)在xi点上的泰勒展开式比较，使其两式右端直到h4的系数相等，经过较复杂的解方程过程便可得到关于ci，ai，bij的一组特解a2=a3=b21=b32=1/2b31=b41=b42=0a4=b43=1c1=c4=1/6c2=c3=1/326四阶(经典)Runge-Kutta方法),()21,2()21,2(),()22(61342312143211KyhxhfKKyhxhfKKyhxhfKyxhfKKKKKyyiiiiiiiiii2020/2/1027例1.使用高阶R-K方法计算初值问题1)0(5.002yxyy.1.0h取解:(1)使用三阶R-K方法时0i1.0201hyK1103.0