龙格库塔法RKF45Matlab实现.docxVIP

龙格库塔法RKF45的Matlab实现??

2007-08-1614:03:32|??分类：?\oMatLab/Maple/MathematicaMatLab/Maple/Mat|字号?订阅

4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiu@uestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。

rkf45.m:

function[RtRx]=rkf45(f,tspan,ya,m,tol)

%Input:

%?????????-f??functioncolumnvector

%?????????-tspan[a,b]leftrightpointof[a,b]

%?????????-ya?initialvaluecolumnvector

%?????????-m???initialguessfornumberofsteps

%?????????-tol?tolerance

%Output:

%?????????-Rt?solution:vectorofabscissas

%?????????-Rx?solution:vectorofordinates

%ProgrambyJohn.Mathews,improvedbyliuliu@uestc

iflength(tspan)~=2

???error(lengthofvectortspanmustbe2.);

end

if~isnumeric(tspan)

???error(TSPANshouldbeavectorofintegrationsteps.);

end

if~isnumeric(ya)

???error(Yashouldbeavectorofinitialconditions.);

end

h=diff(tspan);

ifany(sign(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=-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...y6

?tj=T(j);

?yj=Y(j,:);

?y1=yj;

?k1=h*feval(f,tj,y1);

?y2=yj+b2*k1;????????????????????????????

?ifbigabs(max(y2))return,end

?k2=h*feval(f,tj+a2*h,y2);

?y3=yj+b3*k1+c3*k2;?????????????????????ifbigabs(max(y3))return,end

?k3=h*feval(f,tj+a3*h,y3);

?y4=yj+b4*k1+c4*k2+d4*k3;??????????????ifbi