MATLAB求解矩阵特征值的六种方法
关于这个特征值的求解一共六种方法
幂法
反幂法
QR方法
对称QR方法
jacobi方法
二分法
接下来就着重讲解这些算法的是如何使用的
幂法
算法如下，
输入:
矩阵A、非零矢量x0、maxit(2000）、tol（1.0e-7)
输出：
模的最大特征量a、模的最大特征量对应的特征向量x
function [a,x,n] = pmethod(A,x0,maxit,tol)
if nargin == 3
tol = 1.0e-6;
elseif nargin == 2
maxit = 1000;
tol = 1.0e-6;
end
a0 = 0;
y = A * x0;
a = max(abs(y));
x = y / a;
n = 1;
while (abs(a - a0) > tol) && (n < maxit)
a0 = a;
x0 = x;
y = A * x0;
a = max(abs(y));
x = y / a;
n = n + 1;
end
if (n > maxit)
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
else
disp(['幂法迭代',num2str(n),'步收敛']);
end
反幂法应用幂法于矩阵的反求矩阵的特征值
输入：矩阵A、非零矢量x0、maxit(2000）、tol（1.0e-7)
输出：模的最小特征量a、模的最小特征量对应的特征向量x
function [a,x,n] = vpmethod(A,x0,maxit,tol)
if nargin == 3
tol = 1.0e-6;
elseif nargin == 2
maxit = 1000;
tol = 1.0e-6;
end
[L,U,P] = lu(A);
a0 = 0;
y = utri(U,ltri(L,P*x0));
a = max(abs(y));
x = y / a;
n = 1;
while (abs(a - a0) > tol) && (n < maxit)
a0 = a;
x0 = x;
y = ltri(L,P*x0);
y = utri(U,y);
a = max(abs(y));
x = y / a;
n = n + 1;
end
a = 1 / a;
if (n >= maxit)
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
else
disp(['反幂法迭代',num2str(n),'步收敛']);
end
这里还用到上次的ltri和utri函数，使用时把他们当作函数来使用就可以
%ltri
function y = ltri(L,b)
n=size(b,1);
y=zeros(n,1);
for j = 1:n-1
y(j) = b(j)/L(j,j);
b(j+1:n) = b(j+1:n) - y(j) * L(j+1:n,j);
end
y(n) = b(n)/L(n,n);
%utri
function x = utri(U,y)
n=size(y,1);
x=zeros(n,1);
for j = n:-1:2
x(j) = y(j)/U(j,j);
y(1:j-1) = y(1:j-1) - x(j) * U(1:j-1,j);
end
x(1) = y(1)/U(1,1);
QR方法利用正交相似变换
输入:
矩阵A、maxit(2000）、tol（1.0e-7)
输出：
所有特征值、所有特征值所对应的特征向量
function [a,x] = qrmd(A,maxit,tol)
%a(i,1)为第i个特征值，x(:,i)为第i个特征值对应的特征向量
if nargin == 2
tol = 1.0e-6;
elseif nargin == 1
maxit = 1000;
tol = 1.0e-6;
end
A0 = A;
a0 =diag(A);
[Q,R] = qr(A);
A = R*Q;
a = diag(A);
n = 1;
while (max(abs(a-a0)) > tol) && (n <= maxit)
a0 = a;
[Q,R] = qr(A);
A = R*Q;
a = diag(A);
n = n + 1;
end
if (n > maxit)
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
else
disp(['QR方法迭代',num2str(n),'步收敛']);
n = size(a,1);
x = zeros(n);
for i = 1:n
x0 = ones(n,1);
[b,x0] = vpmethod(A0-a(i,1)*eye(n),x0,maxit,tol);
x(:,i) = x0;
a(i,:) = a(i,:) + b;
end
end
这里也是运用到了一些常用的计算函数，直接复制在上面就可以
%ltri
function y = ltri(L,b)
n=size(b,1);
y=zeros(n,1);
for j = 1:n-1
y(j) = b(j)/L(j,j);
b(j+1:n) = b(j+1:n) - y(j) * L(j+1:n,j);
end
y(n) = b(n)/L(n,n);
%utri
function x = utri(U,y)
n=size(y,1);
x=zeros(n,1);
for j = n:-1:2
x(j) = y(j)/U(j,j);
y(1:j-1) = y(1:j-1) - x(j) * U(1:j-1,j);
end
x(1) = y(1)/U(1,1);
%vpmethod
function [a,x,n] = vpmethod(A,x0,maxit,tol)
if nargin == 3
tol = 1.0e-6;
elseif nargin == 2
maxit = 1000;
tol = 1.0e-6;
end
[L,U,P] = lu(A);
a0 = 0;
y = utri(U,ltri(L,P*x0));
a = max(abs(y));
x = y / a;
n = 1;
while (abs(a - a0) > tol) && (n < maxit)
a0 = a;
x0 = x;
y = ltri(L,P*x0);
y = utri(U,y);
a = max(abs(y));
x = y / a;
n = n + 1;
end
a = 1 / a;
if (n >= maxit)
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
else
disp(['反幂法迭代',num2str(n),'步收敛']);
end
对称QR方法
输入：
矩阵A、maxit(2000）、tol（1.0e-7)
输出：
所有特征值、所有特征值所对应的特征向量
function [a,x] = wilkinsonqr(A,maxit,tol)
%a(i,1)为第i个特征值，x(:,i)为第i个特征值对应的特征向量,A为实对称矩阵
if nargin == 2
tol = 1.0e-6;
elseif nargin == 1
maxit = 1000;
tol = 1.0e-6;
end
n = size(A,1);
A0 = A;
A = hess(A);
a0 =diag(A);
d = (A(n-1,n-1) - A(n,n)) / 2;
u = A(n,n) + min(abs([d + sqrt(d^2 + A(n,n-1)^2),d - sqrt(d^2 + A(n,n-1)^2)]));
[Q,R] = qr(A - u * eye(n));
A = R*Q + u * eye(n);
a = diag(A);
m = 1;
while (max(abs(a-a0)) > tol) && (m <= maxit)
a0 = a;
d = (A(n-1,n-1) - A(n,n)) / 2;
u = A(n,n) + min(abs([d + sqrt(d^2 + A(n,n-1)^2),d - sqrt(d^2 + A(n,n-1)^2)]));
[Q,R] = qr(A - u * eye(n));
A = R*Q + u * eye(n);
a = diag(A);
m = m + 1;
end
if (m >= maxit)
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
else
disp(['隐式对称QR方法迭代',num2str(m),'步收敛']);
n = size(a,1);
x = zeros(n);
for i = 1:n
x0 = ones(n,1);
[b,x0] = vpmethod(A0-a(i,1)*eye(n),x0,maxit,tol);
x(:,i) = x0;
end
这里是要用到已经学过的ltir、utri、 vpmethod函数，直接复制上面就好
jacobi方法最古老的方法
输入：
矩阵A、maxit(2000）、tol（1.0e-7)
输出：
模的最大特征量a、模的最大特征量对应的特征向量x
function [a,V] = jacobi_eig(A,maxit,tol)
%a的元素为A的特征值，V的第j列为对应于特征值a(j,1)的特征向量
if nargin == 2
tol = 1.0e-6;
elseif nargin == 1
maxit = 1000;
tol = 1.0e-6;
end
A0 = A;
n = size(A,1);
V = eye(n);
b = A0 - diag(diag(A0));
w = sqrt(sum(sum(b .* b)));
[r,c] = find(abs(b) > w);
k = 1;
while (w > tol) && (k < maxit)
if (~isempty(r))
for i = 1:size(r,1)
if (r(i,1) > c(i,1))
if (abs(A0(r(i,1),c(i,1))) < tol)
u = 0;
else
u = acot((A0(r(i,1),r(i,1)) - A0(c(i,1),c(i,1)))...
/ (2 * A0(r(i,1),c(i,1)))) / 2;
end
V1 = eye(n);
V1(r(i,1),r(i,1)) = cos(u);
V1(c(i,1),c(i,1)) = cos(u);
V1(r(i,1),c(i,1)) = -sin(u);
V1(c(i,1),r(i,1)) = sin(u);
A0 = V1' * A0 * V1;
V = V * V1;
k = k + 1;
end
end
end
b = A0 - diag(diag(A0));
w = w / n;
[r,c] = find(abs(b) > w);
end
if (k <= maxit)
disp(['Jacobi方法迭代',num2str(k),'步收敛！']);
else
disp(['迭代步数大于',num2str(maxit),'收敛太慢了！']);
end
a = diag(A0);
二分法由两部分模块组成
输入：
矩阵A、maxit(2000）、tol（1.0e-7)
输出：
所有特征值、所有特征值所对应的特征向量
第一部分
function n = switch_num(T,u)
x = diag(T);
y = [0;diag(T,1)];
m = size(T,1);
q = x(1,1) - u;
n = 0;
for k = 1:m
if (q < 0)
n = n + 1;
end
if (k < m)
if (q == 0)
q = abs(y(k+1,1)) * (1.0e-3);
end
q = x(k+1,1) - u - y(k+1,1)^2 / q;
end
end
第二部分
输入：
A为非零矩阵,b0,b1为所求特征值所在的区间端点、tol（1.0e-7)
输出：所有特征值、所有特征值所对应的特征向量
function [a,x] = dichotomy_eig(A,b0,b1,tol)
%A为非零矩阵,b0,b1为所求特征值所在的区间端点，返回值a为区间内的所有特征值所构成的向量
A = hess(A);
mx = max(sum(abs(A)));
if nargin == 3
tol = 1.0e-6;
elseif nargin == 2
b1 = mx;
tol = 1.0e-6;
elseif nargin ==1
b0 = -mx;
b1 = mx;
tol = 1.0e-6;
end
n = switch_num(A,b1) - switch_num(A,b0);
n1 = switch_num(A,b0) - switch_num(A,-mx);
a = ones(n,1)*b0;
for i = 1:n
b0 = a(i,1);
b1 = mx;
while (abs(b1 - b0) > tol)
a(i,1) = (b0 + b1) / 2;
s = switch_num(A,a(i,1));
if (s >= i+n1)
b1 = a(i,1);
else
b0 = a(i,1);
end
end
end
m = size(A,1);
x = zeros(m,n);
maxit = 1000;
for i = 1:n
x0 = ones(m,1);
[b,x0] = vpmethod(A-a(i,1)*eye(m),x0,maxit,tol);
x(:,i) = x0;
end
这里还是会用到已经学过的utri 、ltri 、vpmethod函数，直接复制在上面就好
参考文献
MATLAB从入门到实践–谢汉龙著
