yacobi算法求特征值matlab,求matlab中求特征值特征向量eig的源程序或解此题

匿名用户

1级

2008-12-14 回答

我觉得这个问题可以有以下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

[xnew, vnew] = powermethod(B, itermax, errmax) ; % call power method to obtain the eigenvalue

x(num1) = xnew ;

end

% calculate eigenvalue use The Inverse iteration method

u = 0.1 ; % shift value

for num1 = 1 : N

C = inv(A - (x(num1)-u)*eye(N)) ;

[xnew, vnew] = powermethod(C, itermax, errmax) ; % call power method to obtain the eigenvector

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 ;

[xnew, i] = max(abs(vnew)) ; % get eigenvalue

xnew = vnew(i) ;

vnew = vnew/xnew ; % normlize

errtemp = abs((xnew-xold)/xnew) ; % calculate the error

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