第6章特征值和特征向量修改2

1、1 1 引引 言言第第6 6章章 矩阵特征值数值计算矩阵特征值数值计算物理、力学和工程技术中很多问题在数学上都归结物理、力学和工程技术中很多问题在数学上都归结为求矩阵的为求矩阵的特征值问题。例如，振动问题特征值问题。例如，振动问题(大型桥梁大型桥梁或建筑物的振动、机械的振动、电磁震荡等或建筑物的振动、机械的振动、电磁震荡等)，物理，物理学中的某些临界值的确定学中的某些临界值的确定及一些稳定性分析和相关及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量分析在数学上都可转化为求矩阵特征值与特征向量的问题的问题 下面先复习一些矩阵的特征值和特征向量的基础下面先复习一些矩阵的特征值和特

2、征向量的基础知识知识. 定义定义1 已知已知n阶矩阵阶矩阵A=(aij)，则，则)2()(det)det()(12211212222111211的项的项次数次数 naaaaaaaaaaaaAInnnnnnnnnn 称为称为A的的特征多项式特征多项式. 一般有一般有n个根个根(实的或复的，复根按重数计算实的或复的，复根按重数计算)称为称为A的的特征值特征值. 用用(A)表示表示A的所有特征值的集合的所有特征值的集合. A的特征方程的特征方程)1 . 1(0)det()( AI 设设为为A的特征值，相应的齐次方程组的特征值，相应的齐次方程组 注：注：当当A为实矩阵时，为实矩阵时， ()=0为实系数

3、为实系数n次代数次代数方程，其复根是共轭成对出现方程，其复根是共轭成对出现.)2 . 1(0)( xAI 的的非零解非零解x称为矩阵称为矩阵A的对应于的对应于的的特征向量特征向量. 210131012A 例例1 求求A的特征值及特征向量，其中的特征值及特征向量，其中. 0)4)(2)(1(8147210131012)det()(23 AI 解解 矩阵矩阵A的特征方程为的特征方程为求得矩阵求得矩阵A的特征值为：的特征值为：. 4, 2, 1 对应于各特征值矩阵对应于各特征值矩阵A的特征向量分别为的特征向量分别为：.121,101,111321 xxx 定理定理1 设设为为ARnn的特征值的特征值

4、, 且且Ax=x (x 0)，则有，则有 - -p为为A- -pI的特征值，即的特征值，即(A- -pI)x=(- -p)x ； c为为cA的特征值的特征值(c0为常数为常数)； 下面叙述有关特征值的一些下面叙述有关特征值的一些结论结论： k为为Ak的特征值，即的特征值，即Akx=kx ； 设设A为非奇异矩阵，那么为非奇异矩阵，那么0 , 且且- -1为为A- -1的的特征值，即特征值，即A- -1x=- -1x . 定理定理2 设设i(i=1,2,n)为为n阶矩阵阶矩阵A=(aij)的特征的特征值，则有值，则有)(Atraniiinii 11 称为称为A的的迹迹； .nA21 定理定理3 设

5、设ARnn，则有，则有. )()(AAT 定理定理4 设设A 为分块上三角矩阵，即为分块上三角矩阵，即, mmmmAAAAAAA22211211其中每个对角块其中每个对角块Aii均为方阵，则均为方阵，则. )()(iiniAA1 定理定理5 设设A与与B为相似矩阵（即存在非奇异矩阵为相似矩阵（即存在非奇异矩阵P使使B=P- -1AP），则），则 定理定理5说明，一个矩阵说明，一个矩阵A经过相似变换，其特征经过相似变换，其特征值不变值不变. A与与B有相同的特征值；有相同的特征值； 如果如果y是是B的特征向量，则的特征向量，则Py是是A的特征向量的特征向量. 定义定义2 如果实矩阵如果实矩阵A有

6、一个重数为有一个重数为k的特征值的特征值, 且对应于且对应于的的A的线性无关的特征向量个数的线性无关的特征向量个数1)或下溢成为或下溢成为0 (| 1|2|n|，则对任意非零初，则对任意非零初始向量始向量v0=u0(a1 0)，有幂法计算公式为，有幂法计算公式为则有则有 ,)max(lim11xxukk .lim1 kk 即即：当当k充分大时充分大时1 k)max(11)(uuyk 证明：证明： 设有一向量设有一向量v 0，将其，将其规范化规范化得得向量向量为为其中其中max(v)表示表示v的的绝对值最大的分量绝对值最大的分量. 即如果有即如果有,)max(vvu ,max1iniqvv 则则

7、max(v)=vq，且，且q为所有为所有绝对值最大的分量中的最绝对值最大的分量中的最小下标小下标.在定理在定理6.3的条件下幂法可这样进行：任取一初的条件下幂法可这样进行：任取一初始向量始向量v0 0(a1 0)，构造，构造规范化规范化向量序列为向量序列为 .)max()max(,)max(.)max()max(,)max()max()max(,00010102022220021200111001vAvAvvuvAvAAuvvAvAvvuAvvAAuvAvAvvvuAvAuvkkkkkkkkk，由由 得得)8 . 2(.)/(211110 niikiikkxaxavA 01 1221(0),n

8、nva xa xa xa设 niikiikniikiikkkkxaxaxaxavAvAu211112111100)/(max)/()max( ).()max()/(max)/(1121112111 kxxxaxaxaxaniikiiniikii 收敛速度收敛速度由比值由比值r=|2/1|确定确定. 总结上述结论，有总结上述结论，有同理，可得到同理，可得到 niikiikniikiikkkkxaxaxaxavAvAv211111121111010)/(max)/()max( ).()/(max)/(maxmax12111121111 kxaxaxaxavniikiiniikiik )9 . 2(

9、 .), 2 , 1( ,max, 0100 向向量量的的规规范范化化kkkkkkk/vukvAuvuv 定理定理设设ARnn有有n个线性无关的特征向量，主个线性无关的特征向量，主特征值特征值1满足满足|1|2|n|，则对任意非零初始向，则对任意非零初始向量量v0=u0(a1 0)，有幂法计算公式为，有幂法计算公式为则有则有 ,)max(lim11xxukk .lim1 kk例：用幂法计算矩阵例：用幂法计算矩阵210121012A模最大的特征值及其对应的特征向量。模最大的特征值及其对应的特征向量。则得）（解：取初始向量,1,11T)0()0( vu)0()1(Auv11121012101210

10、11m)max()1(v1)1(uTmv1 , 0 , 11)1()1()2(Auv1012101210122221k2k2m)max()2(v2)2(uTmv1 , 1, 12)2()2()3(Auv1112101210123433m)max()3(v4)3(uTmv75. 0, 1 ,75. 03)3(3k4k)3()4(Auv75. 0175. 02101210125 . 25 . 35 . 24m)max()4(v5 . 3)4(uTmv75, 1 ,754)4(38例例 用幂法求矩阵用幂法求矩阵 210120012A的按模最大的特征值和相应的特征向量的按模最大的特征值和相应的特征向量

11、. 取取 x(0)=(0, 0, 1)T, 要求误差不超过要求误差不超过10 3.解解 ,1 , 0 , 000Txy , 2)max( ,2 , 1, 0)1(101 xAyxT Txy)1, 5 . 0, 0(1)1()1( , 5 . 2)max( ,5 . 2, 2, 5 . 0)2(2)1()2( xAyxT 3938910 0006049. 099909241. 29996973. 2 Txy)1, 8 . 0, 2 . 0(2)2()2( , 8 . 2)max( ,8 . 2, 6 . 2, 2 . 1)3(3)2()3( xAyxT .9996973. 2,9990924.

12、2,9972799. 2,9918619. 2,9756097. 2,9285714. 2987654 .9996973. 291 6.1.2 6.1.2 幂法的其他复杂情况幂法的其他复杂情况1.1.我们假设了我们假设了A具有完全的特征向量系具有完全的特征向量系, ,即即A具有具有n个个线性无关的特征向量线性无关的特征向量. .当当A不具有不具有n个线性无关的特个线性无关的特征向量时征向量时, ,幂法不适用幂法不适用, ,但事前往往无法判断这一点但事前往往无法判断这一点. .因此在运用幂法时因此在运用幂法时, ,发现不收敛或收敛很慢情况发现不收敛或收敛很慢情况, ,要考虑此种可能要考虑此种可能

13、. .2.2.我们假设在我们假设在(6.3)(6.3)中中1 10,0,这在选择这在选择u0 0时时, ,也无法也无法判断判断, ,但这往往不影响幂法的成功使用但这往往不影响幂法的成功使用. .因为若选因为若选u0 0, ,使使1 1=0,=0,由于舍入误差的影响由于舍入误差的影响, ,在迭代某一步会在迭代某一步会产生产生uk k, ,它在它在x1 1方向上的分量不为零方向上的分量不为零, ,这时以后的迭这时以后的迭代仍会收敛代仍会收敛. .A3.如果矩阵 的特征值不满足假设，n21此时幂法的收敛性分析就变得复杂，它可能产生此时幂法的收敛性分析就变得复杂，它可能产生如下情况：如下情况：，且11

14、21rr（1），且3121（2），且3121（3）112(1)111111 111111(1)111 11 11, ()(), ()mkkkkmnmmmmnnmkkmmmmmAnxuuuukxuuuuA（ ） 是 重根，即矩阵 仍有 个线性无关的特征向量。此时有 显然，只要不全为零，当 充分大时，就有因也是矩阵 相1(1)12( )(1) kimkikxxx应于 的特征向量，故有为相应的特征向量，即对这种情况幂法仍然有效。1213(1)1111311 122331122k 2(2)11 12231( )k11 122312, ( 1)()()( 1)()limlim( 1)()kkkkknnn

15、knkikiijjikknkkjkiiijiAnxuuuuxxxxxxxx （ ）且矩阵 有 个线性无关的特征向量。21(1)(2)2(2)( )11,2( )(1)1111 122( )11 122/ ( 1) ( 1)kkkkiiikikkkkkkxkxxxxxuuxuu 由上式可知，是个摆动序列，当 充分大时，有又由(3）若）若 ，则连续迭代两次则连续迭代两次，计算出计算出x(k+1)，x(k+2)，210)()1()2( kjkjkjqxpxx然后对然后对j = 1, 2, , n 解方程解方程求出求出 、 后后，由公式由公式pq2122 pqip 2222 pqip 解出主特征值解出

16、主特征值 1、 2。此时收敛速度取决于此时收敛速度取决于 的程度的程度。113r)(2)1(kkxx)(1)1(kkxx向量向量 、 分别为相应于分别为相应于 1， 2的特征向量的近似值的特征向量的近似值。1求大型稀疏矩阵按模最大的特征值,1为为单单根根或或重重根根矩矩阵阵的的特特征征值值分分布布为为 nm ,11m1 即即幂法适用范围：幂法适用范围：.个个线线性性无无关关的的特特征征向向量量矩矩阵阵有有 n6. 2 幂法的加速方法1、原点平移法、原点平移法 由前面讨论知道，应用幂法计算由前面讨论知道，应用幂法计算A的主特征值的的主特征值的收敛速度主要由比值收敛速度主要由比值 r=|2/1|(

17、收敛因子收敛因子)来决定，但来决定，但当当r 接近于接近于1时，收敛可能很慢时，收敛可能很慢. 这时，一个补救办法这时，一个补救办法是采用加速收敛的方法是采用加速收敛的方法.其中其中p为参数，设为参数，设A的特征值为的特征值为 i，则对矩阵，则对矩阵B的特征的特征值为值为 i- -p ，而且，而且A, B的特征向量相同的特征向量相同. 引进矩阵引进矩阵 B=A- -pI .不妨设不妨设 1 2 n ， 1为为A按模最大的特征按模最大的特征值值， 2为为A按模次大的特征值按模次大的特征值，且且 | 2 | | n |。 1 2 nOp = ( 2 + n ) / 2令令 B = A pE ，则有

18、则有 | AE | = | (B+pE)E| = | B (p)E | A p = B A = B + p 。选择选择p使得使得 1 p为为B按模最大的特征值按模最大的特征值， 2 p为为B按按模次大的特征值模次大的特征值，且且思路：思路：1212max ppini那么，对矩阵那么，对矩阵B=A- -pI应用幂法求其主特征值应用幂法求其主特征值 1- -p, 收收敛速度将会加快敛速度将会加快. 这种通过求这种通过求B=A- -pI的主特征值和特的主特征值和特征向量，而得到征向量，而得到A的主特征值和特征向量的方法叫的主特征值和特征向量的方法叫原原点平移法点平移法. 对于对于A的特征值的某种分布

19、，它是十分有的特征值的某种分布，它是十分有效的效的.例例4 设设AR44有特征值有特征值),4 , 3 , 2 , 1(15 jji 比值比值r=|2/1|0.9. 做变换做变换B=A- -12I (p=12),则则B的特征值为的特征值为. 1, 0, 1, 24321 应用幂法计算应用幂法计算B的主特征值的主特征值1的收敛速度的比值为的收敛速度的比值为. 9 . 021121212 pp 虽然常常能够选择有利的虽然常常能够选择有利的p值值, 使幂法得到加速使幂法得到加速, 但设计一个自动选择适当参数但设计一个自动选择适当参数p的过程是困难的的过程是困难的. 下面考虑当下面考虑当A的特征值是实

20、数时，怎样选择的特征值是实数时，怎样选择p使使采用幂法计算采用幂法计算1得到加速得到加速.ppn 1且使且使收敛速度的比值收敛速度的比值.min,max112 ppppn 设设A的特征值都是实数，且满足的特征值都是实数，且满足)10. 2(,121nn 则对实数则对实数p，使矩阵，使矩阵A- -pI的主特征值为的主特征值为 1- -p或或 n- -p时时，当，当我们计算我们计算 1及及x1时，首先应选取时，首先应选取p使使显然，当显然，当 2- -p=- -( n- -p )时，即时，即 P=( 2+ n)/2P* 时时为最小值，这时为最小值，这时收敛速度的比值收敛速度的比值为为.221211

21、2nnnpppp 当希望计算当希望计算 n时，应选取时，应选取 p=( 1+ n-1)/2P* 使得应用幂法计算使得应用幂法计算 n得到加速得到加速. 当当A的特征值都是实数，满足的特征值都是实数，满足且且 2, n能初步估计出来，我们就能确定能初步估计出来，我们就能确定P*的近似值的近似值.nn 121 例例2 用原点平移加速法求用原点平移加速法求例例1中矩阵中矩阵A的主特征值的主特征值与其对应的特征向量与其对应的特征向量.5 . 36345 . 510235 . 4 B,1634310232 A对对B应用幂法，仍应用幂法，仍取取 v0=(0,0,1)T , 则则 .875. 0 , 1,

22、5 . 01, 4,5 . 3 , 4 , 211111TTvuu ), 2 , 1(max1 k/vuvBuvkkkkkkk 解解 取取p=- -2.5, 做平移变换做平移变换B=A- -pI，则，则迭代迭代5步的计算结果步的计算结果见下表见下表k1 2, 4, 3.54 0.5, 1, 0.8752 7, 14, 10.5625140.5, 1, 0.75453 6.76, 13.5179, 10.140613.5179 0.5, 1, 0.75074 6.7503, 13.5007, 10.125613.5007 0.5, 1, 0.75005 6.7500, 13.5000, 10.1

23、25013.5000 0.5, 1, 0.7500TkuTkv可得到可得到B的主特征值为的主特征值为 1 13.5000, 主特征向量为主特征向量为 v1 (0.5 ,1.0, 0.7500)T ,因此，因此，A的主特征值为的主特征值为 1 = 1 +p 11.0000, 主特征向量仍为主特征向量仍为 x1 =(0.5,1,0.7500)T .k 原点位移的加速方法，是一个矩阵变换方法原点位移的加速方法，是一个矩阵变换方法. 这这种变换容易计算，又不破坏矩阵种变换容易计算，又不破坏矩阵A的稀疏性，但的稀疏性，但p的的选择依赖对选择依赖对A的特征值分布的大致了解的特征值分布的大致了解.0(0)(

24、1)( )00(4)4 140 5 1302.9,102.8(1,1,1) ,()6.9140 510.10100.1(3.1000568,2.214326, 0.968TkkAAxxAI xAIx -4例：,用原点移位法求矩阵 的按模最大的特征值，要求误差不超过10 。解：取按进行计算 4(5)54547661) 3.1000568(3.0999984,2.2142846, 0.9687501) 3.09999840.000058410 x 1123212010 3.09999842.95.9999984 (3.0999984,2.2142846, 0.9687501) . 6,3,2.8,

25、1,20.11 3.131TAxAA所以，矩阵 的按模最大的特征值为相应的特征向量为 不难求出， 的特征值为若对 直接用幂法，则比值 而用原点移位法，则有因此收敛速度明显加快。 1211222112121 lim0 ()22kkkkkkkkkkkkkkkkkkkkkkkaaaacaaaaaakaaaaaaaaaaaaaaaaaaaaaAitkena如果序列线性收敛到 ，即则当 充分大时，有序列比更快地收敛到 ，这就是加速法。将这一方法用于幂法所产生的序列 k，可加快幂法的收敛速度。2、埃特肯、埃特肯(Aitken)加速法加速法，收敛速度很慢，当1r幂法的加速幂法的加速速度的主特征值时，其收敛使

26、用幂法求矩阵A,12的大小取决于r加速技术。此时，可以采用Aitken，）任取（01)0(v)()2(kv)(kukkmv)(km, 2 , 1k)max()(kv,)max()0(1)0(vAvAkkkmkkkkkkmmmmmm1221222)()1( kAu加速时，用当Aitkenk3例：用幂法计算矩阵例：用幂法计算矩阵1634310232A的主特征值及对应的特征向量。并对计算主特征值的迭代的主特征值及对应的特征向量。并对计算主特征值的迭代进行进行Aitken加速加速。，）（解：取T)0(,0,10v)0()1(Auv1m)max()1(v4)1(uTmv25. 0, 1, 5 . 01)

27、1(1k1634310232100142，)0()0(uv)1()2(Auv2m)max()2(v9)2(uTmv8611. 0, 1, 5 . 02)2(2k163431023225. 015 . 075. 795 . 4)2()3(Auv3m)max()3(v4444.11)3(uTmv7306. 0, 1, 5 . 03)3(3k16343102328611. 015 . 03611. 84444.1172. 51m12322332)(mmmmmm4924444.11)94444.11(4444.1123380. 24444.117824.13)3()4(Auv4m)max()3(v92

28、24.10)4(uTmv7535. 0, 1, 5 . 04)4(4k16343102327306. 015 . 02306. 89224.104612. 52m23423442)(mmmmmm94444.1129224.10)4444.119224.10(9224.1020641.111417. 09224.10k01(0,0,1.0000)14(0.5000,1.000,0.2500)29(0.5000,1.000,0.8611)311.4444(0.5000,1.000,0.7306)13.7824410.9224(0.5000,1.000,0.7535)11.0641511.0140(

29、0.5000,1.000,0.7493)11.0001610.9927(0.5000,1.000,0.7501)711.0004(0.5000,1.000,0.7500)811.0000(0.5000,1.000,0.7500)Tku)(kmkm 定义定义4 设设ARnn为对称矩阵，对于任一非零为对称矩阵，对于任一非零向量向量x，称，称,),(),()(xxxAxxR 为对应于向量为对应于向量x的的瑞利瑞利(Rayleigh)商商.其中其中(x, x)=xTx为为内积内积. 定理定理 设设ARnn为对称矩阵为对称矩阵(其特征值次序记其特征值次序记为为12n)，则，则 1. (对任何非零对任何非

30、零xRn)；1 ),(),(xxxAxn 2. ；),(),(maxxxxAxxRxn01 3. .),(),(minxxxAxxRxnn0 3、Rayleigh商加速法商加速法 证明证明 只证只证1，关于，关于2, 3自己作练习自己作练习. 由于由于A为实对称矩阵，可将为实对称矩阵，可将 1,2,n 对应的特对应的特征向量征向量 x1,x2,xn 正交规范化，则有正交规范化，则有 (xi, xj)=ij，设，设x 0为为Rn中任一向量，则有中任一向量，则有. 0,1221 niiniiixxx 于是于是.),(),(11212 niiniiinxxxAx从而从而1成立成立. 结论结论1说明说

31、明瑞利商瑞利商必位于必位于n和和1之间之间.下面我们将瑞利商应用到用幂法计算实对称矩阵下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征的主特征值的加速上来值的加速上来 定理定理 设设ARnn为为对称矩阵对称矩阵，特征值满足，特征值满足对应的特征向量对应的特征向量xi满足满足(xi, xj)=ij (单位正交向量单位正交向量) ，应用幂法公式应用幂法公式(2.9)计算计算A的主特征值的主特征值 1，则规范化，则规范化向量向量uk的的瑞利商瑞利商给出给出 1的较好的近似值为的较好的近似值为,121nn kkkkkkuuuAuuR2121, 由此可见，由此可见，R(uk)比比k更快的收敛于更快的

32、收敛于 1. 证明证明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(2121122112200001 knjkjjnjkjjkkkkkkkkkaauAuAuAuAuuuAuuR 得得 幂法的幂法的瑞利商加速迭代公式瑞利商加速迭代公式可以写为可以写为 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A为为n阶实对称矩阵阶实对称矩阵.,11kkux 对给定的误差限对给定的误差限 ，当，当| kk- -1| 时，取近似值时，取近似值三、反幂法 反幂法是计算矩阵按模

33、最小的特征值及特征向量反幂法是计算矩阵按模最小的特征值及特征向量的方法，也是修正特征值、求相应特征向量的最有效的方法，也是修正特征值、求相应特征向量的最有效的方法。的方法。11111 , 1 An nuAAuuuA uA uuAAAAA设 为阶非奇异矩阵，为 的特征值与相应的特征向量，即此式表明，的特征值是 的特征值的倒数，而相应的特征向量不变。因此，若对矩阵用幂法，即可计算出的按模最大的特征值，其倒数恰为 的按模最小的特征值。这就是反幂法的基本思想。 ), 2 , 1(max11 k/vuvuAvkkkkkkk 反幂法的迭代公式反幂法的迭代公式为为knknux ,1 对给定的误差对给定的误差

34、 ，当，当|kk- -1| n|0，则对任意非零初始向量则对任意非零初始向量u0(an 0) ，由反幂法计算公，由反幂法计算公式构造的向量序列式构造的向量序列vk，uk满足满足 ,)max(limnnkkxxu .1)max(limnkkv 1(1)( )(1)1( )(1) kkkkkAAAxxxA xxA因为的计算比较麻烦，而且往往不能保持矩阵 的一些好性质（如稀疏性），因此，反幂法在实际运算时以求解方程组 代替幂法迭代求得，每迭代一次要解一个线性方程组。由于矩阵在迭代过程中不变，故可对 先进行三角分解，每次迭代只要解两个三角形方程组。反幂法计算的主要步( )( )( )( )( )1(

35、)(1)1.2.max, 3. kkkkkriri nkkAALUxrxxxyLzyUxz 骤：对 进行三角分解求整数 ，使得计算解方程组(0)210 021012(0,0,1) . 2.930.9310 2.9300.931010.93100 0100 1/0.931TAAxAIAI例：，用反幂法求矩阵 接近2.93的特征值，并求相应的特征向量，取解：对作三角分解得40.931000.931000.93 1/0.9333.0000954,310(1, 0.9992431,0.9991478)(1,-1,1)0.001.TTur按算法迭代 次，与准确值 的误差小于，与准确值比较，残差 在反幂法

36、中也可以用在反幂法中也可以用原点平移法原点平移法加速迭代过程加速迭代过程，或或求其它特征值与其对应的特征向量求其它特征值与其对应的特征向量. 如果矩阵如果矩阵(A- -pI)- -1存在，显然其特征值为存在，显然其特征值为,1,1,121pppn 对应的特征向量仍然是对应的特征向量仍然是x1,x2,xn，现对矩阵，现对矩阵(A- -pI)- -1应用幂法，得到反幂法的迭代公式应用幂法，得到反幂法的迭代公式)12. 2()., 2 , 1()max(/.,)(, 01100 kvuupIAvvukkkkk 初初始始向向量量 如果如果p是是A的特征值的特征值 j的一个近似值，且设的一个近似值，且设 j与其它与其它特征值是分离的，即特征值是分离的，即),(jippij 就是说就是说1/( j- -p)是矩阵是矩阵 (A- -pI)-