求汽轮机叶片的动频程序【采用传递矩阵法】

华电机械

求各位matlab高手给看看我的程序，帮忙求下trr1=0 这个方程【其中只有一个未知数f】，最后结果就是要得到频率。。。求助。
       在此感激啊！！
clc
clear
syms E I a1 a2 m l s omiga f pi tr1 tr2  trr trr1 ;

%等截面叶片参数
l=0.328;
b=0.028;
t=0.003;
A=b*t;
R=0.15;
%汽轮机直叶片材料参数（忽略叶片的质量）
u=0.3;
rou=7850;
E=2.17e11;
I=b*t^3/12;
m=rou*A;
s=rou*A*(R+l)*l*314.15926^2/9.8;  %叶片的离心力
omiga=2*3.14*f;
a1=sqrt((sqrt(s^2+4*m*omiga^2*E*I)-s)/(2*E*I));
a2=sqrt((sqrt(s^2+4*m*omiga^2*E*I)+s)/(2*E*I));

%计算总传递矩阵
tr1=[sin(a1*l) cos(a1*l) sinh(a2*l) cosh(a2*l);
    a1*cos(a1*l) -a1*sin(a1*l) a2*cosh(a2*l) a2*sinh(a2*l);
    -E*I*a1^2*sin(a1*l) -E*I*a1^2*cos(a1*l) E*I*a2^2*sinh(a2*l) E*I*a2^2*cosh(a2*l);
    (-E*I*a1^3-s*a1)*cos(a1*l)  (E*I*a1^3+s*a1)*sin(a1*l) (E*I*a2^3-s*a2)*cosh(a2*l)  (E*I*a2^3+s*a2)*sinh(a2*l)];

%文献的矩阵
tr2=(a1^2+a2^2).*[0 a2^2/a1 0 -1/(E*I*a1);a2^2 0 -1/(E*I) 0; 0 a1^2/a2 0 1/(E*I*a2); a1^2 0 1/(E*I) 0];
             tr=tr1.*tr2;   %总传递矩阵

trr=[0, a2*cosh(a2*l)*(a1^2 + a2^2);
           -(cosh(a2*l)*(- E*I*a2^3 + s*a2)*(a1^2 + a2^2))/(E*I), 0];
      trr1=det(trr);   %特征方程
       vpa(trr1,2)