1、3.3 龙格龙格-库塔库塔(Runge-Kutta)方法方法3.3.1 龙格龙格-库塔方法的设计思想库塔方法的设计思想3.3.2 二阶龙格二阶龙格-库塔方法库塔方法3.3.3 三阶龙格三阶龙格-库塔方法库塔方法3.3.4 四阶龙格四阶龙格-库塔方法库塔方法3.3.5 变步长的龙格变步长的龙格-库塔方法库塔方法3.3.1 龙格龙格-库塔方法的设计库塔方法的设计思想思想 对对一阶常微分方程初值问题一阶常微分方程初值问题 00( , )()yf x yy xy ， 欧拉公式欧拉公式1(,)nnnnyyhf xy， 写成写成 111(,)nnnnyyhkkf xy， 调用调用 1 次函数次函数。有。有
2、 1 阶精度阶精度。 3.3 龙格龙格- -库塔库塔（Runge-Kutta)法法对对一阶常微分方程初值问题一阶常微分方程初值问题 00( , )()yf x yy xy 改进欧拉法，预报改进欧拉法，预报-校正公式校正公式 n+111n+1(,) (,)(,)2nnnnnnnnyyh f xyhyyf xyf xy 写成写成 112121()2(,),(,)nnnnnnhyykkkf xykf xh yhk 调用调用 2 次函数次函数。有。有 2 阶精度阶精度。 龙格龙格- -库塔库塔方方法的设计思想法的设计思想将改进欧拉法（预报将改进欧拉法（预报-校正公式）改写成平均化的形式校正公式）改写成
3、平均化的形式 1121211()2(,),(,)nnnnnnhyykkkf xykf xyhk 用用 nx与与 1nx两个点的斜率值两个点的斜率值 1k与与 2k取算术平均作为取算术平均作为平均斜率平均斜率*121()2kkk，而，而 1nx处的斜率值处的斜率值 2k 则通过已知信息则通过已知信息 ny 来预测。来预测。 龙格龙格- -库塔库塔方方法的设计思想法的设计思想利用利用( , )f x y在某些点的在某些点的函数函数值的线性组合值的线性组合构造公构造公式，得到式，得到1()ny x的近似值的近似值1ny，增加调用，增加调用( , )f x y的次数，的次数，可提高精度的阶数。可提高精
4、度的阶数。 设法在设法在1,nnxx这一步内这一步内多预报几个点多预报几个点的斜率， 然后的斜率， 然后将其将其加权平均加权平均，作为平均斜率，则可构造出具有更高精，作为平均斜率，则可构造出具有更高精度的计算公式。度的计算公式。 龙格龙格- -库塔库塔方方法的设计思想法的设计思想3.3.2 二阶龙格二阶龙格-库塔库塔方方法法 考察考察1,nnxx内一点内一点npnxxph，其中，其中01p 用用,nnpxx点的斜率点的斜率12,k k加权平均，得平均斜率，即加权平均，得平均斜率，即 112(1)nnyyhkk 式中式中为待定系数。取为待定系数。取1(,)nnkf xy，2(,)npnpkf x
5、y，先用欧拉法，先用欧拉法提供提供()npy x的预报值的预报值npy，有，有1npnyyphk，因此，因此 21(,)(,)npnpnnkf xyf xph yphk 这样有计算公式：这样有计算公式： 112121(1)(,),(,)nnnnnnyyhkkkf xykf xph yphk 其中含有两个待定参数，选取合适的值，使公式有高的精度。其中含有两个待定参数，选取合适的值，使公式有高的精度。 仍假定仍假定()nnyy x，则有，则有：1(,)()nnnkf xyy x， 并并将将 k2在在(xn,yn)处处进行进行二元二元泰勒泰勒展开，展开，有：有： 21(,)nnkf xph yphk
6、 2(,)(,)(,)(,)()nnxnnnnynnf xyph fxyf xyfxyO h 2()()()nny xphyxO h 这时，有这时，有 )()()( )()(3n2nn21n1nhOxyphxhyxykk1hyy 根据根据局部局部截断截断误差误差的的定义：定义： 由由 231()()()()nnnnyy xhy xph yxO h 和和1()ny x的泰的泰勒展开式勒展开式 2311()()()()()2nnnny xy xhy xh yxO h 比较，比较，要要使使局部局部截断截断误差误差 )()(31n1nhOyxy 从而有构造从而有构造二阶二阶龙格龙格-库塔方法所需条件库
7、塔方法所需条件为：为： 12p 2/1p为不定方程，有无穷多个解为不定方程，有无穷多个解，对应，对应一族一族二阶龙格二阶龙格-库塔格式库塔格式。 ),(),()1(121211phkyphxfkyxfkkkhyynnnnnn 如：如： （1）取）取11,2p，得，得 112121()2(,)(,)nnnnnnhyykkkf xykf xh yhk 这就是这就是改进的欧拉改进的欧拉格格式式。 （2）取）取1,12p，得，得 12121(,)(,)22nnnnnnyyhkkf xyhhkf xyk 该公式称为该公式称为变形的欧拉格式变形的欧拉格式（或（或中点中点格格式式） 。 （3）取）取43,3
8、2p，得，得 ),(),()(1nn2nn121n1nhk32yh32xfkyxfkk3kh41yy 这就是这就是休恩休恩(Heun)公公式式。 3.3.3 三三阶龙格阶龙格-库塔库塔方方法法 为进一步提高精度，在为进一步提高精度，在1,nnxx内内除除点点pnx外再增加一点外再增加一点qhxxnqn，其中，其中1 qp，用用三个点三个点qnpnnxxx,的斜率的斜率321,kkk加权平均，得平均斜率，即加权平均，得平均斜率，即 )1(3211kkkhyynn 式中式中,为待定系数。为待定系数。 取取1(,)nnkf xy，),(12phkyxfknpn，用，用21,kk的加权平均的加权平均得
9、出得出,qnnxx上的平均斜率，从而得到上的平均斜率，从而得到)(qnxy的预报值的预报值qny，有，有)1(21kkqhyynqn 再通过计算再通过计算 f 得到得到),(3qnqnyxfk，运用泰勒展开方法选择参数，运用泰勒展开方法选择参数,qp，使，使)1(3211kkkhyynn具有具有三阶三阶精精度，得到一族度，得到一族三阶龙格三阶龙格-库塔格式库塔格式。 如当如当2,61,64, 1,21qp时，其中的一种为：时，其中的一种为： )2(,()2,(),()4(62113121321121kkhyxfkkhyxfkyxfkkkkhyynnnnnnnn 三阶龙格三阶龙格-库塔格式库塔格
10、式 如当如当1,43, 0,32,31qp时，其中的一种为：时，其中的一种为： )32,()31,(),()3(4231213113231hkyxfkhkyxfkyxfkkkhyynnnnnnnn 称为称为三阶休恩格式三阶休恩格式。 3.3.4 3.3.4 四阶龙格四阶龙格- -库塔方法库塔方法仿照上述的讨论，用四个不同点上的函数值的线性组合仿照上述的讨论，用四个不同点上的函数值的线性组合, ,可导出一族可导出一族四阶龙格库塔格式四阶龙格库塔格式: :四阶龙格库塔公式每步要四阶龙格库塔公式每步要四次四次计算函数值，具有计算函数值，具有四阶四阶精精度，局部截断误差是度，局部截断误差是O(h5)
11、. ),(),(),(),()(3332321313422212123111121443322111hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn四阶龙格库塔格式四阶龙格库塔格式下列下列经典格式经典格式是其中常用的一种是其中常用的一种: :经典龙格库塔公式每步要经典龙格库塔公式每步要四次四次计算函数值，具有计算函数值，具有四阶四阶精度，局部截断误差是精度，局部截断误差是O(h5) . ),()2,21()2,21(),()22(6342312143211hkyhxfkkhyhxfkkhyhxfkyxfkkkkkhyynnnnnnnnnn四阶经典
12、龙格库塔法的框图四阶经典龙格库塔法的框图 112341213243(22)6(,),(,)22(,),(,)22nnnnnnnnnnhyykkkkhhkf xykf xykhhkf xykkf xh yhk框图框图3.3.5 变步长的龙格变步长的龙格- -库塔方法库塔方法对每一步对每一步: : 步长越小步长越小, ,截断误差越小截断误差越小; ;但步长越小但步长越小, ,完成的步数越多完成的步数越多, ,计算量越大计算量越大, ,可能造成舍入误差的严重积累可能造成舍入误差的严重积累. .选择步长的原则选择步长的原则: :1)1) 如何衡量和检验计算结果的精度如何衡量和检验计算结果的精度; ;2
13、)2) 如何依据所判定的精度来调整步长如何依据所判定的精度来调整步长. .四阶经典龙格四阶经典龙格- -库塔方法步长折半前后的库塔方法步长折半前后的局部截断误差局部截断误差误差事后估计)(1)2(1)2(11)(11)2(115)2(115)(11151)(161)()(22)()(hnhnhnnhnnhnnhnnhnnyyyxyyxyyxyhCyxyChyxy可以通过检查步长折半前后两次计算结果的偏差可以通过检查步长折半前后两次计算结果的偏差)(1)2(1hnhnyy来判断所选取的步长是否合适来判断所选取的步长是否合适.1) 对于给定的精度对于给定的精度, ,如果如果 ,则反复将步长折则反复
14、将步长折半进行计算半进行计算,直到直到 为止为止,这时取步长折半后这时取步长折半后的新值作为结果的新值作为结果;2) 如果如果为止为止,这时取步长加倍前的老值作为结果这时取步长加倍前的老值作为结果.变步长方法表面上看为了选择步长变步长方法表面上看为了选择步长,每一步的计每一步的计算量增加了算量增加了,但总体考虑往往是合算的但总体考虑往往是合算的.写出经典龙格-库塔方法求解初值问题 1)0(1yxyy 的计算公式，并取步长2 . 0h，计算)4 . 0(y的近似值。 解解 由题意知1,xyyxf，所以有例例1 11234121132222,16(,)()()122220.20.2110.90.9
15、122(,)()()122220.20.20.90.911022nnnnnnnnnnnnnnnnnnnnnnnnhyykkkkkfxyyxhhhhkf xykykxyyxxyxhhhhkf xykykxyyxx 433.910.911,10.20.910.9110.210.8180.8181nnnnnnnnnnnnyxkfxh yhkyhkxhyyxxyx 代入，有 2 . 0181267. 0818733. 022643211nnnnxykkkkhyy 由于00 x，1)0(0 yy，取步长2 . 0h，取4 . 02, 2 . 00201hxxhxx，有 018733. 1)2 . 0(1
16、 yy, 070324. 1)4 . 0(2 yy 例例2用经典用经典Runge-Kutta方法求解方法求解取取 h = 0.1计算结果见下表计算结果见下表 , 1)0(, 10, 1 yxxyy解解: Euler方法方法 预估校正方法预估校正方法 经典经典R-K方法方法 0.01.0000000.01.0000000.01.000000000.00.11.0000004.81031.0050001.61041.004837508.21080.21.0100008.71031.0190252.91041.018730901.51070.31.0290001.21021.0412184.0104
17、1.040818422.01070.41.0561001.41021.0708004.81041.070320292.41070.51.0904901.61021.1070765.51041.106530932.71070.61.1314411.71021.1494045.91041.148811932.91070.71.1782971.81021.1972106.21041.196585623.11070.81.2304671.91021.2499756.51041.249329293.21070.91.2874201.91021.3072286.61041.306569993.31071.
18、01.3486781.91021.3685146.61041.36789773.3107kykykykxkkyxy )(kkyxy )(kkyxy )(整体截断误差整体截断误差（四级（四级4阶）阶）（二级（二级2阶）阶）（一级（一级1阶）阶） 说明说明:（1）通过比较可以发现在相同步长下经典通过比较可以发现在相同步长下经典Runge-Kutta （2）经典）经典Runge-Kutta方法的计算量是方法的计算量是Euler方法的方法的4倍倍,预估预估注注:（2）中）中不与梯形方法比较，因梯形方法与这三种不同类。不与梯形方法比较，因梯形方法与这三种不同类。方法的结果比方法的结果比Euler方法、梯形方法、预估校正方法、梯形方法、预估校正Euler方法好的方法好的 多，即更准确，误差较小。多，即更准确，误差较小。4h2h的步长取的步长取 ，预估校正，预估校正Euler方法取方法取 ，它们的计算量大致，它们的计算量大致校正方法的校正方法的2倍。倍。相同，此时，经典相同，此时，经典Runge-Kutta方法仍比方法仍比Euler方法、预估校正方法、预估校正 若经典若经典Runge-Kutta方法步长取方法步长取h，Euler方法方法Euler方法好的多。方法好的多。
