|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 煜宸0922 于 2011-5-7 10:14 编辑
function dxdt=fun(t,x)
k1=2;k2=3;k3=0.05;k4=20;k5=1;k6=1;k7=3;
dxdt=zeros(4,1);
dxdt(1)=x(3);
dxdt(2)=x(4);
dxdt(3)=(1-k1)*sin(k2*t+k4)-2*k3*x(3)+2*k3*x(4)-x(1)+x(2);
dxdt(4)=(k1*sin(k2*t+k4)+2*k3*x(3)-2*k3*(1+k5)*x(4)+x(1)-(1+k6)*x(2))/k7;
end
这是所要求数值解的方程组程序,不用ode45,请教各位怎样编写龙格库塔法解这个方程组的程序?万急,求救!!!
我编的程序老是出错,请各位大侠指正,不胜感激,程序如下:
function R=rk41(f,a,b,xa,xb,M)
h=0.01;
x=zeros(1,M+1);
t=a:h:b;
x(3)=xa;
x(4)=xb;
k11=x(3);
k21=x(4);
k31=feval(f,t,x(1),x(3));
k41=feval(f,t,x(2),x(4));
k12=x(3)+h/2*k31;
k22=x(4)+h/2*k41;
k32=feval(f,t+h/2,x(1)+h/2*k11,x(3)+h/2*k31);
k42=feval(f,t+h/2,x(2)+h/2*k21,x(4)+h/2*k41);
k13=x(3)+h/2*k32;
k23=x(4)+h/2*k42;
k33=feval(f,t+h/2,x(1)+h/2*k12,x(3)+h/2*k32);
k43=feval(f,t+h/2,x(2)+h/2*k22,x(4)+h/2*k42);
k14=x(3)+h*k33;
k24=x(4)+h*k43;
k34=feval(f,t+h,x(1)+h*k13,x(3)+h*k33);
k44=feval(f,t+h,x(2)+h*k23,x(4)+h*k43);
x(12)=x(1)+h/6*(k11+2*k12+2*k13+k14);
x(22)=x(2)+h/6*(k21+2*k22+2*k23+k24);
x(32)=x(3)+h/6*(k31+2*k32+2*k33+k34);
x(42)=x(4)+h/6*(k41+2*k42+2*k43+k44);
end
这程序是想只求解一个点,等正确了再循环,可还是错的,在下matlab刚入门,十分幼稚的错误在所难免,只求大侠们给予纠错帮助。临表涕零,不胜感激。
|
|