欢迎来到天天文库
浏览记录
ID:60847769
大小:10.80 KB
页数:4页
时间:2020-12-22
《 龙格库塔法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; ifbig9、; y4=yj+b4*k1+c4*k2+d4*k3; ifbig10、*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
9、; y4=yj+b4*k1+c4*k2+d4*k3; ifbig10、*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
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
此文档下载收益归作者所有