关于矩阵太大导致求解刚度和质量矩阵时出错的问题

南派的三叔

我编了一段简单的程序，求解一个质量矩阵和一个刚度矩阵。两个矩阵都能求解，但是用刚度矩阵除以质量矩阵就报错。有人热心指正小弟感激不尽。
syms x y db dv dc hb hv hc b L;
y=x/L;
%求Nub与其转置地积分q
Nub=[0,0,1-y,0,0,0,y,0];
T=[Nub]'*[Nub];
q=int(T,0,1);
%质量矩阵[Mub]
[Mub]=db*hb*b*q;
Nw=[1-3*y^2+2*y^3,(y-2*y^2+y^3)*L,0,0,3*y^2-2*y^3,(y^3-y^2)*L,0,0];
%对Nw求导
x=sym('x');
dNw=diff(Nw)
%求Nw与其转置的积分qq
TT=Nw'*Nw;
qq=int(TT,0,1);
%质量矩阵 [Mwb]
[Mwb]=db*hb*b*qq;
[Nb]=[0,0,0,1-y,0,0,0,y];
[Nuc]=-(hv+hb/2+hc/2)*dNw+[Nub]+hv*[Nb];
[Nuv]=-(hb+hc)/2*dNw+[Nub]+hv/2+[Nb];
%质量矩阵[Muv]&[Mwv]
qqq=int([Nuv]'*[Nuv],0,1);
[Muv]=dv*hv*b*qqq;
[Mwv]=dv*hv*b*qq;
%求Nuc的转置与其的积分qqqq
TTT=[Nuc]'*[Nuc];
qqqq=int(TTT,0,1);
%质量矩阵[Muc]&[Mwc]
[Muc]=dc*hc*b*qqqq;
[Mwc]=dc*hc*b*qq;
%得到总的质量矩阵
[M]=[Mub]+[Mwb]+[Muv]+[Mwv]+[Muc]+[Mwc];




syms Ib Ic C11d C11e G
%求刚度矩阵
syms Eb C11E C11D Ic;
%Nub求导
x=sym('x');
dNub=diff([Nub]);
%dNub的积分
kkkk=[dNub]'*[dNub];
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 编辑 ]

sigma665

报什么错误

sogooda

试试LR分解之后再进行运算

lq12131010

为什么要除啊  ..........................................................................................

南派的三叔

除了之后求特征值。报错说矩阵太大。

sogooda

看起来你这个好像是一个动力学问题啊，能不能考虑一下自由度缩减？

gh688

原帖由 南派的三叔 于 2008-4-18 15:58 发表

                               
登录/注册后可看大图

除了之后求特征值。报错说矩阵太大。

这和矩阵大小没关系，肯定是报错M矩阵奇异了，不会是矩阵太大。不要用除号，误差太大了。

sigma665

原帖由 gh688 于 2008-4-18 22:50 发表

                               
登录/注册后可看大图

这和矩阵大小没关系，肯定是报错M矩阵奇异了，不会是矩阵太大。不要用除号，误差太大了。

请问怎么才提高精度？

gh688

原帖由 sigma665 于 2008-4-19 13:11 发表

                               
登录/注册后可看大图

请问怎么才提高精度？

呵呵，我看他求用K/M，其实在matlab中就是按照K*inv(M)计算的，这种问题求逆的话不太好啊，傻瓜式办法直接eig(-M,K)要好一点，不用求逆，个人观点

sigma665

我是问有什么方法可以代替/
我也用到那个，做出来的结果不对，是不是这个问题

gh688

原帖由 sigma665 于 2008-4-19 14:02 发表

                               
登录/注册后可看大图

我是问有什么方法可以代替/
我也用到那个，做出来的结果不对，是不是这个问题

如果只是求特征值这一步的话，那上面的命令可以替代，如果程序中别的计算方面要用/，还真的不知道有别的东西能代替它

sigma665

多谢了