2022年第章矩阵特征值的数值解法 .pdf

《2022年第章矩阵特征值的数值解法 .pdf》由会员分享，可在线阅读，更多相关《2022年第章矩阵特征值的数值解法 .pdf（40页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1 第 9 章矩阵特征值的数值解法9.1 引言矩阵特征值问题有广泛的应用背景. 例如动力系统和结构系统中的振动问题、电力系统的静态稳定分析上、工程设计中的某些临界值的确定等，都归结为矩阵特征值问题. 数学中诸如方阵的对角化及解微分方程组等问题，都要用到特征值的理论. 本章介绍n阶实矩阵n nRA的特征值与特征向量的数值解法. 定义 9.1.1已知n阶实矩阵()n nijaRA，如果存在常数和非零向量x，使Axx或()0AI x(9.1.1) 那么称为A的特征值 (eigenvalue)，x为A的相应于的特征向量 (eigenvector). 多项式111212122212( )det()nnn

2、nnnnaaaaaapaaaAI(9.1.2) 称为 特征多项式 (characteristic polynomial) ，det()0AI(9.1.3) 称为 特征方程 (characteristic equation). 注式(9.1.3)是以为未知量的一元n次代数方程，( )det()npAI是的n次多项式 . 显然，A的特征值就是特征方程(9.1.3)的根. 特征方程 (9.1.3)在复数范围内恒有解，其个数为方程的次数(重根按重数计算)，因此n阶矩阵A在复数范围内有n个特征值 . 除特殊情况(如2,3n或A为上 (下)三角矩阵 )外，一般不通过直接求解特征方程(9.1.3)来求A的特

3、征值 , 原因是这样的算法往往不稳定. 在计算上常用的方法是幂法与反幂法和相似变换方法 . 本章只介绍求矩阵特征值与特征向量的这两种基本方法. 为此将一些特征值和特征向量的性质列在此处. 定理 9.1.2设n阶方阵()ijn naA的特征值为12,n，那么(1) 121122nnnaaa; (2) 12detnA. 定理 9.1.3如果是方阵A的特征值，那么(1) k是kA的特征值，其中k是正整数 ; (2) 当A是非奇异阵时，1是1A的特征值 . (3) ( )np是()npA的特征值，其中( )npx是多项式2012( )nnnpxaa xa xa x. 定义9.1.4设,A B都是n阶方

4、阵 . 若有n阶非奇异阵P，使得1PAPB，则称矩阵A与B相似 (similar)，1PAP称为对A进行 相似变换 (similarity transformation) ，P称名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 40 页 - - - - - - - - - 2 为相似变换矩阵 (similarity transformation matrix). 定理 9.1.5若矩阵A与B相似，则A与B的特征值相同 . 定理 9.1.6如果A是n阶正交矩阵，那么(1) 1

5、TAA，且det1A或1; (2) 若yAx，则22yx, 即TTxxyy. 定理 9.1.7设A是任意n阶实对称矩阵，则(1) A的特征值都是实数；(2) A有n个线性无关的特征向量. 定理 9.1.8设A是任意n阶实对称矩阵，则必存在n阶正交矩阵P，使得1TPAPP AP，其中12diag(,)n是以A的n个特征值12,n为对角元素的对角矩阵. 定理9.1.9 (圆盘定理 ) 矩阵()ijn naA的任意一个特征值至少位于复平面上的几个圆盘1,1,2,niiiijjjiDz zaain，中的一个圆盘上。9.2 幂法与反幂法9.2.1幂法及其加速9.2.1.1 幂法幂法是计算矩阵按模最大特征

6、值(largest eigenvalue in magnitude) 及相应特征向量的迭代法. 该方法稍加修改，也可用来确定其他特征值. 幂法的一个很有用的特性是：它不仅可以求特征值， 而且可以求相应的特征向量. 实际上，幂法经常用来求通过其他方法确定的特征值的特征向量 . 下面探讨幂法的具体过程. 设矩阵n nRA的n个特征值满足1230n, (9.2.1) 且有相应的n个线性无关的特征向量12,nxxx，则12,nxxx构成n维向量空间nR的一组基 , 因此1,1,2,nniiiiinRRzzx. 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - -

7、 - - - - - 名师精心整理 - - - - - - - 第 2 页，共 40 页 - - - - - - - - - 3 在nR中选取某个满足10的非零向量01niiizx. 用矩阵A同时左乘上式两边，得011nniiiiiiiAzAxx. 再用矩阵A左乘上式两边，得2201niiiiA zx. 这样继续下去，一般地有01,1,2,.nkkiiiikA zx(9.2.2) 记10,1,2,kkkkzAzA z，则由式 (9.2.2)得0111121, 1,2,.knnkkkikiiiiiiikzA zxxx(9.2.3) 由假设 (9.2.1)，结合式 (9.2.3)，得111lim,

8、kkkzx(9.2.4) 于是对充分大的k有111.kkzx(9.2.5) 式(9.2.4)表明随着k的增大，序列1kkz越来越接近A的对应于特征值1的特征向量1x的1倍, 由此可确定对应于1的特征向量1x. 当k充分大时，可得1x的近似值 . 上述收敛速度取决于比值21. 事实上，由式 (9.2.3)知，321122331111kkkknnnkzxxxx. (9.2.6) 再由式 (9.2.1)得321111n. (9.2.7)结合式 (9.2.6)和式(9.2.7)知，序列1kkz收敛速度取决于比值21. 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - -

9、 - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 40 页 - - - - - - - - - 4 下面计算1. 由式 (9.2.3)知1111011121,knkkikkiiiAAzzzxx当 k 充分大时 , 11111kkzx. 结合式 (9.2.5)，得11kkzz. 这表明两个相邻向量大体上只差一个常数倍，这个倍数就是A的按模最大特征值1. 记T(1)(2)( ),nkkkkzzzz, 则有()11( )lim,1,2,jkjkkzjnz，(9.2.8) 即两个相邻的迭代向量所有对应分量的比值收敛到1. 定义 9.2.1上述由已知非零向量0z及矩

10、阵A的乘幂kA构造向量序列kz来计算A的按模最大特征值1及相应特征向量的方法称为幂法 (power method)，其收敛速度由比值21来确定，越小，收敛越快. 注 由幂法的迭代过程(9.2.3)容易看出，如果11（或11） ，那么迭代向量kz的各个非零的分量将随着k而趋于无穷 （或趋于零） ，这样在计算机上实现时就可能上溢(或下溢 ). 为了克服这个缺点，需将每步迭代向量进行规范化: 记T(1)(2)()1,nkkkkkyyyyAz. 若存在ky的某个分量0()jky，满足0()( )1maxjjkkj nyy，则记0()max()jkkyy. 将ky规范化max()kkkzyy, 这样就把

11、kz的分量全部控制在 1,1中. 例如，设T( 2, 3, 0,5, 1)ky，因为ky的所有分量中，绝对值最大的的是5，所以max()5ky，故Tmax()(0.4,0.6, 0, 1,0.2)kkkzyy. 综上所述，得到下列算法：算法 9.2.1 (幂法 ) 设A是n阶实矩阵，取初始向量0nRz，通常取T0(1,1,1)z，其迭代过程是：对1,2,k，有1,max(),.kkkkkkkmmyAzyzy(9.2.9) 定理 9.2.1对式 (9.2.9)中的序列kz和km有名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师

12、精心整理 - - - - - - - 第 4 页，共 40 页 - - - - - - - - - 5 11limmaxkkxzx, 1limkkm，(9.2.10) 其收敛速度由21确定. 证明由迭代过程 (9.2.9)知211111kkkkkkkkkkkmmm mm mzyA zA yA z0001max()kkkkiimA zA zA z, (9.2.11) 其中01niiizx. 若10，则由 (9.2.3)知：011121knkkiiiiA zxx, 代入式(9.2.11)得11211121maxkniiiikkniiiixxzxx，(9.2.12) 故11limmaxkkxzx.

13、(9.2.13) 而212011120maxmaxmaxkkkkkkkkAAAyA zA zyAzyzz111211111121m a xknkiiiiknkiiiiaaaaxxxx, (9.2.14) 于是1121111121maxmax()maxkniiiikkkniiiiaamaaxxyxx, (9.2.15) 故1limkkm. (9.2.16) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 40 页 - - - - - - - - - 6 由式 (9.2.6)

14、和式 (9.2.7)知：上述收敛速度由21确定 . 证毕！例 9.2.1用幂法求方阵123213335A按模最大特征值及相应的特征向量，要求3110kkmm. 解选取初始向量T0(1,1,1)z，按式 (9.2.9)迭代，结果见表9.2.1. 表 9.2.1 kkzkm1kkmm1 T(0.545455, 0.545455, 1)112 T(0.560441, 0.560441, 1)8.27272.72733 T(0.559787, 0.559787, 1)8.362729104 T(0.559818, 0.559818, 1)8.358734105 T(0.559817, 0.559817

15、, 1)8.358930.2 10因此，所求按模最大特征值18.3589，相应特征向量T1(0.559817, 0.559817, 1)x. 事 实 上 ，A的 按 模 最 大 特 征 值14198.3588989， 相 应 特 征 向 量T1(0.55981649, 0.55981649, 1)x, 故所得结果具有较高的精度. 9.2.1.2 幂法的加速从上面的讨论可知，由幂法求按模最大特征值，可归结为求数列km的极限值，其收敛速度由21确定 . 当21接近1时, 收敛速度相当缓慢. 为了提高收敛速度, 可以利用外推法进行加速. 因为序列km的收敛速度由21确定， 所以若km收敛， 当k充分

16、大时， 则有121kkmO，或改写为121kkmc, 其中c是与k无关的常数 . 由此可得11211kkmm, (9.2.17) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 40 页 - - - - - - - - - 7 这表明幂法是线性收敛的. 由式(9.2.17)知1121111kkkkmmmm. 由上式解出1，并记为2km，即2221122121()22kkkkkkkkkkkkkmmmmmmmmmmmmm，(9.2.18) 这就是计算按模最大特征值的加速公式.

17、 将上面的分析归结为如下算法：算法9.2.2 (幂法的加速) 设A是n阶实矩阵，给定非零初始向量0nRz，通常取T0(1,1,1)z. 对1,2k，用迭代式1,max(),kkkkkkkmmyAzyzy求出12,mm及12,zz. 再对3,4,k，迭代过程为1212212,max(),(),2.kkkkkkkkkkkkkkmmmmmmmmmyAzyzy(9.2.19) 当1kkmm(0是预先给定的精度) 时，迭代结束，并计算kz；否则继续迭代，直至满足迭代停止条件1kkmm. 有关加速收敛技术，读者请参考文献11. 8.2.2 反幂法及其加速反幂法是计算矩阵按模最小特征值及相应特征向量的迭代法

18、, 其基本思想是：设矩阵n nRA非奇异， 用其逆矩阵1A代替A, 矩阵A的按模最小特征值就是矩阵1A的按模最大特征值 . 这样用1A代替A做幂法，即可求出1A的按模最大特征值， 也就是矩阵A的按模最小特征值；这种方法称为反幂法 (inverse power method). 因为矩阵A非奇异，所以由iiiAxx可知：11iiiAxx. 这说明：如果A的特征值满足1210n-n, (9.2.20) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 40 页 - - - -

19、- - - - - 8 那么1A的特征值满足1211111nn-, (9.2.21) 且对A的应于特征值i的特征向量ix也是1A的对应于特征值1i的特征向量 . 由上述分析知：对1A应用幂法求按模最大的特征值1n及相应的特征向量nx，就是求A的按模最小的特征值n及相应的特征向量nx. 算法 9.2.3 (反幂法 ) 任取初始非零向量0nRz，通常取T0(1,1,1)z. 为了避免求1A，对1,2,k，将迭代过程 (9.2.9)改写为 : 1,max(),.kkkkkkkmmAyzyzy(9.2.22) 仿定理 9.2.1 的证明，可得：定理 9.2.2对式 (9.2.22)中的序列kz和km有

20、1lim,limmax()nkkkknnmxzx, (9.2.23) 其收敛速度由1nn确定. 注 按(9.2.22)进行计算，每次迭代都需要解一个方程组1kkAyz. 若利用三角分解法求解方程组，即ALU，其中L是下三角矩阵，U是上三角矩阵，这样每次迭代只需解两个三角方程组11,.kkLzUyv9.2.3 原点平移法为了提高收敛速度，下面介绍加速收敛的原点平移法. 设矩阵pBAI，其中p是一个待定的常数，A与B除主对角线上的元素外，其他元素相同. 设A的特征值为12,n，则B的特征值为12,nppp,且A与B的特征向量相同. 9.2.3.1 原点平移法的幂法设1是A的按模最大的特征值，选择p

21、，使12,3,4,ipppin，及名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 40 页 - - - - - - - - - 9 2211pp. (9.2.24) 对B应用幂法，得算法 9.2.4 对1,2,k，有1(),max(),kkkkkkkpmmyAI zyzy(9.2.25) 且111lim,lim,maxkkkkmpxzx(9.2.26) 其收敛速度由21() ()pp确定. 由式 (9.2.24)知：在计算B的按模最大特征值1p的过程 (9.2.25)中，

22、收敛速度得到加速；算法9.2.4 又称为 原点平移下的幂法(power method with shift) . 9.2.3.2 原点平移下的反幂法设n是A的按模最小的特征值，选择p，使1,1,2,2nnipppin. (9.2.27) 及11nnnnpp. (9.2.28) 若矩阵pBAI可逆，则1B的特征值为12111,nppp且有1111,1,2,2nniinppp. (9.2.29) 对B应用反幂法，得：算法 9.2.5 对1,2,k，有名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - -

23、 - 第 9 页，共 40 页 - - - - - - - - - 10 1(),max(),kkkkkkkpmmAI yzyzy(9.2.30) 且111lim,lim,maxkkkknmpxzx(9.2.31) 其收敛速度由1() ()nnpp确定 . 由式 (9.2.28)知：在计算1B的按模最大特征值1np的过程 (9.2.30)中，收敛速度得到加速 . 算法 9.2.5 又称为 原点平移下的反幂法(inverse power method with shift) . 定义 9.2.8 原点平移下的幂法与原点平移下的反幂法统称为原点平移法 .注 有的资料上的原点平移法专指原点平移下的反

24、幂法；而有的资料上的反幂法指的就是原点平移下的反幂法. 原点平移下的反幂法( 算法 9.2.7)的主要应用是：已知矩阵的近似特征值后，求矩阵的特征向量 . 其主要思想是：若已知A的特征值m的近似特征值为m，则mAI的特征值就是(1,2,)imim，其中(1,2,)iim是A的特征值 . 而按模最小的特征值是mm，相应的特征向量与A的特征向量相同. 利用公式 (9.2.30)可求出,kkmz，并且有1lim,limmax()mkkkkmmmmxzx，其收敛速度由0011(min)mmimimimim确定 . 于是当k充分大时，可取mkxz，1mmkm, 只要近似值m适当好，收敛速度就很快，往往迭

25、代几次就能得到满意的结果. 例 9.2.2已知特征值的近似值0.3589，用原点平移下的反幂法求方阵123213335A得对应特征值的特征向量 . 解取0.3589p，对矩阵名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 40 页 - - - - - - - - - 11 1.3589230.358921.35893335.3589ApIAI. 迭代公式 (9.2.30)中的ky是通过解方程组1()kkpAIyz求得的 . 为了节省工作量，可先将pAI进行 LU 分解

26、. 在 LU 分解中尽量避免较小的rru当除数，通常可以先对矩阵pAI的行进行调换后再分解 . 为此，我们可用001010100P乘pAI后再进行LU 分解，即60011.358923335.3589()01021.3589321.35893100335.35891.358923100335.35890.66671000.64110.5726,0.453011003.0710pP AILU1(1.2679 )kkP AIyPz， 即1kkLUyPz. 令11,.kkkkLvPzUyv选取0z，使1T10(1,1,1)UyL Pz，得T1(290929.45, 290927.56,325732.

27、90)y，11max()325732.90my，T111( 0.8932,0.8931, 1)mzy. 由121UyL Pz得T2( 845418.49,845418.49, 946558.42)y，22max()946558.42my，T222( 0.8932,0.8932, 1)mzy. 由于1z与2z的对应分量几乎相等，故A的特征值为2110.35890.35889894354117946558.42m, 相应的特征向量为T2( 0.8932,0.8932, 1)xz. 而矩阵A的一个特征值为4190.358898943540674，相应的特征向量为T( 0.89315,0.89315,

28、 1)，由此可见得到的结果具有较高的精度. 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 40 页 - - - - - - - - - 12 9.3 QR 算法上一节我们介绍了求矩阵特征值的幂法和反幂法. 幂法主要用来求矩阵的按模最大特征值，而反幂法主要用于求特征向量. 本节将介绍幂法的推广和变形 QR 算法，它是求一般中小型矩阵全部特征值的最有效的方法之一，其基本思想就是利用矩阵的QR 分解 . 矩阵n nRA的QR分解就是：用Householder 变换将矩阵A分

29、解成正交矩阵Q与上三角矩阵R的乘积，即AQR. 下面首先介绍Householder 变换 .9.3.1 豪斯霍尔德变换定义9.3.1设()ijn nbB是n阶方阵，若当1ij时，0ijb，则称矩阵B为上Hessenberg矩阵 (Hessenberg matrix)，又称为 准上三角矩阵 ，它的一般形式为1112121222323,1nnnn nnnbbbbbbbbbbB. (9.3.1) 下面讨论如何将矩阵A用正交相似变换化成(9.3.1)的形式 . 为此先介绍一个对称正交矩阵 Householder 矩阵 . 定义 9.3.2设向量nRu的欧氏长度21u，I为n阶单位矩阵，则称n阶方阵T(

30、 )2HH uIuu(9.3.2) 为 Householder 矩阵 (Householder matrix) . 对任何nRx，称由( )HH u确定的变换yHx(9.3.3) 为 镜 面 反 射 变 换 (specular reflection transformation) ， 或Householder (Householder transformation) . 注Householder 变换，最初由A.C Aitken在 1932 年提出 . Alston Scott Householder在1958 年指出了这一变换在数值线性代数上的意义. 这一变换将一个向量变换为由一个超平面反射

31、的镜像，是一种线性变换. 运用线性代数的知识，很容易证明：定理9.3.3(9.3.2)式定义的矩阵H是对称正交矩阵；对任何nRx，由线性变换Hyx得到y的欧氏长度满足22yx. 反之，有下列结论：定理 9.3.4设,nRx y，xy. 若22xy，则一定存在由单位向量确定的镜面反射矩阵( )H u，使得Hxy. 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页，共 40 页 - - - - - - - - - 13 证令2xyuxy，显然21u. 构造单位向量u确定的镜面反射

32、矩阵T( )2HH uIuu，TTTT2222()()2()()22xy xyxy x xy xHxIuuxIxxxyxy. 又因为22xy，即TTx xy y，所以海森伯格 (Karl Adolf Hessenberg1904 年 9 月 8 日 1959 年 2 月 22 日) 是德国数学家和工程师. 豪斯霍尔德 (Alston ScottHouseholder 1904 年 5 月 5 日 1993 年 7 月 4 日) 是美国数学家， 在数学生物学和数值分析等领域卓有建树. 2TTT2TTTTTTTTTT() ()()()2(),xyxyxyxyxyx xx yy xy yx xx y

33、x yx xx xy x于是TT222()()()xyx xy xHxxxxyyxy. 证毕！由定理 9.3.4 得：算法 9.3.1若T12(,)nx xxx，其中2,nxx不全为零，则由12T11212T1T22sgn(),(1,0,0),(),( )2nxxR12其中=xuxeeuuuHH uIIuuu(9.3.4) 确定的镜面反射矩阵H，使得1Hxe，其中1,0,sgn( )0,0,1,0.aaaa例9.3.1设T( 1,2,2)x， 按 式 (9.3.4) 的 方 法 构造 镜 面 反 射 矩 阵H， 使 得T( ,0, 0)*Hx(*表示某非零元素). 解22212sgn()( 1

34、) ( 1)2( 2)3xx；TTT1( 1,2,2)(3, 0, 0)( 4, 2,2)uxe, 其中T1(1,0,0)e；名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 13 页，共 40 页 - - - - - - - - - 14 212()3 3( 1)12x12u，则所求镜面反射矩阵为1T10041 32 32 310102( 4, 2,2)2 32 31 31200122 31 32 3HIuu，且1 32 32 3132 32 31 3202 31 32 320Hx.

35、 可以证明：定理 9.3.6对任意n阶方阵()ijn naA，存在正交矩阵Q，使得TBQ AQ为形如式 (9.3.1)的上 Hessenberg阵. 证记(1)(1)(1)1112111121(1)(1)(1)21222212221(1)(1)(1)1212nnnnnnnnnnnnaaaaaaaaaaaaaaaaaaAA，(1)21(1)311(1)1naaax. 由式 (9.3.4)可构造1n阶对称正交矩阵1H：1 22121121122T1111 11211112121T1111sgn()sgn(),(1,0,0),(),niinaaaaR12其中=xuxeeuHIu u(9.3.5) 使

36、得111 1H xe. 记111IQH，且1n nRQ，1I是 1 阶单位矩阵 . 显然1Q是对称正交矩阵. 用1Q对A作相似变换，得名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页，共 40 页 - - - - - - - - - 15 (1)(2)(2)11121(2)(2)12221(2)(2)111112323(2)(2)200nnnnnnaaaaaaaaa记Q AQQ AQA. (9.3.6) 记(2)(2)(2)T2232422(,)nnaaaRx，同理可构造2n

37、阶对称正交矩阵2H，使得2221H xe(T21(1,0, 0)nRe). 记222IQH，其中2I为 2 阶单位矩阵，则2Q仍是对称正交矩阵，用2Q对2A作相似变换，得(1)(2)(3)(3)1112131(2)(3)(3)122232(3)(3)123332222223(3)(3)434(3)(3)300000nnnnnnnaaaaaaaaaaaaa记Q A QQ A QA. (9.3.7) 依此类推，经过k步对称正交相似变换，得1111111kkkkkkQAQQAQ(1)(2)(3)(1)( )( )()1112131,111,11(2)(3)(1)( )()( )122232,122,

38、12(3)(1)( )( )()2333,133,133(1)( )()()1,11,1,11,000000kkkkkkknkkkkkkknkkkkkkknkkkkkkkkkkknaaaaaaaaaaaaaaaaaaaaaa( )()()1,1( )()()1,1,11,( )()(),100000000000kkkkkkkk kknkkkkkkkknkkknkn knnaaaaaaaaa记A. (9.3.8) 重复上述过程，则有(1)(2)(3)(1)( )1112131,11(2)(3)(1)( )122232,12(3)(1)( )12333,1322222213(1)( )1,11,(

39、 )1nnnnnnnnnnnnnnnnnnnnnnnnnnnnnaaaaaaaaaaaaaaa记QAQQAQA. (9.3.9) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 40 页 - - - - - - - - - 16 由式 (9.3.6)-(9.3.9) ，可得122223332231132nnnnnnnnnnnnnAQAQQQAQQQQQ AQQQ. 若记1nBA，122nQQ QQ，则Q为正交矩阵，且有TBQ AQ. 证毕！注 1由定理 9.3.6 知：

40、因为任意n阶方阵A与n阶上 Hessenberg矩阵B相似，所以求A的阵特征值的问题，就可转化为求上Hessenberg矩阵B的特征值的问题. 注 2若A是对称矩阵，则B也是对称矩阵 . 再由B的形式 (9.3.1)知，此时B一定是对称三对角阵. 于是，求对称矩阵A的阵特征值的问题，便可转化为求对称三对角阵B的特征值问题 . 例 9.3.2设矩阵1212221111112111A试用镜面反射变换求正交矩阵Q，使TQ AQ为上 Hessenberg矩阵 . 解第 1 步记1AA，(1)(1)(1)TT1213141(,)(2, 1, 2)aaax，利用式 (9.3.4)构造三阶镜面反射阵1H：2

41、22112sgn(2)2123x；TTT111 1(2, 1, 2)(3, 0, 0)(5, 1, 2)uxe, 其中T1(1, 0, 0)e；2(1)1111212()3(32)15a12u；则所求镜面反射矩阵为1T111110050.66670.33330.666710101 (5, 1, 2)0.33330.93330.1333 ,1500120.66670.13330.7333HIu u111100000.66670.33330.6667,00.33330.93330.133300.66670.13330.7333IQHT211111130032.33330.46670.066700.

42、46671.57331.346700.06671.34670.0933AQ AQQ AQ. 第 2 步记(2)(2)TT23242(,)(0.4667,0.0667)aax, 利用式 (9.3.4)构造 2 阶镜面反射阵2H：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 40 页 - - - - - - - - - 17 22222sgn(0.4667)0.4667( 0.0667)0.4714x；TTT2221(0.4667,0.0667)(0.4714, 0)(0

43、.9381,0.0667)uxe, 其中T1(1, 0)e；2(2)2222322()0.4714(0.47140.4667)0.4422a12u；则所求镜面反射矩阵为1T2222100.93810.9901 0.14151(0.9381,0.0667),010.06670.14150.98990.4422HIu u22210000100,000.99010.1415000.14150.9899IQHT3222222130032.33330.4714000.47141.57331.5000001.50000.5000AQ A QQ A Q. 由此得正交矩阵12100000.66670.2357

44、0.707100.33330.94290.000100.66670.23570.7070QQ Q，使得T3130032.33330.4714000.47141.16671.5000001.50000.5000Q AQA为上 Hessenberg矩阵。9.3.2 QR 算法QR 算法的基本思想是：利用QR分解得到一系列与A相似的矩阵kA，在一定的条件下，当k时，kA收敛到一个以A的特征值(1,2, )iin为主对角线元素的上三 角 矩 阵 . 首 先 介 绍 用QR分 解(QR decomposition 或QR factorization) ; 即 用Householder 变换将矩阵A分解成

45、正交矩阵Q与上三角矩阵R的乘积，即AQR. 9.3.2.1 QR 分解算法 9.3.2 (QR 分解) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页，共 40 页 - - - - - - - - - 18 第 1 步记A的第 1 列为(1)(1)(1)T111211(,)naaax，(1)1()ijn naAA. 利用式 (9.3.4)：1 22(1)(1)11111T111 11(1)111111T1111sgn(),(1,0,0),(),niinaaaR其中=uxee

46、HIu u构造出的1H是n阶对称正交矩阵，使得111 1H xe，从而有(2)(2)1121(2)(2)222211(2)(2)20.0nnnnnaaaaaaAH A第 2 步记(2)(2)(2)T222322(,)naaax，同理可构造出1n阶对称正交矩阵2H，使得2222H xe，其中1 2(2)222222sgn()niiaa，T12(1,0, 0)nRe. 若记221,HH它仍是对称正交矩阵，于是有(2)(2)(2)112131(3)(3)2232(3)(3)322333(3)(3)30.0000nnnnnnaaaaaaaaaAH A如此继续下去，直到完成第1n步后，得到上三角矩阵(2

47、)(2)(2)112131(3)(3)223211( )1,.nnnnnnnnnaaaaaaAHA于是有111221211121.nnnnnnnnnnAHAHHAHHH AHHH A令nRA，121nQH HH，其中Q是对称正交矩阵，则RQA. 因为Q是对称正交名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 18 页，共 40 页 - - - - - - - - - 19 矩阵，所以得AQR. 注 1若A非奇异，则上三角矩阵R也非奇异，从而R的主对角线元素不为零. 注 2若要求R的主

48、对角线元素取正数，则A的QR分解是唯一的. 例 9.3.3求矩阵101214230A的 QR 分解AQR，并使矩阵R的主对角线上的元素都是正数。解对A运用算法 9.3.2. 第 1 步记1AA, T1(1, 2,2)x，则2221sign(1) 12( 2)3, 11111()3(31)12a=, TTTT111 11(1, 2,2)(3, 0, 0)(4, 2,2) ,(1, 0, 0)uxee, 1T11111004122110102(4, 2,2)2211230012212HIu u, 21134 37 305 310 307 32 3AH A. 第 2 步记T2(5 3, 7 3)x，

49、222sign(5 3) (5 3)(7 3)2.86744, (2)22222()2.86744(2.867445 3)13.0013a=, TTT2221(5 3, 7 3)(2.86744, 0)(4.53411, 2.33333)uxe, 1T2222104.534111(4.53411, 2.33333)012.3333313.0013HIu u0.581240.813730.813730.58124, 记名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 19 页，共 40

50、页 - - - - - - - - - 20 221001000.581240.81373000.813730.58124HH，则32231.333332.3333302.867442.47995002.32494AH A. 为了使R的主对角线上的元素都是正数，取3100010001H，显然3H是正交阵，且43331.333332.3333302.867442.47995002.32494AH A. 令431.333332.3333302.867442.47995002.32494RA，1230.333330.154990.929980.666670.658740.348740.666670.