幂法_反幂法求解矩阵最大最小特征值和对应的特征向量

数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介：当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。

矩阵A 需要满足的条件为： (1) 的特征值为A i n λλλλ,0||...||||21≥≥≥>(2) 存在n 个线性无关的特征向量，设为n x x x ,...,,211.1计算过程：i ni i i u xx αα,1)0()0(∑==，有对任意向量不全为0，则有1111112211211111111011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k ni ik i i ni i i k )(k (k))(k αλλλλλα++++=+=+++≈⎥⎦⎤⎢⎣⎡+++======∑∑Λ 可见，当||12λλ越小时，收敛越快；且当k 充分大时，有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k(k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。

2 算法实现.,, 3,,1 , ).5()5(,,,,||).4();max(,).3()(max(;0,1).2(,).1()()()(停机否则输出失败信息转置若转否则输出若计算最大迭代次数，误差限，初始向量输入矩阵βλβεβλβλε←+←<<-←←=←←k k N k y x Ay x x abs x y k N x A k k k3 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 % 判断第一次迭代后是否满足要求t=max(x);return;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值t=x(index); % 是原值，而非其绝对值。

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

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

设置初始向量为x0=ones(15,1)，结果显示如下可见，结果正确。

得到了15阶Hilb 矩阵的按模最大特征值和对应的特征向量。

二.反幂法1.反幂法简介及其理论在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。

其基本理论如下，与幂法基本相同：x x A x A x x Ax λλλ1)(11==⇒=--则，由，可知，A 和A -1的特征值互为倒数，求A 按模最小特征值即求A -1的按模最大特征值,取倒数即为A 的按模最小特征值所以算法基本相同，区别就是在计算)1()()1()(1-)1()()1()1(;,+++++===k k k k k k k k x LU A y Ax y A x Ay x x 分解，来计算做对具体计算时，变换为而是时，不是令2.算法实现., );4( ,))(max(, ,1,).7(),7(,,,1,||).6( ),max().5(),,().4().3(,))(max(,0 ,1).2(,,,,).1(000停机否则输出失败信息转置若否则转停机输出若解方程组作三角分解置最大迭代次数误差限初始向量输入矩阵x abs xy k k N k y x z Ux y Lz y LUx LUA x abs xy k N x A ←←+←<<-←←====←←←λλλελλμλμλε3 matlab 程序代码function [s,y]=invpower(A,x0,eps,n) % s 为按模最小特征值，y 是对应特征向量k=1; r=0; % r 相当于0λ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 % 判断第一次迭代后是否满足终止条件returnendwhile abs(u-r)>eps && k<n % 终止条件.k=k+1;r=u;y=x./max(abs(x));z=L\y;x=U\z;u=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值s=1/x(index); % 是原值，而非其绝对值。

end4 举例验证同幂法一样，选取一个矩阵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 A A cond T 2121,,)(2λλλλ=,而如果A 为对称矩阵，如Hilb 矩阵，A A T 的最大最小特征值，分别为A 的最大最小特征值的平方。

所以cond(A) 为A 的最大最小特征值得比值。

对于本例中的15阶Hilb 矩阵来说，利用上面计算结果得其条件数(选择第二种条件数）为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017 在同一数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明Hilb 矩阵是病态矩阵。

四.Aitken 商加速法1.简介与原理{}{}.,ˆˆ:2)( 2)(,,;0lim,122112212221211加速法这种方法称为逼近用有充分大时当线性收敛即且收敛与若Aitken a aaa a a a a a a x x x x x x y a a aa a a a a k a c a a aa a a k k kk k k k k nn n n n n n k k k k k kk k k =+---≈⇒+---=--≈--≠=--+++++++++++++∞→同幂法和反幂法计算最大和最小特征值类似，如果计算最大特征值，则迭代格式为)()1(k k Ay x =+；计算最小特征值时，迭代格式为)1()()(1)1(,+-+==k k k k Ax y y A x 即。

2.算法实现计算按模最大特征值算法如下：.,,),3()max(,,1 , , ,,).6(),6(,,,).5(,2)().4(,)max( , ).3(,)max( ,0.1 ,0 ,0 ,1).2(,,,),().1(01201001220102010停机输出失败信息否则转置若否则转停机输出若计算置计算置最大迭代次数误差限初始向量输入y x xk k N k y x Ay x x xy k N x a A ij ⇒⇒+⇒⇒⇒<<-+---=⇒=======λλααααλελλααααααλαλααε类似幂法和反幂法可以写出按模最小特征值算法，此处不再赘述。

3.matlab 程序代码function [r,y]=aitken(A,x0,eps,n) % r 按模最大特征值,y 为对应特征向量 k=1;a0=0; % a 相当于 0α a1=1; % a1 相当于 1α r0=1; % 相当于2中的 0λy=x0./max(abs(x0)); % 规范化初始向量 x=A*y;a2=max(abs(x)); % a2相当于2αr=a0-(a1-a0)^2/(a2-2*a1+a0); % 相当于 λif (a2-2*a1+a0)==0 % 若上式中分母为0，则迭代失败，返回 disp "初始向量迭代失败"return;endif abs(r-r0)<eps% 判断第一次迭代后是否满足要求，如满足，则返回结果returnendwhile abs(r-r0)>eps && k<n % 终止条件k=k+1;a0=a1;a1=a2;r0=r;y=x./max(abs(x));x=A*y; % 迭代格式a2=max(abs(x));if (a2-2*a1+a0)==0 % 若分母为0，则迭代失败，返回return;endr=a0-(a1-a0)^2/(a2-2*a1+a0);[m,index]=max(abs(eig(A))); % 以下代码保证取出来的按模最大特征值aa=eig(A); % 是原值，而非其绝对值。

if aa(index)>0 ||aa(index)==0r=r;e lser=-r;endendend类似可得按模最小特征值和特征向量的代码如下：与上面类似，所不同的只是迭代格式不同.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);i f (a2-2*a1+a0)==0disp "初始向量迭代失败"return;endif abs(r-r0)<eps% 判断第一次迭代后是否满足要求，如满足，则返回结果returnendwhile abs(r-r0)>eps && k<nk=k+1;a=b;b=c;r0=r;y=x./max(abs(x));z=L\y;x=U\z;a2=max(abs(x));if (a2-2*a1+a0)==0return;endr=a0-(a1-a0)^2/(a2-2*a1+a0);end[m,index]=min(abs(eig(A))); % 以下代码保证取出来的按模最大特征值aa=eig(A); % 是原值，而非其绝对值。