龙格库塔法求助（万急）
本帖最后由 煜宸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刚入门，十分幼稚的错误在所难免，只求大侠们给予纠错帮助。临表涕零，不胜感激。
给个出错信息！ 回复 2 # meiyongyuandeze 的帖子
??? Error using ==> fun
Too many input arguments.
Error in ==> rk41 at 9
k31=feval(f,t,x(1),x(3));
回复 3 # 煜宸0922 的帖子
首先你定义的fun函数只有两个输入量（t,x），所以k31=feval(f,t,x(1),x(3));这句是错的，
回复 4 # meiyongyuandeze 的帖子
呵呵，哪我该怎么改一下才行？k31那一句要用的是方程组的第三个方程，大侠指点下吧！ 本帖最后由 meiyongyuandeze 于 2011-5-7 11:02 编辑
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;
clc
clear
k=0.001;
T=0:0.001:10;
X(:,1)=';
for j=1:length(T)
t(j,1)=T(j);
k1=k*feval('fun',T(j),X(:,j));
k2=k*feval('fun',T(j)+k/2,X(:,j)+k1/2);
k3=k*feval('fun',T(j)+k/2,X(:,j)+k2/2);
k4=k*feval('fun',T(j)+k,X(:,j)+k3);
X(:,(j+1))=X(:,j)+(k1+2*k2+2*k3+k4)/6;
end
plot(T,X(1,2:end))比较忙，实在是不想找错误，就帮你直接编写了，结果：
回复 5 # 煜宸0922 的帖子
帮你写了下程序，你看下吧，图不知道是不是你想要的! 回复 7 # meiyongyuandeze 的帖子
十分感谢大侠，我是要用这个画相图最后根据初值和其他参数画分岔图和混沌。我本来觉得，要每个方程都有四个斜率值，四个方程得十六个斜率值，而对每个方程算对应的斜率我就不会调了，十分感谢大侠百忙之中给予的帮助，以后又问题还得向大侠请教呢，希望大侠到时不吝赐教。 回复 8 # 煜宸0922 的帖子
对了，大侠，你能解释下你画图的那个命令里面的东西吗？ 回复 7 # meiyongyuandeze 的帖子
对了，大侠，如果精度不够，要采用变步长的方法，程序该怎么调一下？ 回复 10 # 煜宸0922 的帖子
那改动还是很大的，你可以参考数值计算相关的书来编写！ 回复 11 # meiyongyuandeze 的帖子
麻烦大侠给推荐些书 回复 2 # meiyongyuandeze 的帖子
大侠，还有个小问题，要是加入了边界条件，该怎么变一下程序？边界条件是：x(1)=b,x(1)=-b,x'(a1)=-Rx'(b1)，这里面，x'(a1)表示碰撞前的速度，x'(b1)表示碰撞后的速度。请大侠指点一二。 回复 13 # 煜宸0922 的帖子
应该是在function dxdt=fun(t,x)中加入判定语句！
回复 14 # meiyongyuandeze 的帖子
求大侠详细解释！万急！！！！
页:
[1]
2