由于原问题较简单,可以直接通过定义来求解特征值及特征向量。
|A-xI|=0, 即化为简单的1元3次方程,求解的x=[-1 -2 -3]。
然后根据(A-xI)v=0,分别将以上三个值代入求解,记得3组特征向量。从而将a对角化。
这个方法没有牵涉到特殊的矩阵运算,只有简单的解方程。matlab以及c实现都非常简单。
这个方法比较适合小型问题的求解。以下是基于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
这个方法比较适合中大型问题的求解。但是需要预处理。即该方法只适用于对称矩阵的特征值求解。所以需要先将原矩阵化为对称矩阵,例如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