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