《龙格-库塔方法》PPT课件

1、9.3 龙格龙格-库塔方法库塔方法 9.3.1 显式龙格显式龙格-库塔法的普通方式库塔法的普通方式 上节给出了显式单步法的表达式),(1hyxhyynnnn其部分截断误差为),(),()()(11pnhOhyxhxyhxyT对欧拉法 ，即方法为 阶.)(21hOTn1p).,(,(),(21nnnnnnnnyxhfyhxfyxfhyy3.1 假设用改良欧拉法，它可表示为 此时增量函数 ).,(,(),(21),(nnnnnnnnyxhfyhxfyxfhyx3.2与欧拉法的 相比，添加了计算一个右函数 的值，可望 . ),(),(nnnnyxfhyxf2p 假设要使得到的公式阶数 更大， 就必需

2、包含更多的 值. pf,)(,()()(1d1nnxxnnxxyxfxyxy3.3 从方程 等价的积分方式2.4，即 ),(yxfy 假设要使公式阶数提高，就必需使右端积分的数值求积公式精度提高，必然要添加求积节点. 为此可将3.3的右端用求积公式表示为11nnxxriininihxyhxfchxxyxf.)(,()(,(d点数 越多，精度越高，r上式右端相当于增量函数 ，),(hyx为得到便于计算的显式方法，可类似于改良欧拉法，将公式表示为 ),(1hyxhyynnnn3.4其中 ,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK3.5,2ri这

3、里 均为常数. ijiic, (3.4)和(3.5)称为 级显式龙格-库塔Runge-Kutta法，r简称R-K方法. 当 时，就是欧拉法，),(),(, 1nnnnyxfhyxr此时方法的阶为 . 1p 当 时，改良欧拉法(3.1)，(3.2)也是其中的一种.2r 下面将证明阶 . 2p 要使公式(3.4),(3.5)具有更高的阶 ,就要添加点数 . pr 下面就 推导R-K方法. 2r,),(1riiinnKchyx),(1nnyxfK ),(11ijjijniniKhyhxfK),(1hyxhyynnnn3.43.5,2ri 9.3.2 9.3.2 二阶显式二阶显式R-KR-K方法方法

4、对 的R-K方法，计算公式如下 2r).,(),(),(12122122111hKyhxfKyxfKKcKchyynnnnnn3.6这里 均为待定常数.21221,cc 希望适中选取这些系数，使公式阶数 尽量高. p 根据部分截断误差的定义，3.6的部分截断误差为 )()(11nnnxyxyT3.7),(),(21221nnnnnhfyhxfcyxfch这里 . ),(),(nnnnnyxffxyy 为得到 的阶 ，要将上式各项在 处做泰勒展开，1nTp),(nnyx由于 是二元函数，故要用到二元泰勒展开，),(yxf),(!32)(4321hOyhyhyhyxynnnnn 其中 ,),(nn

5、nnfyxfy各项展开式为 ,),(),()(,(nnnynnxnnnfyxfyxfxyxfxy dd),(),(2),(2nnyynnnxynnnxxnyxffyxffyxfy 3.8);,(),(),(nnynnnxnnyyxffyxfyxf),(212nnnfhyhxf将以上结果代入部分截断误差公式那么有 )(),(),(),(),(232122121hOhfyxfhyxffcfchfyxfyxfhfhTnnnynnxnnnnnynnxnn).(),(),(2212hOfhyxfhyxffnnnynnxn).(),()21(),()21()1(3221222221hOhfyxfchyxf

6、chfccnnnynnxn要使公式3.6具有 阶，必需使 2p,0121cc即 .1,21,212121222cccc3.9的解是不独一的. 令 ，那么得 02 ac.21,12121aac这样得到的公式称为二阶R-K方法， 如取 ，那么2/1a.1,2/121221cc这就是改良欧拉法3.1.,02122c3.9,021212c假设取,1a那么 . 2/1, 1,021221cc,21hKyynn称为中点公式， 3.10也可表示为 ).,(2,2(1nnnnnnyxfhyhxfhyy得计算公式 ),(1nnyxfK 3.10).2,2(12KhyhxfKnn相当于数值积分的中矩形公式. 的R

7、-K公式(3.6)的部分误差不能够提高到 . 2r)(4hO 把 多展开一项，从3.8的 看到展开式中 的项是不能经过选择参数消掉的.2Kny ffffyxy 实践上要使 的项为零，需添加3个方程，要确定4个参数 ，这是不能够的. 3h21221,cc 故 的显式R-K方法的阶只能是 ，而不能得到三阶公式. 2r2p 9.3.3 9.3.3 三阶与四阶显式三阶与四阶显式R-KR-K方法方法 要得到三阶显式R-K方法，必需 . 3r),(),(),(),(232131321212213322111hKhKyhxfKhKyhxfKyxfKKcKcKchyynnnnnnnn3.11其中 及 均为待定

8、参数.321,ccc32313212,.)()(33221111KcKcKchxyxyTnnn此时3.4,3.5的公式表示为 公式3.11的部分截断误差为 只需将 按二元函数泰勒展开，使 ，32, KK)(41hOTn可得待定参数满足方程 .61,31,21, 13223233222332232313212321cccccccc3.12这是8个未知数6个方程的方程组，解也不是独一的. 所以这是一簇公式. 满足条件3.12的公式3.11统称为三阶R-K公式. 一个常见的公式为 ).2,(),2,2(),(),4(62131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnn

9、nnn此公式称为库塔三阶方法. 继续上述过程，经过较复杂的数学演算，可以导出各种四阶龙格-库塔公式，以下经典公式是其中常用的一个： 可以证明其截断误差为 . )(5hO 四阶龙格-库塔方法的每一步需求计算4次函数值 ，f).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn3.13 例3.1)0(),10(2yxyxyy 解设取步长 ，从 到 用四阶龙格-2.0h0 x1x库塔方法求解初值问题 这里 ，经典的四阶龙格-库塔公式为yxyyxf2),(),22(643211KKKKhyynn,2

10、1nnnyxyK,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的计算结果，显然龙格-库塔方法的精度高. 但由于这里放大了步长 ，所以表9-3和表9-2所耗费的计算量几乎一样. )2.0(h 龙格-库塔方法的推导基于泰勒展开方法，因此它要求所求的解具有较好的光滑性质. 反之，假设解的光

11、滑性差，那么，运用四阶龙格-库塔方法求得的数值解，其精度能够反而不如改良的欧拉方法. 四阶龙格-库塔方法每一步要4次计算函数 ,其计算量f比每一步只需2次计算函数 的二阶龙格-库塔方法-f改良的欧拉方法大一倍. 9.3.4 9.3.4 变步长的龙格变步长的龙格- -库塔方法库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的减少，在一定求解范围内所要完成的步数就添加了. 步数的添加不但引起计算量的增大，而且能够导致舍入误差的严重积累. 因此同积分的数值计算一样，微分方程的数值解法也有个选择步长的问题. 在选择步长时，需求思索两个问题： 1. 怎样衡量和检验计算结果的精度？ 2. 如何

12、根据所获得的精度处置步长？ 调查经典的四阶龙格-库塔公式: ).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn3.13 从节点 出发，先以 为步长求出一个近似值 ，nxh)(1hny,)(5)(11chyxyhnn3.14然后将步长折半，即取 为步长从 跨两步到 ，2hnx1nx再求得一个近似值 ，每跨一步的截断误差是)2(1hny,52 hc.)()(521122hcyxyhnn3.15比较3.14式和3.15式我们看到，步长折半后，由于公式的部分截断误差为 ，故有 )(5hO因此有 误差大约减少到 .即有161.161)()()(11)2(11hnnhnnyxyyxy由此易得以下事后估计式 .151)()(1)2(1)2(11hnhnhnnyyyxy这样，可以经过检查步长、折半前后两次计算结果的偏向 )(1)2(1hnhn