关于矩阵太大导致求解刚度和质量矩阵时出错的问题
我编了一段简单的程序，求解一个质量矩阵和一个刚度矩阵。两个矩阵都能求解，但是用刚度矩阵除以质量矩阵就报错。有人热心指正小弟感激不尽。syms x y db dv dc hb hv hc b L;
y=x/L;
%求Nub与其转置地积分q
Nub=;
T='*;
q=int(T,0,1);
%质量矩阵
=db*hb*b*q;
Nw=;
%对Nw求导
x=sym('x');
dNw=diff(Nw)
%求Nw与其转置的积分qq
TT=Nw'*Nw;
qq=int(TT,0,1);
%质量矩阵
=db*hb*b*qq;
=;
=-(hv+hb/2+hc/2)*dNw++hv*;
=-(hb+hc)/2*dNw++hv/2+;
%质量矩阵&
qqq=int('*,0,1);
=dv*hv*b*qqq;
=dv*hv*b*qq;
%求Nuc的转置与其的积分qqqq
TTT='*;
qqqq=int(TTT,0,1);
%质量矩阵&
=dc*hc*b*qqqq;
=dc*hc*b*qq;
%得到总的质量矩阵
=+++++;
syms Ib Ic C11d C11e G
%求刚度矩阵
syms Eb C11E C11D Ic;
%Nub求导
x=sym('x');
dNub=diff();
%dNub的积分
kkkk='*;
p=int(kkkk,x,0,1);
Kub=Eb*b*hb*p;
a=diff(Nw);
aa=diff(a);
pp=int(aa'*aa,x,0,1);
Kwb=Eb*Ib*pp;
dNuc=diff(Nuc);
p=int(dNuc'*dNuc,x,0,1);
Kuc=C11e*b*hc*p;
Kwc=C11d*Ic*pp;
p=int(Nb'*Nb);
Kbv=G*b*hv*p;
K=Kub+Kwb+Kuc+Kwc+Kbv;
K/M
[ 本帖最后由 eight 于 2008-4-18 10:39 编辑 ]
回复 楼主 的帖子
报什么错误 试试LR分解之后再进行运算 为什么要除啊.......................................................................................... 除了之后求特征值。报错说矩阵太大。 看起来你这个好像是一个动力学问题啊，能不能考虑一下自由度缩减？ 原帖由 南派的三叔 于 2008-4-18 15:58 发表 http://www.chinavib.com/forum/images/common/back.gif除了之后求特征值。报错说矩阵太大。
这和矩阵大小没关系，肯定是报错M矩阵奇异了，不会是矩阵太大。不要用除号，误差太大了。 原帖由 gh688 于 2008-4-18 22:50 发表 http://www.chinavib.com/forum/images/common/back.gif
这和矩阵大小没关系，肯定是报错M矩阵奇异了，不会是矩阵太大。不要用除号，误差太大了。
请问怎么才提高精度？ 原帖由 sigma665 于 2008-4-19 13:11 发表 http://www.chinavib.com/forum/images/common/back.gif
请问怎么才提高精度？
呵呵，我看他求用K/M，其实在matlab中就是按照K*inv(M)计算的，这种问题求逆的话不太好啊，傻瓜式办法直接eig(-M,K)要好一点，不用求逆，个人观点
回复 9楼 的帖子
我是问有什么方法可以代替/我也用到那个，做出来的结果不对，是不是这个问题 原帖由 sigma665 于 2008-4-19 14:02 发表 http://www.chinavib.com/forum/images/common/back.gif
我是问有什么方法可以代替/
我也用到那个，做出来的结果不对，是不是这个问题
如果只是求特征值这一步的话，那上面的命令可以替代，如果程序中别的计算方面要用/，还真的不知道有别的东西能代替它
回复 11楼 的帖子
多谢了
页:
[1]