第四讲龙格-库塔方法

龙格-库塔方法3.2Runge-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.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)，要求选择参数，使公式的精度阶p尽量高，由局部截断误差定义11122211()()[(,())(,)]nnnnnnnTyxyxhcfxyxcfxahybhk(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方法只能达到二阶，而不可能达到三阶.因为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方法程序如下functiony=DELGKT3_kuta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1)y(i-1)]);K2=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h/2]);K3=Funval(f,varvec,[x(i-1)+hy(i-1)-h*K1+K2*2*h]);y(i)=y(i-1)+h*(K1+4*K2+K3)/6;%满足c1+c2+c3=1,(1/64/61/6)endformatshort;3.2.3四阶R-K方法及步长的自动选择利用二元函数Taylor展开式可以确定(3.2.4)中r=4,p=4的R-K方法，其迭代公式为111223344()nnyyhckckckck其中1(,)nnkfxy，2221(,(,))nnnnkfxahybhfxy，而33311322(,)nnkfxahybhkbhk44411422433(,)nnkfxahybhkbhkbhk共计13个参数待定，Taylor展开分析局部截断误差，使得精度达到四阶，即误差为5()Oh。于是，r=4，p=4的13个参数(c4不能为0)引出了多种方案和挑战，如：参数优化使阶数增加到5阶，得到四阶五阶R-K方法，matlab中有程序ode45；四级四阶R-K方法的步长自动选取；结合新算法的应用算法构造；适应于新的领域实现求解；……经典的四阶R-K方法是：(3.2.12)其中也需满足12341cccc，这里为(1/62/62/61/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)精度较高，但要从理论上给出误差的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h？通常是通过不同步长在计算机上的计算结果近似估计.设在处的值,当时，的近似为()1hny，于是由四阶R-K方法有若以为步长，计算两步到，则有于是得即或(3.2.13)它给出了误差的近似估计.如果(ε为给定精度)，则认为以为步长的计算结果满足精度要求，若，则还可放大步长.因此(3.2.13)提供了自动选择步长的方法.附：经典的4级(也是4阶)R-K方法的程序：functiony=DELGKT4_Rungkuta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1)y(i-1)]);K2=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h/2]);K3=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K2*h/2]);K4=Funval(f,varvec,[x(i-1)+hy(i-1)+h*K3]);y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformatshort;讲解：求初值问题（1.2.1）的单步法主要是指Runge-Kutta法，本节主要讨论显式R－K方法，建立具体的计算公式使用的是Taylor展开，形如（3.2.4）的显式R－K方法，当r＝1时就是Euler法，因此只要讨论的计算公式，在r确定后如何推导公式都是一样的，只是r越大计算越复杂，为了掌握了解公式来源，只要以r＝2为例推导计算公式即可.因此本节重点就是用Taylor展开求出r＝2的显式R-K方法的计算公式，由于方法的局部截断误差为（3.2.6），的右端有的项，要对它做Taylor展开，就要用到二元函数的Taylor展开，按照二元函数Taylor级数（3.2.14）将它用到（3.2.6）的的展开式中，即可得到按升幂整理出的结果，对r＝2的公式只能得到＝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方法11122()nnyyhckck，其中1(,)nnkfxy，2221(,(,))nnnnkfxahybhfxy．四参数满足举例：中点公式四个参数2122111,0,2ccab，即1,(,)22nnnnnnhhyyhfxyfxy，程序如下functiony=DELGKT2_mid(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1v1=Funval(f,varvec,[x(i-1)y(i-1)]);t=y(i-1)+h*v1/2;v2=Funval(f,varvec,[x(i-1)+h/2t]);y(i)=y(i-1)+h*v2;endformatshort;休恩（suen）法c1=1/4,c2=3/4,a2=b21=2/3functiony=DELGKT2_suen(f,h,a,b,y0,varvec)formatlong;N=uint16((b-a)/h);y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1)y(i-1)]);K2=Funval(f,varvec,[x(i-1)+2*h/3y(i-1)+K1*2*h/3]);y(i)=y(i-1)+h*(K1+3*K2)/4;endformatshort;改进的Euler法c1=1/2,c2=1/2,a2=b21=1functiony=DEModifEuler(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1v1=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;endformatshort;实验对比参数取法不同得到不同的二阶R-K方法精确程度各异：这个已经可以做(同学们已做？！)注意：二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.2．三阶三级R-K方法其中1(,)nnkfxy，2221(,(,))nnnnkfxahybhfxy，而共计8个参数待设定。Taylor展开分析局部截断误差，使得精度达到三级，即误差为4()Oh。于是，r=3，p=3的8个参数所满足的关系式(c3不能为0)为Kutta法functiony=DELGKT3_kutta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1)y(i-1)]);K2=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h/2]);K3=Funval(f,varvec,[x(i-1)+hy(i-1)-h*K1+K2*2*h]);y(i)=y(i-1)+h*(K1+4*K2+K3)/6;endformatshort;Suen法functiony=DELGKT3_suen(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=Funval(f,varvec,[x(i-1