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 == 3tol = 1.0e-6;
elseif nargin == 2maxit = 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),'收敛太慢了！']);
elsedisp(['幂法迭代',num2str(n),'步收敛']);
end

反幂法应用幂法于矩阵的反求矩阵的特征值
输入：矩阵A、非零矢量x0、maxit(2000）、tol（1.0e-7)
输出：模的最小特征量a、模的最小特征量对应的特征向量x

function [a,x,n] = vpmethod(A,x0,maxit,tol)
if nargin == 3tol = 1.0e-6;
elseif nargin == 2maxit = 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),'收敛太慢了！']);
elsedisp(['反幂法迭代',num2str(n),'步收敛']);
end

这里还用到上次的ltri和utri函数，使用时把他们当作函数来使用就可以

%ltri
function y = ltri(L,b)
n=size(b,1);
y=zeros(n,1);
for j = 1:n-1y(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:2x(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 == 2tol = 1.0e-6;
elseif nargin == 1maxit = 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),'收敛太慢了！']);
elsedisp(['QR方法迭代',num2str(n),'步收敛']);n = size(a,1);x = zeros(n);for i = 1:nx0 = 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-1y(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:2x(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 == 3tol = 1.0e-6;
elseif nargin == 2maxit = 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),'收敛太慢了！']);
elsedisp(['反幂法迭代',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 == 2tol = 1.0e-6;
elseif nargin == 1maxit = 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),'收敛太慢了！']);
elsedisp(['隐式对称QR方法迭代',num2str(m),'步收敛']);n = size(a,1);x = zeros(n);for i = 1:nx0 = 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 == 2tol = 1.0e-6;
elseif nargin == 1maxit = 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;elseu = acot((A0(r(i,1),r(i,1)) - A0(c(i,1),c(i,1))).../ (2 * A0(r(i,1),c(i,1)))) / 2;endV1 = 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;endendendb = A0 - diag(diag(A0));w = w / n;[r,c] = find(abs(b) > w);
end
if (k <= maxit)disp(['Jacobi方法迭代',num2str(k),'步收敛！']);
elsedisp(['迭代步数大于',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:mif (q < 0)n = n + 1;endif (k < m)if (q == 0)q = abs(y(k+1,1)) *  (1.0e-3); endq = 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 == 3tol = 1.0e-6;
elseif nargin == 2b1 = mx;tol = 1.0e-6;
elseif nargin ==1b0 = -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:nb0 = 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);elseb0 = a(i,1);endend
end
m = size(A,1);
x = zeros(m,n);
maxit = 1000;
for i = 1:nx0 = ones(m,1);[b,x0] = vpmethod(A-a(i,1)*eye(m),x0,maxit,tol);x(:,i) = x0;
end

这里还是会用到已经学过的utri 、ltri 、vpmethod函数，直接复制在上面就好

参考文献
MATLAB从入门到实践–谢汉龙著

MATLAB求解矩阵特征值的六种方法相关推荐

使用MTL库求解矩阵特征值和特征向量
关于矩阵的特征值和特征向量求解,大部分的数学运算库都进行了提供,下面是使用MTL库的接口进行封装. #include <mtl/matrix.h> #include <mtl/mtl ...
MATLAB求解线性方程组的八种方法
MATLAB求解线性方程组的八种方法求解线性方程分为两种方法–直接法和迭代法常见的方法一共有8种直接法 Gauss消去法 Cholesky分解法迭代法 Jacobi迭代法 Gauss-Seid ...
2021-01-07 matlab数值分析矩阵特征值与特征向量的计算改进乘幂法反幂法
matlab数值分析矩阵特征值与特征向量的计算 1改进乘幂法 function [t,y]=eigIPower(A,v0,ep) [tv,ti]=max(abs(v0)); lam0=v0(ti) ...
【NA】矩阵特征值的雅可比方法
前面在求解非线性方程时介绍过Jacobi迭代法,是采用同步更新策略的迭代法.这里的雅可比方法,人名是同一个,但区别很大.其它以Jacobi命名的数学概念还有Jacobi矩阵,其元素是 nnn 元函数偏 ...
关于Matlab中矩阵元素的表示方法
首先一点要说明的是,在matlab中,矩阵中的元素序号是按照"先行后列"的顺序排列的. 设如下随机矩阵: A=rand(4,6) A = 0.6551 0.9597 0 ...
关于matlab中矩阵取值的方法
在matlab中,取出矩阵中某一个值的方法如下: 1.对于二维数组: a(i, j) % 表示取出二维数组a的第 i 行,第 j 列的数据 a(:, j) % 表示取出二维数组a的第 j 列的所有数据 ...
matlab中矩阵取值的方法
在matlab中,取出矩阵中某一个值的方法如下: 1.对于二维数组: a(i, j) % 表示取出二维数组a的第 i 行,第 j 列的数据 a(:, j) % 表示取出二维数组a的第 j 列的所有数据 ...
matlab求微分方程的系数,如何利用matlab求解矩阵系数的二阶微分方程
M=[2,0;0 1 ]; %质量矩阵 K=[6 -2;-2 4]; %刚度矩阵 a=0;b=0; C=a*K+b*M; ...
matlab求矩阵特征值和特征向量、行列式
如果A为方阵,满足AX=λX的λ称为A的特征值,X称为A的特征向量. 计算A的特征值用eig(A). 例:A=[1 2 3;4 5 6;7 8 9]; Z=eig(A) Z = 16.1168 ...
matlab二阶节的系数,如何利用matlab求解矩阵系数的二阶微分方程
M=[2,0;0 1 ]; %质量矩阵 K=[6 -2;-2 4]; %刚度矩阵 a=0;b=0; C=a*K+b*M; ...

MATLAB求解矩阵特征值的六种方法

MATLAB求解矩阵特征值的六种方法

MATLAB求解矩阵特征值的六种方法相关推荐

最新文章

热门文章