[求助] 直接传递矩阵法求双转子系统临界转速报错Conversion to logical from sym

rou=7754;
E=2.06e11;
%L1=1.20;内转轴的总长
%L2=0.448;外转轴的总长
K1=1.236e6; K2=2.472e6;K3=4.025e6;K4=1.342e6;%弹性支撑刚度
L1=[0.11 0.04 0.11 0.11 0.11 0.11 0.11 0.130 0.068 0.106 0.146 0.05]; %内转轴的分段
L2=[0.09 0.05 0.11 0.06 0.06 0.06 0.018];%外转轴的分段
r1=0.01;r2=0.018;R2=0.0225;%内外轴半径
rp1=0.09;h1=0.03;rp2=0.115;h2=0.018;%两轮盘半径及厚度
J1=pi*(r1^4)/4;J2=pi*(R2^4-r2^4)/4;%内外轴截面惯性矩
Ip1=rou*pi*rp1^4/2;Id1=Ip1/2;Ip2=rou*pi*rp2^4/2;Id2=Ip2/2;%内外轴轮盘极转动惯量及直径转动惯量
M1=zeros(1,length(L1)+1);M2=zeros(1,length(L2)+1);
for jk1=1:length(L1)
m1(jk1)=rou*pi*r1^2*L1(jk1)/2 ;%内转轴的质量分布
end
for jk2=1:length(L2)
m2(jk2)=rou*pi*(R2^2-r1^2)*L2(jk2)/2 ;%外转轴的质量分布
end
M1(:,2:length(L1))=m1(:,1:length(L1)-1)+m1(:,2:length(L1));
M1(1)=m1(1);
M1(length(L1)+1)=m1(length(L1));
M2(:,2:length(L2))=m2(:,1:length(L2)-1)+m2(:,2:length(L2));
M2(1)=m2(1);
M2(length(L2)+1)=m2(length(L2));
M1(10)=M1(10)+rou*pi*(rp1^2-r1^2)*h1;%低压轴轮盘的等效质量
M2(4)=M1(4)+rou*pi*(rp2^2-R2^2)*h2;%高压轴轮盘的等效质量

k=0;
Tit=['第一阶频率的振型和弯矩图';'第二阶频率的振型和弯矩图';'第三阶频率的振型和弯矩图'];
syms w;%for w=0:0.01:2000;
for i1=1:length(L1)
T1(:,:,i1)=[1 L1(i1) L1(i1)^2/(2*E*J1) L1(i1)^3/(6*E*J1);0 1 L1(i1)/(E*J1) L1(i1)^2/(2*E*J1);0 0 1 L1(i1);0 0 0 1]*[1 0 0 0;0 1 0 0;0 0 1 0;M1(i1)*w^2 0 0 1];
T2(:,:,i1)=[1 0 0 0;0 1 0 0;0 0 1 0;M1(i1+1)*w^2 0 0 1]*[1 L1(i1) L1(i1)^2/(2*E*J1) L1(i1)^3/(6*E*J1);0 1 L1(i1)/(E*J1) L1(i1)^2/(2*E*J1);0 0 1 L1(i1);0 0 0 1];
T3=[1 L1(3) L1(3)^2/(2*E*J1) L1(3)^3/(6*E*J1);0 1 L1(3)/(E*J1) L1(3)^2/(2*E*J1);0 0 1 L1(3);0 0 0 1]*[1 0 0 0;0 1 0 0;0 0 1 0;(M1(3)*w^2-K1) 0 0 1];
T4=[1 0 0 0;0 1 0 0;0 0 1 0;(M1(12)*w^2-K2) 0 0 1]*[1 L1(11) L1(11)^2/(2*E*J1) L1(11)^3/(6*E*J1);0 1 L1(11)/(E*J1) L1(11)^2/(2*E*J1);0 0 1 L1(11);0 0 0 1];
T5=[1 0 0 0;0 1 0 0;0 (Ip1*w^2/1.2-Id1*w^2) 1 0;(M1(10)*w^2-K4) 0 0 1];%Hn1R(n1L)
end
H1=T1(:,:,4);
for i2=5:9;
H1=T1(:,:,i2)*H1;
end
H10L=H1*T3*T1(:,:,2)*T1(:,:,1);
H10R=T2(:,:,12)*T4*T2(:,:,10);
Hd1=H10R*T5*H10L;%Hn2R(1L)
Hd2=H10R*T5;%Hn2R(n1L)
for i3=1:length(L2)
T6(:,:,i3)=[1 L2(i3) L2(i3)^2/(2*E*J2) L2(i3)^3/(6*E*J2);0 1 L2(i3)/(E*J2) L2(i3)^2/(2*E*J2);0 0 1 L2(i3);0 0 0 1]*[1 0 0 0;0 1 0 0;0 0 1 0;M2(i3)*w^2 0 0 1];
T7=[1 L2(2) L2(2)^2/(2*E*J2) L2(2)^3/(6*E*J2);0 1 L2(2)/(E*J2) L2(2)^2/(2*E*J2);0 0 1 L2(2);0 0 0 1]*[1 0 0 0;0 1 0 0;0 0 1 0;(M2(2)*w^2-K3) 0 0 1];
T8=[1 L2(4) L2(4)^2/(2*E*J2) L2(4)^3/(6*E*J2);0 1 L2(4)/(E*J2) L2(4)^2/(2*E*J2);0 0 1 L2(4);0 0 0 1]*[1 0 0 0;0 1 0 0;0 (Ip2*w^2-Id2*w^2) 1 0;M2(4)*w^2 0 0 1];
T9(:,:,i3)=[1 0 0 0;0 1 0 0;0 0 1 0;M2(i3+1)*w^2 0 0 1]*[1 L2(i3) L2(i3)^2/(2*E*J2) L2(i3)^3/(6*E*J2);0 1 L2(i3)/(E*J2) L2(i3)^2/(2*E*J2);0 0 1 L2(i3);0 0 0 1];
T10=[1 0 0 0;0 1 0 0;0 0 1 0;(M2(7)*w^2-K4) 0 0 1];
end
H20L=T6(:,:,6)*T6(:,:,5)*T8*T6(:,:,3)*T7*T6(:,:,1);
syms x;syms y;
Zn3L=[x,y,0,0]';
Zn4L=H20L*Zn3L;
syms x10;
Zn5R=T9(:,:,7)*T10*(Zn4L+[0 0 0 K4*x10]');

%%
X1=Zn5R(3,:);
X2=Zn5R(4,:);
S=solve(X1,X2,x,y);
U=S.x/x10;
V=S.y/x10;
a=H20L(1,1)*U+H20L(1,2)*V;
G1=Hd1(3,1)+a*Hd2(3,4)*K4*H10L(1,1);
G2=Hd1(3,2)+a*Hd2(3,4)*K4*H10L(1,2);
G3=Hd1(4,1)+a*Hd2(4,4)*K4*H10L(1,1);
G4=Hd1(4,2)+a*Hd2(4,4)*K4*H10L(1,2);
F=G1*G4-G2*G3 ; %剩余量
for w=0:0.01:2000;
if F*(-1)^k<0 %求解临界转速
k=k+1;
        wi(k)=w ;    %固有圆频率
        w=wi(k);
         ni(k)=wi(k)*30/pi ; %临界转速

end
end
ni=ni'  
wi=wi'
f=ni/60