Runge-Kutta法

17.2Runge-Kutta法§龙格-库塔(Runge-Kutta)方法简称R-K法，是一种应用较广的高精度的单步法。所谓单步法就是在计算yi时只用到前一步信息yi-1的方法。本节介绍R-K法的构造原理、常用公式。2一、Runge-Kutta方法的构造原理0)(),(yaybxayxfy对于常微分方程的初值问题的解),(xyy1[,],iixx在区间上使用微分中值定理有11()()()()iiiiyxyxyxx1(,)iixx1()()()iiyxyxhy即----------(3)1iihxx其中，为步长31()[,]iiKyyxxx可以认为是在区间上的平均斜率引入记号()Ky[,()]fy1()[,]iiyxxxK只要使用适当的方法求出在区间上平均斜率的近似值就可得到相应的Runge-Kutta方法1nxnxxy)(xyyK4二、低阶Runge-Kutta方法1nxnxxy)(xyy如下图111.()()[,]iiiyxxyxxxK如果以在处的斜率作为在上的平均斜率即1()iKyx11[,()]iifxyx则(3)式化为111(,)iiiiyyhfxy11(,)iifxy即Euler方法Euler方法也称为一阶Runge-Kutta方法2()()iehOh由于----(4)KK51nxnxxy)(xyy11212.()()[,]iiiiyxxxKKyxxx如果以在和处的斜率和的算术平均值作为在上的平均斜率11()iKyx11(,)iifxy2()iKyx[,()]iifxyx(,)iifxy11(,)iifxyhK(由(4)式)令221KKK则(3)式化为K1K2K611211211()()222iiihyyhKKyKK111(,)iiKfxy211(,)iiKfxyhK)(00xyy-----------(5)即改进Euler公式,也称为二阶Runge-Kutta法3()()iehOh7三、高阶Runge-Kutta方法112[,]iimxxKKK如果上多预估几个点上的斜率值、、1()[,]iiyxxx用它们的加权平均值作为在上的平均斜率1nxnx21hxnxy)(xyyK1K2KmK111221222()(,)(,)(,)iimmiiiimimimyyhKKKKfxyKfxhyhKfxhyh01,()jjjjjyhyxh其中是的预估值。8RKjjjh上述公式就是公式。式中，，都是待定参数，可以通过是公式的局部截断误差关于的阶数尽量高来确定。下面以二阶R-K为例，说明R-K公式的构造过程。1112212221()(,)(,)iiiiiiyyhKKKfxyKfxhyhK122RK上述公式就是二阶公式。式中，，都是待定参数。-----------(6)1226为确定，，，需分析式（）的局部截断误差。91112212221212221222(,)(,)(,)(,)(,)(,)(,)()(,)(,)()(,)(,)iiiiiiiiiiiiiiiiiiiiiiiiiiiyyhKhKyhfxyhfxhyhKfxyfxyyhfxyhfxyhfxyohxyfxyfxyyhfxyhfxyxy3()oh()iiyyx假定1122322231222()(,)(,)(,)(,)()()()()()()()()()()()iiiiiiiiiiiiiyhfxfxfxhfxohxyyxhyxhyyxyxyxohxxxyy另一方面231()()()()()2iiiihyxyxhyxyxoh1011231222()()11()()()()2niiiiehyxyhyxhyxoh6式（）的局部截断误差:122211222RK1110121212有无穷多组解，从而可以得到许多具体的二阶-公式如：，，，111123(4)6iihyyKKK111(,)iiKfxy2111(,)22iihhKfxyK31121(,(2))iiKfxhyhKK-----------(7))(00xyy(7)式称为三阶Runge-Kutta方法用类似的方法也可以构造三阶R-K公式4()()()iiiyxyehOh因而三阶R-K方法(7)具有3阶精度1211234(22)6iihyyKKKK111(,)iiKfxy2111(,)22iihhKfxyK4113(,)iiKfxhyhK)(00xyy类似于(7)式,还可构造四阶(经典)Runge=Kutta方法3112(,)22iihhKfxyK5()()iehOh-----------(8)因而方法(8)有4阶精度13例1.使用高阶R-K方法计算初值问题1)0(5.002yxyy.1.0h取解:(1)使用三阶R-K方法时1n201yK12102)21.0(KyK1025.121203))2(1.0(KKyK2555.1)4(61.032101KKKyy1111.1RK.m14其余结果如下:(2)如果使用四阶R-K方法时1n201yK12102)21.0(KyK1025.1nxnk1k2k3yn1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993152203)21.0(KyK1133.12304)1.0(KyK2351.1)22(61.0432101KKKKyy1111.1其余结果如下:ixiK1K2K3K4yi1.00000.10001.00001.10251.11331.23511.11112.00000.20001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38922.78051.66675.00000.50002.77773.26003.34764.00572.00001600.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精Euler精00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精Euler精精精-精精精精00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精Euler精精精-精精精精精精Euler精00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精Euler精精精-精精精精精精Euler精3精RK精00.050.10.150.20.250.30.350.40.450.511.21.41.61.822.2精精精Euler精精精-精精精精精精Euler精3精RK精4精RK精170.29650.2970.29750.2980.29850.2990.29950.30.30051.4181.421.4221.4241.4261.428精精精精精Euler精3精RK精4精RK精