龙格库塔法求助（万急）

煜宸0922

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刚入门，十分幼稚的错误在所难免，只求大侠们给予纠错帮助。临表涕零，不胜感激。

meiyongyuandeze

给个出错信息！

煜宸0922

回复 2 # meiyongyuandeze 的帖子

??? Error using ==> fun
Too many input arguments.

Error in ==> rk41 at 9
    k31=feval(f,t,x(1),x(3));

meiyongyuandeze

回复 3 # 煜宸0922 的帖子

首先你定义的fun函数只有两个输入量（t,x），所以k31=feval(f,t,x(1),x(3));这句是错的，

煜宸0922

回复 4 # meiyongyuandeze 的帖子

呵呵，哪我该怎么改一下才行？k31那一句要用的是方程组的第三个方程，大侠指点下吧！

meiyongyuandeze

1. function dxdt=fun(t,x)
   
2. k1=2;k2=3;k3=0.05;k4=20;k5=1;k6=1;k7=3;
   
3. dxdt=zeros(4,1);
   
4. dxdt(1)=x(3);
   
5. dxdt(2)=x(4);
   
6. dxdt(3)=(1-k1)*sin(k2*t+k4)-2*k3*x(3)+2*k3*x(4)-x(1)+x(2);
   
7. dxdt(4)=(k1*sin(k2*t+k4)+2*k3*x(3)-2*k3*(1+k5)*x(4)+x(1)-(1+k6)*x(2))/k7;

复制代码

1. clc
   
2. clear
   
3. k=0.001;
   
4. T=0:0.001:10;
   
5. X(:,1)=[0.005 0 0 0]';
   
6. for j=1:length(T)
   
7. t(j,1)=T(j);
   
8. k1=k*feval('fun',T(j),X(:,j));
   
9. k2=k*feval('fun',T(j)+k/2,X(:,j)+k1/2);
   
10. k3=k*feval('fun',T(j)+k/2,X(:,j)+k2/2);
    
11. k4=k*feval('fun',T(j)+k,X(:,j)+k3);
    
12. X(:,(j+1))=X(:,j)+(k1+2*k2+2*k3+k4)/6;
    
13. end
    
14. plot(T,X(1,2:end))

复制代码

比较忙，实在是不想找错误，就帮你直接编写了，结果：

meiyongyuandeze

回复 5 # 煜宸0922 的帖子

帮你写了下程序，你看下吧，图不知道是不是你想要的!

煜宸0922

回复 7 # meiyongyuandeze 的帖子

十分感谢大侠，我是要用这个画相图最后根据初值和其他参数画分岔图和混沌。我本来觉得，要每个方程都有四个斜率值，四个方程得十六个斜率值，而对每个方程算对应的斜率我就不会调了，十分感谢大侠百忙之中给予的帮助，以后又问题还得向大侠请教呢，希望大侠到时不吝赐教。

煜宸0922

回复 8 # 煜宸0922 的帖子

对了，大侠，你能解释下你画图的那个命令里面的东西吗？

煜宸0922

回复 7 # meiyongyuandeze 的帖子

对了，大侠，如果精度不够，要采用变步长的方法，程序该怎么调一下？

meiyongyuandeze

回复 10 # 煜宸0922 的帖子

那改动还是很大的，你可以参考数值计算相关的书来编写！

煜宸0922

回复 11 # meiyongyuandeze 的帖子

麻烦大侠给推荐些书

煜宸0922

回复 2 # meiyongyuandeze 的帖子

大侠，还有个小问题，要是加入了边界条件，该怎么变一下程序？边界条件是：x(1)=b,x(1)=-b,x'(a1)=-Rx'(b1)，这里面，x'(a1)表示碰撞前的速度，x'(b1)表示碰撞后的速度。请大侠指点一二。

meiyongyuandeze

回复 13 # 煜宸0922 的帖子

应该是在function dxdt=fun(t,x)中加入判定语句！

煜宸0922

回复 14 # meiyongyuandeze 的帖子

求大侠详细解释！万急！！！！