幂法-反幂法求解矩阵最大最小特征值及其对应的特征向量.doc

数值计算解矩阵的按模最大最小特征值及对应的特征向量
一.幂法
1. 幂法简介：
当矩阵A满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。矩阵A需要满足的条件为：
(1)
(2) 存在n个线性无关的特征向量，设为
1.1计算过程：
不全为0，则有

可见，当越小时，收敛越快；且当k充分大时，有，对应的特征向量即是。

2 算法实现

3 matlab程序代码
function [t,y]=lpowerA,x0,eps,N) % t 为所求特征值，y是对应特征向量
k=1;

z=0; % z 相当于
y=x0./max(abs(x0)); % 规范化初始向量
x=A*y; % 迭代格式
b=max(x); % b 相当于
if abs(z-b)eps && k);
x=A*y;
b=max(x);
end
[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值
t=x(index); % 是原值，而非其绝对值。
end

4 举例验证
选取一个矩阵A，代入程序，得到结果，并与eig(A)的得到结果比较，再计算 A*y-t*y，验证y是否是对应的特征向量。结果如下：

结果正确，表明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与eig(A)的得到结果比较，再计算 A*y-t*y，验证y是否是对应的特征向量。设置初始向量为x0=ones(15,1)，结果显示如下

可见，结果正确。得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。
二.反幂法
1.反幂法简介及其理论
在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。其基本理论如下，与幂法基本相同：
，可知，A和A-1的特征值互为倒数，求A按模最小特征值即求A-1的按模最大特征值,取倒数即为A的按模最小特征值所以算法基本相同，区别就是在计算
2. 算法实现

3 matlab程序代码

function [s,y]=invpower(A,x0,eps,n) % s 为按模最小特征值，y是对应特征向量
k=1; r=0; % r相当于
y=x0./max(abs(x0)); % 规范化初始向量
[L,U]=lu(A);
z=L\y;
x=U\z;
u=max(x);
s=1/u; % 按模最小为A-1按模最大的倒数.
if abs(u-r)eps && k
x=U\z;
u=max(x);
end
[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值
s=1/x(index); % 是原值，而非其绝对值。
end

4 举例验证
同幂法一样，选取一个矩阵A，代入程序，得到结果，并与eig(A)的得到结果比较，再计算 A*y-t*y，验证y是否是对应的特征向量。

可见结果正确，然后利用此程序计算15阶Hilb矩阵，eig(A)的得到结果比较，再计算 A*y-s*y，验证y是否是对应的特征向量。设置初始向量为x0=ones(15,1)，结果显示如下

可见，结果真确。得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。
三. 计算条件数
矩阵A的条件数等于A的范数与A的逆的范数的乘积，即cond(A)=‖A‖·‖A^(-1)‖，对应矩阵的3种范数，可以定义3种条件数。 函数 cond(A,1)、cond(A)或cond(A inf)是判断矩阵病态与否的一种度量，条件数越大表明矩阵的病态程度越大.
,
而如果A为对称矩阵，如Hilb矩阵，的最大最小特征值，分别为A的最大最小特征值的平方。所以cond(A) 为A的最大最小特征值得比值。对于本例中的15阶Hilb矩阵来说，利用上面计算结果得其条件数(选择第二种条件数）为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017 在同一数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明Hilb矩阵是病态矩阵。
四. Aitken商加速法
1. 简介与原理

同幂法和反幂法计算最大和最小特征值类似，如果计算最大特征值，则迭代格式为；计算最小特征值时，迭代格式为。
2. 算法实现
计算按模最大特征值算法如下：

类似幂法和反幂法可以写出按模最小特征值算法，此处不再赘述。
3. matlab 程序代码
function [r,y]=aitken(A,x0,eps,n) % r按模最大特征值,y为对应特征向量
k=1;
a0=0; % a 相当于
a1=1; % a1 相当于
r0=1; % 相当于2中的
y=x0./max(abs(x0)); % 规范化初始向量
x=A*y;
a2=max(abs(x)); % a2相当于
r=a0-(a1-a0)^2/(a2-2*a1+a0); % 相当于
if (a2-2*a1+a0)==0 % 若上式中分母为0，则迭代失败，返回
disp "初始向量迭代失败"
return;
end
if abs(r-r0)要求，如满足，则返回结果
return
end
while abs(r-r0)>eps && km,index]=max(abs(eig(A))); % 以下代码保证取出来的按模最大特征值
aa=eig(A); % 是原值，而非其绝对值。
if aa(index)>0 ||aa(index)==0
r=r;
else
r=-r;
end
end
end
类似可得按模最小特征值和特征向量的代码如下：与上面类似，所不同的只是迭代格式不同.
function [r,y]=invaitken(A,x0,eps,n)
k=1;
a0=0;
a1=1;
r0=1;
y=x0./max(abs(x0));
[L,U]=lu(A); % 迭代格式的不同
z=L\y;
x=U\z;
a2=max(abs(x));
r=a0-(a1-a0)^2/(a2-2*a1+a0);
if (a2-2*a1+a0)==0
disp "初始向量迭代失败"
return;
end
if abs(r-r0)eps && k
b=c;
r0=r;
y=x./max(abs(x));
z=L\y;
x=U\z;
a2=max(abs(x));
if (a2-2*a1+a0)==0
return;
end
r=a0-(a1-a0)^2/(a2-2*a1+a0);
end
[m,index]=min(abs(eig(A))); % 以下代码保证取出来的按模最大特征值 aa=eig(A); % 是原值，而非其绝对值。
if aa(index)>0 ||aa(index)==0
r=1/r;
else
r=-1/r;
end
end

4. 计算Hilb矩阵特征值
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones(15,1)，结果如下，并将结果与幂法和反幂法得到结果比较

这与幂法得到的特征值和特征向量一致，表明算法和代码正确；同理，最小特征值结果如下：

这与反幂法得到的结果一致，表明结果正确。
五，对称矩阵的Rayleigh商加速法
1. 简介与原理

原理如下：

2. 算法实现

3. Matlab程序代码
function [r,y]=rayleigh(A,x0,eps,n) % r 是特征值，y是特征向量
k=1; r0=0;
y=x0./max(abs(x0));
x=A*y; % 迭代格式计算新的x
r=dot(y,x)/dot(y,y); % Reyleigh商
if abs(r-r0)eps && k r=dot(y,x)/dot(y,y);
end
end
类似得计算按模最小特征值的Rayleigh商加速法，如下：
function [r,y]=invrayleigh(A,x0,eps,n)
k=1;r0=0;
y=x0./max(abs(x0));
[L,U]=lu(A); % 迭代格式不同
z=L\y;
x=U\z;
r=dot(y,x)/dot(y,y);
if abs(r-r0)eps && k r0=r;
y=x./max(abs(x));
z=L\y;
x=U\z;
r=dot(y,x)/dot(y,y);
end
r=1/r;
end
4. 计算Hilb矩阵特征值
此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为ones(15,1)，结果如下，并将结果与幂法和反幂法得到结果比较

这与幂法得到结果一致，表明算法和代码正确。
同理，最小特征值如下：

与反幂法得到结果一致，表明代码和算法正确。