93-龙格-库塔方法

19.3龙格-库塔方法29.3.1显式龙格-库塔法的一般形式上节给出了显式单步法的表达式),,,(1hyxhyynnnn其局部截断误差为),(),,()()(11pnhOhyxhxyhxyT对欧拉法，即方法为阶.)(21hOTn1p))].,(,(),([21nnnnnnnnyxhfyhxfyxfhyy（3.1）若用改进欧拉法，它可表示为3此时增量函数))].,(,(),([21),,(nnnnnnnnyxhfyhxfyxfhyx（3.2）与欧拉法的相比，增加了计算一个右函数的值，可望.),(),,(nnnnyxfhyxf2p若要使得到的公式阶数更大，就必须包含更多的值.pf,))(,()()(1d1nnxxnnxxyxfxyxy（3.3）从方程等价的积分形式（2.4），即),(yxfy4若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，必然要增加求积节点.为此可将（3.3）的右端用求积公式表示为11nnxxriininihxyhxfchxxyxf.))(,())(,(d点数越多，精度越高，r上式右端相当于增量函数，),,(hyx为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为),,,(1hyxhyynnnn（3.4）其中5,),,(1riiinnKchyx),,(1nnyxfK),,(11ijjijniniKhyhxfK（3.5）,,,2ri这里均为常数.ijiic,,(3.4)和(3.5)称为级显式龙格-库塔（Runge-Kutta）法，r简称R-K方法.当时，就是欧拉法，),(),,(,1nnnnyxfhyxr此时方法的阶为.1p当时，改进欧拉法(3.1)，(3.2)也是其中的一种.2r6下面将证明阶.2p要使公式(3.4),(3.5)具有更高的阶,就要增加点数.pr下面就推导R-K方法.2r,),,(1riiinnKchyx),,(1nnyxfK),(11ijjijniniKhyhxfK),,,(1hyxhyynnnn（3.4）（3.5）,,,2ri79.3.2二阶显式R-K方法对的R-K方法，计算公式如下2r).,(),,(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6）这里均为待定常数.21221,,,cc希望适当选取这些系数，使公式阶数尽量高.p根据局部截断误差的定义，（3.6）的局部截断误差为)()(11nnnxyxyT（3.7）)],,(),([21221nnnnnhfyhxfcyxfch8这里.),(),(nnnnnyxffxyy为得到的阶，要将上式各项在处做泰勒展开，1nTp),(nnyx由于是二元函数，故要用到二元泰勒展开，),(yxf),(!32)(4321hOyhyhyhyxynnnnn其中,),(nnnnfyxfy各项展开式为,),(),())(,(nnnynnxnnnfyxfyxfxyxfxydd),(),(2),(2nnyynnnxynnnxxnyxffyxffyxfy（3.8）)];,(),()[,(nnynnnxnnyyxffyxfyxf9),(212nnnfhyhxf将以上结果代入局部截断误差公式则有)(]),(),(([]),(),([232122121hOhfyxfhyxffcfchfyxfyxfhfhTnnnynnxnnnnnynnxnn).(),(),(2212hOfhyxfhyxffnnnynnxn).(),()21(),()21()1(3221222221hOhfyxfchyxfchfccnnnynnxn要使公式（3.6）具有阶，必须使2p10,0121cc即.1,21,212121222cccc（3.9）的解是不惟一的.令，则得02ac.21,12121aac这样得到的公式称为二阶R-K方法，如取，则2/1a.1,2/121221cc这就是改进欧拉法（3.1）.,02122c（3.9）,021212c11若取,1a则.2/1,1,021221cc,21hKyynn称为中点公式，（3.10）也可表示为)).,(2,2(1nnnnnnyxfhyhxfhyy得计算公式),,(1nnyxfK（3.10）).2,2(12KhyhxfKnn相当于数值积分的中矩形公式.的R-K公式(3.6)的局部误差不可能提高到.2r)(4hO12把多展开一项，从（3.8）的的项是不能通过选择参数消掉的.2Knyffffyxy实际上要使的项为零，需增加3个方程，要确定4个参数，这是不可能的.3h21221,,,cc故的显式R-K方法的阶只能是，而不能得到三阶公式.2r2p139.3.3三阶与四阶显式R-K方法要得到三阶显式R-K方法，必须.3r),,(),,(),,(),(232131321212213322111hKhKyhxfKhKyhxfKyxfKKcKcKchyynnnnnnnn（3.11）其中及均为待定参数.321,,ccc32313212,,,,].[)()(33221111KcKcKchxyxyTnnn此时（3.4）,（3.5）的公式表示为公式（3.11）的局部截断误差为14只要将按二元函数泰勒展开，使，32,KK)(41hOTn可得待定参数满足方程.61,31,21,,,13223233222332232313212321cccccccc（3.12）15这是8个未知数6个方程的方程组，解也不是惟一的.所以这是一簇公式.满足条件（3.12）的公式（3.11）统称为三阶R-K公式.一个常见的公式为).2,(),2,2(),,(),4(62131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn此公式称为库塔三阶方法.16继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为.)(5hO四阶龙格-库塔方法的每一步需要计算4次函数值，f).,(),2,2(),2,2(),,(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13）17例3.1)0(),10(2yxyxyy解设取步长，从到用四阶龙格-2.0h0x1x库塔方法求解初值问题这里，经典的四阶龙格-库塔公式为yxyyxf2),(),22(643211KKKKhyynn,21nnnyxyK18,222223KhyhxKhyKnnn,222112KhyhxKhyKnnn.)(2334hKyhxhKyKnnn表9-3列出了计算结果，同时列出了相应的精确解.7321.17321.10.16125.16125.18.04832.14833.16.03416.13417.14.01832.11832.12.0)(nnnxyyx计算结果3表9比较例3和例2的计算结果，显然龙格-库塔方法的精度高.19但由于这里放大了步长，所以表9-3和表9-2所耗费的计算量几乎相同.)2.0(h龙格-库塔方法的推导基于泰勒展开方法，因而它要求所求的解具有较好的光滑性质.反之，如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法.四阶龙格-库塔方法每一步要4次计算函数,其计算量f比每一步只要2次计算函数的二阶龙格-库塔方法----f改进的欧拉方法大一倍.209.3.4变步长的龙格-库塔方法单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了.步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累.因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题.在选择步长时，需要考虑两个问题：1.怎样衡量和检验计算结果的精度？212.如何依据所获得的精度处理步长？考察经典的四阶龙格-库塔公式:).,(),2,2(),2,2(),,(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13）从节点出发，先以为步长求出一个近似值，nxh)(1hny22,)(5)(11chyxyhnn（3.14）然后将步长折半，即取为步长从跨两步到，2hnx1nx再求得一个近似值，每跨一步的截断误差是)2(1hny,52hc.)()(521122hcyxyhnn（3.15）比较（3.14）式和（3.15）式我们看到，步长折半后，由于公式的局部截断误差为，故有)(5hO因此有误差大约减少到.即有16123.161)()()(11)2(11hnnhnnyxyyxy由此易得下列事后估计式].[151)()(1)2(1)2(11hnhnhnnyyyxy这样，可以通过检查步长、折半前后两次计算结果的偏差)(1)2(1hnhnyy来判定所选的步长是否合适.具体地说，将区分以下两种情况处理：241.对于给定的精度，如果，反复将步长折半进行计算，直至为止.这时取最终得到的作为结果；)2(1hny2.如果，反复将步长加倍，直到为止，这种通过加倍或折半处理步长的方法称为变步长方法.这时再将步长折半一次，就得到所要的结果.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.