eigfem_s.m

来自「结构力学中的有限元例子,包含了7个分类文件夹」· M 代码 · 共 50 行

function [ E_Vec, eig_Values] = eigFEM_s (K, M, N)

% eigFEM_s
% ---------------------------------------------------------------------
% N smallest eigenvalues and eigenvectors for stiffness & mass: K, M
% INPUT:
%       K           - stiffness matrix
%       M           - mass matrix
%       N           - number of smallest eig-values to compute
% OUTPUT:
%       E_Vect     - matrix of eigen vectors
%       eig_Values - eigen values
%
% 2004
% ---------------------------------------------------------------------
NEG = N;               % No of smallest eig-values to be computed
% eigen values and vectors calculation -------------------------------
opts.maxit = 40;   % maximum number of iterations
opts.tol = 1e+1;   % convergence
% opts.v0 = ones(size(M,1),1).*6.2;             % starting vector
[eig_vect, eig_val]=eigs(K,M,NEG,'SM',opts);  % NEG smallest eig-vales
eig_val=sqrt(real(eig_val));             % omega eigen-values

eig_V=diag(eig_val);                     % freq
eig_vect=real(eig_vect);                 % shapes

D = size(eig_vect,1);                    % # of dof's 

% ordering eigen vectors matrix --------------------------------------
NOM=0; DIN=realmax; S=0;

for i=1:NEG
   for k=1:NEG
      if abs(eig_V(k))>=abs(NOM)
         if abs(eig_V(k))<abs(DIN)
             NOM=eig_V(k); S=k;
         end;
      end;
   end;
   eig_Values(i)=NOM; DIN=NOM;    % ordered eigen-values:    eig_Values
   for p=1:D
       E_Vec(p,i)=eig_vect(p,S);  % ordered eigen-vectors:   E_Vec(p,mode#) 
   end;
   NOM=0;
end;
disp('Periods: (s)  ');
disp('-----------------------');
disp(2*pi./eig_Values);
disp('-----------------------');