几种求解多自由度系统的特征值与特征向量的数值方法

博大广阔

%use standard eigenvalue problem to quit eigenvalue and eigenvector
%the eigenvalue pronlem was stated is :
%     K*X=w^2*M*X      %目的是求W和X
%which can be    rewritten in the form of a standard forms:
1）当D动力矩阵不对称的时候，使用choleski decomposition 分解方法，在求解相应的特征值与特征向量
   D*X=r*X   where D=inv(K)*M   and  r=1/w^2 %直接求解要求D是对称的
%assuming the matrix K is symmetric and positive definite,we can us
%choleski decomposition ans K express is:
%  K=U.'*U where U is upper triangular matrix
% so easy obtain r*U.'*U*X=M*X :premultipling this equation by inv(U.')
%and define a new vector Y=U*X we can get:
%        D*Y=r*Y     D=inv(U.')*M*inv(U)  and can slve Y and r
% then we can only apply invere transformation and find the desired
% eigenvectors.X=inv(U)*Y
function [r,X]=standard_eigenvalue(K,M)
if nargin<1;error('wrong inpiut');end
if size(K)~=size(M);error('wrong matrix');end
if K~=K.';error('K unsymmetric and cant usse this metohd');end
U=chol(K);
D=inv(U.')*M*inv(U);
[V,rr]=eig(D);
X=inv(U)*V;
r=rr;
end

2）雅克比方法求，适用于动力矩阵是对称的（不对称的可以先分解）可求解，所有的系统所有的固有频率和阵型
%jacobi method is also an interative method whitch produces all the
%eigenvalues and eigenvalues and eigenvectors of D simultaneously.
%this method is based on a theorom in linear algebra stating that a real
%symmetric matrix D such has only real eigenvalues ans that exist a real
%orthogonal matrix R such that R.'*D*R is diagonal.
%and the diagonal elments is the eigenvaous and the columns is eigenvectors
function [r,X]=jacobi_engenvalue(K,M,D)
%输入参数的判断
   if nargin<1;error('wrong inpiut');end
if nargin==2; if size(K)~=size(M);error('wrong matrix');end
D=inv(K)*M;
end
if D~=D.';error('D unsymmetric and cant use this metohd');end
%迭代主体
maxvalue=10^4;R=diag(ones(1,size(D,2)));
  while abs(maxvalue)>=0.001
      [Rn,D,maxvalue]=othogonal_matrixr(D);
      R=R*Rn;      %阵型R=R1*R2*R3
  end
  %输出固有频率和阵型
r=D;
X=R;
end

function [Rn,D,maxvalue]=othogonal_matrixr(D)
%每次迭代，更新矩阵D（固有频率），计算新的阵型，最大值用于截止判断
st=2;max=D(1,2);R=diag(ones(1,size(D,2)));
for row=1:1:size(D,2)-1
  for colum=st:1:size(D,2)
       if max<=D(row,colum)        
         max=D(row,colum);
         rr=row;
         cc=colum;
     end
  end
  st=st+1;
end
c=rr;
s=cc;
q=0.5*atan((2*D(c,s))/(D(c,c)-D(s,s)));
R(c,c)=cos(q);R(c,s)=-sin(q);
R(s,c)=sin(q);R(s,s)=cos(q);
Rn=R;
D=Rn.'*D*Rn;
maxvalue=max;
end
3）矩阵迭代法%matrix iteration method
%this method assumes that the  frequencies are distrinv and well
%seprated that w1<w2<w3...
%base princile:according to the expansion theorem:
%X1=c1*u1+c2*u2+c3*u3...... =C*Uwhitch ui is eigenvector
%conbine with D*X=r*X
%we obtain: X(r+1)=D*Xr=[c1/w1^2*r,c2/w2^2r,c3/w3^2r....]*U
%since the natrual frequences are assumed to be w1<w2<w3...
%so we can neglect the second to end
%X(r+1)=(c1/w1^2*r )*u1
%so :w1^2=Xr/X(r+1);  终止条件Xr=X(r+1); X(r+1)就是对应的一阶阵型。because D*Xr=r*Xr
%为计算第R+1阶的特征向量和特征值，应用正交条件消去前R阶
%Ci=ui*M*X/ui*M*ui 当ui是第i个阵型的正则化后的阵型ui*M*ui=1；X是假设的阵型，Ci是第i个阵型对该假设阵型的贡献量
function [R,XX]=matrix_interation(K,M,number,accuracy)
  if nargin<3;warndlg('wrong inpiut');end
if nargin<4; accuracy=0.001;end
D=K*M;  
for n=1:1:number                    %number is the need solved digenvector
     Xold=ones(size(D,2),1);erro=10^4;
      while erro>=accuracy
   Xnew=D*Xold;
  r=Xnew(1);
   Xnew=Xnew/Xnew(1);
   erro=norm((Xnew-Xold),2);
   Xold=Xnew;
   end
Xs=Xold/sqrt(Xold.'*M*Xold);
D=D-r*Xs*Xs.'*M;                 %消去前面几项,更新D矩阵
%%%保存各个阵型和固有频率
R(n)=r;
XX(:,n)=Xs;
end
R=sqrt(1./R);
XX=XX;
end







4）子空间跌代法%子空间迭代法
%其实质是不断改善李籽法中的假设阵型质量，使得李兹法的假设阵型所张成的子空间不断接近前R阶主阵型
function [R,XX]=zmatrix_interation(K,M,number,accuracy)
if nargin<3;warndlg('wrong inpiut');end
if nargin<4; accuracy=0.001;end
D=inv(K)*M;  
%%%
那个
%%%

erro=10^4;
  while erro>=accuracy  
  X1=D*Xold
  a=max(X1);
  for k=1:1:size(X1,2)   
  X1(:,k)=X1(:,k)./a(k);
  end
K1=X1.'*K*X1;
M1=X1.'*M*X1;
[V,rr]=eig(inv(K1)*M1);
Xnew=Xold*V;
a=Xnew-Xold;
erro=norm(a(:,1),2);
Xold=Xnew;
   end
R=rr;
XX=Xold;
end

！！！若有误请指正