龙格库塔法RKF45Matlab实现.docx

《 龙格库塔法RKF45Matlab实现.docx》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、龙格库塔法RKF45的Matlab实现  2007-08-1614:03:32

2、  分类： MatLab/Maple/Mat

3、字号 订阅　　4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiu@uestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。rkf45.m:function[RtRx]=rkf45(f,tspan,ya,m,tol)%Input

4、:%         -f  functioncolumnvector%         -tspan[a,b]left&rightpointof[a,b]%         -ya initialvaluecolumnvector%         -m   initialguessfornumberofsteps%         -tol tolerance%Output:%         -Rt solution:vectorofabscissas%         -Rx solution:vectorofordinates%

5、ProgrambyJohn.Mathews,improvedbyliuliu@uestciflength(tspan)~=2   error('lengthofvectortspanmustbe2.');endif~isnumeric(tspan)   error('TSPANshouldbeavectorofintegrationsteps.');endif~isnumeric(ya)   error('Yashouldbeavectorofinitialconditions.');endh=diff(tspan);ifany(sign

6、(h(1))*h<=0)   error('EntriesofTSPANarenotinorder.');end a=tspan(1);b=tspan(2);ya=ya(:);a2=1/4;b2=1/4;a3=3/8;b3=3/32;c3=9/32;a4=12/13;b4=1932/2197;c4=-7200/2197;d4=7296/2197;a5=1;b5=439/216;c5=-8;d5=3680/513;e5=-845/4104;a6=1/2;b6=-8/27;c6=2;d6=-3544/2565;e6=1859/4104;f6=

7、-11/40;r1=1/360;r3=-128/4275;r4=-2197/75240;r5=1/50;r6=2/55;n1=25/216;n3=1408/2565;n4=2197/4104;n5=-1/5;big=1e15;h=(b-a)/m;hmin=h/64;%步长自适应范围下限hmax=64*h;%步长自适应范围上限max1=200;%迭代次数上限Y(1,:)=ya;T(1)=a;j=1;%tj=T(1);br=b-0.00001*abs(b);while(T(j)br),h=b-T(j);end

8、   %caculatevaluesofk1...k6,y1...y6 tj=T(j); yj=Y(j,:); y1=yj; k1=h*feval(f,tj,y1); y2=yj+b2*k1;                             ifbig

9、; y4=yj+b4*k1+c4*k2+d4*k3;              ifbig

10、*h,y6); err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6); ynew=yj+n1*k1+n3*k3+n4*k4+n5*k5; %errorandstepsi