matlab中三种求解特征值和特征向量的方法

weixin

　　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 [x, v] = 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

           [xnew, vnew] = 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

           [xnew, vnew] = powermethod(C, itermax, errmax) ;

           v(:, num1) = vnew ;

           end

           function [x, v] = powermethod(A, itermax, errmax)

           N = size(A, 1) ;

           xold = 0 ;

           vold = ones(N, 1) ;

           for num2 = 1 : itermax

           vnew = A * vold ;

          % get eigenvalue

           [xnew, i] = 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 [v,d,history,historyend,numrot]=jacobi(a_in,itermax)

           % [v,d,history,historyend,numrot]=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

　　来源：新浪夜帝的博客，博主：夜帝，略有修改