MATLAB入门学习-#7-幂法&幂法的加速&反幂法编程练习
这几种求特征值和特征向量的方法,理论讲解可以去看数值分析这门课的课本或者下面这个链接: 【图文】幂法及反幂法_百度文库
另外还有一个讲的很好的幂法相关的文档: 科学网—图说幂法求特征值和特征向量 - 康建的博文
下面就给出我写的程序了,统一用的是先用一个函数检验是否具有使用某一方法的条件,若有条件再调用计算函数,因此每种方法都有两个函数。程序还是有问题的,比如对误差限epsilon,设置好以后会多算一次,就是说程序的顺序还是有问题,最近时间很紧张,我就不再去修改了,就把我第一次写的放上来,把中间一些注意的地方记录到就好了。
1.幂法(power method)
function [lambda,beta]=power_method(A,times,x0,epsilon)
%这个是一个检验A能不能用幂法的函数
%[lambda,beta]=power_method(A,times,x0,epsilon)
%lambda是最后得出的特征值
%beta是最后得到的特征值对应的特征向量
%A是待求矩阵
%times是最大迭代次数
%x0是初始向量
%epsilon是最大误差
format long;
b=eig(A); %求A的特征值所组成的向量
c=size(b);
d=c(1)*c(2); % 矩阵元素数量,即特征值的总个数
e=length(unique(b)); % 有几个代表值
if(d==e)
fprintf('特征值互异,是对角化矩阵\n'); %QUESTION1:fprintf
[lambda,beta]=power_method_cal(A,times,x0,epsilon) %再去调用计算函数,代码在下面
else
error('特征值有重复,不是对角化矩阵\n');
end
function [alpha,x]=power_method_cal(A,times,x0,epsilon)
%[alpha,x]=power_method_cal(A,times,x0,epsilon)
format long;
k=1;
u=0;
[m,n]=size(x0);
y=zeros(m,n);
y=x0;
while k <= times
x=A*y;
alpha=max(x); %用max()求一个矩阵中最大的元素
y=x./alpha;
fprintf('第%d次迭代\n',k); %QUESTON3:%d %f %s
fprintf('alpha=%.8f\n',alpha)
disp(x);
disp(y);
if abs(alpha-u)<epsilon
break
else
k=k+1;
u=alpha;
end %QUESTION2:end
end
QUESTION1:matlab disp、sprintf、fprintft函数 - hengyaha的博客 - CSDN博客
QUESTION2:matlab中,if,while,for这些语句后面都要用end来结束
QUESTION3:
%3d就是说按照长度为3的整型输出,比如10,输出就是“10”,“”代表空格。
%6.2f就是小数点后保留2位,输出总长度为6,比如3.14159,输出后就是“_ _ _3.14”(前面三个空格)
%c就是输出字符串
%s就是输出字符串,和%c是一样的
2.幂法的加速——原点平移法
function [lambda,beta]=origin_displacement_method(A,times,x0,epsilon,lambda0)
%[lambda,beta]=origin_displacement_method(A,times,x0,epsilon,lambda0)
%lambda是最后得出的特征值
%beta是最后得到的特征值对应的特征向量
%A是待求矩阵
%times是最大迭代次数
%x0是初始向量
%epsilon是最大误差
%lambda0就是lambda0
b=eig(A);
c=size(b);
d=c(1)*c(2); % 矩阵元素数量
e=length(unique(b)); % 有几个代表值
if(d==e)
fprintf('特征值互异,是对角化矩阵\n');
[lambda,beta]=origin_displacement_method_cal(A,times,x0,epsilon,lambda0)
else
error('特征值有重复,不是对角化矩阵\n');
end
function [alpha,x]=origin_displacement_method_cal(A,times,x0,epsilon,lambda0)
%[lambda,beta]=origin_displacement_method_cal(A,times,x0,epsilon,lambda0)
k=1;
u=0;
[m,n]=size(x0);
y=zeros(m,n);
y=x0;
I=eye(m);
A=A-lambda0*I;
while k <= times
x=A*y;
alpha=max(x);
y=x./alpha;
fprintf('第%d次迭代\n',k);
fprintf('alpha=%.8f\n',alpha)
fprintf('特征值为%.8f\n',alpha+lambda0)
disp(x);
disp(y);
if abs(alpha-u)<epsilon
break
else
k=k+1;
u=alpha;
end
end
3.幂法的加速——埃特肯(Aitken)
function [lambda,beta]=power_method_aitken(A,x0,times,epsilon)
%[lambda,beta]=power_method_aitken(A,x0,times,epsilon)
%lambda是最后得出的特征值
%beta是最后得到的特征值对应的特征向量
%A是待求矩阵
%x0是初始向量
%times是最大迭代次数
%epsilon是最大误差
format long;
b=eig(A);
c=size(b);
d=c(1)*c(2); % 矩阵元素数量
e=length(unique(b)); % 有几个代表值
if(d==e)
fprintf('特征值互异,是对角化矩阵\n');
[lambda,beta]=power_method_aitken_cal(A,x0,times,epsilon);
else
error('特征值有重复,不是对角化矩阵\n');
end
function [lambda,x]=power_method_aitken_cal(A,x0,times,epsilon)
%[alpha,x]=power_method_aitken_cal(A,x0,times,epsilon)
format long;
k=1;
u=1;
[m,n]=size(x0);
y=zeros(m,n);
x=x0;
alpha0=0;
alpha1=0;
while k <= times
alpha=max(x);
y=x./alpha;
x=A*y;
alpha2=max(x);
lambda=alpha0-((alpha1-alpha0)^2)/(alpha2-2*alpha1+alpha0);
fprintf('第%d次迭代\n',k);
fprintf('alpha0=%.8f\n',alpha0);
fprintf('alpha1=%.8f\n',alpha1);
fprintf('alpha2=%.8f\n',alpha2);
fprintf('lambda=%.8f\n',lambda);
disp(x);
disp(y);
if abs(lambda-u)<epsilon
break
else
k=k+1;
u=lambda;
alpha0=alpha1;
alpha1=alpha2;
end
end
4.反幂法(inverse power method)
function [lambda,beta]=inverse_power_method(A,x0,lambda_star,times,epsilon)
%[lambda,beta]=inverse_power_method(A,x0,lambda_star,times,epsilon)
%lambda是最后得出的特征值
%beta是最后得到的特征值对应的特征向量
%A是待求矩阵
%x0是初始向量
%lambda_star是特征值lambda的近似值
%times是最大迭代次数
%epsilon是最大误差
format long;
b=eig(A);
c=size(b);
d=c(1)*c(2); % 矩阵元素数量
e=length(unique(b)); % 有几个代表值
if(d==e)
fprintf('特征值互异,是对角化矩阵\n');
[lambda,beta]=inverse_power_method_cal(A,x0,lambda_star,times,epsilon);
else
error('特征值有重复,不是对角化矩阵\n');
end
function [lambda,x]=inverse_power_method_cal(A,x0,lambda_star,times,epsilon)
%[lambda,beta]=inverse_power_method_cal(A,x0,lambda_star,times,epsilon)
format long;
k=1;
u=1;
[m,n]=size(x0);
y=zeros(m,n);
x=x0;
I=eye(m);
H=A-lambda_star.*I;
[L,U]=lu(H); %LU分解函数(matlab自带)
while k <= times
beta_b=max(x);
y=x./beta_b;
z=L\y;
x=U\z;
beta_b=max(x);
fprintf('第%d次迭代\n',k);
lambda=lambda_star+1/beta_b;
fprintf('lambda=%.8f\n',lambda);
disp(x);
disp(z);
if abs(1/beta_b-1/u)<epsilon
break
else
k=k+1;
u=beta_b;
end
end