龙格库塔方法

Runge-Kuttua方法和matlab原理主讲人：龙格－库塔法（Runge-Kutta）数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。四阶Runge-Kutta方法这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。当四个斜率取平均时，中点的斜率有更大的权值：误差分析：注意上述公式对于标量或者向量函数(y可以是向量)都适用。四阶R-K方法的每一步需要计算四次函数值f，可以证明其局部截断误差为O(h5).R-K(高阶)方法不唯一,选择不同的参数能得到不同的R-K公式注意的问题R-K方法的推导是基于Taylor展开法，因而要求解具有较好的光滑性，如果光滑性较差精度可能不如改进Euler方法,最好采用低阶算法而将步长h取小。Runge-Kutta法的主要运算在于计算Ki的值，即计算f的值。计算量与可达到的最高精度阶数的关系：753可达到的最高精度642每步须算Ki的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n四阶Runge-Kutta方法的MATLAB实现原理：四阶R-K方法实现开始00,,,xyhN输入10100200100200312342232210226;(,),(,)(,),(,)()xxhhhKfxyKfxyKhhKfxyKfxhyhKhyyKKKK输出x1,y1?nN10101;,nNxxyy结束YNfunctionff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间c=(b-a)/h+1;i1=1;%c为迭代步数；i1为迭代步数累加值y=y0;z=zeros(c,6);%z生成c行，6列的零矩阵存放结果；%每行存放c次迭代结果，每列分别存放k1~k4及y的结果不断迭代运算：forx=a:h:bifi1=ck1=feval(yy,x,y);k2=feval(yy,x+h/2,y+(h*k1)/2);k3=feval(yy,x+h/2,y+(h*k2)/2);k4=feval(yy,x+h,y+h*k3);y=y+(h/6)*(k1+2*k2+2*k3+k4);z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y;i1=i1+1;endend1)0()10(2)1.0()43(yxyxyyhKR步长解初值问题公式阶阶、求例4解yxyyxf/2),(.2,,2,2,,,462131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn例题4xnYn|yn-y(xn)|R-K3误差y(xn)0.11.09590.00051.095440.45e-41.09540.21.18410.00091.183220.17e-41.18320.31.26620.00131.264910.15e-41.26490.41.34340.00181.341650.48e-41.34160.51.41640.00221.414220.25e-41.41420.61.48600.00281.483260.55e-41.48320.71.55250.00331.549210.14e-41.54920.81.61650.00401.6124780.21e-41.61250.91.67820.00491.673350.54e-41.67331.01.73790.00581.732090.06e-41.7321改进Euler法一步需要计算两个函数值(h=0.1)四阶Runge-Kutta方法一步需要计算四个函数值（h=0.2）总计算量大致相当，但四阶Runge-Kutta方法精度更高五、变步长Runge-Kutta方法从每一步看，步长越小，截断误差越小；但随着步长的缩小，在一定求解范围内所要完成的步数就会增加，步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累，因此需要选择步长•如何衡量和检验计算结果的精度•如何依据所判定的精度来处理步长实施方案以经典四阶Runge-Kutta方法为例155151()()()()[,]hnnhnnnxhyChCyxxxn+1设从节点出发，先以为步长求出一个近似值，记为因为经典格式的局部截断误差为O(h),因此有y(x)-y其中与在内的值有关2112()+hnnnhxxy将步长折半，即取为步长从跨两步到，求得一个近似值5521222()(),()hnhhCCn+1每跨一步截断误差是有y(x)-y21111116116()()()()hnnhnnyxyyxy步长折半后，误差大约减少为原来的，即有221111115()()()()hhhnnnnyxyyy事后估计式为可以通过检查步长折半前后两次计算结果的偏差211()()=hhnnyy来判断所选取的步长是否合适,（1）对于给定的精度，如果则反复将步长折半进行计算直至为止，这时取折半以后的“新值”作为结果；（2）如果，则反复将步长加倍，直至为止，这时取步长加倍前的“老值”作为结果变步长方法Thanks！