悬臂梁的模态分析【传递矩阵法】

悬臂梁受定常轴向力的作用，将悬臂梁分成50段，每一段悬臂梁的传递矩阵方程：Z(i+1)=[tr1]i*[tr2]i* Z(i);  其中[tr1]i=[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)*cos(a1*l)  (E*I*a1^3)*sin(a1*l) (E*I*a2^3)*cosh(a2*l)  (E*I*a2^3)*sinh(a2*l)];[tr2]i=(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];那么总传递矩阵是不是 trr=(tr1.*tr2)^50;  不知理解正确否，请各位大虾指点。。如果正确，我按照这种思路编了相应的程序：clcclearsyms E I a1 a2 m l s omiga f pi tr1 tr2  trr trr1 ;%文献悬臂梁的参数E=2.1856e11; I=5e-8;l=0.44;m=3/l;s=1e4;omiga=2*pi*f;%传递矩阵中的各个参数a1=((sqrt(s^2+4*m*omiga^2*E*I)-s)/(2*E*I))^(1/2);a2=((sqrt(s^2+4*m*omiga^2*E*I)+s)/(2*E*I))^(1/2);%由边界条件得到的总传递矩阵（这里只是trr=tr1.*tr2，个人感觉trr的50次与一次得到的特征方程的解相差不大）trr=[0, a2*cosh(a2*l)*(a1^16 + a2^16);            -(cosh(a2*l)*(- E*I*a2^3 + s*a2)*(a1^16 + a2^16))/(E*I), 0];      trr1=det(trr);       solve(vpa(trr1,2))得到的结果为 84.071，与文献中116.267相差很大，不知道是怎么回事，，望有能力的大虾们指点迷津，非常感激