各位大侠帮看看程序吧~~！

zizhisuole

要用MATLAB做，总是出错~希望大家帮忙看看，我是个菜鸟，谢谢！程序如下

%主程序
clear
clc
clf
Ld=0.6;
r=0.1;
lx=0.5;
p=1;
K1=0.025;
K2=0.03;
m2=10;
Rx=-40;
W=6.0;  %参数
[T,Y]=ode45('dy1',[0 2*pi/w],[1.001 0 0 0],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
[n1,mm1]=size(Y);
    j(1,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
for i=2:1000
    i
[T,Y]=ode45('dy1',[0 2*pi/w],[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
[n1,mm1]=size(Y);
j(i,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
end;
plot(j(:,1),j(:,2),'.');
xlabel('x1');ylabel('x2');
%子程序
function dy=odel(t,y,flag,w,lx,r,Ld,p,K1,K2,Rx,m2)
g=9.81;
k1=K1/(m2*lx*lx);
k2=K2/(m2*lx*lx);
kg=g/lx;
dy=[0,0,0,0];
xg=Ld-r*cos(w*t);
yg=r*sin(w*t);
xgt=r*w*sin(w*t);
ygt=r*w*cos(w*t);
xgtt=r*w^2*cos(w*t);
ygtt=-r*w^2*sin(w*t);
d=sqrt(xg*xg+yg*yg);
k=sqrt(4*lx^2-d^2);
cf1=(xg-yg*k/d)/(2*lx);
cf2=d^2/(2*lx^2)-1;
sf1=(yg+xg*k/d)/(2*lx);
sf2=d*k/(2*lx^2);
cf21=(xg+yg*k/d)/(2*lx);
sf21=(-yg+xg*k/d)/(2*lx);
f1t=-(xgt-ygt*k/d+4*lx^2*yg*(xg*xgt+yg*ygt)/(k*d^3))/(yg+xg*k/d);
f2t=-2*(xg*xgt+yg*ygt)/(d*k);
h1=d*d/(yg^4-xg^4+2*xg*yg*d*k+4*lx*lx*xg^2);
h2=(xg*ygtt-xgt*ygt)*4*lx*lx/d/d-(xg*ygtt-2*xgt*ygt+xgtt*yg)+(xgt^2-ygt^2-xg*xgtt+yg*ygtt)*k/d;
h3=(xg*xgt+yg*ygt)*((xgt*yg-xg*ygt)/d/d-(xg*xgt+yg*ygt)/d/k)-(xgt^2+xg*xgtt+ygt^2+yg*ygtt)*(xg*yg/d/d+yg^2/d/k);
h4=-16*lx*lx*(xg*yg*(d^2-2*lx^2)/d/d-yg^2*(3*lx*lx-d*d)/d/k)*(xg*xgt+yg*ygt)^2/(d^4*k^2);
f1tt=h1*(h2+4*lx*lx*h3/d/d+h4);
f2tt=-2*(2*(xg*xgt+yg*ygt)^2*(d*d-2*lx*lx)/d^2/k^2+xgt^2+xg*xgtt+ygt^2+yg*ygtt)/d/k;
a21=3*(2*kg*sf1*p+4*kg*sf1+4*Rx+3*kg*sf21*cf2)/(4*p+12-9*cf2^2);
a22=-6*(2*k1-2*f2t*sf2+2*f1t*sf2+3*f1t*sf2*cf2)/(4*p+12-9*cf2^2);
a23=(-6*f2t^2*cf2+12*f2t*cf2*f1t+6*f1tt*sf2-6*sf2*f2tt-6*f1t^2*cf2+4*Rx-9*f1t^2*cf2^2-9*kg*sf21*cf2-9*cf2*f1tt*sf2+6*cf2*Rx)/(4*p+12-9*cf2^2);
a24=-6*(2*f2t*sf2-2*f1t*sf2+3*k2*cf2)/(4*p+12-9*cf2^2);
a41=3*(2*kg*sf1*p+4*kg*sf1+6*kg*sf21+4*Rx+3*cf2*kg*sf1*p+6*cf2*kg*sf1+3*kg*sf21*cf2+6*cf2*Rx+2*kg*sf21*p)/(4*p+12-9*cf2^2);
a42=-6*(2*k1-2*f2t*sf2+3*cf2*k1-3*cf2*f2t*sf2+2*f1t*sf2*p+8*f1t*sf2+6*f1t*sf2*cf2)/(4*p+12-9*cf2^2);
a43=(-18*kg*sf21-6*f2t^2*cf2+12*f2t*cf2*f1t-12*f1tt*sf2-6*sf2*f2tt-9*kg*sf21*cf2-9*f2t^2*cf2^2+18*f2t*cf2^2*f1t-9*cf2*sf2*f2tt-6*p*f1t^2*cf2-6*kg*sf21*p-6*p*f1tt*sf2+4*p*Rx-24*f1t^2*cf2+16*Rx-18*f1t^2*cf2^2+12*cf2*Rx)/(4*p+12-9*cf2^2);
a44=-6*(6*k2+2*f2t*sf2-2*f1t*sf2+3*k2*cf2+3*cf2*f2t*sf2-3*f1t*sf2*cf2+2*k2*p)/(4*p+12-9*cf2^2);
dy=[y(2);a21*y(1)+a22*y(2)+a23*y(3)+a24*y(4);y(4);a41*y(1)+a42*y(2)+a43*y(3)+a44*y(4)];

zizhisuole

有个问题是调用子程序时出现
Input argument 'K1' is undefined.
Error in ==> d:\MATLAB6p5\work\odel.m
On line 3  ==> k1=K1/(m2*lx*lx);
是怎么回事啊，求解~！

meiyongyuandeze

首先有两个错误
1. 你程序中w写成了W。
2. 调用ODE45的语句中的函数名不对，应该是odel不是dy1
这样修该下没问题。

1. clear
   
2. clc
   
3. clf
   
4. Ld=0.6;
   
5. r=0.1;
   
6. lx=0.5;
   
7. p=1;
   
8. K1=0.025;
   
9. K2=0.03;
   
10. m2=10;
    
11. Rx=-40;
    
12. w=6.0;  %参数
    
13. [T,Y]=ode45('odel',[0 2*pi/w],[1.001 0 0 0],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
    
14. [n1,mm1]=size(Y);
    
15.     j(1,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
    
16. for i=2:1000
    
17.     i
    
18. [T,Y]=ode45('odel',[0 2*pi/w],[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)],[],w,lx,r,Ld,p,K1,K2,Rx,m2);
    
19. [n1,mm1]=size(Y);
    
20. j(i,:)=[Y(n1,1) Y(n1,2) Y(n1,3) Y(n1,4)];
    
21. end;
    
22. plot(j(:,1),j(:,2),'.');
    
23. xlabel('x1');ylabel('x2');
    
24. 
    

复制代码

附一个计算图：

zizhisuole

回复 3 # meiyongyuandeze 的帖子

谢谢你~~应经改好了@！！！