龙格-库塔方法(Runge-Kutta).docxVIP

龙格-库塔方法（Runge-Kutta） Runge-Kutta 法 显式 Runge-Kutta 法的一般形式 上节已给出与初值问题(1.2.1)等价的积分形式 (3.2.1) 只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题 (1.2.1)的数值方法，若用显式单步法 (3.2.2) 当 ，即数值求积用左矩形公式，它就是 Euler 法 (3.1.2)，方法只有一阶精度，若取 (3.2.3) 就是改进 Euler 法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数 f 的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的 f 值，根据数值积分公式，可将(3.2.1)右端积分表示为 注意，右端 f 中 还不能直接得到，需要像改进 Euler 法(3.1.11) 一样，用前面已算得的 f 值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为 (3.2.4) 其中 这里 均为待定常数，公式(3.2.2)，(3.2.4) 称为 r 级的显式 Runge-Kutta 法，简称 R-K 方法.它每步计算 r 个 f 值(即 )，而k 由前面(i-1)个已算出的 表示，故公式是显式的.例 i 如当 r=2 时，公式可表示为 (3.2.5) 其中 一个二级显式 R-K 方法.参数 二、三级显式 R-K 方法 .改进 Euler 法(3.1.11)就是取不同的值，可得到不同公式. 对 r=2 的显式 R-K 方法(3.2.5)，要求选择参数阶 p 尽量高，由局部截断误差定义 ，使公式的精度 T ? y(x ) ? y(x ) ? h[c f (x , y(x )) ? c f (x a h, y b hk )] (3.2.6) n?1 n?1 n 1 n n 2 n 2 n 21 1 令 ，对(3.2.6)式在 处按 Taylor 公式展开，由于 将上述结果代入(3.2.6)得 要使公式(3.2.5)具有的阶 p=2,即 ，必须 (3.2.7) 即 由（3.2.5）式可知 由此三式求 ，于是有解 的解不唯一.因 r=2， (3.2.8) 它表明使(3.2.5)具有二阶的方法很多，只要 都可得到二阶精度 R-K 方法.若取则得 ，则 ，则得改进 Euler 法(3.1.11)，若取 ， ，此时(3.2.5)为 (3.2.9) 其中 称为中点公式.改进的Euler 法(3.1.11)及中点公式(3.2.9) 是两个常用的二级 R-K 方法，注意二级 R-K 方法只能达到二阶，而不可能达到三阶.因为 r=2 只有 4 个参数，要达到 p=3 则在(3.2.6)的展开式中要增加 3 项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求 4 个待定参数，验证得 出是无解的.当然 r=2,p=2 的 R-K 方法(3.2.5)当 他公式，但系数较复杂，一般不再给出. 对 r=3 的情形，要计算三个 k 值，即 其中 取其他数时，也可得到其 将 按二元函数在 处按 Taylor 公式展开，然后代入局部截断误差表达式，可得 可得三阶方法，其系数共有 8 个，所应满足的方程为 这是 8 个未知数 6 个方程的方程组，解也是不唯一的，通常三级三阶 R-K 方法是下面的三级 Kutta 方法： .一种常见的 (3.2.11) 附：R-K 的三级 Kutta 方法程序如下 1/6) function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]); y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足 c1+c2+c3=1,(1/6 4/6 end format short; 四阶 R-K 方法及步长的自动选择 利用二元函数 Taylor 展开式可以确定(3.2.4)中 r=4,p=4 的 R-K 方法，其迭代公式为 y ? y h(c k ? c k c k ? c k ) n?1 n 1 1 2 2 3 3 4 4 其中k ? f (x , y ) ， k ? f (x a h, y b hf (x , y )) ，而