声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2163|回复: 0

[编程技巧] matlab中三种求解特征值和特征向量的方法

[复制链接]
发表于 2017-10-10 10:13 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 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 [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

  来源:新浪夜帝的博客,博主:夜帝,略有修改

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-29 15:30 , Processed in 0.229454 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表