利用传递矩阵法计算叶片的动频。。。求指点

华电机械

其实这个问题我的动机就是 要计算汽轮机叶片的动频，，参考了文献【利用Euler梁模型计算叶片静频和动频的传递矩阵法】，根据文献编制了相应求解的程序，可是求不出，求指教。。
现在不知道是程序中叶片转动的离心力计算有问题，还是传递矩阵有问题，结果f 就是求不出，，求指点啊！由于叶片分了110段，感觉得110个矩阵相乘得到它的总传递矩阵。
%等截面叶片参数
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*pi*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));
for i=1:50;
%传递矩阵1
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)];
%传递矩阵2
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)
求高手指点。。。。