第二节 龙格-库塔方法

1、.基本思想：基本思想：利用利用 在某些特殊点上的在某些特殊点上的函数值的函数值的线性线性组合组合来构造高阶单步法的来构造高阶单步法的平均斜率平均斜率。( , )f x y第二节第二节 龙格龙格- -库塔法库塔法什么叫平均斜率？什么叫平均斜率？对差商对差商 应用微分中值定理，有，应用微分中值定理，有，1()()iiy xy xh 1()()()iiiy xy xhy xh 利用微分方程利用微分方程 ，有，有( , )yf x y 1()()(, ()iiiiy xy xhf xh y xh 这里的这里的 称为称为平均斜率平均斜率。(, ()iif xh y xh .可将改进的欧拉格式改写成可将改

2、进的欧拉格式改写成112121112(),(,),(,).iiiiiiyKKKyhfxyKyhfxK 的算术平均值作为平均斜率。的算术平均值作为平均斜率。该公式可以看作是用该公式可以看作是用 和和 两个点处的斜率两个点处的斜率 和和1 iixx 12 KK由改进型欧拉公式我们可以猜想，如果在由改进型欧拉公式我们可以猜想，如果在1,iix x 内多预测几个点的斜率，再对他们进行加权平均，内多预测几个点的斜率，再对他们进行加权平均，可能得到精度更好的平均斜率！可能得到精度更好的平均斜率！.下面以下面以2阶龙格阶龙格- -库塔方法库塔方法为例来阐述这种思想为例来阐述这种思想考察区间考察区间 上的一点

3、上的一点 ，1,iix x 01, ipixxphp 用用 和和 的斜率的斜率 和和 的加权平均作为平均的加权平均作为平均 iipxx 12 KK斜率斜率 的近似值：的近似值：*K1122*KKK 即取即取11122()iiyyhKK 其中其中 和和 是待定常数。若取是待定常数。若取 ，则，则12 1(,)iiKf xy 问题在于如何确定问题在于如何确定 处的斜率处的斜率 和常数和常数 和和 。ipx 2K12 .仿照改进的欧拉方法，用欧拉方法预测仿照改进的欧拉方法，用欧拉方法预测 的值，的值，()ipy x 1ipiyyphK 并用它来估计斜率并用它来估计斜率 ：2K2(,)ipipKf x

4、y 于是得到如下形式的算法：于是得到如下形式的算法：111221211(),(,),(,).iiiiiiyyhKKKf xyKf xyphK 通过适当选取参数通过适当选取参数 和和 的值，使得公式具有的值，使得公式具有12, p 2阶精度！阶精度！.由泰勒公式展开，要使公式具有由泰勒公式展开，要使公式具有 2 阶精度阶精度，只需，只需12211 2,p方程组有方程组有无穷无穷多解：多解：二级二级方法有无穷多种方法有无穷多种常见的常见的3种二级方法：种二级方法： 中点法（修正的中点法（修正的Euler法法）1221012,cca 取取122(,(,)nnnnnnhhyyhf xyf xy 二阶二

5、阶龙格库塔方法龙格库塔方法122112,cca 取取12 (,)(,(,)nnnnnnnnhyyf xyf xh yhf xy 3()O h.三级方法：三级方法：N = 3 类似于类似于N = 2的推导方法，可得到的推导方法，可得到1231;ccc 223312;c ac a 22223313;c ac a 323216c a b 4()O h常见的常见的2种三阶方法：种三阶方法： 库塔库塔三阶方法三阶方法112346()nnhyykkk 1(,);nnkf xy 2122(,)nnhhkf xyk3122(,)nnkf xh yhkhk . 四级方法：四级方法：N = 45()O h局部截断

6、误差局部截断误差常见的常见的2种四阶方法：种四阶方法：经典经典龙格龙格- -库塔库塔方法方法11234226()nnhyykkkk 1(,)nnkf xy 2122(,)nnhhkf xyk3222(,)nnhhkf xyk43(,)nnkf xh yhk.解：解：201( )dyxydxyy 01( , )x 例例2：用用经典的经典的龙格龙格- -库塔库塔方法方法求解下列初值问题求解下列初值问题 。0 1 .h 经典的四阶经典的四阶龙格龙格- -库塔公式：库塔公式：11234226()nnhyykkkk 12;nnnxkyy 2112222()nnnhxhkykhyk 4332()nnnxh

7、kyhkyhk 3222222();nnnhxhkykhyk . 0.1 0.2 0.3 0.4 0.5 1.0954 1.1832 1.2649 1.3416 1.4142nxny 0.6 0.7 0.8 0.9 1.0 1.4832 1.5492 1.6125 1.6733 1.7321nxny同保留同保留5位的位的精确值精确值完全一致：完全一致： 0.1 0.2 0.3 0.4 0.5 1.0954 1.1832 1.2649 1.3416 1.4142nxny 0.6 0.7 0.8 0.9 1.0 1.4832 1.5492 1.6125 1.6733 1.7321nxny21yx.

8、二、二、高阶高阶和隐式和隐式Runge-Kutta方法方法注注：对于显式对于显式N级级R-K方法，最多只能得到方法，最多只能得到N级方法；级方法； N 1,2,3,4 5,6,7 8,9 10,11, N N-1 N-2()p N已经证明已经证明N级级R-K方法的方法的阶阶 具有下列关系：具有下列关系：()p N2N 若要得到若要得到N阶以上方法，则使用阶以上方法，则使用N级隐式级隐式R-K方法方法 N级隐式级隐式R-K方法的一般形式：方法的一般形式：11Nnniiiyyhc k 11(,);,.,Nininijjjkf xa h yhb kiN N级隐式级隐式R-K法法可以达到可以达到2N阶

9、阶.(1)(1)一一级级二阶二阶的隐式的隐式中点中点方法：方法：11nnyyhk 1122(,)nnhkhkf xy (2)(2)二二级级四阶四阶的隐式的隐式R-K方法：方法：1212()nnh kkyy 112131326446() ,()nnhkkf xh yhk221131326464() ,()nnhkkf xh yhk .三、三、变步长变步长方法方法基本基本思想思想：根据精度：根据精度自动自动地选择地选择步长步长对于对于经典经典Runge-Kutta方法：方法：0 1 2, , ,n Step1:设从设从 出发，以出发，以 为步长，经过为步长，经过一步一步计算得到计算得到nxh511

10、( )()hnny xyCh Step2:取取 为步长，再从为步长，再从 出发，经过出发，经过两步两步计算得到计算得到nx2h521122()()( )hnnhy xyC.21111116( / )( )()()hnnhnny xyy xy 2111116( / )( )( ()()hhnnnny xyy xy 221111115( / )( / )( )()hhhnnnny xyyy 211( / )( )|hhnnyy 记记如果如果 ，则将步长，则将步长折半折半进行计算，直到进行计算，直到 为止为止 此时取此时取 为最终结果；为最终结果；21( / )hny 如果如果 ，则将步长，则将步长

11、加倍加倍进行计算，直到进行计算，直到 为止为止 此时将步长此时将步长折半折半一次计算，得到的为最终结果。一次计算，得到的为最终结果。.一、一、收敛性收敛性 /* *Convergence* */3 单步法的单步法的收敛性收敛性、相容性相容性和和绝对稳定性绝对稳定性1(, )nnnnyyhxy h 对于初值问题对于初值问题 的一种的一种000( , )();dyf x ydxy xyxx ( ) 1Def单步法单步法 产生的近似解，如果产生的近似解，如果 对于任一对于任一固定固定的的 ，均有，均有 ，0nxxnh 则称该单步法是则称该单步法是收敛收敛的。的。0lim()nnhyy x 类似地可以

12、定义类似地可以定义隐式隐式单步法、多步法（单步法、多步法（4）的）的收敛性收敛性.3 1 .Th设初值问题（设初值问题（*）对应的下列）对应的下列单步法单步法是是 阶的，阶的，1(, )nnnnyyhxy h p且函数且函数 满足对满足对 的的Lipschitz条件，即存在常数条件，即存在常数y 0L 121212| ( , )( , )|,x y hx y hL yyy y 则该则该单步法单步法是收敛的，且是收敛的，且()()pnny xyO h 证明：证明：()nnney xy 记记由由截断截断误差的定义误差的定义11()()(, (), )nnnnny xy xhxy xhT 11 (,

13、 (), )(, )nnnnnnneehxy xhxyhT .因为因为单步法单步法是是 阶的：阶的：p000,hhh 满足满足11|pnTCh 11| |pnnneehL eCh |ne 其中其中11,phLCh 212|nnneee 3231|()ne 2101|(.)nnnee 10001| exp ()|exp ()pnnneL xxeCh LL xx 101exp ()pnCh LL xx ()pneO h 00()h.二、二、相容性相容性 /* *Consistency* */()( )( , ( ), )y xhy xhx y x h 0 ( )( ).( ) ( , ( ), )

14、.y xhy xy xhx y x 0( )( , ( ), ).h y xx y x 1()pO h 100( )( , ( ), )py xx y x 对于对于 阶方法：阶方法：1(, )nnnnyyhxy h p() 若方法（若方法（*）的）的增量增量函数满足：函数满足：2Def0( , , )( , )x yf x y 则称该方法与初值问题（则称该方法与初值问题（*）相容相容。.设方法（设方法（*）与初值问题（）与初值问题（*）相容相容，且，且 满足满足L-条件，条件， 则该方法（则该方法（*）是）是收敛收敛的，即当的，即当 固定，固定， 时时nxx 0h()nnyy x1(, )nn

15、nnyyxyhh 0(, (), )()nnnxy xy x ()(, ()nnny xf xy x 再由再由相容性相容性得：得：上式说明：当上式说明：当 时，方法（时，方法（*）趋于）趋于原微分方程原微分方程0h 本章讨论的数值方法都是与原初值问题本章讨论的数值方法都是与原初值问题相容相容的的 .三、三、绝对稳定性绝对稳定性 /* *Absolute Stibility* */计算过程中产生的计算过程中产生的舍入误差舍入误差对计算结果的影响对计算结果的影响首先以首先以Euler公式为例，来讨论一下公式为例，来讨论一下舍入误差舍入误差的传播的传播:1(,)nnnnyyhf xy 设设实际实际计

16、算得到的点计算得到的点 的的近似近似函数值为函数值为 ，nnnyy nxny其中其中 为为精确值精确值， 为误差为误差n 1(,)nnnnyyhf xy 111nnnyy 11 (,)(,)(, )nnnnnnynnh f xyf xyhfx 如果如果 ，则误差是，则误差是不增不增的，故可认为是的，故可认为是稳定稳定的的11|yhf .例如：例如：对于初值问题对于初值问题0()yyy xa 精确解精确解为为0 x xyae 而而实际求解实际求解的初值问题为的初值问题为0()yyy xaa 精确解精确解为为0()x xyaa e 在在 处的误差为处的误差为nx0nxxae 可见误差随着可见误差随

17、着 的增加呈的增加呈指数函数指数函数增长增长nx如果初值问题为如果初值问题为0()yyy xa 精确解精确解为为0 xxyae .实际求解实际求解的初值问题为的初值问题为0()yyy xaa 精确解精确解为为0()xxyaa e 在在 处的误差为处的误差为nx0nxxae 可见误差随着可见误差随着 的增加呈的增加呈指数函数指数函数递减递减nx当当 时，微分方程是时，微分方程是不稳定不稳定的；的；0yf 而而 时，微分方程是时，微分方程是稳定稳定的。的。0yf 上面讨论的上面讨论的稳定性稳定性，与，与数值方法数值方法和方程中和方程中 有关有关f.实验实验方程：方程：0,Re( )yyC 1(,

18、)nnnnyyhxy h 对单步法对单步法 应用应用实验实验方程，方程，1()nnyEh y 3Defh 如果如果 ，当，当 时，则称该时，则称该1()Eh 单步法是单步法是绝对稳定绝对稳定的，在复平面上复变量的，在复平面上复变量 满足满足1()Eh 的区域，称为该单步法的绝对稳定的区域，称为该单步法的绝对稳定域域，它与它与实轴实轴的的交集交集称为绝对稳定称为绝对稳定区间区间。1111()()pnnnTy xyO h 若单步法是若单步法是 阶的，则阶的，则p由由实验实验方程可得：方程可得：1()exp()nny xyh 1exp()()()pnnyhEh yO h ()exp()Ehh .11()()nnnnyyhyhy 例例3：分别求分别求Euler法和法和经典的经典的R-K法的法的绝对稳定绝对稳定区间区间。1()Ehh 解：解： Euler公式：公式：1(,)nnn