MATLAB入门学习-#7-幂法&幂法的加速&反幂法编程练习
1.幂法(power method)2.幂法的加速——原点平移法3.幂法的加速——埃特肯（Aitken）4.反幂法(inverse power method)
这几种求特征值和特征向量的方法，理论讲解可以去看数值分析这门课的课本或者下面这个链接：
【图文】幂法及反幂法_百度文库
另外还有一个讲的很好的幂法相关的文档：
科学网—图说幂法求特征值和特征向量 - 康建的博文
下面就给出我写的程序了，统一用的是先用一个函数检验是否具有使用某一方法的条件，若有条件再调用计算函数，因此每种方法都有两个函数。程序还是有问题的，比如对误差限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的博客 - 博客 QUESTION2:matlab中，if,while,for这些语句后面都要用end来结束 QUESTION3:
=就是说按照长度为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
