第5章 矩阵特征值与特征向量的计算n阶方阵A的特征值是特征方程det(A-E)=0 的根. Gerschgorin圆盘定理 设矩阵A=(aij)nn ,记复平面上以aii为圆心,以ri=为半径的n个圆盘为Ri=aiiri,i=1,2,n A的特征向量是齐次线性方程组(A-E)x=0 的非零解. 则 (1)A的任一特征值至少位于其中一个圆盘内;(2)在m个圆盘相互连通(而与其余n-m个圆盘互不连通 )的区域内,恰有A的m个特征值(重特征值按重数记).试讨论A的特征值的分布.解 由A确定的3个圆盘分别为所以315-22n , 这时,(5.1)式可写成若a10, 则对充分大的k有 因而有 或取 而特征向量 x1 v(k).乘幂法的收敛速度取决于|2/1|的大小. 求矩阵A的按模最大的特征值解 取v(0)=(1,0)T ,计算v(k)=Av(k-1), 结果如下例2kv1(k)v2(k)v1(k)/v1(k-1)v2(k)/v2(k-1)010 10.250.2 20.102500.0833330.410.41665 30.0422920.0343890.412600.41267 40.0174510.0141900.412630.41263可取0.41263 ,x1(0.017451,0.014190)T .对非零向量v,用max(v)表示v的按绝对值最大的分量,称向量u=v/max(v)为向量v的规范化向量. 例如, 设向量 v=(2,1,-5,-1)T,则max(v)=-5,u=(-0.4,-0.2,1,0.2)T.可见规范化向量u总满足u=1.乘幂法的规范化计算公式为:任取初始向量u(0)=v(0) 0,计算可得所以其收敛速度由比值|2/1|来确定. 又由于所以因此,当k充分大时可取: 1 k , x1 u(k) .如用规范化乘幂法解例2,仍取u(0)=v(0)=(1,0)T,则有故可取 1 0.412627, x1 (1,0.813138)T.k01234k0.250.410.4126020.412627 u1(k)11111 u2(k)00.80.8130080.8131360.813138用乘幂法求A的按模最大的特征值和相应特征向量.例3 设解 取初值u(0)=v(0)=(1,1,1)T,计算得k ku(k) 0 1 2 3 10 11 1210 7.2 6.5 6.003352 6.001675 6.000837(1,1,1)T (1,0.8,0.1)T (1,0.75,-0.111)T (1,0.730769,-0.188034)T . (1,0.714405,-0.249579)T (1,0.714346,-0.249790)T (1,0.714316,-0.249895)T可取 1 6.000837, x1 (1,0.714316,-0.249895)T.实际上,A的3个特征值分别为1=6,2=3,3=2.2. 设1=2=r,且 |1r+1n ,这时,(5.1)式可写成若a1,a2,ar不全为零, 则对充分大的k有 由于a1x1+a2x2+arxr 是对应1的特征向量, 若仍记为x1 ,则有: v(k) 1kx1 ,故前面的结论仍然成立.3. 设1=-2,且 |1=|2|3 n ,这时,(5.1)式可写成则对充分大的k有 v(2i)12i(a1x1+a2x2) , v(2i+1)12i+1(a1x1-a2x2)于是有x1v(k+1)+1v(k) , x2v(k+1)-1v(k)的按模最大特征值和相应的特征向量。例4 用乘幂法求矩阵解 取初始向量u(0)=v(0)=(1,1,2)T ,计算可得k ku(k) 0 1 2 3 4 5 6 7 8 9 10 1111 3.553628 4.679204 3.461124 4.635465 3.452655 4.632116 3.454315 4.631929 3.454291 4.631920(1,1,2)T (0.454545, 0.909091, 1)T (0.537222, 0.972116, 1)T (0.465201, 0.994041, 1)T (0.539392, 0.998269, 1)T (0.465721, 0.999627, 1)T (0.539487, 0.999892, 1)T (0.465890, 0.999975, 1)T (0.539495, 0.999993, 1)T (0.465893, 0.999999, 1)T (0.539495, 1, 1)T (0.465893, 1, 1)T1.2 加速技术由于所以,乘幂法收敛速度取决于比值|2/1|,当|2/1|1时,收敛是很慢的.1.Aitken 加速方法 由(5.2)式可知12 133.454288 4.631924(0.539495, 1, 1)T (0.465893, 1, 1)Tx2=13u(13)-1u(12)=(0 , 0.631924 , 0.631924)T.x1=13u(13)+1u(12)=(4.315961, 8.631924, 8.631924)T,实际上,A的特征值为1=4,2=-4,3=1.可见,序列k线性收敛于1 .会达到加速收敛的目的.构造Aitken序列如把Aitken加速方法用于例3,则有 ku(k)10 7.2 6.5 6.003352 6.001675 6.000837(1,1,1)T (1,0.8,0.1)T (1,0.75,-0.111)T (1,0.730769,-0.188034)T . (1,0.714405,-0.249579)T (1,0.714346,-0.249790)T (1,0.714316,-0.249895)Tk 0 1 2 3 10 11 126.266667 6.000017 6.000003 6.0000002.原点位移法作矩阵B=A-pE, 则B的特征值为mi=i-p(i=1,2,n),而且对应的特征向量相同.则对B应用乘幂法可达到加速收敛的目的。解 由于A的特征值为1=6,2=3,3=2,故取p=2.5,则B的特征值为m1=3.5,m2=0.5,m3=-0.5,则如果选取p,使m1仍然是B的按模最大特征值，且满足取初始向量u(0)=v(0)=(1,1,1)T,由规范化计算公式:例5 用原点位移法求例3中矩阵A的按模最大的特征值和特征向量.计算可得kkU(k) 0 1 2 3 4 5 67.5 3.766662 3.535396 3.505002 3.500718 3.500102(1,1,1,)T (1,0.733333,-0.2)T (1,0.716814,-0.238938)T (1,0.714643,-0.249061)T (1,0.714337,-0.249777)T (1,0.714293,-0.249981)T (1,0.714287,-0.249995)T这是因为|2/1|=1/2,而|m2/m1|=1/7,故对B应用乘幂法远比对A应用乘幂法收敛的快.反幂法是求矩阵按模最小的特征值和相应特征向量的方法.取,16+2.5=6.000102,x1u(6)=(1,0.714287,0.249995)T1.3 反幂法设A是n阶非奇异矩阵, 其特征值为|1| |2| |n-1| |n| 0对应的特征向量为x1,x2,xn, 则有A-1的特征值为对应的特征向量为xn,xn-1,x1.要想求n和xn只需对A-1应用 乘幂法,任取初始向量u(0)=v(0)0, 作也可将上式改写成式(5.3)称为反幂法. 显然有每一步求v(k)需要求解线性方程组, 可采用LU分解法求解.反幂法还可结合原点位移法应用.设已求得矩阵A的特征值i的某个近似值对B应用反幂法可求出精度更高的i和xi.设已求得例3中矩阵A的特征值的近似值16.003,和相应的特征向量x1(1,0.714405,-0.249579)T, 试用带原点位移的反幂法求1和x1的更精确的值.作原点位移,令B=A- E, 则B的特征值为 例6解 取p=6.003, 作矩阵B=A-6.003E,则取初始向量u(0)=(1,0.714405,-0.249579)T,对B用反幂法计算可得:可见收敛速度非常快,这是因为B的3个特征值为1=-4.003, 2=-3.003,3=-0.003,|3/2|0.000999很小.Jacobi方法是求实对称矩阵全部特征值和特征向量的一种矩阵变换方法。2 Jacobi 方法实对称矩阵A具有下列性质： （1）A的特征值均为实数； （2）存在正交矩阵R，使RTAR=diag(1,2,n),而 R的第i个列向量恰为i的特征向量;直接求正交矩阵R是困难的 . Jacobi提出用一系列所谓平面旋转矩阵逐次将A约化为对角矩阵.平面解析几何中的平面坐标旋转变换表示平面上坐标轴旋转角的变换.（3）若记A1=RTAR,则A1仍为对称矩阵. 2.1 平面旋转矩阵 在三维空间直角坐标系中,ox1y1平面绕着oz1轴旋转角的坐标变换为一般地, 在n维向量空间Rn中, 沿着xpyq平面旋转角的变换矩阵为称Rpq()为平面旋转矩阵. Rpq()具有下列性质: 设实对称矩阵A=(aij)nn ,记B=RpqT()ARpq()=(bij)nn则它们元素之间有如下关系:（1）Rpq()为正交矩阵,即Rpq-1()=RpqT(); （2）如果A为对称矩阵, 则RpqT()ARpq()也为对称矩阵, 且与A有相同的特征值. （3）RpqT()A仅改变A的第p行与第q行元素,ARpq()仅改变A的第p列与第q列元素.所以有从而有(5.5)、(5.6)式可得如果apq0, 适当选取角, 使只需角满足从而如果取|apq|=若记于是则上式可记为由式(5.7),令t=tan,则t满足方程t2+2t-1=0经典Jacobi算法是对A(0)=A施行一系列平面旋转变换:为保证|/4,取绝对值较小的根,有于是2.2 Jacobi 方法A(1)=R1TA(0)R1 ,A(2)=R2TA(1)R2 , A(k)=RkTA(k-1)Rk ,每一步变换选择A(k-1)=(aij(k-1)nn 的非对角线元素中绝对值最大者apq(k-1)(称为主元素)作为歼灭对象, 构造平面旋转矩阵Rk=Rpq(), 经变换得到A(k)=(aij(k)nn ,且apq(k)=0,这时由(5.8)式有从而由此递推得到当k充分大时,或者(A(k),或者另外,由于A(k)=RkTA(k-1)Rk=RkTRk-1TR1TAR1R2Rk=RTAR 是给定的精度要求,则A的特征值可取为iaii(k),i=1,2,n.的全部特征值.解 记 A(0)=A,取p=1,q=2,apq(0)=a12(0)=2,于是有因此,R=R1R2Rk 的