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

1、龙格-库塔方法（Runge-Kutta）3.2 Runge-Kutta法3.2.1 显式 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

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值(即)，而ki由前面(i-1)个已算出的表示，故公式是显式的.例如当r=2时，公式可表示为 (3.2.5)其中.改进Euler法(3.1.11)就是一个二级显式R-K方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K方法对r=2的显式R-K方法(3.2.5)，要求选择

3、参数，使公式的精度阶p尽量高，由局部截断误差定义 (3.2.6)令，对(3.2.6)式在处按Taylor公式展开，由于 将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须 (3.2.7)即由此三式求的解不唯一.因r=2，由（3.2.5）式可知，于是有解 (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方法只能达到二阶，而

4、不可能达到三阶.因为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个，所应满足的方程为 (3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法： (3.2.11)附：R-

5、K的三级Kutta方法程序如下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

6、*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short;3.2.3 四阶R-K方法及步长的自动选择利用二元函数Taylor展开式可以确定(3.2.4)中r=4,p=4的R-K方法，其迭代公式为其中，而共计13个参数待定，Taylor展开分析局部截断误差，使得精度达到四阶，即误差为。于是，r=4，p=4的13个参数(c4不能为0)引出了多种方案和挑战，如：参数优化使阶数增加到5阶，得到四阶五阶R-K方法，matlab中有程序ode45；四级四阶R-K方法的步长自动选取；结合新算法的应用算法构造；适应于新的领域实现求解；经典的四阶R-

7、K方法是： (3.2.12)其中也需满足，这里为(1/6 2/6 2/6 1/6).它的局部截断误差，故p=4，这是最常用的四阶R-K方法，数学库中都有用此方法求解初值问题的软件.这种方法的优点是精度较高，缺点是每步要算4个右端函数值，计算量较大.例3.3 用经典四阶R-K方法解例1.1的初值问题，仍取h=0.1，计算到，并与改进Euler法、梯形法在处比较其误差大小.解 用四阶R-K方法公式(3.2.12)，此处，于是当n=0时于是，按公式(3.2.12)可算出此方法误差：改进Euler法误差：梯形法误差：可见四阶R-K方法的精度比二阶方法高得多.用四阶R-K方法求解初值问题(1.2.1)精

8、度较高，但要从理论上给出误差的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h？通常是通过不同步长在计算机上的计算结果近似估计.设在处的值,当时，的近似为，于是由四阶R-K方法有若以为步长，计算两步到，则有 于是得即 或 (3.2.13)它给出了误差的近似估计.如果(为给定精度)，则认为以为步长的计算结果满足精度要求，若，则还可放大步长.因此(3.2.13)提供了自动选择步长的方法.附：经典的4级(也是4阶)R-K方法的程序：function y = DELGKT4_Rungkuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y

9、 = 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/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformat short;讲解：求初值问

10、题（1.2.1）的单步法主要是指Runge-Kutta法，本节主要讨论显式RK方法，建立具体的计算公式使用的是Taylor展开，形如（3.2.4）的显式RK方法，当r1时就是Euler法，因此只要讨论的计算公式，在r确定后如何推导公式都是一样的，只是r越大计算越复杂，为了掌握了解公式来源，只要以r2为例推导计算公式即可.因此本节重点就是用Taylor展开求出r2的显式R-K方法的计算公式，由于方法的局部截断误差为（3.2.6），的右端有的项，要对它做Taylor展开，就要用到二元函数的Taylor展开，按照二元函数Taylor级数（3.2.14）将它用到（3.2.6）的的展开式中，即可得到按升

11、幂整理出的结果，对r2的公式只能得到2阶的公式，即，于是2级R-K方法（3.2.5）的系数必须满足（3.2.7）给出的方程，它的解由（3.2.8）给出，只要，求出的公式都是r=2的2阶R-K方法.而常用的就是得到的改进Euler法（3.1.11）和得到的中点公式（3.2.9）.我们知道，显式单步法的数值公式为：如果取定右边的则称该数值公式为显式Runge-Kutta方法（一）显式Runge-Kutta方法 1. 二阶二级R-K方法，其中，四参数满足举例：中点公式四个参数，即，程序如下function y = DELGKT2_mid(f, h,a,b,y0,varvec)format long;

12、N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 v1 = Funval(f,varvec,x(i-1) y(i-1); t = y(i-1) + h*v1/2; v2 = Funval(f,varvec,x(i-1)+h/2 t); y(i) = y(i-1)+h*v2;endformat short;休恩（suen）法 c1=1/4, c2=3/4, a2=b21=2/3function y = DELGKT2_suen(f, h,a,b,y0,varvec)format long;N

13、 = uint16(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)+2*h/3 y(i-1)+K1*2*h/3); y(i) = y(i-1)+h*(K1+3*K2)/4;endformat short;改进的Euler法 c1=1/2, c2=1/2, a2=b21=1function y = DEModifEuler(f, h,a,b,y0,varvec)forma

14、t 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 v1 = Funval(f,varvec,x(i-1) y(i-1); t = y(i-1) + h*v1; v2 = Funval(f,varvec,x(i) t); y(i) = y(i-1)+h*(v1+v2)/2;endformat short;实验对比参数取法不同得到不同的二阶R-K方法精确程度各异：这个已经可以做(同学们已做？！)注意：二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=

15、3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.2三阶三级R-K方法其中，而共计8个参数待设定。Taylor展开分析局部截断误差，使得精度达到三级，即误差为。于是，r=3，p=3的8个参数所满足的关系式(c3不能为0)为Kutta法function y = DELGKT3_kutta(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 =

16、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;endformat short;Suen法function y = DELGKT3_suen(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 = findsy

17、m(f);for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/3 y(i-1)+K1*h/3); K3 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)+K2*2*h/3); y(i) = y(i-1)+h*(K1+3*K3)/4;endformat short;3. 四阶四级R-K方法举例：经典四级Rung-Kutta法程序function y = DELGKT4_Rungkutta(f, h,a,b,y0,varvec)format long;N = (b-a)

18、/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/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformat short;基尔(

19、jer、Gear)法function y = DELGKT4_jer(f, h,a,b,y0,varvec)format long;N = (b-a)/h;C2 = sqrt(2);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/2 y(i-1)+K1*h*(C2-1)/2+K2*h*(2

20、-C2)/2);K4 = Funval(f,varvec,x(i-1)+h y(i-1)-K2*h*C2/2+K3*h*(2+C2)/2); y(i) = y(i-1)+h*(K1+(2-C2)*K2+(2+C2)*K3+K4)/6;endformat short;变型Rung-Kutta法function y = DELGKT4_qt(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,

21、x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/3 y(i-1)+K1*h/3); K3 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)-K1*h/3+K2*h); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*(K1-K2+K3); y(i) = y(i-1)+h*(K1+3*K2+3*K3+K4)/8;endformat short;说明：四阶四级R-K方法的步长自动选取问题用四阶R-K方法求解初值问题(1.2.1)精度已经足够高了，但要从理论上给出误差的估计式则比较困难.那么应如何判断计

22、算结果的精度以及如何选择合适的步长h？通常是通过不同步长在计算机上的计算结果近似估计.设在处的值,当时，的近似为，于是由四阶R-K方法有若以为步长，计算两步到，则有 于是得即 或 (3.2.13)它给出了误差的近似估计.如果 (为给定精度，如1E-3)，则认为以为步长的计算结果满足精度要求，若，则还可放大步长.因此(3.2.13)提供了自动选择步长的方法.（二）隐式(半隐式)R-K方法举例一种半隐式举例，该格式可以通过除法化为显式单步法格式，受到除法复杂度限制，容易造成耗时，不稳，相对精度降低。罗赛布诺克半隐式法function y = DELSBRK(f, h,a,b,y0,varvec)f

23、ormat long;a1 = 1.40824829;a2 = 0.59175171;c1 = 0.17378667;c2 = c1;w1 = -0.41315432;w2 = 1.41315432;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);dy = diff(f, varvec(2);for i=2:N+1 f1 = Funval(f,varvec,x(i-1) y(i-1); dy1 = Funval(dy,varvec,x(i-1) y(i-1); k1 = h*f1/(1-h*a1*dy1); d

24、y2 = Funval(dy,varvec,x(i-1)+c1*h y(i-1)+c2*k1); f2 = Funval(f,varvec,x(i-1)+c1*h y(i-1)+c2*k1); k2 = h*f2/(1-h*a2*dy2); y(i) = y(i-1)+w1*k1+w2*k2;endformat short;3.3 单步法的收敛性与绝对稳定性3.3.1 单步法的收敛性定义3.1 设y(x)是初值问题(1.2.1)的精确解，是单步法(3.2.2)在处产生的近似解，若则称方法(3.2.2)产生的数值解收敛于.实际上，定义中是一固定点，当h0时n，n不是固定的.因显然方法收敛，则在固

25、定点处的整体误差，当p1时.下面定理给出方法(3.2.2)收敛的条件.定理3.1设初值问题(1.2.1)的单步法(3.2.2)是p阶方法(p1)，且函数对y满足Lipschitz条件，即存在常数L0，使对，均有 则方法(3.2.2)收敛，且.定理证明略.3.3.2 绝对稳定性用单步法(3.2.2)求数值解，由于原始数据及计算过程舍入误差影响，实际得到的不是而是，其中是误差，再计算下一步得到以Euler法为例，若令,则(3.3.1)如果，则从计算到误差不增长，它是稳定的.但如果条件不满足就不稳定.例3.4 y=-100y,y(0)=1,精确解为，用Euler法求解得若取h=0.025，则,当，而

26、，显然计算是不稳定的.如果用后退Euler法(3.1.5)解此例，仍取h=0.025,则,即显然当，计算是稳定的.由此看到稳定性与方法有关，也与有关，在此例中.在研究方法的稳定性时，通常不必对一般的f(x,y)进行讨论，而只针对模型方程 (3.3.2)这里可能为复数.规定是因为时微分方程(3.3.2)本身是不稳定的，而讨论数值方法(3.2.2)的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题(1.2.1)，若将f(x,y)在处线性展开，可得于是方程(1.2.1)可近似表示为它表明用模型方程(3.3.2)是合理的，至于模型方程(3.3.2)中所以用复数是因为初值问题(1.2.1

27、)如果是方程组，即，则是(mm)阶矩阵，其特征值可能是复数.当然对单个方程，就是实数，此时只要规定0即可.用单步法(3.2.2)解模型方程(3.3.2)可得到 (3.3.3)其中依赖所选方法，如用Euler法，则(3.3.4)此时由(3.3.1)看到误差方程也为，与(3.3.4)是一样的.因此对一般单步法(3.2.2)误差方程也与(3.3.3)一致.下面再考虑二阶R-K方法有对四阶R-K方法，可得定义3.2将单步法(3.2.2)用于解模型方程(3.3.2)，若得到(3.3.3)中的 则称方法是绝对稳定的.在复平面上复变量满足 的区域，称为方法(3.2.2)的绝对稳定域，它与实轴的交点称为绝对稳定区间.例如对Euler法， 在复平面上是以(-1，0)为圆心，以1为半径的单位圆域内部，当为实数时，则得绝对稳定区间为，因0，故有.在例3.4中 时方法稳定，而例中取h=0.025故不稳定.对后退Euler法(3.1.5)，因0，故，其绝对稳定域是以(1，0)为圆心的单位圆外部，绝对稳定区间为，即对任何h0方法都是绝对稳定的