怎么实现矩阵相乘的啊！总共50个矩阵连乘

华电机械

编程代码:
%计算总传递矩阵
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;
for i=1:50;
tr（i）=tr(i).*tr(i);
end
想要把tr连乘50次,怎么实现啊！ 求高手指点。。

ChaChing

水平有限, 实在不清楚为何一定得用符号运算!?

华电机械

ChaChing 发表于 2013-4-26 21:58

                               
登录/注册后可看大图

水平有限, 实在不清楚为何一定得用符号运算!?

那用什么方法能实现呢？求指点。。

ChaChing

1.LZ未给齐完整代码!? 这边复製贴上matlab即有报错
2.水平有限, 不清楚为何一定得用符号运算, 不能使用数值运算!?
3.a连乘50次, 不就是a^50!?

华电机械

ChaChing 发表于 2013-4-27 10:37

                               
登录/注册后可看大图

1.LZ未给齐完整代码!? 这边复製贴上matlab即有报错
2.水平有限, 不清楚为何一定得用符号运算, 不能使用数值 ...

1.恩，那个矩阵那样不能运行。。
2 . 请问数值运算是怎样的啊
3. 由于矩阵比较大，用a^50运算好像是运行不了，以下试了一下，的确不行，，求指教。。
syms E I a1 a2 m l s omiga f pi tr1 tr2  trr  ;
tr =[ 0, (a2^2*cos(a1*l)*(a1^2 + a2^2))/a1, 0, -(cosh(a2*l)*(a1^2 + a2^2))/(E*I*a1);
      a1*a2^2*cos(a1*l)*(a1^2 + a2^2),  0,   -(a2*cosh(a2*l)*(a1^2 + a2^2))/(E*I) ,0;
      0, -(E*I*a1^4*cos(a1*l)*(a1^2 + a2^2))/a2, 0, a2*cosh(a2*l)*(a1^2 + a2^2);
-a1^2*cos(a1*l)*(E*I*a1^3 + s*a1)*(a1^2 + a2^2), 0, -(cosh(a2*l)*(- E*I*a2^3 + s*a2)*(a1^2 + a2^2))/(E*I),0];
trr=tr^50

ChaChing

华电机械 发表于 2013-4-27 11:07

                               
登录/注册后可看大图

1.恩，那个矩阵那样不能运行。。
2 . 请问数值运算是怎样的啊
3. 由于矩阵比较大，用a^50运算好像是运行 ...

1.养成好习惯,给齐完整代码! 不然别人不好复製状况
2.矩阵不大, 才4*4, 数值运算a^50不超过一秒, 但符号运算才a^5就busy了, 所以个人才会有疑问, 为何一定得用符号运算!?
3."请问数值运算是怎样的啊!" 难道术语不同?

华电机械

ChaChing 发表于 2013-4-27 11:18

                               
登录/注册后可看大图

1.养成好习惯,给齐完整代码! 不然别人不好复製状况
2.矩阵不大, 才4*4, 数值运算a^50不超过一秒, 但符号 ...

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

华电机械

ChaChing 发表于 2013-4-27 11:18

                               
登录/注册后可看大图

1.养成好习惯,给齐完整代码! 不然别人不好复製状况
2.矩阵不大, 才4*4, 数值运算a^50不超过一秒, 但符号 ...

您能给我好好看看上面的东西吗？ 这几天的确弄的头都大了，在此，非常的感激您，，求帮忙！！现在感觉个人的能力很有限啊！

华电机械

华电机械 发表于 2013-4-27 11:28

                               
登录/注册后可看大图

其实这个问题我的动机就是 要计算汽轮机叶片的动频啊，，参考了文献【利用Euler梁模型计算叶片静频和动频 ...

f是频率，恩，好像是少了个end，但结果还是不行

321forever

华电机械 发表于 2013-4-27 15:49

                               
登录/注册后可看大图

f是频率，恩，好像是少了个end，但结果还是不行

你给的程序里面 f 没有赋值，缺的end要放在哪里

华电机械

321forever 发表于 2013-4-27 19:23

                               
登录/注册后可看大图

你给的程序里面 f 没有赋值，缺的end要放在哪里

现在的问题是要求频率f 啊，也就是求方程trr1=0的解 ，里面就f 一个未知数。。。怎么求，知道吗？

ChaChing

华电机械 发表于 2013-4-27 11:28

                               
登录/注册后可看大图

其实这个问题我的动机就是 要计算汽轮机叶片的动频啊，，参考了文献【利用Euler梁模型计算叶片静频和动频 ...

水平/时间有限, 这已经不是编程问题了! ~
建议看看陈教授的书(Ch6.7) http://forum.chinavib.com/thread-99956-1-1.html
我相信一定还有许多其他资料可参考, 花些苦功google下吧
毕竟我没真的玩过

华电机械

ChaChing 发表于 2013-4-29 09:13

                               
登录/注册后可看大图

水平/时间有限, 这已经不是编程问题了! ~
建议看看陈教授的书(Ch6.7) http://forum.chinavib.com/thread ...

好的，，谢谢您啊。。

mxlzhenzhu

1、在这里定义符号变量；
s=rou*A*(R+l)*l*314.15926^2/9.8;
syms f;
omiga=2*pi*f;
2、最后添加一个end
for &&&

vpa()
end

3、没看到特征方程啊我？

4、外行弱弱地问一下，有通用有限元了，干嘛还要用传递矩阵法呢？

谢谢。

最后乱修改了一下，不晓得你的本意，你看看如何吧。 delete_this_file.txt (1.02 KB, 下载次数: 2)

2013-4-30 16:51 上传

点击文件名下载附件

后缀直接改为.m就可以运行了。

华电机械

mxlzhenzhu 发表于 2013-4-30 16:45

                               
登录/注册后可看大图

1、在这里定义符号变量；
s=rou*A*(R+l)*l*314.15926^2/9.8;
syms f;

1.运行了一下你改的那 个，得到的答案是f=0，明显不对啊，可能是我的理解有错了。
2 特征方程就是trr1=0，，解这个方程的解，也就是频率了
3 由于鄙人正在学传递矩阵法，还没怎么接触有限元法，可能有限元法求叶片的频率跟简单，但现在就是想试试传递矩阵法应该也能求吧。
4 谢谢你的帮助，非常感激。。。