zizhisuole 发表于 2011-5-28 09:04

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

要用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;%参数
=ode45('dy1',,,[],w,lx,r,Ld,p,K1,K2,Rx,m2);
=size(Y);
    j(1,:)=;
for i=2:1000
    i
=ode45('dy1',,,[],w,lx,r,Ld,p,K1,K2,Rx,m2);
=size(Y);
j(i,:)=;
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=;
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=;

zizhisuole 发表于 2011-5-28 09:44

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

meiyongyuandeze 发表于 2011-5-28 11:15

首先有两个错误
1. 你程序中w写成了W。
2. 调用ODE45的语句中的函数名不对,应该是odel不是dy1
这样修该下没问题。
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;%参数
=ode45('odel',,,[],w,lx,r,Ld,p,K1,K2,Rx,m2);
=size(Y);
    j(1,:)=;
for i=2:1000
    i
=ode45('odel',,,[],w,lx,r,Ld,p,K1,K2,Rx,m2);
=size(Y);
j(i,:)=;
end;
plot(j(:,1),j(:,2),'.');
xlabel('x1');ylabel('x2');

附一个计算图:


zizhisuole 发表于 2011-5-30 13:45

回复 3 # meiyongyuandeze 的帖子

谢谢你~~应经改好了@!!!
页: [1]
查看完整版本: 各位大侠帮看看程序吧~~!