Source: Sina Blog by Ye Di, modified slightly
How to solve eigenvalues and eigenvectors in MATLAB? There are three methods to address this problem.1. Definition Method
Since the original problem is relatively simple, we can directly solve for the eigenvalues and eigenvectors using definitions.
|A–xI|=0, which simplifies to a simple cubic equation, solving for x gives x=[-1 -2 -3].
Then, based on (A–xI)v=0, substitute the above three values to solve for the eigenvectors. Thus, we can diagonalize A.
This method does not involve special matrix operations, just simple equation solving. Implementations in MATLAB and C are quite straightforward.
2. Power Method
This method is more suitable for solving small problems. Below is the solution to this problem based on the power method. It can directly obtain the eigenvalues and eigenvectors. There are no very complex matrix operations, and it can be implemented with simple MATLAB or C programs. For more details, refer to this link.
The program is original and is likely hard to find on other websites.
function [x, v] = findeigen(A)% Usage:
% compute the subsequent eigenvalue and eigenvector
% Input:
% A original matrix
% x0 initial eigenvalue
% v0 initial eigenvector
% Output:
% x final eigenvalue
% v final eigenvector
% 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 using 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 using 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) ;
% normalize
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
This method is more suitable for medium to large problems. However, it requires preprocessing. This method is only applicable to the eigenvalue problem of symmetric matrices. Therefore, the original matrix needs to be transformed into a symmetric matrix, for example, through Hermitian Transformation, and then solve the problem. Here is the Jacobi method:
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 Rutishauser’s modifications of the classical
% Jacobi rotation method with threshold 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 eigenvector matrix
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
MATLAB Calculation of the Natural Characteristics of Vibration Systems
Disclaimer: This WeChat reposted article is for non-commercial educational and research purposes and does not imply endorsement of its viewpoints or verification of its content authenticity. Copyright belongs to the original author; if there are copyright issues with the reposted article, please contact us immediately, and we will make changes or delete related articles to protect your rights!