第8章矩阵特征值问题的数值解法

1、上页上页下页下页第第8章章 矩阵特征值问题的数值解法矩阵特征值问题的数值解法 8.1 引言引言 8.2 幂法及反幂法幂法及反幂法 8.3 雅可比方法雅可比方法 8.4 QR算法算法上页上页下页下页8.1 引引 言言 工程技术中有多种振动问题，如桥梁或建筑物的工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题征向量的问题. 下面先复习一些矩阵的特征值和特征向量的基础下面先复习一些矩阵的特征值和特征向量的基础知识

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

3、 ()=0为实系数为实系数n次代数次代数方程，其复根是共轭成对出现方程，其复根是共轭成对出现.)2 . 1(0)( xAI 的的非零解非零解x称为矩阵称为矩阵A的对应于的对应于的的特征向量特征向量.上页上页下页下页 定理定理1 设设为为ARnn的特征值的特征值, 且且Ax=x (x 0)，则有则有 - -p为为A- -pI的特征值，即的特征值，即(A- -pI)x=(- -p)x ； c为的为的cA特征值特征值(c0为常数为常数)； 下面叙述有关特征值的一些下面叙述有关特征值的一些结论结论： k为为Ak的特征值，即的特征值，即Akx=kx ； 设设A为非奇异矩阵，那么为非奇异矩阵，那么0 ,

4、且且- -1为为A- -1的特的特征值，即征值，即A- -1x=- -1x .上页上页下页下页 定理定理2 设设i(i=1,2,n)为为n阶矩阵阶矩阵A=(aij)的特征值，的特征值，则有则有)(Atraniiinii 11 称为称为A的的迹迹； .nA21 定理定理3 设设ARnn，则有，则有. )()(AAT 上页上页下页下页 定义定义2 设设n阶矩阵阶矩阵A=(aij)，令，令 下面讨论下面讨论矩阵特征值界的估计矩阵特征值界的估计.).,(niarnijjiji211 ； 集合集合 称称为复平面上以为复平面上以aii为圆心，以为圆心，以ri为半径的为半径的n阶矩阵阶矩阵A的的n个个Ger

5、schgorin圆盘圆盘. ),(,|niCzrazzDiiii21 上页上页下页下页 定理定理4 (Gerschgorin圆盘定理圆盘定理) 特别地，如果特别地，如果A的一个圆盘的一个圆盘Di是与其它圆盘分离是与其它圆盘分离(即即孤立圆盘孤立圆盘)，则，则Di中精确地包含中精确地包含A的一个特征值的一个特征值.).,(niaranijjijiii211 设设n阶矩阵阶矩阵A(aij)，则，则A的每一个特征值必属的每一个特征值必属于下面某个圆盘之中于下面某个圆盘之中 如果如果A有有m个个圆盘圆盘组成一个连通的并集组成一个连通的并集S，且，且S与余下与余下n- -m个个圆盘圆盘是分离的，则是分离

6、的，则S内恰包含内恰包含A的的m个个特征值特征值.或者说或者说 A的特征值都在的特征值都在n个圆盘的个圆盘的并集并集中中.上页上页下页下页 证明证明 只就给出证明只就给出证明. 设设为为A的特征值，即的特征值，即.1knkjjkjkkraa Ax=x，其中，其中x=(x1,x2, xn)T 0.或或 记记 ，考虑，考虑Ax=x的第的第k个个方程，即方程，即0max1 xxxinik,1knjjkjxxa ,)( kjjkjkkkxaxa 于是于是即即, kjkjkkjjkjkkkaxxaxa 上页上页下页下页这说明，这说明，A的每一个特征值必位于的每一个特征值必位于A的一个圆盘的一个圆盘中，并

7、且相应的特征值中，并且相应的特征值一定位于第一定位于第k个圆盘中个圆盘中(其中其中k是对应特征向量是对应特征向量x绝对值最大的分量的下标绝对值最大的分量的下标).上页上页下页下页.411101014 A 例例2 估计矩阵估计矩阵A的特征值范围，其中的特征值范围，其中 解解 矩阵矩阵A的的3个圆盘为个圆盘为. 24:, 2:, 14:321 DDD 由定理由定理8，可知，可知A的的3个特征值位于个特征值位于3个圆盘的并个圆盘的并集中，由于集中，由于D1是孤立圆盘，所以是孤立圆盘，所以D1内恰好包含内恰好包含A的的一个特征值一个特征值1(为实特征值为实特征值)，即，即. 531 A的其它两个特征值

8、的其它两个特征值2, 3包含在包含在D2, D3的并集中的并集中.上页上页下页下页 定义定义3 设设ARnn为对称矩阵，对于任一非零为对称矩阵，对于任一非零向量向量x，称，称,),(),()(xxxAxxR 为对应于向量为对应于向量x的的瑞利瑞利(Rayleigh)商商. 定理定理5 设设ARnn为对称矩阵为对称矩阵(其特征值次序记其特征值次序记为为12n)，则，则 1. (对任何非零对任何非零xRn)；1 ),(),(xxxAxn 2. ；),(),(maxxxxAxxRxn01 3. .),(),(minxxxAxxRxnn0 上页上页下页下页 证明证明 只证只证1，关于，关于2, 3自己

9、作练习自己作练习. 由于由于A为实对称矩阵，可将为实对称矩阵，可将 1,2,n 对应的特对应的特征向量征向量 x1,x2,xn 正交规范化，则有正交规范化，则有 (xi, xj)=ij，设，设x 0为为Rn中任一向量，则有中任一向量，则有. 0,1221 niiniiixxx 于是于是.),(),(11212 niiniiinxxxAx从而从而1成立成立. 结论结论1说明说明瑞利商瑞利商必位于必位于n和和1之间之间.上页上页下页下页 关于计算矩阵关于计算矩阵A的特征值问题，当的特征值问题，当n2,3时，我时，我们还可按行列式展开的办法求们还可按行列式展开的办法求 ()=0的根的根. 但当但当n

10、较较大时，如果按展开行列式的办法，首先求出大时，如果按展开行列式的办法，首先求出 ()的系的系数，再求数，再求 ()的根，工作量就非常大，用这种办法求的根，工作量就非常大，用这种办法求矩阵的特征值是不切实际的，由此需要研究矩阵的特征值是不切实际的，由此需要研究求求A的特的特征值及特征向量的数值解法征值及特征向量的数值解法. 本章将介绍一些计算机上常用的本章将介绍一些计算机上常用的两类方法两类方法，一，一类是类是幂法及反幂法幂法及反幂法（迭代法），另一类是（迭代法），另一类是正交相似正交相似变换的方法变换的方法（变换法）（变换法）.上页上页下页下页 幂法与反幂法都是求幂法与反幂法都是求实矩阵实矩

11、阵的特征值和特征的特征值和特征向量的向量的向量迭代法向量迭代法，所不同的是幂法是计算矩阵，所不同的是幂法是计算矩阵的的主特征值主特征值(矩阵矩阵按模最大的特征值按模最大的特征值称为称为主特征值主特征值，其模就是该矩阵的其模就是该矩阵的谱半径谱半径)和和相应特征向量相应特征向量的一种的一种向量迭代法，而反幂法则是计算非奇异向量迭代法，而反幂法则是计算非奇异(可逆可逆)矩阵矩阵按模最小的特征值按模最小的特征值和和相应特征向量相应特征向量的一种向量迭的一种向量迭代法代法. 下面分别介绍幂法与反幂法下面分别介绍幂法与反幂法.8.2 幂法及反幂法幂法及反幂法上页上页下页下页现讨论求现讨论求1及及x1的方

12、法的方法.), 2 , 1(nixAxiii 设实矩阵设实矩阵A=(aij)有一个有一个完全的特征向量组完全的特征向量组，即，即A有有n个线性无关的特征向量个线性无关的特征向量，设矩阵，设矩阵A的特征值为的特征值为1,2,n, 相应的特征向量为相应的特征向量为x1,x2,xn. 已知已知A的主的主特征值特征值1是实根，且满足条件是实根，且满足条件)1 . 2(|,|21n 8.2.1 幂法（又称幂法（又称乘幂法乘幂法）上页上页下页下页)3 . 2(),0(122110 axaxaxavnn设设 幂法的幂法的基本思想基本思想是是: 任取任取非零非零的初始向量的初始向量v0 , 由矩由矩阵阵A构造

13、一向量序列构造一向量序列vk称为迭代向量，由假设，称为迭代向量，由假设，v0可唯一表示为可唯一表示为)2 . 2(.,.,011021201 vAAvvvAAvvAvvkkk上页上页下页下页)3 . 2(),0(122110 axaxaxavnn设设于是于是).()/(1112111122211101kkniikiiknknnkkkkkxaxaxaxaxaxavAAvv 其中其中.)/(21 niikiikxa 由假设由假设 故故 从而从而), 3 , 2( 1/1nii , 0lim kk )4 . 2(.lim111xavkkk 为为1的特征向量的特征向量.上页上页下页下页所以当所以当k充

14、分大时，有充分大时，有)5 . 2(,111xavkk 即为矩阵即为矩阵A的的对应特征值对应特征值 1 的一个近似特征向量的一个近似特征向量.用用(vk)i 表示表示vk的第的第i个分量，则当个分量，则当k充分大时，有充分大时，有 )7 . 2(.11 ikikvv即为即为A的的主特征值主特征值 1的近似值的近似值.)6 . 2(,111111kkkkvxaAvv 由于由于上页上页下页下页 迭代公式实质上是由矩阵迭代公式实质上是由矩阵A的乘幂的乘幂 Ak与非零向与非零向量量v0相乘来构造向量序列相乘来构造向量序列vk=Akv0，从而计算主特，从而计算主特征值征值1及其对应的特征向量，这就是及其

15、对应的特征向量，这就是幂法幂法的思想的思想. ).(11 kvvikik 的收敛速度由比值的收敛速度由比值,12 r来确定，来确定，r越小收敛越快，但当越小收敛越快，但当r1 1时收敛可能很慢时收敛可能很慢.上页上页下页下页 定理定理6 设设ARnn有有n个线性无关的特征向量，个线性无关的特征向量，主特征值主特征值1满足条件满足条件 |1|2|n|，则对任何非零向量则对任何非零向量v0(a1 0)，幂法的算式成立，幂法的算式成立.又设又设A有有n个线性无关的特征向量，个线性无关的特征向量，1对应的对应的r个线性个线性无关的特征向量为无关的特征向量为x1,x2,xr，则由，则由(2.2)式有式有

16、 如果如果A的主特征值为实的重根的主特征值为实的重根, 即即1=2=r, 且且 |r|r+1|n|，,)/(11110 nriikiiriiikkkxaxavAv 上页上页下页下页).0(lim111 riiiriiikkkxaxav设设 为为A的特征向量，这说明当的特征向量，这说明当A的主特征值是实的重根的主特征值是实的重根时，定理时，定理6的结论还是正确的的结论还是正确的. 应用应用幂法幂法计算计算A的主特征值的主特征值1及其对应的特征向及其对应的特征向量时，如果量时，如果|1|1(或或|1|2|n|，则对任意非零初始，则对任意非零初始向量向量v0=u0(a1 0)，有幂法计算公式为，有幂

17、法计算公式为则有则有 ,)max(lim11xxukk .lim1 kk上页上页下页下页例例1 用幂法计算矩阵用幂法计算矩阵 1634310232A的主特征值与其对应的特征向量的主特征值与其对应的特征向量.解解取取 v0=u0=(0,0,1)T , 则则 TTvuv25. 0 , 1, 5 . 01, 4,1 , 4 , 211111 ), 2 , 1(max1 k/vuvAuvkkkkkkk 上页上页下页下页直到直到k=8 时的计算结果见下表时的计算结果见下表k1 2, 4, 1, 4 0.5, 1, 0.252 4.5, 9, 7.75 90.5, 1, 0.86113 5.7222, 1

18、1.4444, 8.36111.44440.5, 1, 0.73604 5.4621, 10.9223, 8.2306 10.92230.5, 1, 0.75365 5.5075, 11.0142, 8.2576 11.01420.5, 1, 0.74946 5.4987, 10.9974, 8.2494 10.99740.5, 1, 0.75017 5.5002, 11.0005, 8.2501 11.00050.5, 1, 0.75008 5.5000, 11.0000, 8.2500 11.00000.5, 1, 0.7500TkvTku Tx7500. 0,0 . 1,5 . 0,00

19、00.1111 从而从而k 上页上页下页下页8.2.2 幂法的加速方法幂法的加速方法1、原点平移法、原点平移法 由前面讨论知道，应用幂法计算由前面讨论知道，应用幂法计算A的主特征值的的主特征值的收敛速度主要由比值收敛速度主要由比值 r=|2/1|来决定，但当来决定，但当r 接近于接近于1时，收敛可能很慢时，收敛可能很慢. 这时，一个补救办法是采用加速这时，一个补救办法是采用加速收敛的方法收敛的方法.其中其中p为参数，设为参数，设A的特征值为的特征值为 i，则对矩阵，则对矩阵B的特征的特征值为值为 i- -p ，而且，而且A, B的特征向量相同的特征向量相同. 引进矩阵引进矩阵 B=A- -pI

20、 .上页上页下页下页 如果要计算如果要计算A的主特征值的主特征值 1, 只要只要选择合适的数选择合适的数p，使使 1- -p为矩阵为矩阵B=A- -pI 的主特征值，且的主特征值，且 1212max ppini那么，对矩阵那么，对矩阵B=A- -pI应用幂法求其主特征值应用幂法求其主特征值 1- -p, 收收敛速度将会加快敛速度将会加快. 这种通过求这种通过求B=A- -pI的主特征值和特的主特征值和特征向量，而得到征向量，而得到A的主特征值和特征向量的方法叫的主特征值和特征向量的方法叫原原点平移法点平移法. 对于对于A的特征值的某种分布，它是十分有的特征值的某种分布，它是十分有效的效的.上页

21、上页下页下页例例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的过程是困难的的过程是困难的.上页上页下页下页 定理定理 设设ARnn为为对称矩阵对称矩阵，特征值满足，特

22、征值满足对应的特征向量对应的特征向量xi满足满足(xi, xj)=ij (单位正交向量单位正交向量) ，应用幂法公式应用幂法公式(2.9)计算计算A的主特征值的主特征值 1，则规范化，则规范化向量向量uk的的瑞利商瑞利商给出给出 1的较好的近似值为的较好的近似值为,121nn kkkkkkuuuAuuR2121, 由此可见，由此可见，R(uk)比比k更快的收敛于更快的收敛于 1.2、瑞利商加速、瑞利商加速上页上页下页下页 证明证明 由由(2.8)式及式及,)max(,)max(00100uAuAAuvuAuAukkkkkkk )11. 2(.),(),(),(),()(212112211220

23、0001 knjkjjnjkjjkkkkkkkkkaauAuAuAuAuuuAuuR 得得上页上页下页下页 幂法的幂法的瑞利商加速迭代公式瑞利商加速迭代公式可以写为可以写为 kkkkkkkkkkvukuuuvAuv /), 2 , 1(),(),(1111其中其中A为为n阶实对称矩阵阶实对称矩阵.,11kkux 对给定的误差限对给定的误差限 ，当，当| kk- -1| 时，取近似值时，取近似值上页上页下页下页反幂法反幂法 反幂法是用于求非奇异矩阵反幂法是用于求非奇异矩阵A的的按模最小的特征按模最小的特征值和对应特征向量值和对应特征向量的方法的方法. 而结合原点平移法的反幂而结合原点平移法的反幂

24、法则可以求矩阵法则可以求矩阵A的任何一个的任何一个具有先了解的特征值和具有先了解的特征值和对应的特征向量对应的特征向量。设矩阵设矩阵A非奇异非奇异,其特征值其特征值 i (i=1,2,n) ,满足满足0121 nn 其相应的特征向量其相应的特征向量x1,x2,xn线性无关，则线性无关，则 A- -1 的特征的特征值为值为1/ i ，对应的特征向量仍为，对应的特征向量仍为xi (i=1,2,n).iiiiiixxAxAx11 上页上页下页下页此时此时, A- -1的特征值满足的特征值满足11111 nn因此因此, 对对A- -1应用幂法应用幂法,可求出其可求出其主特征值主特征值1/ n k 和和

25、特征向量特征向量 xn uk .从而求得从而求得A的的按模最小按模最小特征值特征值 n 1/k 和对应的和对应的特征向量特征向量 xn uk ，这种求这种求A- -1的方法就称的方法就称为为反幂法反幂法.上页上页下页下页为了避免求为了避免求A- -1, 可通过解线性方程组可通过解线性方程组Avk=uk- -1得到得到vk，采用采用LU分解法，即先对分解法，即先对A进行进行LU分解分解A=LU, 此时此时反幂反幂法的迭代公式法的迭代公式为为 , 2 , 1/max,1 kvuvvzUvzuLzkkkkkkkkkkk 求求出出解解求求出出解解 ), 2 , 1(max11 k/vuvuAvkkkk

26、kkk knknux ,1 反幂法的迭代公式反幂法的迭代公式为为上页上页下页下页对给定的误差对给定的误差 ，当，当|kk- -1| n|0，则对任意非零初始向量则对任意非零初始向量u0(an 0) ，由反幂法计算公，由反幂法计算公式构造的向量序列式构造的向量序列vk，uk满足满足 ,)max(limnnkkxxu .1)max(limnkkv 上页上页下页下页 在反幂法中也可以用在反幂法中也可以用原点平移法原点平移法加速迭代过程加速迭代过程，或或求其它特征值与其对应的特征向量求其它特征值与其对应的特征向量. 如果矩阵如果矩阵(A- -pI)- -1存在，显然其特征值为存在，显然其特征值为,1,

27、1,121pppn 对应的特征向量仍然是对应的特征向量仍然是x1,x2,xn. 如果如果p是是A的特征值的特征值 j的一个近似值，且设的一个近似值，且设 j与其它与其它特征值是分离的，即特征值是分离的，即),(jippij 就是说就是说1/( j- -p)是矩阵是矩阵 (A- -pI)- -1的主特征值，可用反幂的主特征值，可用反幂法应用于矩阵法应用于矩阵(A- -pI)- -1计算计算 j及对应的特征向量及对应的特征向量.上页上页下页下页现对矩阵现对矩阵(A- -pI)- -1应用反幂法，得到如下迭代公式应用反幂法，得到如下迭代公式 00110,(),(1,2,).(2.12)max/ ma

28、x()kkkkkkkuvvApIukvuv 初初始始向向量量上页上页下页下页 定理定理9 设设ARnn有有n个线性无关的特征向量，个线性无关的特征向量， 矩阵矩阵A的特征值及对应的特征向量分别记为的特征值及对应的特征向量分别记为 i 及及xi (i=1,2,n)，而，而p为为 j的近似值，的近似值，(A- -pI)- -1存在，且存在，且 ,)max(limjjkkxxu jkjkkvppv )max(1,1)max(lim即即则对任意非零初始向量则对任意非零初始向量u0(aj 0) ，由反幂法计算公式，由反幂法计算公式(2.12)构造的向量序列构造的向量序列vk，uk满足满足. )( |ji

29、ppij . |min/ |pprijij 且收敛速度为且收敛速度为上页上页下页下页 由该定理知，对由该定理知，对A- -pI(其中其中p j)应用反幂法，可应用反幂法，可用来计算特征向量用来计算特征向量xj，只要选择，只要选择p是是 j的一个较好的近的一个较好的近似且特征值分离情况较好，一般似且特征值分离情况较好，一般r很小，常常只要迭很小，常常只要迭代一二次就可完成特征向量的计算代一二次就可完成特征向量的计算.反幂法迭代公式中的反幂法迭代公式中的vk是通过解方程组是通过解方程组.)(1 kkuvpIA求得的求得的, 为了节省工作量为了节省工作量, 可以先将可以先将A- -pI进行三角分解进

30、行三角分解.)(LUpIAP 于是求于是求vk相对于解两个三角形方程组相对于解两个三角形方程组.,1kkkkyUvPuLy 上页上页下页下页实验表明实验表明, 按下述方法选择按下述方法选择u0是较好的是较好的: 选选u0使使)13. 2()1 , 1 , 1(011 PuLUv用回代求解用回代求解(2.13)即得即得v1,然后再按公式然后再按公式(2.12)进行迭代进行迭代.上页上页下页下页例例5 求矩阵求矩阵A最接近于最接近于1.9的特征值和相应的特征向量的特征值和相应的特征向量 3 1- 2-1- 4 3 2- 3 7 A取取Ty) 1 , 1 , 1 (0作迭代，结果如表：作迭代，结果如

31、表：上页上页下页下页上页上页下页下页1222,1,(1) , diag(,) (1,2, ),2(),(), Tniijn nnTijn nijiji ji jAUU AUDinAUAaUBU AUbab 任任意意实实对对称称矩矩阵阵 可可通通过过正正交交相相似似变变换换化化成成对对角角阵阵 即即存存在在 正正交交矩矩阵阵使使得得其其中中是是 的的特特征征值值中中各各列列即即为为相相应应的的特特征征向向量量。（ ）在在正正交交相相似似变变换换下下，矩矩阵阵元元素素的的平平方方和和不不变变。设设为为正正交交矩矩阵阵，记记则则理理论论基基础础：1n Jacobi方方法法用用来来求求实实对对称称矩矩

32、阵阵的的全全部部特特征征值值及及相相应应特特征征向向量量。8.3 Jacobi方法方法上页上页下页下页,A通通过过一一次次正正交交变变换换 将将 中中一一对对非非零零的的非非对对角角元元素素化化成成零零 并并且且使使得得非非对对角角元元素素的的平平方方和和减减少少。反反复复进进行行上上述述过过程程，使使变变换换后后的的矩矩阵阵的的非非对对角角元元素素的的平平方方和和趋趋于于零零，从从而而使使该该矩矩阵阵近近似似为为对对角角矩矩阵阵，得得到到全全部部特特征征值值和和特特征征向向量量。Jacobi 方法的基本思路方法的基本思路上页上页下页下页矩阵的旋转变换矩阵的旋转变换cos ,1sin ,pqp

33、pqqpqqppqUuuunuU 其其中中，的的主主对对角角线线元元素素中中其其余余为为 ；而而其其非非主主对对角角线线元元素素中中维维空空间间中中的的二二维维坐坐其其余余为为标标旋旋0 0。称称为为转转矩矩阵阵。( )pqU 坐坐标标旋旋转转矩矩阵阵是是正正交交矩矩阵阵. .11cossin1( )1sincos11pqU 上页上页下页下页(1)(1)(1)22(1)22(1)(1)0, ()cossinsin2 sincossin2 cossin (, ) TpqqppqpqijppppqqpqqqppqqpqpiippiqiAaaAU AUaaaaaaaaaaaaaip qa 设设 为为

34、实实对对称称矩矩阵阵，且且若若记记则则有有(1)(1)(1)(1)(1)(1)sincos ( , )1()sin2cos2 2()cot2 2qiiqpiqiijjiijpqqpqqpppqppqqpqaaaaaai jp qaaaaaaaa 如如果果取取 使使得得(1)(1)(/4)0,.pqqppqqpaaAaa 则则有有得得到到一一个个使使 中中非非零零的的非非对对角角元元素素变变成成零零的的正正交交相相似似变变换换上页上页下页下页 (1)(2)( )(1)2(1)(1)2(1)(1)(1)(1)(1)(1), ( ),()()( , ),cossin ,kijijijijijjiij

35、piippiqiqiiqAAAAAE AaE Aaaaai jp q aaaaaa 对对重重复复上上述述过过程程得得到到一一个个矩矩阵阵序序列列。可可证证，虽虽然然这这种种变变换换不不一一定定能能使使矩矩阵阵中中非非对对角角元元中中零零元元素素的的个个数数单单调调增增加加，但但可可保保证证非非对对角角元元的的平平方方和和递递减减。以以 与与为为例例：设设，则则由由(1)(1)(1)(1)2(1)2(1)2(1)2(1)2,2222,sincos ,(, )=0()()()2()() +2() 2()( )2( ) piqipqqpijijpiqipqijijip qi jp qijpiqipq

36、ijip qi jp qaaip qaaE AaaaaaaaaE AaE A ，上上式式表表明明，在在上上述述旋旋转转变变换换下下，非非对对角角元元的的平平方方和和严严格格单单调调递递减减，因因而而对对角角元元的的平平方方和和单单调调递递增增。上页上页下页下页 ( )(0)2( )( )Jacobi, lim ()=lim0,Jacobikkkijkki jAnAAAAE Aa 设设 为为 阶阶实实对对称称阵阵，对对 用用法法得得到到序序列列其其中中则则即即法法收收敛敛。Jacobi 迭代法收敛定理迭代法收敛定理 上页上页下页下页说说明明o1 Jacobi;定定理理表表明明，法法是是收收敛敛的

37、的Jacobi法法是是求求中中小小型型稠稠密密实实对对称称矩矩阵阵的的全全部部特特征征值值与与特特征征向向量量的的较较好好方方法法。o2 JacobiAnAn当当 的的阶阶数数 不不太太高高时时，算算法法的的收收敛敛速速度度很很快快；但但当当 的的 阶阶数数 变变得得较较大大时时，其其收收敛敛速速度度将将会会变变慢慢，即即法法 为为适适合合计计算算中中等等规规模模的的实实对对称称矩矩阵阵的的特特征征值值问问题题；o3 对对中中等等规规模模问问题题，具具有有较较好好的的数数值值稳稳定定性性；求求得得的的结结果果 的的精精度度也也很很高高，得得到到的的特特征征向向量量正正交交性性很很好好；o4 不

38、不足足之之处处：运运算算量量大大，不不能能保保持持矩矩阵阵的的特特殊殊形形状状 （如如稀稀疏疏性性）。上页上页下页下页 00(0)12T(1)(0)(0)(0)0:1,2,cot20,4cos0.7071068,sin0.70710680.70710680.70710680 ( )0.70710680.70710680 ,001100.7071068030.70710680.70710680.70710682kpqUUAUA U 210121012AA 用用J Ja ac co ob bi i方方法法求求 的的特特征征值值。解：解： 例例6： 上页上页下页下页11(1)230.70710678

39、0.47765830.88807380.45970081 0 00 0.888071:2,3,cot2,cos,sin 38 -0.45970080 0.4597008 ( ) 0.88807kpqUU T(2)(1)(1)(1)38 3.0000000 0.3250576 0.6279630 0.3250576 0.6339746 -0.0000000 0.6279630 0 2.3,U660254UAA 上页上页下页下页22(2)230.50478660.55166340.85165390.52410460.8516539 0 -0.5241046 0 1 2:1,3,cot2,cos,s

40、in 00.5241046 0 0.8516( )53 9kpqUU T(3)(2)(2)(2) 3.3864461 0.2768366 -0.0000000 0.2768366 0.6339746 -0.1703642-0.0000000 -0.1703642 1.9795,UU793AA 上页上页下页下页66(6)236:1,2,cot2-693.88568,-0.00072057929,cos0.9999997,sin0.0007205792 0.9999997 0.0007205792 0 ( ) -0.0007205792 0.9999997 0 0 kpqUU T(7)(6)(6)

41、(6), 0 1 3.4142136 0 0.0000000UU -0.0000000 0.5857864 -0.0000242 0.0000000 -0.0000242 1.9999999(AAE (7)1.1668492e-009A 上页上页下页下页123123(0)(1)(6)3.4142136,0.5857864,1.9999999;223.4142136,220.5857864,20.0000001 0.5000000 0.5000121 -0.7070982-0.7071068 0.7071068 0.0000121 0.5000000 QUUU所所以以，准准确确特特征征值值最最大

42、大误误差差不不超超过过；又又 0.4999879 0.7071153 0.50000000.5000121 -0.7070982-0.70710680.7071068 0.0000121 0.50000000.4999879 0.7071153 所所以以，对对应应的的特特征征向向量量分分别别为为，上页上页下页下页8.4 QR算法算法8.4.1. 化矩阵为化矩阵为Hessenberg形形8.4.2 QR 算法及其收敛性算法及其收敛性（1）镜面反射矩阵）镜面反射矩阵 （2）化一般实矩阵为上）化一般实矩阵为上Hessenberg矩阵矩阵 （3）化对称矩阵为三对角矩阵）化对称矩阵为三对角矩阵上页上页下

43、页下页 当当A为实矩阵，如果限制用正交相似变换，由于为实矩阵，如果限制用正交相似变换，由于A有复的特征值，有复的特征值， A不能用正交相似变换约化为上三不能用正交相似变换约化为上三角阵角阵. 用正交相似变换能约化到什么程度呢？用正交相似变换能约化到什么程度呢？(Schur定理定理) 设设ARnn，则存在酉矩阵，则存在酉矩阵U使使其中其中rii(i=1,2,n)为为 A的特征值的特征值. 下面给出理论上有关通过酉相似变换及正交变下面给出理论上有关通过酉相似变换及正交变换可以约化一般矩阵换可以约化一般矩阵A到什么程度的问题到什么程度的问题.),(22211211上三角阵上三角阵RrrrrrrAUU

44、nnnnH 上页上页下页下页其中其中Rii(i=1,2,m)为一阶或二阶方阵，且每个一阶为一阶或二阶方阵，且每个一阶Rii是是A的实特征值，每个二阶对角块的实特征值，每个二阶对角块Rii的两个特征值的两个特征值是是 A的两个共轭复特征值的两个共轭复特征值. (实实Schur分解分解) 设设ARnn，则存在正交矩阵，则存在正交矩阵Q使使,22211211 mmmmTRRRRRRAQQ上页上页下页下页63（1）镜面反射矩阵）镜面反射矩阵。或称为或定义：矩阵r r变变换换H Ho ou us se eh ho ol ld de e镜镜面面反反射射变变换换,RuuuuuEHuuuuEHnTTT),2(

45、222。，使，变换存在、为对称正交变换；即、：TnTexeHxH,xRxHHHH)0 , 0 , 1 (, rHouseholde0,2 ,1 1211性性质质8.4.1. 化矩阵为化矩阵为Hessenberg形形镜面反射矩阵的意义是镜面反射矩阵的意义是“成批成批”消去向量的非零元素消去向量的非零元素。上页上页下页下页.,11221yxeyxT满足变成求变换矩阵将, 1121623. 5,1623. 3)(, 1121121TTexuxxsignx9387.00613.01225.03162.00613.09387.01225.03162.01225.01225.0755.06325.0316

46、2.03162.06325.06325.02uuuuEHTT例：例：解解:THxy0001623. 3#上页上页下页下页 (1)用用初等反射矩阵作正交相似变换初等反射矩阵作正交相似变换约化一般约化一般实矩阵实矩阵A为为上上Hessenberg矩阵矩阵.化矩阵为上化矩阵为上Hessenberg矩阵矩阵讨论讨论两个问题两个问题 (2)用用初等反射矩阵作正交相似变换初等反射矩阵作正交相似变换约化对称约化对称矩阵矩阵A为为对称三对角矩阵对称三对角矩阵. 于是，于是，求原矩阵特征值问题求原矩阵特征值问题，就，就转化为转化为求上求上Hessenberg矩阵矩阵或或对称三对角矩阵的特征值对称三对角矩阵的特征

47、值问题问题.上页上页下页下页（2）化一般实矩阵为上）化一般实矩阵为上Hessenberg矩阵矩阵 设设ARnn，下面来说明，可选择初等反射矩，下面来说明，可选择初等反射矩阵阵U1,U2,Un- -2使使A经正交相似变换约化为一个上经正交相似变换约化为一个上Hessenberg阵阵.(1) 设设,)1(221)1(1211212222111211 AcAaaaaaaaaaaAnnnnnn上页上页下页下页21/212112111 111121sgn()(),(3.1)().niiaaucea 其中其中c1=(a21,an1)TRn- -1 ,不妨设不妨设c10，否则这一步不，否则这一步不需要约化需

48、要约化. 于是于是, 可选择初等反射阵可选择初等反射阵11111TRIu u 1 11 1R ce 令令,111 RU使使 其中其中上页上页下页下页则则 1)1(221111)1(12111112RARcRRAaUAUA(2)(2)(2)1112131(2)(2)(2)122232(2)(2)1112(2)(2)(2)32333(2)122(2)(2)(2)23,000nnnnnnnaaaaaaaAAaaacAaaa 其中其中.,),()2()2()2(222)2(2)2(322 nnnTnRARaac上页上页下页下页(2) 第第k步约化：重复上述过程，设对步约化：重复上述过程，设对A已完成第

49、已完成第1步步,第第k- -1步正交相似变换，即有步正交相似变换，即有,111 kkkkUAUA或或,11111 kkkUUAUUA且且 )()(1,)()(, 1)(1, 1)(, 1)()(1,)(1)(2)(1,2)(2)1(1,2)2(221)(1)(1, 1)(1)1(1, 1)2(12)1(11knnkknknkknkkkkkkkkknkkkkkkkknkkkkkkknkkkkkkkaaaaaaaaaaaaaaaaaaaaA 上页上页下页下页,0)(22)(12)(11knkAcAAknkkkkk 其中其中 为为k阶上阶上Hessenberg阵，阵，设设ck0, 于是可选择初等反射

50、阵于是可选择初等反射阵Rk使使 其中，其中，Rk计算公式为计算公式为1,kkkR ce )(11)()(, 1,),(kknTknkkkkkARaac .)()()(22knknkRA 上页上页下页下页( )( )2 1/21,11( )1,1sgn()() ),(3.2)(),.nkkkkkiki kkkkkkkkkkTkkkkaauceaRIu u 令令, kkRIU则则)3 . 3(00)1(221)1(12)1(11)(22)(12)1(111 kkkkkkkkkkkkkkkkAcAARARcRRAAUAUA上页上页下页下页221122 nnUUAUUUU其中其中 为为k+1阶上阶上H

51、essenberg阵，第阵，第k步约化只需步约化只需计算计算 及及 (当当A为对称矩阵时，只需要为对称矩阵时，只需要计算计算 ).)1(11 kAkkRA)(12kkkRAR)(22kkkRAR)(22(3) 重复上述过程，则有重复上述过程，则有.1)(1)1(1, 12)3(332)2(22111 nnnnnnnnnAaaaaa 上页上页下页下页 定理定理11 (Household约化矩阵为上约化矩阵为上Hessenberg阵阵) 设设ARnn则存在初等反射矩阵则存在初等反射矩阵U1,U2,Un- -2 使使.00221122HAUUUUAUUUUTnn 为为上上Hessenberg矩阵矩阵

52、.总结上述结论，有总结上述结论，有上页上页下页下页例例7 用用Household方法方法将矩阵将矩阵 7242327341AA约化为约化为上上Hessenberg阵阵. 解解 选取初等反射阵选取初等反射阵R1使使 , 其中其中c1=(2,4)T.1111ecR (1) 计算计算R1:)() 1 , 5 . 0(, 4) 4 , 2max(11规规范范化化Tcc .,472136. 4,809017. 1)5 . 0(,)1,618034. 1(,118034. 125. 11111111111TTuuIRecu 上页上页下页下页.1111ecR 则有则有(2) 约化计算约化计算:,111 RU

53、则得到则得到上上Hessenberg阵阵为为.200000. 2399999. 00400000. 0799999. 7472136. 4447214. 0602631. 74112HAUUA 上页上页下页下页（3）化对称矩阵为对称三对角矩阵）化对称矩阵为对称三对角矩阵推论推论 (Household约化对称矩阵为对称三对角矩阵约化对称矩阵为对称三对角矩阵) 设设ARnn为对称矩阵，则存在初等反射矩阵为对称矩阵，则存在初等反射矩阵U1,U2,Un- -2使使.111222111221122CcbbcbbcbbcUUAUUUUnnnnnn 为为对称三对角矩阵对称三对角矩阵.上页上页下页下页8.4.

54、2 QR 算法及其收敛性算法及其收敛性QR算法算法是一种变换方法，是计算一般矩阵是一种变换方法，是计算一般矩阵(中小中小型矩阵型矩阵)全部特征值全部特征值问题的问题的最有效方法之一最有效方法之一. (1) 上上Hessenberg矩阵矩阵的的全部特征值全部特征值问题；问题； (2) 计算计算对称三对角矩阵对称三对角矩阵的的全部特征值全部特征值问题，问题， 目前目前QR算法算法主要主要用来计算：用来计算：且且QR算法算法具有收敛快，算法稳定等特点具有收敛快，算法稳定等特点.上页上页下页下页 对于一般矩阵对于一般矩阵ARnn (或对称矩阵或对称矩阵)，则首先用，则首先用Household方法将方法

55、将A化为上化为上Hessenberg阵阵B(或对称三或对称三对角阵对角阵)，然后再用，然后再用QR算法算法计算计算B的全部特征值的全部特征值. 设设ARnn，且，且A对进行对进行QR分解，即分解，即A=QR，其中其中R为上三角阵为上三角阵, Q为正交阵为正交阵, 上页上页下页下页于是可得到一新矩阵于是可得到一新矩阵B=RQ=QTAQ.显然，显然，B是由是由A经过正交相似变换得到，因此经过正交相似变换得到，因此B与与A的的特征值相同特征值相同. 再对再对B进行进行QR分解，又可得一新矩阵分解，又可得一新矩阵,重复这一过程可得到矩阵序列：重复这一过程可得到矩阵序列： 设设A=A1 将将A1进行进行

56、QR分解分解A1=Q1R1 作矩阵作矩阵A2=R1Q1=Q1TR1Q1 QR算法，就是利用矩阵的算法，就是利用矩阵的QR分解，按上述递分解，按上述递推法则构造矩阵序列推法则构造矩阵序列Ak的过程的过程. 只要只要A为非奇异矩为非奇异矩阵，则由阵，则由QR算法就完全确定算法就完全确定Ak.上页上页下页下页 定理定理13 (基本基本QR算法算法) 设设A=A1Rnn, 构造构造QR算法算法:)1 . 4(), 2 , 1(),(1 kQRARIQQRQAkkkkkTkkkk为为上上三三角角阵阵其其中中则则有有记记,1221RRRRQQQQkkkk ;,111kTkkkkAQQAAA 即即相相似似于

57、于）（;)()()2(1211211kTkkTkkQAQQQQAQQQA .)3(kkkkRQAA 的的分分解解式式为为上页上页下页下页.)()(11111111111111111kkkkkkTkkkkkkkkkkkkkkAAARQARQAQQRAQRRAQQRRRQQQRQ 于是于是 证明证明 (1),(2)显然，现证显然，现证(3). 用归纳法，显然，当用归纳法，显然，当k=1时有时有 , 设设Ak- -1有分解式有分解式1111RQRQA ,111 kkkRQA上页上页下页下页 定理定理14 (QR算法的收敛性算法的收敛性) 设设A=(aij)Rnn, (1) 如果如果A的特征值满足的特征值满足:; 021 n (2) A有标准型有标准型A=XDX- -1, 其中其中D=diag( 1, 2, n),且设且设X- -1有三角分解有三角分解X- -1=LU(L为单位下三