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

矩阵A 需要满足的条件为：⑴I 1 I I 2|n |- 0, i 为A 的特征值(2)存在n 个线性无关的特征向量，设为 X i ,X 2，…，X n 1.1计算过程：n对任意向量x (0)，有x (0)八：-M —不全为0,则有i 4X (k 岀)=Ax (k)== A k 岀乂。

)nn A k 1aq a 扌15i =1i =1■k 12 可见，当1—1越小时，收敛越快；且当k 充分大时，有 ?"12算法实现⑶.计算x Ay,… max(x);⑷若| •一十：；,输出-,y,否则,转(5)(5)若N ,置k 「k 1^-,转3,否则输出失败信息，停机.3 matlab 程序代码(冲1%叫x(k 1)[x(k)k二ux(k)>(k+1)1，对应的特征向量即是x(1).输入矩阵A ，初始向量X ，误差限 最大迭代次数N(k)0; y(k)max(abs(x (k))k=1;z=0;y=x0./max(abs(x0)); x=A*y; % z相当于■%规范化初始向量%迭代格式b=max(x); % b相当于：if abs(z-b)<epst=max(x);return;%判断第一次迭代后是否满足要求endwhile abs(z-b)>eps && k<N k=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,i ndex]=max(abs(x)); %这两步保证取出来的按模最大特征值t=x(i ndex); end%是原值，而非其绝对值。

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

结果如下：» A=[l 1 0. 5;1 1 0. 25;0. 5, 0. 25, 2]A =1.0000 E 00000. 50001.0000 1. 00000, 25000. 50000. 2500 2.0000» x0=[l;l;ll; eps=le-5; 220; >> y]=lpower (A, xO, eps, X)^0. 01661.4801 2.5365» A#y-t*yans -1. 0e-004 +-0.1603 -0, 1684结果正确，表明算法和代码正确，然后利用此程序计算 15阶Hilb 矩阵，与eig(A)的得到结果比较，再计算 A*y-t*y ，验证y 是否是对应的特征向量。

设置初始向量为x0=ones(15,1),结果显示如下» A 二hilb(15); » x0=ones (15,1):】；eps-le^6:» N=30;>> Lt, y] -1 power (A, xO, eps, X}t 二L 84592. 5365 0. 7482 0. 6497 L 0000ans =可见，结果正确。

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

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

其基本理论如下，与幕法基本相同：1由Ax V x= x = A‘(’X),则A’x x，可知，A和A-1的特征值互为倒数，Z求A按模最小特征值即求A-1的按模最大特征值，取倒数即为A的按模最小特征值所以算法基本相同，区别就是在计算x(k1)时，不是令X(k1)= Ay(k),而是X(k J A-1y(k)具体计算时，变换为Ax(k-1^ y(k);对A做LU分解，来计算x(k 1)2■算法实现(1) .输入矩阵A,初始向量x,误差限:,最大迭代次数N, (2) •置k 「1,，° “0, y ,max(abs(x))(3) .作三角分解A = LU(4) .解方程组 LUx = y(Lz = y,Ux = z), (5) . ：— max(x), - ■1(6) 若- ’0卜：；,输出—,y,停机,否则转(7), (7).若k ：： N,置k ・ k 1,•，y ・否则输出失败信息，停机.3 matlab 程序代码fun ctio n [s,y]=i nvpower(A,xO,eps, n) % s 为按模最小特征值，y 是对应特征向量 k=1; r=0;% r 相当于-oy=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 %判断第一次迭代后是否满足终止条件return endwhile 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,i ndex]=max(abs(x)); %这两步保证取出来的按模最大特征值 s=1/x(i ndex); %是原值，而非其绝对值。

endA ，代入程序，得到结果，并与eig(A)的得到结果比ma 莎面，转⑷; 4举例验证同幕法一样，选取一个矩阵较，再计算A*y-t*y，验证y是否是对应的特征向量y 二L 0000 -0. 9517 -0.1300可见结果正确，然后利用此程序计算 15阶Hilb 矩阵，eig(A)的得到结果比 较，再计算 A*y-s*y ，验证y 是否是对应的特征向量。

设置初始向量为x0=ones(15,1),结果显示如下» A=liilb(15): x0=ones (157 1); eps=le~6; n=30; >> [s, y] = invpower (A ? eps^ n)-5.9673e-018» A 二[1 J, 0. 5； 1? 1, 0. 25；0. 5, 0. 25, 2]: xO=[l;l ：j]:eps=le~6;n=30;_s, y] =invpower (A, xOj eps, n)s =-0. 0166» eig(A) ans 二-0. 01661. 0e-016 宕 -0. 1388 -0. 0694 -0. 19081. 48012. 5365ans =可见，结果真确。

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

三•计算条件数矩阵A 的条件数等于A 的范数与A 的逆的范数的乘积，即 cond(A)= II A A(-1) II,对应矩阵的3种范数，可以定义3种条件数。

函数cond(A,1)、cond(A)或cond(A inf)是判断矩阵病态与否的一种度量，条件数越 大表明矩阵的病态程度越大•这里我们选择矩阵的2范数，即cond(A) j 分别为矩阵A TA 的最大和最小特征值而如果A 为对称矩阵，如Hilb 矩阵，A TA 的最大最小特征值，分别为 A 的 最大最小特征值的平方。

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

对于本 例中的15阶Hilb 矩阵来说，利用上面计算结果得其条件数(选择第二种条件数) 为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017在同一0.0000 -0.0000 0.0000 -0.0006 0.0050 -0. 0245 0.0699 -0. 0968 -Q.0421 0. 4497 -0.9084 1.0000 -0. 6527 0.2379» eig(A)ans =-0. 0000 0. 0000 0. 0000 0, 0000 0. 0000 0+ 0000 0+ 0000 0, 0000 0” 00000. 0000 0, 0004 0. 0056 0. 0572 0. 4266>> A*y-s*y ans -L 0e-017 *0. 3903 -0. 5638 0. 0000 -0. 5208 -0. 4741 0. 3540 0. 4537 -0. 2746 CL 9724 -0. 5556 -0.6288 0.6618 0. 0659数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明 Hilb 矩阵是病态矩阵。

四.Aitken 商加速法1.简介与原理若：a k收敛与a,且lim 一二c = 0;即：a k：线性收敛， k护 a k_ a当k 充分大时，有亟口玉心ak -aakd1 - a、， 、， (Xn 雄 一 xn 出)y n =X n 2 -X n 七 一 2X n* + X n(a k 1 一 a k )2 O a kak 2 - 2ak -1 ' a k用a k逼近a,这种方法称为Aitken 加速法.同幕法和反幕法计算最大和最小特征值类似， 如果计算最大特征值，则迭代格式 为x (k ° = Ay (k)；计算最小特征值时，迭代格式为x(k1^A-1y (k),即y (k)= Ax ("2.算法实现计算按模最大特征值算法如下：(1).输入A-©),初始向量x,误差限；,最大迭代次数N,⑵.置k=1,X 0=0, -1=0, ° =1.0, y = ⑶.计算x = Ay,置max(x)=匚2,(4).计算丸= 2(G -Of ) ° ( 1 0丿 0 )_2口 乜(5).若丸一町 ＜名,输出&, y 停机,否则转(6),(6).若k :: N,置-：：1 = 0,二2 = ■■ 1，八=八0,否则,输出失败信息，停机.类似幕法和反幕法可以写出按模最小特征值算法，此处不再赘述3.matlab 程序代码 function [r,y]=aitken(A,xO,eps,n) % r 按模最大特征值$为对应特征向量k=1;k 1= k,,xmax(x)y 转(3),a0=0; % a相当于>0a仁1; % a1相当于:'1r0=1; %相当于2中的‘0y=x0./max(abs(x0)); %规范化初始向量x=A*y;a2=max(abs(x)); % a2相当于：2r=a0-(a1-a0)A2/(a2-2*a1+a0); % 相当于'if (a2-2*a1+a0)==0 %若上式中分母为0,则迭代失败，返回disp "初始向量迭代失败"return;endif abs(r-r0)<eps %判断第一次迭代后是否满足要求，如满足，则返回结果return endwhile 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)A2/(a2-2*a1+a0);[m,index]=max(abs(eig(A))); %以下代码保证取出来的按模最大特征值aa=eig(A); %是原值，而非其绝对值。