数值分析教教案15

3.5龙格-库塔法龙格-库塔方法的根本原理是：用方程中函数在前一节点上取值的线性组合构造一个表示的近似值公式，从而防止了求时用的高阶导数。该方法有多种推导方法，这里用数值积分法推导。为此，先将微分方程略加变形得出：，对该式两边在相邻两节点和之间求积分，移项得出：(3-15)采用不同的方法计算式(3-15)中定积分，便可得出数值积分的不同近似结果。如果用矩形求积公式计算，代入式(3-15)就可得出，这个结果和由Taylor公式得出的Euler公式结果完全一样。二阶龙格-库塔公式假设用梯形求积公式计算式(3-15)中的积分，那么有：令，上式中里的用欧拉公式代换，那么可得出，把和代入上式那么有：令：，那么得出式(3-15)的一个近似结果：(3-16)这就是二阶龙-格库塔公式，它的局部截断差为。式(3-16)右边的函数是和的线性组合，而和是把和的值代入函数得出的。这样计算时不再像用Taylor公式那样求的导数。三阶龙格-库塔公式假设用抛物线(Simpson)求积公式计算式(3-15)中的积分，那么有：式中，而和都是未知的，可以用欧拉格式估算和。类似二阶龙格-库塔公式的推导，并令：，，把它们代入式(3-15)那么得出三阶龙格-库塔公式，它的局部截断误差为。〔3-17〕三阶龙格-库塔公式的MATLAB实现三阶龙格-库塔法计算公式〔还有休恩法〕为：三阶龙格-库塔公式的MATLAB程序代码如下所示：functiony=DELGKT3_kuta(f,h,a,b,y0,varvec)%f：一阶常微分方程的一般表达式的右端函数%h：积分步长%a：自变量取值下限%b：自变量取值上限%y0：函数初值%varvec：常微分方程的变量组%subs(S)表示：用数值替代所有的符号变量。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;DELGKT3_kuta函数运行时需要调用以下函数：functionfv=Funval(f,varvec,varval)var=findsym(f);iflength(var)<4ifvar(1)==varvec(1)fv=subs(f,varvec(1),varval(1));elsefv=subs(f,varvec(2),varval(2));endelsefv=subs(f,varvec,varval);end【例3-9】三阶龙格-库塔求解一阶常微分方程应用实例。用三阶龙格-库塔法求下面常微分方程的数值解。在编辑窗口输入以下程序段，然后执行该程序：symsxy;z=2*x-3*y+2;yy=DELGKT3_kuta(z,0.1,0,1,1,[xy])程序执行后得到如下结果：yy=1.00000.92250.88240.87000.87820.90150.93600.97891.02801.08161.1387四阶龙格-库塔公式假设用三次多项式代替(3-15)式中积分里的去计算定积分，那么有：〔参见表2-1，取n=3〕，式中，未知的和用欧拉公式代换后，一并代入式〔3-15〕，并令：，，那么可得出4阶龙格-库塔公式：〔3-18〕这里用3/8辛甫森求积法〔积分区间被分为3的倍数份〕得出的结果，它的局部截断误差为。更多使用的是基于1/3辛甫森求积法〔求积区间被分为偶数份〕得出的4阶龙格-库塔公式，它的推导方法和截断误差与式(3-18)相同：〔3-19〕其中,,，还可以构造出更高阶的龙格-库塔公式，它的一般形式是：(3-20)其中,,把式(3-20)中的与式(3-6)中的取成相同的值，然后进行拟合，让它们同幂次项和系数相等，确定出参数和，从而得出阶龙格-库塔公式，根据Taylor级数理论可知，它的总体截断误差是。四阶龙格-库塔公式的MATLAB实现四阶龙格-库塔法计算公式〔还有变形，基尔法〕为：四阶龙格-库塔公式的MATLAB程序代码如下所示：functiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec)%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;【例3-10】四阶龙格-库塔求解一阶常微分方程应用实例。用四阶龙格-库塔法求下面常微分方程的数值解。在编辑窗口输入以下程序段，然后执行该程序：z=1+log(x+1)yy=DELGKT4_lungkuta(z,0.1,0,1,1,[x,y])程序执行后得到如下结果：yy=1.00001.10481.21881.34111.47111.60821.75201.90212.05802.21952.38633.6阿达姆斯法泰勒公式和龙格-库塔公式在计算值时，只用前一点的和近似值，在给定初始后可以逐步往后推算。这种被称为一步法。此方法的优点是只要知道初值就可进行计算。但是，要提高精度就需要增加中间函数值〔如等值〕的计算，这就增大了计算量。下面以常用的4阶阿达姆斯法为例介绍多步法，这种方法在计算的值时，除了利用外还要用到已计算出的等，虽然也增多了预知函数值，但都是已经算出的，不再新增计算量。假设微分方程(3-4)解函数在的函数的值。由于式(3-15)中被积函数说明是的函数，那么可用这四点构造出一个三次Lagrange多项式近似地代替：于是得出：把代入上式积分。这里，仅计算第4项的积分，令，，那么：和都是常量，，于是得出：同样