求救！各位帮我看看我哪里做错了

libbyzifeng

这是我求桁架固有频率的程序 ，求不出结果，而且结构质量矩阵和刚度矩阵也不对，请各位帮我看看哪里出错了，谢谢！
E=2.1e11;
A=1e-4;
density=7.3e3;
node_number=5;
element_number=7;
nc=[0,0;1,0;2,0;0,1;1,1];
en=[1,2;2,3;1,4;2,4;2,5;3,5;4,5];
ed(1:node_number,1:2)=1;
constraint=[1,1;1,2;3,2];
for loopi=1:length(constraint);
       ed(constraint(loopi,1),constraint(loopi,2))=0;
end
dof=0;
for loopi=1:node_number
      for loopj=1:2
         if ed(loopi,loopj)~=0
           dof=dof+1;
           ed(loopi,loopj)=dof;
         end
      end
end
ek=E*A*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];
em=(density*A)/2*eye(4);
k(1:dof,1:dof)=0;
m=k;
theta(1:7)=0;
el(1:7)=0;
e2s(1:4)=0;
for loopi=1:element_number
     for zi=1:2
        e2s((zi-1)*2+1)=ed(en(loopi,zi),1);
        e2s((zi-1)*2+2)=ed(en(loopi,zi),1);
     end
el(loopi)=sqrt((nc(en(loopi,1),1)-nc(en(loopi,2),1))^2+(nc(en(loopi,1),2)-nc(en(loopi,2),2))^2);
theta(loopi)=asin((nc(en(loopi,1),2)-nc(en(loopi,2),2))/el(loopi));
lmd=[cos(theta(loopi)) sin(theta(loopi));-sin(theta(loopi)) cos(theta(loopi))];
t=[lmd zeros(2);zeros(2) lmd];
dk=t'*ek*t/el(loopi);
dm=t'*em*t*el(loopi);
    for jx=1:4
     for jy=1:4
       if(e2s(jx)*e2s(jy)~=0)
         k(e2s(jx),e2s(jy))=k(e2s(jx),e2s(jy))+dk(jx,jy);
         m(e2s(jx),e2s(jy))=m(e2s(jx),e2s(jy))+dm(jx,jy);
       end
     end
   end
end
[v,d]=eig(k,m);
frequency=sqrt(diag(d))/(2*pi);
[frequency,indexf]=sort(frequency);
d=d(:,indexf);

libbyzifeng

各位帮我看看吧~~~