matlab 龙格库塔法 变步长龙格库塔法.doc

河北科技大学硕士学位研究生20122013学年第二学期Matlab 语言及应用结课论文学 院 信息科学与工程学院专 业 电 路 与 系 统姓 名 张 利 超学 号 S2012014011经典龙格库塔法及变步长龙格库塔法1.经典龙格库塔法及变步长龙格库塔法 matlab 代码a.经典龙格库塔法文件 Rungkutta4.mfunction RRungkutta4f,a,b,N,yah b-a/N;x zeros1,N1;y zeros1,N1;x ahb;y1 ya;for i1Nk1 ff,xi,yi; k2 ff,xih/2,yih/2*k1; k3 ff,xih/2,yih/2*k2; k4 ff,xih,yih*k3; yi1 yih/6*k12*k22*k3k4;endb.变步长龙格库塔法文件change_step_RK.mfunction change_step_RKfun; p212p-1; while xendAbsTol; x1x2; y1y2; hh/2; x2,y2RK_ffun,xn,yn,h/2; end end xa,yaRK_ffun,h,xn,yn; xn1xa; yn1ya; nn1; end plotx,y,k; function xa,yaRK_ffun,h,x,y; k1funx,y; k2funxh/2,yh*k1/2; k3funxh/2,yh*k2/2; k4funxh,yh*k3; xaxh; yayh*k1k2*22*k3k4/6; 2.利用两种方法求解初值问题0AbsTol; x1x2; y1y2; hh/2; x2,y2RK_ffun,xn,yn,h/2; end end xa,yaRK_ffun,h,xn,yn; xn1xa; yn1ya; nn1; end fprintf i xi yin; for i1n fprintf2d 4.6f 4.6fn,i,xi,yi; end plotx,y,rfunction xa,yaRK_ffun,h,x,y k1funx,y; k2funxh/2,yh*k1/2; k3funxh/2,yh*k2/2; k4funxh,yh*k3; xaxh; yayh*k1k2*22*k3k4/6; 运行结果i xi yi1 0.000000 0.0000002 0.002500 0.0025033 0.007500 0.0075284 0.010000 0.0100505 0.015000 0.0151126 0.017500 0.0176537 0.022500 0.0227538 0.025000 0.0253129 0.030000 0.03045010 0.032500 0.03302811 0.037500 0.03820312 0.040000 0.04080013 0.045000 0.04601214 0.047500 0.04862815 0.052500 0.05387816 0.055000 0.05651217 0.060000 0.06179918 0.062500 0.06445219 0.067500 0.06977720 0.070000 0.0724493.结论从运行结果来看，变步长龙格库塔法产生的误差|y(i1)-y(i)| 比经典龙格库塔法的小，但是运算量比较大。单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累，经过推算得以验证。经典的四阶龙格-库塔公式从节点 出发，先以 h 为步长求出一个近似值 ，nx hny1由于公式的局部截断误差为 ，故有5o然后将步长折半，即取 为步长从 跨两步到 ，再求得一个近似值 ，2hnx1n 21hny每跨一步的截断误差是 ，因此有5o比较(3.14)式和(3.15)式我们看到，步长折半后，误差大约减少到 ，既有16由此易得下列事后估计式 .15 12121 hnhnhnn yyxy 这样，可以通过检查步长，折半前后两次计算结果的偏差来判定所选的步长是否合适. 具体地说，将区分以下两种情况处理1. 对于给定的精度 ，如果 ，反复将步长折半进行计算，直至 为止.这时取最终得到的 作为结果； 21hny2. 如果 ，反复将步长加倍，直到 为止，这时再将步长折半一次，就得到所要的结果.表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的.