《matlab 龙格库塔法 变步长龙格库塔法》由会员共享，可在线阅读。 更多相关《matlab 龙格库塔法 变步长龙格库塔法(7页珍藏版)》请在金锺文库搜索。
1、 河北科技大学硕士学位研究生20122013学年第二学期Matlab语言与应用结科论文学院：信息科学与工程学院专业：电路与系统名称： jsddy学号： S2012014011经典龙格库塔法与变长龙格库塔法1 .经典龙格库塔法与变长龙格库塔法Matlab代码a .经典龙格库塔法ngkutta4. mfunctionr=rung kutta4(f，a，b，n，ya ) h=(B-a )/N； x=Zeros(1，N 1 )； y=Zeros(1，N 1 )； x=a:h:b； y(1)=ya； forI1:NK1=feval(f，x ) I，y ) I； k。
2、2=feval(f，x ) I ) h/2、y ) I ) h/2 ) *k1 )； k3=feval(f，x ) I ) h/2，y ) I ) h/2 ) *k2 )； k4=feval(f，x ) I ) h，y ) I ) h*k3 )； y(I1 )=y (I ) ) h/6 ) ) k12*k3k4)； endb.step strange _ cooker方法文档： change _ step _ rk.mfunctionchange _ step _ rk (fun )； p21=2p-1； whilex(end ) AbsTol； x1=x2； y1=y2； h=h/2； x2，y2=rk_f(fun，x ) n，y ) n，h/2 )； 结束xa。
3，ya=rk_f(fun，h，x(n )，y ) n； x(n1 )=xa； y(n1 )=ya； n=n 1； endplot(x，y，k )； function xa，ya=rk_f(fun，h，x，y )； k1=fun(x，y )； k2=fun(xh/2，y h*k1/2 )； k3=fun(xh/2，y h*k2/2 )； k4=fun(xh，y h*k3 )； xa=x h； ya=yh*(k1k2*22*k3k4)/6； 2 .用两种方法求解初值问题0AbsTol的x1=x2； y1=y2； h=h/2； x2，y2=rk_f(fun，x ) n，y ) n，h/2 )； 结束。
4、end xa，ya=rk_f(fun，h，x ) n，y ) n )； x(n1 )=xa； y(n1 )=ya； n=n 1； endfprintf(Ix ) I ) y ) I ) n； fori=1: n fprintf (-% 4.6f % 4.6 fn，I，x ) I，y ) I； endplot(x，y，r ) function xa，ya=rk_f ) fun，h，x，y ) k1=fun(x ) x，y )； k2=fun(xh/2，y h*k1/2 )； k3=fun(xh/2，y h*k2/2 )； k4=fun(xh，y h*k3 )； xa=x h； ya=yh*(k1k2*22*k3k4)/6；
5、 运行结果Ix(I ) Ix(I ) ) 1.000000.0000020.00250330.007500.00752840.0100505.015000.0151126.017000 17500 129.030000.03045010.032500.033028110.037500.038203120.040000.040800130.04601214.0475000.04040
6， 053878160.055000.056512170.061799180.062500.06445219.067500.06977720.070000.0724493 .结论从运行结果来看，应当步长越小，截断误差越小，但随着步长的减小，在一定求解范围内需要完成的步数增加。 估计并验证了步骤数的增加不仅会导致计算量的增大，还可能导致舍入误差的严重积累。 经典的四阶龙格库塔公式从节点出发，首先以h为步长求近似值，nx hny1。
7、由于公式局部切割误差为，有5o，将步骤对折，即步骤跨两步，再求近似值。 由于2hnx1n 21hny的每步切割误差为，有5o的比较(3.14 )式和(3.15 )式，可见在将舞步对折后，误差约减少。 由此，容易得到以下的事后推算式. 15 () ) ) )1)2)1)如1hnhnnnyyxy那样，检查步长，根据对半分割前后2次的计算结果的偏差，判定选择的步长是否合适从表面上看，为了选择步长，每步的计算量会增加，但从整体上考虑往往很合算。