用Matlab求特征值（自振频率）时遇到的问题

zishendemi

大家好！

最近我在用matlab编有限元结构动力分析方面的程序，但在求特征值时遇到些问题。

我要求一个结构前10阶自振频率，这个问题归根结底其实就是特征值问题。先介绍一下大体背景：

[K-(w^2)M]V=0

V就是特征向量，而w^2就是特征值，w在工程上的物理意义是自振频率。

K和M都是130×130的矩阵，我已经建好。我用
W=sqrt(eigs(K,M,10,'sa'))
求前10个最小的特征值，但返回的10个特征值都是0。
当我求前30个最小的特征值时，前19个都是0，从20个开始就不是0了。
当我用W=sqrt(eig(K,M))
求全部特征值时，返回：第一个特征值是复数，其余的全是非零实数。
这就非常奇怪了：1.求全部特征值时返回的特征值没有0，而求前n个最小的特征值时却全返回0，为什么？；2.求全部特征值时为什么返回了一个复数特征值？自振频率为复数在工程上没有意义。

ps:已经设了format long，不会是因为小数点位数不够约等于0。.

还望高手能够指点迷津！非常感谢！

huangjinysf

可以把矩阵贴上来，大家看看啊！

ChaChing

回复 1 # zishendemi 的帖子

好奇想试试, 同上!:@)
mat档若不能上传,可先更名可上传格式再註明!

zishendemi

谢谢楼上两位！代码如下：
clear;
Rou=17000;
E=2.1*10^11;
I=0.548;
N=2*1000*[146.4
146.4
284.4
410.9
528.2
653.3
772.1
884.8
995.8
1105.7
1203.3
1293.8
1384.2
1462.4
1534.1
1594.1
1643.2
1621.1
1593.8
1562.9
1527.8
1487.9
1436.9
1376.1
1310.9
1240.0
1165.1
1090.8
1013.8
907.2
799.7
688.7
580.0
688.7
799.7
907.2
1013.8
1090.8
1165.1
1240.0
1310.9
1376.1
1436.9
1487.9
1527.8
1562.9
1593.8
1621.1
1643.2
1594.1
1534.1
1462.4
1384.2
1293.8
1203.3
1105.7
995.8
884.8
772.1
653.3
528.2
410.9
284.4
146.4
146.4 ];
L=[1.5
1.5
3
3
9
9
9
9
9
9
9
9
9
9
9
9
47
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
15
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
47
9
9
9
9
9
9
9
9
9
9
9
9
3
3
1.5
1.5];
cable=[1.30E+07
1.17E+07
1.19E+07
1.05E+07
1.66E+07
1.79E+07
1.95E+07
2.14E+07
2.41E+07
2.59E+07
2.90E+07
3.09E+07
3.45E+07
3.49E+07
3.97E+07
7.74E+07
5.00E+07
3.30E+07
2.55E+07
2.07E+07
1.67E+07
1.75E+07
1.94E+07
1.68E+07
1.52E+07
1.09E+07
8.70E+06
7.74E+06
1.05E+07
9.11E+06
6.80E+06
5.53E+06
5.53E+06
6.80E+06
9.11E+06
1.05E+07
7.74E+06
8.70E+06
1.09E+07
1.52E+07
1.68E+07
1.94E+07
1.75E+07
1.67E+07
2.07E+07
2.55E+07
3.30E+07
5.00E+07
7.74E+07
3.97E+07
3.49E+07
3.45E+07
3.09E+07
2.90E+07
2.59E+07
2.41E+07
2.14E+07
1.95E+07
1.79E+07
1.66E+07
1.05E+07
1.19E+07
1.17E+07
1.30E+07
];
M=zeros(132);
K=zeros(132);
for j=1:65   %i=4:2:124
k=BeamElementStiffness(E,I,L(j),N(j));
K=BeamAssemble(K,k,j,j+1);
m=BeamElementMass(Rou,L(j));
M=BeamMassAssemble(M,m,j,j+1);
end;
M(3,:)=[];
M(:,3)=[];
M(129,:)=[];
M(:,129)=[];
K(3,:)=[];
K(:,3)=[];
K(129,:)=[];
K(:,129)=[];
cable=zeros(64); %暂时不把cable加到总刚K中，即一个简支梁
a=2;
for i=4:2:126
    K(i,i)=K(i,i)+cable(a);
    a=a+1;
end;
K(1,1)=K(1,1)+cable(1);
K(129,129)=K(129,129)+cable(64);
W=sqrt(eigs(K,M,10,'sa'));
W

zishendemi

补充一下：
BeamElementStiffness
BeamAssemble
BeamElementMass
BeamMassAssemble
这四个是我调用的函数，我把它们传上来。

zishendemi

自己顶！

ChaChing

回复 4 # zishendemi 的帖子

1.本想LZ仅上传K,M矩阵即可, LZ竟传齐了
2.试跑下eigs(K,M,10,'sa'), 出现一堆讯息, 尤其iteration 301次, 直觉没收敛
3.help eigs, 看了下说明(因不常用, 记不得)
4.试下[V2,D2,flag] =eigs(K,M,10,'sa'); 发现flag不为0, 没收敛!
(If flag is 0 then all the eigenvalues converged; otherwise not all converged.)
5.改用[V2,D2,flag] =eigs(K,M,10,'sm'); 发现flag为1, 收敛!
6.所以将'sa'改为'sm'即可!
7.至於为何, 有点懒又没多少时间查资料了! 就等LZ或大牛们挖掘说明了