用matlab求解简支梁的故有频率的程序

jiaozi20040560

%Define all vaiables.
N1=8;
EJ=1;
pA=1;
L=1;
m=pA*L/16;
k=EJ/16*(L^3);
%Define the mass matrix and the rigid matrix when both N1=8.
format long
%when the number of the beam elements is N1:
L1=L/N1;
M1=pA*L1/420*[156 22*L1 54 -13*L1;
                                22*L1 4*L1^2 13*L1 -3*L1^2;
                                54 13*L1 156 -22*L1;
                                -13*L1 -3*L1^2 -22*L1 4*L1^2];
K1=(EJ/L1^3)*[12 6*L1 -12 6*L1;
                            6*L1 4*L1^2 -6*L1 2*L1^2;
                            -12 -6*L1 12 -6*L1;
                            6*L1 2*L1^2 -6*L1 4*L1^2];
%Define matrix S1(4*18)(num=18)
S=zeros(4,20);
for num=1:8
            for r=1:4
                S(r,2*(num-1)+r,num)=1;
            end
end
%transimition of the mass matrix and the rigid matrix:
Me1=S(:,:,1)'*M1*S(:,:,1);
Ke1=S(:,:,1)'*K1*S(:,:,1);
for k=2:8
    Me1=Me1+S(:,:,k)'*M1*S(:,:,k);
    Ke1=Ke1+S(:,:,k)'*K1*S(:,:,k);
end
%the mass matrix and rigid matrix after considering the additional mass and
%string
Me1(19,19)=Me1(19,19)+m;
Me1(20,20)=Me1(20,20)+m;
Ke1(5,5)=Ke1(5,5)+k;
Ke1(19,19)=Ke1(19,19)+k;
Ke1(5,19)=Ke1(5,19)-k;
Ke1(19,5)=Ke1(19,5)-k;
Ke1(13,13)=Ke1(13,13)+k;
Ke1(20,20)=Ke1(20,20)+k;
Ke1(13,20)=Ke1(13,20)-k;
Ke1(20,13)=Ke1(20,13)-k;
Me1(1,:)=[];
Me1(:,1)=[];
Me1(17,:)=[];
Me1(:,17)=[];
Ke1(1,:)=[];
Ke1(:,1)=[];
Ke1(17,:)=[];
Ke1(:,17)=[];
w=sqrt(eig(inv(Me1)*Ke1))
求解的结果不对
正确的结果是w1=1，w2=1，w3=pi^2,w4=(2*pi)^2...........

jiaozi20040560

请各位不吝赐教；振动结构未下图

[ 本帖最后由 jiaozi20040560 于 2009-5-16 10:18 编辑 ]

jiaozi20040560

经过反复实践，正确的程序应该为：
%Define all vaiables.
N1=8;
EJ=1;
pA=1;
L=1;
m=pA*L/16;
k=EJ/16*(L^3);
%Define the mass matrix and the rigid matrix when both N1=8.
format long
%when the number of the beam elements is N1:
L1=L/N1;
M1=pA*L1/420*[156 22*L1 54 -13*L1;
                                22*L1 4*L1^2 13*L1 -3*L1^2;
                                54 13*L1 156 -22*L1;
                                -13*L1 -3*L1^2 -22*L1 4*L1^2];
K1=(EJ/L1^3)*[12 6*L1 -12 6*L1;
                            6*L1 4*L1^2 -6*L1 2*L1^2;
                            -12 -6*L1 12 -6*L1;
                            6*L1 2*L1^2 -6*L1 4*L1^2];
%Define matrix S1(4*18)(num=18)
S=zeros(4,20);
for num=1:8
            for r=1:4
                S(r,2*(num-1)+r,num)=1;
            end
end
%transimition of the mass matrix and the rigid matrix:
Me1=S(:,:,1)'*M1*S(:,:,1);
Ke1=S(:,:,1)'*K1*S(:,:,1);
for j=2:8
    Me1=Me1+S(:,:,j)'*M1*S(:,:,j);
    Ke1=Ke1+S(:,:,j)'*K1*S(:,:,j);
end
%the mass matrix and rigid matrix after considering the additional mass and
%string
Me1(19,19)=Me1(19,19)+m;
Me1(20,20)=Me1(20,20)+m;
Ke1(5,5)=Ke1(5,5)+k;
Ke1(19,19)=Ke1(19,19)+k;
Ke1(5,19)=Ke1(5,19)-k;
Ke1(19,5)=Ke1(19,5)-k;
Ke1(13,13)=Ke1(13,13)+k;
Ke1(20,20)=Ke1(20,20)+k;
Ke1(13,20)=Ke1(13,20)-k;
Ke1(20,13)=Ke1(20,13)-k;
Me1(17,:)=[];
Me1(:,17)=[];
Ke1(17,:)=[];
Ke1(:,17)=[];
Me1(1,:)=[];
Me1(:,1)=[];
Ke1(1,:)=[];
Ke1(:,1)=[];

%[v,d]=eig(inv(Mlow)*Klow):
w=sqrt(eig(inv(Me1)*Ke1))
计算结果为：
w =

  1.0e+003 *

   3.21277450189085
   3.04282729351354
   2.64008614395003
   2.18124279527326
   1.76223447547116
   1.40774559762493
   1.11437149328907
   0.87286529922085
   0.70108487360661
   0.49867300912581
   0.36179846558380
   0.24902787295142
   0.15853693939645
   0.08894142804986
   0.03949183931220
   0.00987616297255
   0.00099934295905
   0.00099991857824

从结果可以看到：w1=0.99，w2=0.99正好是两个弹簧振子的故有频率，w3=9.876=pi^2,w4=(2pi)^2恰好是简支梁的第一介和第二节故有频率。