clc

clear

a=2/3;

u=0.3;

rou=7800;

E=2.0e11;

G=E/(2*(1+u));

%参数的数组形式

D=[0.0325 0.0325 0.0325 0.0325 0.037 0.035 0.0325 0.0325 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0286 0.0255 0.0252 0.0252];%轴的外径

d=[0.0245 0.016 0.016 0.016 0.016 0.016 0.016 0.014 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012];%轴的内径

L=[0.03 0.0305 0.0325 0.0075 0.012 0.0008 0.01 0.0008 0.0003 0.012 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.013 0.043 0.006 0.024 0.05];%轴段长

A=pi.*(D.^2-d.^2)/4;

m=rou.*A.*L;

bearingL=0.03;

bearingD=0.08;

bearingd=0.0325;

bearingM=rou*pi*(bearingD^2-bearingd^2)/4*bearingL;

zhuanziL=0.02;

zhuanziD=0.07;

zhuanzid=0.06;

zhuanziM=rou*pi.*(zhuanziD.^2-zhuanzid.^2)/4*zhuanziL;

M=m+[0 0 bearingM 0 0 0 0 0 0 0 zhuanziM zhuanziM zhuanziM zhuanziM zhuanziM zhuanziM zhuanziM 0 0 0 0 0];

Jp=rou.*pi.*(D.^2+d.^2).^2.*L/2;

Jd=Jp/2;

I=pi.*(D.^4-d.^4)/64;

V=6*E.*I./(a*G.*A.*L.^2);

K1=1.5e8;

K=[0 0 K1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 K1 0];

k=0;

FF=1;

for w=0:0.01:2000;

for i=1:22;

T(:,:,i)=[1+(L(i)^3)*(1-V(i))*M(i)*w^2-K(i)/6*E*I(i)   L(i)+L(i)^2*[(Jp(i)-Jd(i))*w^2+K(i)]/(2*E*I(i))     L(i)^2/(2*E*I(i))    L(i)^3*(1-V(i))/(6*E*I(i));

(L(i)^2)*(M(i)*w^2-K(i))                          1+(L(i)^2)*[(Jp(i)-Jd(i))*w^2+K(i)]/(E*I(i))         L(i)/(E*I(i))          L(i)^2/(2*E*I(i));

(M(i)*w^2-K(i))*L(i)                               (Jp(i)-Jd(i))*w^2+K(i)                                   1                        L(i);

M(i)*w^2-K(i)                                       0                                                         0                        1                               ];

end

H=T(:,:,1);

for i2=2:22;

H=T(:,:,i2)*H;

end

F=H(3,1)*H(4,2)-H(3,2)*H(4,1)

if F*FF< 0 %求解临界转速

k=k+1;

wi(k)=w;            %角速度 rad/s

f=wi(k)/(2*pi)      %频率   HZ

ni(k)=wi(k)*30/pi;  %转速   转/min

end

FF=F;

end

着急用,感谢大家了.帮我看看,错的话错在哪里,

我运行了下不对!!