龙格库塔方法专业课件

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 1欢迎下载欢迎下载 可修改可修改可将改进的欧拉

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

3、想为例来阐述这种思想考察区间考察区间 上的一点上的一点 ，1,iix x 01, ipixxphp 用用 和和 的斜率的斜率 和和 的加权平均作为平均的加权平均作为平均 iipxx 12 KK斜率斜率 的近似值：的近似值：*K1122*KKK 即取即取11122()iiyyhKK 其中其中 和和 是待定常数。若取是待定常数。若取 ，则，则12 1(,)iiKf xy 问题在于如何确定问题在于如何确定 处的斜率处的斜率 和常数和常数 和和 。ipx 2K12 3欢迎下载欢迎下载 可修改可修改仿照改进的欧拉方法，用欧拉方法预测仿照改进的欧拉方法，用欧拉方法预测 的值，的值，()ipy x 1ipi

4、yyphK 并用它来估计斜率并用它来估计斜率 ：2K2(,)ipipKf xy 于是得到如下形式的算法：于是得到如下形式的算法：111221211(),(,),(,).iiiiiiyyhKKKf xyKf xyphK 通过适当选取参数通过适当选取参数 和和 的值，使得公式具有的值，使得公式具有12, p 2阶精度！阶精度！4欢迎下载欢迎下载 可修改可修改由泰勒公式展开，要使公式具有由泰勒公式展开，要使公式具有 2 阶精度阶精度，只需，只需12211 2,p方程组有方程组有无穷无穷多解：多解：二级二级方法有无穷多种方法有无穷多种常见的常见的3种二级方法：种二级方法： 中点法（修正的中点法（修正的

5、Euler法法）1221012,cca 取取122(,(,)nnnnnnhhyyhf xyf xy 二阶二阶龙格库塔方法龙格库塔方法122112,cca 取取12 (,)(,(,)nnnnnnnnhyyf xyf xh yhf xy 3()O h5欢迎下载欢迎下载 可修改可修改三级方法：三级方法：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

6、xy 2122(,)nnhhkf xyk3122(,)nnkf xh yhkhk 6欢迎下载欢迎下载 可修改可修改 四级方法：四级方法：N = 45()O h局部截断误差局部截断误差常见的常见的2种四阶方法：种四阶方法：经典经典龙格龙格- -库塔库塔方法方法11234226()nnhyykkkk 1(,)nnkf xy 2122(,)nnhhkf xyk3222(,)nnhhkf xyk43(,)nnkf xh yhk7欢迎下载欢迎下载 可修改可修改解：解：201( )dyxydxyy 01( , )x 例例2：用用经典的经典的龙格龙格- -库塔库塔方法方法求解下列初值问题求解下列初值问题 。

7、0 1 .h 经典的四阶经典的四阶龙格龙格- -库塔公式：库塔公式：11234226()nnhyykkkk 12;nnnxkyy2112222()nnnhxhkykhyk 4332()nnnxhkyhkyhk 3222222();nnnhxhkykhyk 8欢迎下载欢迎下载 可修改可修改 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.

8、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.7321nxny21yx9欢迎下载欢迎下载 可修改可修改10欢迎下载欢迎下载 可修改可修改9、 人的价值，在招收诱惑的一瞬间被决定。2022-2-262022-2-26Saturday, February 26, 202210、低头要有勇气，抬头要有低气。2022-2-262022-2-262022-2-262/26/2022 12:11:24 AM11、人总是珍惜为得到。2022-2-262022

9、-2-262022-2-26Feb-2226-Feb-2212、人乱于心，不宽余请。2022-2-262022-2-262022-2-26Saturday, February 26, 202213、生气是拿别人做错的事来惩罚自己。2022-2-262022-2-262022-2-262022-2-262/26/202214、抱最大的希望，作最大的努力。2022年2月26日星期六2022-2-262022-2-262022-2-2615、一个人炫耀什么，说明他内心缺少什么。2022年2月2022-2-262022-2-262022-2-262/26/202216、业余生活要有意义，不要越轨。20

10、22-2-262022-2-26February 26, 202217、一个人即使已登上顶峰，也仍要自强不息。2022-2-262022-2-262022-2-262022-2-26二、二、高阶高阶和隐式和隐式Runge-Kutta方法方法注注：对于显式对于显式N级级R-K方法，最多只能得到方法，最多只能得到N级方法；级方法； N 1,2,3,4 5,6,7 8,910,11, N N-1 N-2()p N已经证明已经证明N级级R-K方法的方法的阶阶 具有下列关系：具有下列关系：()p N2N 若要得到若要得到N阶以上方法，则使用阶以上方法，则使用N级隐式级隐式R-K方法方法 N级隐式级隐式R

11、-K方法的一般形式：方法的一般形式：11Nnniiiyyhc k 11(,);,.,Nininijjjkf xa h yhb kiN N级隐式级隐式R-K法法可以达到可以达到2N阶阶12欢迎下载欢迎下载 可修改可修改9、 人的价值，在招收诱惑的一瞬间被决定。2022-2-262022-2-26Saturday, February 26, 202210、低头要有勇气，抬头要有低气。2022-2-262022-2-262022-2-262/26/2022 12:11:24 AM11、人总是珍惜为得到。2022-2-262022-2-262022-2-26Feb-2226-Feb-2212、人乱于心

12、，不宽余请。2022-2-262022-2-262022-2-26Saturday, February 26, 202213、生气是拿别人做错的事来惩罚自己。2022-2-262022-2-262022-2-262022-2-262/26/202214、抱最大的希望，作最大的努力。2022年2月26日星期六2022-2-262022-2-262022-2-2615、一个人炫耀什么，说明他内心缺少什么。2022年2月2022-2-262022-2-262022-2-262/26/202216、业余生活要有意义，不要越轨。2022-2-262022-2-26February 26, 202217、

13、一个人即使已登上顶峰，也仍要自强不息。2022-2-262022-2-262022-2-262022-2-26(1)(1)一一级级二阶二阶的隐式的隐式中点中点方法：方法：11nnyyhk 1122(,)nnhkhkf xy (2)(2)二二级级四阶四阶的隐式的隐式R-K方法：方法：1212()nnh kkyy 112131326446() ,()nnhkkf xh yhk 221131326464() ,()nnhkkf xh yhk 14欢迎下载欢迎下载 可修改可修改三、三、变步长变步长方法方法基本基本思想思想：根据精度：根据精度自动自动地选择地选择步长步长对于对于经典经典Runge-Kut

14、ta方法：方法：0 1 2, , ,n Step1:设从设从 出发，以出发，以 为步长，经过为步长，经过一步一步计算得到计算得到nxh511( )()hnny xyCh Step2:取取 为步长，再从为步长，再从 出发，经过出发，经过两步两步计算得到计算得到nx2h521122()()( )hnnhy xyC 15欢迎下载欢迎下载 可修改可修改21111116( / )( )()()hnnhnny xyy xy 2111116( / )( )( ()()hhnnnny xyy xy 221111115( / )( / )( )()hhhnnnny xyyy211( / )( )|hhnnyy

15、记记如果如果 ，则将步长，则将步长折半折半进行计算，直到进行计算，直到 为止为止 此时取此时取 为最终结果；为最终结果；21( / )hny 如果如果 ，则将步长，则将步长加倍加倍进行计算，直到进行计算，直到 为止为止 此时将步长此时将步长折半折半一次计算，得到的为最终结果。一次计算，得到的为最终结果。16欢迎下载欢迎下载 可修改可修改一、一、收敛性收敛性 /* *Convergence* */3 单步法的单步法的收敛性收敛性、相容性相容性和和绝对稳定性绝对稳定性1(, )nnnnyyhxy h 对于初值问题对于初值问题 的一的一种种000( , )();dyf x ydxy xyxx ( )

16、 1Def单步法单步法 产生的近似解，如果产生的近似解，如果 对于任一对于任一固定固定的的 ，均有，均有 ，0nxxnh 则称该单步法是则称该单步法是收敛收敛的。的。0lim()nnhyy x 类似地可以定义类似地可以定义隐式隐式单步法、多步法（单步法、多步法（4）的）的收敛性收敛性17欢迎下载欢迎下载 可修改可修改3 1 .Th设初值问题（设初值问题（*）对应的下列）对应的下列单步法单步法是是 阶的，阶的，1(, )nnnnyyhxy h p且函数且函数 满足对满足对 的的Lipschitz条件，即存在常数条件，即存在常数y 0L 121212| ( , )( , )|,x y hx y h

17、L yyy y 则该则该单步法单步法是收敛的，且是收敛的，且()()pnny xyO h 证明：证明：()nnney xy 记记由由截断截断误差的定义误差的定义11()()(, (), )nnnnny xy xhxy xhT 11 (, (), )(, )nnnnnnneehxy xhxyhT 18欢迎下载欢迎下载 可修改可修改因为因为单步法单步法是是 阶的：阶的：p000,hhh 满足满足11|pnTCh 11| |pnnneehL eCh |ne 其中其中11,phLCh 212|nnneee 3231|()ne 2101|(.)nnnee 10001| exp ()|exp ()pnnn

18、eL xxeCh LL xx 101exp ()pnCh LL xx ()pneO h 00()h19欢迎下载欢迎下载 可修改可修改二、二、相容性相容性 /* *Consistency* */()( )( , ( ), )y xhy xhx y x h 0 ( )( ).( ) ( , ( ), ).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( , ,

19、)( , )x yf x y 则称该方法与初值问题（则称该方法与初值问题（*）相容相容。20欢迎下载欢迎下载 可修改可修改设方法（设方法（*）与初值问题（）与初值问题（*）相容相容，且，且 满足满足L-条件，条件， 则该方法（则该方法（*）是）是收敛收敛的，即当的，即当 固定，固定， 时时nxx 0h()nnyy x1(, )nnnnyyxyhh 0(, (), )()nnnxy xy x ()(, ()nnny xf xy x 再由再由相容性相容性得：得：上式说明：当上式说明：当 时，方法（时，方法（*）趋于）趋于原微分方程原微分方程0h 本章讨论的数值方法都是与原初值问题本章讨论的数值方法

20、都是与原初值问题相容相容的的 21欢迎下载欢迎下载 可修改可修改三、三、绝对稳定性绝对稳定性 /* *Absolute Stibility* */计算过程中产生的计算过程中产生的舍入误差舍入误差对计算结果的影响对计算结果的影响首先以首先以Euler公式为例，来讨论一下公式为例，来讨论一下舍入误差舍入误差的传播的传播:1(,)nnnnyyhf xy 设设实际实际计算得到的点计算得到的点 的的近似近似函数值为函数值为 ，nnnyy nxny其中其中 为为精确值精确值， 为误差为误差n 1(,)nnnnyyhf xy 111nnnyy 11 (,)(,)(, )nnnnnnynnh f xyf xy

21、hfx 如果如果 ，则误差是，则误差是不增不增的，故可认为是的，故可认为是稳定稳定的的11|yhf 22欢迎下载欢迎下载 可修改可修改例如：例如：对于初值问题对于初值问题0()yyy xa 精确解精确解为为0 x xyae 而而实际求解实际求解的初值问题为的初值问题为0()yyy xaa 精确解精确解为为0()x xyaa e 在在 处的误差为处的误差为nx0nxxae 可见误差随着可见误差随着 的增加呈的增加呈指数函数指数函数增长增长nx如果初值问题为如果初值问题为0()yyy xa 精确解精确解为为0 xxyae 23欢迎下载欢迎下载 可修改可修改实际求解实际求解的初值问题为的初值问题为0

22、()yyy xaa 精确解精确解为为0()xxyaa e 在在 处的误差为处的误差为nx0nxxae 可见误差随着可见误差随着 的增加呈的增加呈指数函数指数函数递减递减nx当当 时，微分方程是时，微分方程是不稳定不稳定的；的；0yf 而而 时，微分方程是时，微分方程是稳定稳定的。的。0yf 上面讨论的上面讨论的稳定性稳定性，与，与数值方法数值方法和方程中和方程中 有关有关f24欢迎下载欢迎下载 可修改可修改实验实验方程：方程：0,Re( )yyC 1(, )nnnnyyhxyh 对单步法对单步法 应用应用实验实验方程，方程，1()nnyEh y 3Defh 如果如果 ，当，当 时，则称该时，则

23、称该1()Eh 单步法是单步法是绝对稳定绝对稳定的，在复平面上复变量的，在复平面上复变量 满足满足1()Eh 的区域，称为该单步法的绝对稳定的区域，称为该单步法的绝对稳定域域，它与它与实轴实轴的的交集交集称为绝对稳定称为绝对稳定区间区间。1111()()pnnnTy xyO h 若单步法是若单步法是 阶的，则阶的，则p由由实验实验方程可得：方程可得：1()exp()nny xyh 1exp()()()pnnyhEh yO h ()exp()Ehh 25欢迎下载欢迎下载 可修改可修改11()()nnnnyyhyhy 例例3：分别求分别求Euler法和法和经典的经典的R-K法的法的绝对稳定绝对稳定

24、区间区间。1()Ehh 解：解： Euler公式：公式：1(,)nnnnyyhf xy 将其应用于将其应用于实验实验方程方程绝对稳定绝对稳定域：域：11h 1120hh 当当 时，时，R 绝对稳定绝对稳定区间：区间：2 0(, ) 经典的经典的R-K公式：公式：11234226()nnhyykkkk 26欢迎下载欢迎下载 可修改可修改1(,)nnkf xy 211222(,)()nnnhhhkf xykyk 322222(,)()nnnhhhkf xykyk 433(,)()nnnkf xh yhkyhk 22()nhy 23224()nhhy 3nyhk 2341234()()()()!hh

25、hEhh 22334411234()!nnhhhyhy ny 27欢迎下载欢迎下载 可修改可修改当当 时，时，R 绝对稳定绝对稳定区间：区间：2 785 0(., ) 23411234()()()()!hhhEhh 可以证明：可以证明： 存在唯一存在唯一极小值点极小值点2341234( )!tttg tt 1 5960 270.; ( ).tg t 23411234( )!tttg tt 由由 得得28欢迎下载欢迎下载 可修改可修改例例4：求求梯形梯形公式（公式（隐式隐式方法）方法）的的绝对稳定绝对稳定区间区间。解：解： 梯形梯形公式：公式：1112 (,)(,)nnnnnnhyyf xyf xy 112nnnnhyyyy 将其应用于将其应用于实验实验方程方程11212nnhyyh 1212()hEhh 10()Ehh 当当 时，时，R 绝对稳定绝对稳定区间：区间：0(, )29欢迎下载欢迎下载 可修改可修改9、 人的价值，在招收诱惑的一瞬间被决定。22.2.2622.2.26Saturday, February 26, 202210、低头要有勇气，抬头要有低气。*2/26/2022 12:11:24 AM