龙格库塔法RKF45Matlab实现

龙格库塔法RKF45的Matlab实现2007-08-1614:03:32|分类：MatLab/Maple/Mat|字号订阅4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiu@uestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。rkf45.m:function[RtRx]=rkf45(f,tspan,ya,m,tol)%Input:%-ffunctioncolumnvector%-tspan[a,b]left&rightpointof[a,b]%-yainitialvaluecolumnvector%-minitialguessfornumberofsteps%-toltolerance%Output:%-Rtsolution:vectorofabscissas%-Rxsolution:vectorofordinates%ProgrambyJohn.Mathews,improvedbyliuliu@uestciflength(tspan)~=2error('lengthofvectortspanmustbe2.');endif~isnumeric(tspan)error('TSPANshouldbeavectorofintegrationsteps.');endif~isnumeric(ya)error('Yashouldbeavectorofinitialconditions.');endh=diff(tspan);ifany(sign(h(1))*h=0)error('EntriesofTSPANarenotinorder.');enda=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=-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)b),if((T(j)+h)br),h=b-T(j);end%caculatevaluesofk1...k6,y1...y6tj=T(j);yj=Y(j,:);y1=yj;k1=h*feval(f,tj,y1);y2=yj+b2*k1;ifbigabs(max(y2))return,endk2=h*feval(f,tj+a2*h,y2);y3=yj+b3*k1+c3*k2;ifbigabs(max(y3))return,endk3=h*feval(f,tj+a3*h,y3);y4=yj+b4*k1+c4*k2+d4*k3;ifbigabs(max(y4))return,endk4=h*feval(f,tj+a4*h,y4);y5=yj+b5*k1+c5.*k2+d5*k3+e5*k4;ifbigabs(max(y5))return,endk5=h*feval(f,tj+a5*h,y5);y6=yj+b6*k1+c6.*k2+d6*k3+e6*k4+f6*k5;ifbigabs(max(y6))return,endk6=h*feval(f,tj+a6*h,y6);err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);ynew=yj+n1*k1+n3*k3+n4*k4+n5*k5;%errorandstepsizecontrolif((errtol)|(h2*hmin)),Y(j+1,:)=ynew;if((tj+h)br),T(j+1)=b;elseT(j+1)=tj+h;endj=j+1;tj=T(j);endif(max(err)==0),s=0;elses1=0.84*(tol.*h./err).^(0.25);%最佳步长值s=min(s1);endif((s0.75)&(h2*hmin)),h=h/2;endif((s1.50)&(2*hhmax)),h=2*h;endif((bigabs(Y(j,:)))|(max1==j)),return,endend%[RtRx]=[T'Y];Rt=T';Rx=Y;使用方法：首先编写方程（组）文件（注意与ode45不同，这儿方程组为1Xn数组：functiondx=fun(t,x)dx=zeros(1,2);dx(1)=x(1)+x(2)*2+0*t;dx(2)=3*x(1)+x(2)*2+0*t;然后使用：[Rt,Rx]=rkf45(@fun,[0,0.2],[6;4],100,1e-7)