首页/文章/ 详情

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

2年前浏览2888

   

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

【免责声明】本文来自新浪夜帝的博客,博主:夜帝,版权归原作者所有,仅用于学习等,对文中观点判断均保持中立,若您认为文中来源标注与事实不符,若有涉及版权等请告知,将及时修订删除,谢谢大家的关注!





来源:CAE之家
科普求解技术代码&命令MATLAB
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2022-07-09
最近编辑:2年前
CAE之家
硕士 | CAE仿真负责人 个人著作《汽车NVH一本通》
获赞 1137粉丝 5971文章 921课程 20
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈