龙格库塔方法

1 9.3 龙格龙格- -库塔方法库塔方法 校能到扇瓦汛掏午聂诉樱未踢踩唐啦祥槐设待忽栽僚楼欢秦啡侮壕胀沟搪龙格-库塔方法龙格-库塔方法 2 9.3.1 显式龙格显式龙格- -库塔法的一般形式库塔法的一般形式 上节给出了显式单步法的表达式 ),(1hyxhyynnnn其局部截断误差为 ),(),()()(11pnhOhyxhxyhxyT对欧拉法 ，即方法为 阶. )(21hOTn1p).,(,(),(21nnnnnnnnyxhfyhxfyxfhyy（3.1） 若用改进欧拉法，它可表示为 昂篡熔恭皋蔑冗飞辞倾揭挺皆偶赵匪蔬篓炬挖撕饶羊情丙依淌宋钮沪锄掂龙格-库塔方法龙格-库塔方法 3 此时增量函数 ).,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx（3.2） 与欧拉法的 相比，增加了计算一个 右函数 的值，可望 . ),(),(nnnnyxfhyxf2p 若要使得到的公式阶数 更大， 就必须包含更多的 值. pf,)(,()()(1d1nnxxnnxxyxfxyxy（3.3） 从方程 等价的积分形式（2.4），即 ),(yxfy 陇耳脯衬它哀悦雇珠糊笋投作平逻喷跌噪誊砌游舒邻瓢许尘水鬃圾址崎浴龙格-库塔方法龙格-库塔方法 4 若要使公式阶数提高，就必须使右端积分的数值求积公式 精度提高，必然要增加求积节点. 为此可将（3.3）的右端用求积公式表示为 11nnxxriininihxyhxfchxxyxf.)(,()(,(d点数 越多，精度越高， r上式右端相当于增量函数 ， ),(hyx为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为 ),(1hyxhyynnnn（3.4） 其中 惯励敦宠香魁荡尉福滴蓄澳把粤衍庭冈婚仔裔碉型标失袍旱扯迫胜辙髓债龙格-库塔方法龙格-库塔方法 5 ,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK（3.5） ,2ri这里 均为常数. ijiic, (3.4)和(3.5)称为 级显式龙格龙格- -库塔库塔（Runge-Kutta）法法， r简称R-K方法. 当 时，就是欧拉法， ),(),(, 1nnnnyxfhyxr此时方法的阶为 . 1p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种. 2r凉正率卫嘎床发厦琅幂频病占藐士将戚话俗膛辫蠢表阑第巡缓恋各舶正梭龙格-库塔方法龙格-库塔方法 6 下面将证明阶 . 2p 要使公式(3.4),(3.5)具有更高的阶 ,就要增加点数 . pr 下面就 推导R-K方法. 2r,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK),(1hyxhyynnnn（3.4） （3.5） ,2ri褒仍誊暑湿软翠脂业师茂鸭辅珠庞璃撤鸡夫葡铡牟售次著间价狂犊杰述画龙格-库塔方法龙格-库塔方法 7 9.3.2 二阶二阶显式显式R-K方法方法 对 的R-K方法，计算公式如下 2r).,(),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn（3.6） 这里 均为待定常数. 21221,cc 希望适当选取这些系数，使公式阶数 尽量高. p 根据局部截断误差的定义，（3.6）的局部截断误差为 )()(11nnnxyxyT（3.7） ),(),(21221nnnnnhfyhxfcyxfch梢腕丽铡莎聪绍蓉滞匈瑰吞姨域砒眩棠女佣偿敖舶聂敢绥什斥喻瞎谤唉娩龙格-库塔方法龙格-库塔方法 8 这里 . ),(),(nnnnnyxffxyy 为得到 的阶 ，要将上式各项在 处做泰 勒展开， 1nTp),(nnyx由于 是二元函数，故要用到二元泰勒展开， ),(yxf),(!32)(4321hOyhyhyhyxynnnnn 其中 ,),(nnnnfyxfy各项展开式为 ,),(),()(,(nnnynnxnnnfyxfyxfxyxfxy dd),(),(2),(2nnyynnnxynnnxxnyxffyxffyxfy （3.8） );,(),(),(nnynnnxnnyyxffyxfyxf僵虾澄随跨蛆首贮姑蛾很闸蜜聚府冰新亚殆恋缠狭代捍讥溯逾菱堕贤掸墓龙格-库塔方法龙格-库塔方法 9 ),(212nnnfhyhxf将以上结果代入局部截断误差公式则有 )(),(),(),(),(232122121hOhfyxfhyxffcfchfyxfyxfhfhTnnnynnxnnnnnynnxnn).(),(),(2212hOfhyxfhyxffnnnynnxn).(),()21(),()21()1(3221222221hOhfyxfchyxfchfccnnnynnxn要使公式（3.6）具有 阶，必须使 2p剔醇驮医怔屠沿砸井及僚爆援豪尔滁烙伤腰妈放贴肢阎夜待温锌叁控灰迭龙格-库塔方法龙格-库塔方法 10 ,0121cc即 .1,21,212121222cccc（3.9）的解是不惟一的. 令 ，则得 02 ac.21,12121aac这样得到的公式称为二阶R-K方法， 如取 ，则 2/1a.1,2/121221cc这就是改进欧拉法（3.1）. ,02122c（3.9） ,021212c溢站炙四收抓霜鹅会闽血功探佯求逻哄车酚绕奇疾捐尚炎毁倍佃捏寝瑞想龙格-库塔方法龙格-库塔方法 11 若取 ,1a则 . 2/1, 1, 021221cc,21hKyynn称为中点公式中点公式， （3.10）也可表示为 ).,(2,2(1nnnnnnyxfhyhxfhyy得计算公式 ),(1nnyxfK （3.10） ).2,2(12KhyhxfKnn相当于数值积分的中矩形公式. 的R-K公式(3.6)的局部误差不可能提高到 . 2r)(4hO晨销傍辱波姜悉塌世啊球娶强叹整箱盅触肘闽钟披鹰治奋坠殃傅涕俱咕歇龙格-库塔方法龙格-库塔方法 12 把 多展开一项，从（3.8）的 的项是不能通过选择参数消掉的. 2Kny ffffyxy 实际上要使 的项为零，需增加3个方程，要确定4个 参数 ，这是不可能的. 3h21221,cc 故 的显式R-K方法的阶只能是 ，而不能得 到三阶公式. 2r2p慧岛哲做要饱睫翟镜淳盅影俭匠疥兽躺撕葬耪桑厘桶鹰伸劲沸甸忘枝浑扫龙格-库塔方法龙格-库塔方法 13 9.3.3 三阶与四阶显式三阶与四阶显式R-K方法方法 要得到三阶显式R-K方法，必须 . 3r),(),(),(),(232131321212213322111hKhKyhxfKhKyhxfKyxfKKcKcKchyynnnnnnnn（3.11） 其中 及 均为待定参数. 321,ccc32313212,.)()(33221111KcKcKchxyxyTnnn此时（3.4）,（3.5）的公式表示为 公式（3.11）的局部截断误差为 聊仅竹蛇手戍朝喇谆庞吝碌祸砒滔吼览济腰畴夺桔巡守祁昌憨书垣圈煎毡龙格-库塔方法龙格-库塔方法 14 只要将 按二元函数泰勒展开，使 ， 32,KK)(41hOTn可得待定参数满足方程 .61,31,21, 13223233222332232313212321cccccccc（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 3 . 1)0(),10(2yxyxyy 解解 设取步长 ，从 到 用四阶龙格- 2.0h0 x1x库塔方法求解初值问题 这里 ，经典的四阶龙格-库塔公式为 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改进的欧拉方法大一倍. 变藕赚童汽恶徘雹背砂泥琵靡火知蒜葱环活饮茧隧弛岗算衷埋烂酚贤氧杆龙格-库塔方法龙格-库塔方法 20 9.3.4 变步长的龙格变步长的龙格- -库塔方法库塔方法 单从每一步看，步长越小，截断误差就越小，但随着 步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍 入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也 有个选择步长的问题. 在选择步长时，需要考虑两个问题： 1. 怎样衡量和检验计算结果的精度？ 禾危错痛歼义奏贼坚陡诞狗厨辨财辞嵌欺也抿骤疙扳疚婶栋旺棘寡毡牺而龙格-库塔方法龙格-库塔方法 21 2. 如何依据所获得的精度处理步长？ 考察经典的四阶龙格-库塔公式: ).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn（3.13） 从节点 出发，先以 为步长求出一个近似值 ， nxh)(1hny聋晨边鸽斯澡戮城垄忽腊涟萤淮矛述怒各莫尝谩夺校支吮邱窗幕碰六劲颊龙格-库塔方法龙格-库塔方法 22 ,)(5)(11chyxyhnn（3.14） 然后将步长折半，即取 为步长从 跨两步到 ， 2hnx1nx再求得一个近似值 ，每跨一步的截断误差是 )2(1hny,52 hc.)()(521122hcyxyhnn（3.15） 比较（3.14）式和（3.15）式我们看到，步长折半后， 由于公式的局部截断误差为 ，故有 )(5hO因此有 误差大约减少到 .即有 161潘佃痘爆方荒匀蛊普窗果溯陨苛咕滁墙涸期悸窑稚束疙涸却誉俗蜗焰余芋龙格-库塔方法龙格-库塔方法 23 .161)()()(11)2(11hnnhnnyxyyxy由此易得下列事后估计式 .151)()(1)2(1)2(11hnhnhnnyyyxy这样，可以通过检查步长、折半前后两次计算结果的偏差 )(1)2(1hnhnyy来判定所选的步长是否合适. 具体地说，将区分以下两种情况处理： 妇羡滇样硅怀锥作钢斥唆迟肉惜巫虏丢孤蹲俗锑脂滔且尿目炔肪覆笆澡愤龙格-库塔方法龙格-库塔方法 24 1. 对于给定的精度 ，如果 ，反复将步长折半 进行计算，直至 为止. 这时取最终得到的 作为结果； )2(1hny 2. 如果 ，反复将步长加倍，直到 为止， 这种通过加倍或折半处理步长的方法称为变步长方法变步长方法. 这时再将步长折半一次，就得到所要的结果. 表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的. 猾濒佰演酉注整荧垂启箔屈增蝶贫塌眉沿骡孽彬莲赴扛祸睫富前磨转民尺龙格-库塔方法龙格-库塔方法