matlab中三种求解特征值和特征向量的方法
本帖最后由 weixin 于 2017-10-10 10:31 编辑matlab如何求解特征值和特征向量?这个问题可以有以下3种解法。
1. 定义法
由于原问题较简单，可以直接通过定义来求解特征值及特征向量。
|A-xI|=0， 即化为简单的1元3次方程，求解的x=[-1 -2 -3]。
然后根据(A-xI)v=0，分别将以上三个值代入求解，记得3组特征向量。从而将a对角化。
这个方法没有牵涉到特殊的矩阵运算，只有简单的解方程。matlab以及c实现都非常简单。
2. power method
这个方法比较适合小型问题的求解。以下是基于power method对该问题进行求解。可以直接求得特征值和特征向量。没有非常复杂的矩阵操作，可以用简单的matlab或者c程序实现。介绍可以参考http://www.miislita.com/information-retrieval-tutorial/matrix-tutorial-3-eigenvalues-eigenvectors.html
程序为原创，估计很难在其他网站找到。
function = findeigen(A) % Usage: % compute the subsequent eigenvalue and eigenvector % Input: % A orginal matrix % x0 initial eigen value % v0 initial eigen vector % Output: % x final eigen value % v final eigen vector
% Author: % % Date: %
% maximum iteration itermax = 100 ; % minimum error errmax = 1e-8 ;
N = size(A, 1) ; xnew = 0 ; vnew = ones(N, 1) ;
x = zeros(1, N) ; v = zeros(N, N) ;
% calculate eigenvalue use The Deflation Method B = A ; for num1 = 1 : N if num1 > 1 B = B - xnew * vnew * vnew' ; else end % call power method to obtain the eigenvalue = powermethod(B, itermax, errmax) ; x(num1) = xnew ;
end
% calculate eigenvalue use The Inverse iteration method % shift value u = 0.1 ; for num1 = 1 : N C = inv(A - (x(num1)-u)*eye(N)) ; % call power method to obtain the eigenvector = powermethod(C, itermax, errmax) ; v(:, num1) = vnew ; end
function = powermethod(A, itermax, errmax)
N = size(A, 1) ; xold = 0 ; vold = ones(N, 1) ;
for num2 = 1 : itermax vnew = A * vold ; % get eigenvalue = max(abs(vnew)) ; xnew = vnew(i) ; % normlize vnew = vnew/xnew ; % calculate the error errtemp = abs((xnew-xold)/xnew) ; if(errtemp < errmax) x = xnew ; v = vnew ; break ; end xold = xnew ; vold = vnew ; end
3. Jacobi's Method
这个方法比较适合中大型问题的求解。但是需要预处理。即该方法只适用于对称矩阵的特征值求解。所以需要先将原矩阵化为对称矩阵，例如Hermitian Transformation，然后再对该问题求解。附带jocobi方法：
function =jacobi(a_in,itermax) % =jacobi(a_in,itermax) % computes the eigenvalues d and % eigenvectors v of the real symmetric matrix a_in, % using rutishausers modfications of the classical % jacobi rotation method with treshold pivoting. % history(1:historyend) is a column vector of the length of % total sweeps used containing the sum of squares of % strict upper diagonal elements of a. a is not % touched but copied locally % the upper triangle is used only % itermax is the maximum number of total sweeps allowed % numrot is the number of rotations applied in total % check arguments siz=size(a_in); if siz(1) ~= siz(2) error('jacobi : matrix must be square ' ); end if norm(a_in-a_in',inf) ~= 0 error('jacobi ; matrix must be symmetric '); end if ~isreal(a_in) error(' jacobi : valid for real matrices only'); end n=siz(1); v=eye(n); a=a_in; history=zeros(itermax,1); d=diag(a); bw=d; zw=zeros(n,1); iter=0; numrot=0; while iter < itermax iter=iter+1; history(iter)=sqrt(sum(sum(triu(a,1).^2))); historyend=iter; tresh=history(iter)/(4*n); if tresh ==0 return; end for p=1:n for q=p+1:n gapq=10*abs(a(p,q)); termp=gapq+abs(d(p)); termq=gapq+abs(d(q)); if iter>4 & termp==abs(d(p)) & termq==abs(d(q)) % annihilate tiny elements a(p,q)=0; else if abs(a(p,q)) >= tresh %apply rotation h=d(q)-d(p); term=abs(h)+gapq; if term == abs(h) t=a(p,q)/h; else theta=0.5*h/a(p,q); t=1/(abs(theta)+sqrt(1+theta^2)); if theta < 0 t=-t; end end c=1/sqrt(1+t^2); s=t*c; tau=s/(1+c); h=t*a(p,q); zw(p)=zw(p)-h; %accumulate corrections to diagonal elements zw(q)=zw(q)+h; d(p)=d(p)-h; d(q)=d(q)+h; a(p,q)=0; %rotate, use information from the upper triangle of a only %for a pipelined cpu it may be better to work %on full rows and columns instead for j=1:p-1 g=a(j,p); h=a(j,q); a(j,p)=g-s*(h+g*tau); a(j,q)=h+s*(g-h*tau); end for j=p+1:q-1 g=a(p,j); h=a(j,q); a(p,j)=g-s*(h+g*tau); a(j,q)=h+s*(g-h*tau); end for j=q+1:n g=a(p,j); h=a(q,j); a(p,j)=g-s*(h+g*tau); a(q,j)=h+s*(g-h*tau); end % accumulate information in eigenvectormatrix for j=1:n g=v(j,p); h=v(j,q); v(j,p)=g-s*(h+g*tau); v(j,q)=h+s*(g-h*tau); end numrot=numrot+1; end end %if end % for q end % for p bw=bw+zw; d=bw; zw=zeros(n,1); end %while
来源：新浪夜帝的博客，博主：夜帝，略有修改
页:
[1]