华电机械 发表于 2013-4-26 15:32

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

编程代码:
%计算总传递矩阵
tr1=;tr2=(a1^2+a2^2).*;tr=tr1.*tr2;
for i=1:50;
tr(i)=tr(i).*tr(i);
end
想要把tr连乘50次,怎么实现啊! 求高手指点。。



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

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

华电机械 发表于 2013-4-27 10:21

ChaChing 发表于 2013-4-26 21:58 static/image/common/back.gif
水平有限, 实在不清楚为何一定得用符号运算!?

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

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

本帖最后由 ChaChing 于 2013-4-27 10:38 编辑

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

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

ChaChing 发表于 2013-4-27 10:37 static/image/common/back.gif
1.LZ未给齐完整代码!? 这边复製贴上matlab即有报错
2.水平有限, 不清楚为何一定得用符号运算, 不能使用数值 ...

1.恩,那个矩阵那样不能运行。。
2 . 请问数值运算是怎样的啊
3. 由于矩阵比较大,用a^50运算好像是运行不了,以下试了一下,的确不行,,求指教。。
syms E I a1 a2 m l s omiga f pi tr1 tr2trr;
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:18

华电机械 发表于 2013-4-27 11:07 static/image/common/back.gif
1.恩,那个矩阵那样不能运行。。
2 . 请问数值运算是怎样的啊
3. 由于矩阵比较大,用a^50运算好像是运行 ...

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

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

ChaChing 发表于 2013-4-27 11:18 static/image/common/back.gif
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).*;
%总传递矩阵
             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)

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

ChaChing 发表于 2013-4-27 11:18 static/image/common/back.gif
1.养成好习惯,给齐完整代码! 不然别人不好复製状况
2.矩阵不大, 才4*4, 数值运算a^50不超过一秒, 但符号 ...

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

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

华电机械 发表于 2013-4-27 11:28 static/image/common/back.gif
其实这个问题我的动机就是 要计算汽轮机叶片的动频啊,,参考了文献【利用Euler梁模型计算叶片静频和动频 ...

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

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

华电机械 发表于 2013-4-27 15:49 static/image/common/back.gif
f是频率,恩,好像是少了个end,但结果还是不行

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

华电机械 发表于 2013-4-28 08:40

321forever 发表于 2013-4-27 19:23 static/image/common/back.gif
你给的程序里面 f 没有赋值,缺的end要放在哪里

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

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

华电机械 发表于 2013-4-27 11:28 static/image/common/back.gif
其实这个问题我的动机就是 要计算汽轮机叶片的动频啊,,参考了文献【利用Euler梁模型计算叶片静频和动频 ...

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

华电机械 发表于 2013-4-29 09:38

ChaChing 发表于 2013-4-29 09:13 static/image/common/back.gif
水平/时间有限, 这已经不是编程问题了! ~
建议看看陈教授的书(Ch6.7) http://forum.chinavib.com/thread ...

好的,,谢谢您啊。。

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

本帖最后由 mxlzhenzhu 于 2013-4-30 16:51 编辑

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、外行弱弱地问一下,有通用有限元了,干嘛还要用传递矩阵法呢?

谢谢。

最后乱修改了一下,不晓得你的本意,你看看如何吧。
后缀直接改为.m就可以运行了。

华电机械 发表于 2013-4-30 19:04

mxlzhenzhu 发表于 2013-4-30 16:45 static/image/common/back.gif
1、在这里定义符号变量;
s=rou*A*(R+l)*l*314.15926^2/9.8;
syms f;


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