几种求解多自由度系统的特征值与特征向量的数值方法
本帖最后由 博大广阔 于 2011-10-30 13:23 编辑%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 andr=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 =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);
=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 =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
=othogonal_matrixr(D);
R=R*Rn; %阵型R=R1*R2*R3
end
%输出固有频率和阵型
r=D;
X=R;
end
function =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 thefrequencies 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=*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 =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 =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;
=eig(inv(K1)*M1);
Xnew=Xold*V;
a=Xnew-Xold;
erro=norm(a(:,1),2);
Xold=Xnew;
end
R=rr;
XX=Xold;
end
！！！若有误请指正
页:
[1]