《第9章--矩阵特征值的数值解法(共40页).doc》由会员分享，可在线阅读，更多相关《第9章--矩阵特征值的数值解法(共40页).doc（40页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上第9章 矩阵特征值的数值解法9.1 引言矩阵特征值问题有广泛的应用背景. 例如动力系统和结构系统中的振动问题、电力系统的静态稳定分析上、工程设计中的某些临界值的确定等，都归结为矩阵特征值问题. 数学中诸如方阵的对角化及解微分方程组等问题，都要用到特征值的理论. 本章介绍n阶实矩阵的特征值与特征向量的数值解法.定义9.1.1 已知n阶实矩阵，如果存在常数和非零向量x，使 或 (9.1.1)那么称为A的特征值(eigenvalue)，为的相应于的特征向量(eigenvector). 多项式 (9.1.2)称为特征多项式(characteristic polynomial)
2、， (9.1.3)称为特征方程(characteristic equation).注 式(9.1.3)是以为未知量的一元n次代数方程，是的n次多项式. 显然，的特征值就是特征方程(9.1.3)的根. 特征方程(9.1.3)在复数范围内恒有解，其个数为方程的次数(重根按重数计算)，因此n阶矩阵在复数范围内有n个特征值. 除特殊情况 (如或为上(下)三角矩阵)外，一般不通过直接求解特征方程(9.1.3)来求的特征值, 原因是这样的算法往往不稳定. 在计算上常用的方法是幂法与反幂法和相似变换方法. 本章只介绍求矩阵特征值与特征向量的这两种基本方法. 为此将一些特征值和特征向量的性质列在此处.定理9.
3、1.2 设阶方阵的特征值为，那么(1) ;(2) .定理9.1.3 如果是方阵的特征值，那么(1) 是的特征值，其中是正整数;(2) 当是非奇异阵时，是的特征值.(3) 是的特征值，其中是多项式.定义9.1.4 设都是阶方阵. 若有阶非奇异阵，使得，则称矩阵与相似(similar)，称为对进行相似变换(similarity transformation)，称为相似变换矩阵(similarity transformation matrix).定理9.1.5 若矩阵与相似，则与的特征值相同.定理9.1.6 如果是阶正交矩阵，那么(1) ，且或;(2) 若，则, 即.定理9.1.7 设是任意阶实对称
4、矩阵，则(1) 的特征值都是实数；(2) 有个线性无关的特征向量.定理9.1.8 设是任意阶实对称矩阵，则必存在阶正交矩阵，使得，其中是以的个特征值为对角元素的对角矩阵.定理9.1.9 (圆盘定理) 矩阵的任意一个特征值至少位于复平面上的几个圆盘，中的一个圆盘上。9.2 幂法与反幂法9.2.1 幂法及其加速9.2.1.1 幂法幂法是计算矩阵按模最大特征值(largest eigenvalue in magnitude)及相应特征向量的迭代法. 该方法稍加修改，也可用来确定其他特征值. 幂法的一个很有用的特性是：它不仅可以求特征值，而且可以求相应的特征向量. 实际上，幂法经常用来求通过其他方法确
5、定的特征值的特征向量. 下面探讨幂法的具体过程.设矩阵的n个特征值满足, (9.2.1)且有相应的n个线性无关的特征向量，则构成n维向量空间的一组基, 因此. 在中选取某个满足的非零向量.用矩阵同时左乘上式两边，得.再用矩阵左乘上式两边，得.这样继续下去，一般地有 (9.2.2)记，则由式(9.2.2)得 (9.2.3)由假设(9.2.1)，结合式(9.2.3)，得 (9.2.4)于是对充分大的k有 (9.2.5)式(9.2.4)表明随着k的增大，序列越来越接近A的对应于特征值的特征向量的倍, 由此可确定对应于的特征向量. 当k充分大时，可得的近似值. 上述收敛速度取决于比值. 事实上，由式(
6、9.2.3)知，. (9.2.6)再由式(9.2.1)得 . (9.2.7)结合式(9.2.6)和式(9.2.7)知，序列收敛速度取决于比值. 下面计算. 由式(9.2.3)知当k充分大时, . 结合式(9.2.5)，得.这表明两个相邻向量大体上只差一个常数倍，这个倍数就是A的按模最大特征值. 记, 则有， (9.2.8)即两个相邻的迭代向量所有对应分量的比值收敛到.定义9.2.1 上述由已知非零向量及矩阵的乘幂构造向量序列来计算的按模最大特征值及相应特征向量的方法称为幂法(power method)，其收敛速度由比值来确定，越小，收敛越快.注 由幂法的迭代过程(9.2.3)容易看出，如果（或
7、），那么迭代向量的各个非零的分量将随着而趋于无穷（或趋于零），这样在计算机上实现时就可能上溢(或下溢). 为了克服这个缺点，需将每步迭代向量进行规范化: 记.若存在的某个分量，满足，则记. 将规范化,这样就把的分量全部控制在中. 例如，设，因为的所有分量中，绝对值最大的的是，所以，故. 综上所述，得到下列算法：算法9.2.1 (幂法) 设是阶实矩阵，取初始向量，通常取，其迭代过程是：对，有 (9.2.9)定理9.2.1 对式(9.2.9)中的序列和有, ， (9.2.10)其收敛速度由确定. 证明 由迭代过程(9.2.9)知, (9.2.11)其中. 若，则由(9.2.3)知：, 代入式(9.
8、2.11)得， (9.2.12)故. (9.2.13)而 , (9.2.14)于是, (9.2.15)故. (9.2.16)由式(9.2.6)和式(9.2.7)知：上述收敛速度由确定. 证毕！ 例9.2.1 用幂法求方阵按模最大特征值及相应的特征向量，要求.解 选取初始向量，按式(9.2.9)迭代，结果见表9.2.1.表9.2.112345因此，所求按模最大特征值，相应特征向量. 事实上，的按模最大特征值，相应特征向量, 故所得结果具有较高的精度.9.2.1.2 幂法的加速从上面的讨论可知，由幂法求按模最大特征值，可归结为求数列的极限值，其收敛速度由确定. 当接近时, 收敛速度相当缓慢. 为了
9、提高收敛速度, 可以利用外推法进行加速.因为序列的收敛速度由确定，所以若收敛，当充分大时，则有，或改写为,其中是与无关的常数. 由此可得, (9.2.17)这表明幂法是线性收敛的. 由式(9.2.17)知.由上式解出，并记为，即， (9.2.18)这就是计算按模最大特征值的加速公式.将上面的分析归结为如下算法：算法9.2.2 (幂法的加速) 设是阶实矩阵，给定非零初始向量，通常取. 对，用迭代式求出及. 再对，迭代过程为 (9.2.19)当(是预先给定的精度) 时，迭代结束，并计算；否则继续迭代，直至满足迭代停止条件.有关加速收敛技术，读者请参考文献11.8.2.2 反幂法及其加速反幂法是计算
10、矩阵按模最小特征值及相应特征向量的迭代法, 其基本思想是：设矩阵非奇异，用其逆矩阵代替, 矩阵的按模最小特征值就是矩阵的按模最大特征值. 这样用代替做幂法，即可求出的按模最大特征值，也就是矩阵的按模最小特征值；这种方法称为反幂法(inverse power method). 因为矩阵非奇异，所以由可知：. 这说明：如果的特征值满足, (9.2.20)那么的特征值满足, (9.2.21)且对的应于特征值的特征向量也是的对应于特征值的特征向量.由上述分析知：对应用幂法求按模最大的特征值及相应的特征向量，就是求的按模最小的特征值及相应的特征向量. 算法9.2.3 (反幂法) 任取初始非零向量，通常取
11、. 为了避免求，对，将迭代过程(9.2.9)改写为: (9.2.22)仿定理9.2.1的证明，可得：定理9.2.2 对式(9.2.22)中的序列和有, (9.2.23)其收敛速度由确定.注 按(9.2.22)进行计算，每次迭代都需要解一个方程组. 若利用三角分解法求解方程组，即，其中是下三角矩阵，是上三角矩阵，这样每次迭代只需解两个三角方程组9.2.3 原点平移法为了提高收敛速度，下面介绍加速收敛的原点平移法. 设矩阵，其中是一个待定的常数，与除主对角线上的元素外，其他元素相同. 设的特征值为，则的特征值为,且与的特征向量相同. 9.2.3.1 原点平移法的幂法设是的按模最大的特征值，选择，使
12、，及. (9.2.24)对应用幂法，得算法9.2.4 对，有 (9.2.25)且 (9.2.26)其收敛速度由确定. 由式(9.2.24)知：在计算的按模最大特征值的过程(9.2.25)中，收敛速度得到加速；算法9.2.4又称为原点平移下的幂法(power method with shift). 9.2.3.2 原点平移下的反幂法设是的按模最小的特征值，选择，使. (9.2.27)及. (9.2.28)若矩阵可逆，则的特征值为且有. (9.2.29)对应用反幂法，得：算法9.2.5 对，有 (9.2.30)且 (9.2.31)其收敛速度由确定.由式(9.2.28)知：在计算的按模最大特征值的过
13、程(9.2.30)中，收敛速度得到加速. 算法9.2.5又称为原点平移下的反幂法 (inverse power method with shift).定义9.2.8 原点平移下的幂法与原点平移下的反幂法统称为原点平移法. 注 有的资料上的原点平移法专指原点平移下的反幂法；而有的资料上的反幂法指的就是原点平移下的反幂法.原点平移下的反幂法(算法9.2.7)的主要应用是：已知矩阵的近似特征值后，求矩阵的特征向量. 其主要思想是：若已知的特征值的近似特征值为，则的特征值就是，其中是的特征值. 而按模最小的特征值是，相应的特征向量与的特征向量相同. 利用公式(9.2.30)可求出，并且有，其收敛速度由
14、确定. 于是当充分大时，可取， ,只要近似值适当好，收敛速度就很快，往往迭代几次就能得到满意的结果.例9.2.2 已知特征值的近似值，用原点平移下的反幂法求方阵得对应特征值的特征向量.解 取，对矩阵.迭代公式(9.2.30)中的是通过解方程组求得的. 为了节省工作量，可先将进行LU分解.在LU分解中尽量避免较小的当除数，通常可以先对矩阵的行进行调换后再分解. 为此，我们可用乘后再进行LU分解，即， 即 .令 选取，使，得，.由得，.由于与的对应分量几乎相等，故的特征值为,相应的特征向量为. 而矩阵的一个特征值为，相应的特征向量为，由此可见得到的结果具有较高的精度.9.3 QR算法上一节我们介绍
15、了求矩阵特征值的幂法和反幂法. 幂法主要用来求矩阵的按模最大特征值，而反幂法主要用于求特征向量. 本节将介绍幂法的推广和变形QR算法，它是求一般中小型矩阵全部特征值的最有效的方法之一，其基本思想就是利用矩阵的QR分解. 矩阵的分解就是：用Householder变换将矩阵分解成正交矩阵与上三角矩阵的乘积，即. 下面首先介绍Householder变换.9.3.1 豪斯霍尔德变换定义9.3.1 设是阶方阵，若当时，则称矩阵为上Hessenberg矩阵(Hessenberg matrix)，又称为准上三角矩阵，它的一般形式为. (9.3.1)下面讨论如何将矩阵用正交相似变换化成(9.3.1)的形式.
16、为此先介绍一个对称正交矩阵Householder矩阵.定义9.3.2 设向量的欧氏长度，为阶单位矩阵，则称阶方阵 (9.3.2)为Householder矩阵(Householder matrix). 对任何，称由确定的变换 (9.3.3)为镜面反射变换(specular reflection transformation)，或Householder (Householder transformation).注 Householder变换，最初由A.C Aitken在1932年提出. Alston Scott Householder在1958年指出了这一变换在数值线性代数上的意义. 这一变换将一
17、个向量变换为由一个超平面反射的镜像，是一种线性变换. 运用线性代数的知识，很容易证明：定理9.3.3 (9.3.2)式定义的矩阵是对称正交矩阵；对任何，由线性变换得到的欧氏长度满足.反之，有下列结论：定理9.3.4 设，. 若，则一定存在由单位向量确定的镜面反射矩阵，使得.证 令，显然. 构造单位向量确定的镜面反射矩阵，.又因为，即，所以海森伯格(Karl Adolf Hessenberg 1904年9月8日1959年2月22日) 是德国数学家和工程师.豪斯霍尔德(Alston Scott Householder 1904年5月5日1993年7月4日) 是美国数学家，在数学生物学和数值分析等领
18、域卓有建树.于是.证毕！由定理9.3.4得：算法9.3.1 若，其中不全为零，则由 (9.3.4)确定的镜面反射矩阵，使得，其中例9.3.1 设，按式(9.3.4)的方法构造镜面反射矩阵，使得(表示某非零元素).解 ；, 其中；，则所求镜面反射矩阵为，且.可以证明：定理9.3.6 对任意阶方阵，存在正交矩阵，使得为形如式(9.3.1)的上Hessenberg阵.证 记，.由式(9.3.4)可构造阶对称正交矩阵： (9.3.5)使得.记，且，是1阶单位矩阵. 显然是对称正交矩阵. 用对作相似变换，得. (9.3.6)记，同理可构造阶对称正交矩阵，使得 (). 记，其中为2阶单位矩阵，则仍是对称正
19、交矩阵，用对作相似变换，得. (9.3.7)依此类推，经过步对称正交相似变换，得. (9.3.8)重复上述过程，则有. (9.3.9)由式(9.3.6)-(9.3.9)，可得.若记，则为正交矩阵，且有. 证毕！ 注1 由定理9.3.6知：因为任意阶方阵与阶上Hessenberg矩阵相似，所以求的阵特征值的问题，就可转化为求上Hessenberg矩阵的特征值的问题.注2 若是对称矩阵，则也是对称矩阵. 再由的形式(9.3.1)知，此时一定是对称三对角阵. 于是，求对称矩阵的阵特征值的问题，便可转化为求对称三对角阵的特征值问题.例9.3.2 设矩阵试用镜面反射变换求正交矩阵，使为上Hessenbe
20、rg矩阵.解 第1步 记，利用式(9.3.4)构造三阶镜面反射阵：；, 其中；则所求镜面反射矩阵为 .第2步 记, 利用式(9.3.4)构造2阶镜面反射阵：；, 其中；则所求镜面反射矩阵为.由此得正交矩阵，使得 为上Hessenberg矩阵。9.3.2 QR算法QR算法的基本思想是：利用分解得到一系列与相似的矩阵，在一定的条件下，当时，收敛到一个以的特征值为主对角线元素的上三角矩阵. 首先介绍用QR分解 (QR decomposition 或 QR factorization); 即用Householder变换将矩阵分解成正交矩阵与上三角矩阵的乘积，即. 9.3.2.1 QR分解算法9.3.2
21、 (QR分解) 第1步 记的第1列为，. 利用式(9.3.4)：构造出的是阶对称正交矩阵，使得，从而有第2步 记，同理可构造出阶对称正交矩阵，使得，其中，.若记 它仍是对称正交矩阵，于是有如此继续下去，直到完成第步后，得到上三角矩阵于是有令，其中是对称正交矩阵，则. 因为是对称正交矩阵，所以得.注1 若非奇异，则上三角矩阵也非奇异，从而的主对角线元素不为零. 注2 若要求的主对角线元素取正数，则的分解是唯一的. 例9.3.3 求矩阵的QR分解，并使矩阵的主对角线上的元素都是正数。解 对运用算法9.3.2. 第1步 记, ，则, , , ,.第2步 记，, , , ,记，则.为了使的主对角线上的
22、元素都是正数，取，显然是正交阵，且.令，且. 9.3.2.2 QR算法了解了分解后，下面介绍算法(QR algorithm).算法9.3.3 (QR算法) 第1步 令，利用算法9.3.2将进行分解，得，其中为正交矩阵，为上三角矩阵；然后将与逆序相乘，得.因为，所以有，即与相似.第2步 以代替，再作分解，得，其中为正交矩阵，为上三角矩阵；再将与逆序相乘，并记.因为，所以，即与相似.依此类推，可得算法公式：对， (9.3.10)因为, 所以，即与相似. 故序列相似于. 这里，我们不加证明地给出QR算法收敛的充分条件：定理9.3.9 (QR算法的收敛性) 设, 是由QR算法产生的矩阵序列，其中. 若
23、(1) 的特征值满足；(2) , 其中，且有三角分解(是单位下三角矩阵，是上三角矩阵)，则 (9.3.11)即收敛到一个以的特征值为主对角线元素的上三角矩阵. 推论 若矩阵是对称矩阵，且满足定理9.3.9中的条件，则由QR算法产生的矩阵序列收敛到对角矩阵.9.3.3 带原点位移的QR算法QR算法收敛速度是线性的，需要提高收敛速度. 经分析知：定理9.3.9中的速度依赖于比值，当越小，收敛速度越快. 这启发我们把9.2节介绍的原点平移法加速幂法和反幂法的思想运用到QR算法的收敛加速，这就是下面将要介绍的带原点位移的QR算法(QR algorithm with shift). 因为求上Hessen
24、berg矩阵的特征值比一般的方阵简单，所以在下面的算法中首先将所考察的矩阵化成上Hessenberg矩阵，然后再用位移加速方法进行加速.算法9.3.4 (带原点位移的QR算法)第1步 将矩阵化成上Hessenberg矩阵.第2步 对，给定原点位移数列，进行迭代加速 (9.3.12)并称由到的变换为带原点位移的QR变换 (QR transformation with shift)，且记. 其中的选取方法主要有(a) 取；或(b) 取为的右下角主子矩阵的2个特征值中最靠近的那一个.注1 由算法9.3.4产生的与相似，即；且它们都是上Hessenberg矩阵. 注2 计算实践表明，取方法(b)中的算
25、法比取方法(a)中的算法收敛速度更快. 特别是对称矩阵，带方法(b)中位移的QR算法是无条件收敛的，且收敛阶为3.注3 在迭代过程中，当时，由定理9.3.4知：，因此可取作为的近似特征值. 而的速度依赖于比值，越小，收敛速度越快；因此平移值取显然是一个很好选择.注4 在迭代过程中，当时，取的2个特征值作为的特征值的近似值.注5 判别和约等于0的标准可取为其中是预先给定的精度. 注6 当求得矩阵的特征值或后，可以将作降阶处理：对左上角的阶或阶主子矩阵继续运用带原点位移的QR算法, 以求的其他特征值，这样可以大大节省运算量.9.4 Jacobi方法上一节我们介绍了Householder变换将矩阵化
26、为上Hessenberg矩阵的方法. 如果矩阵是实对称矩阵，用平面旋转变换将矩阵化为上Hessenberg矩阵比用Householder变换更好. 本节介绍的Jacobi方法就是这种方法. 它主要用于求实对称矩阵的全部特征值和特征向量. Jacobi方法(Jacobi method)的基本思想是：对实对称矩阵，一定存在正交矩阵，使, (9.4.1)其中，就是矩阵的特征值，而正交矩阵的第列就是对应于的特征向量. 由此可见，Jacobi方法的实质和关键就是找一个正交矩阵，将化为对角矩阵.9.4.1 Jacobi方法定义9.4.1 设是阶实对称矩阵，称阶矩阵 (9.4.2)为旋转矩阵(rotatio
27、n matrix)，或Givens矩阵(Givens matrix)，简记为. 对进行的变换 (9.4.3)称为Givens旋转变换(Givens rotation).注1 Givens矩阵是在阶单位矩阵的第行第列、第行第列、第行第列、第行第列的交叉的位置上分别换上、而成的. 注2 Givens矩阵是正交矩阵，变换(9.4.3)是正交相似变换. 注3 Jacobi方法就是通过一系列Givens旋转变换，把化为对角矩阵，从而求得特征值及相应的 特征向量的方法. 因此，Jacobi方法也称为平面旋转法(plane rotation method).下面介绍具体介绍将阶实对称矩阵化为对角矩阵的Jac
28、obi方法.设是阶实对称矩阵，记. (9.4.4)因为,所以仍是对称矩阵. 通过直接计算可得基文斯(J. W. Givens 日 年日) 是美国数学家，计算机领域的先驱之一. (9.4.5)不难看出，经过作用后，与相比，只有的第行、第列、第行、第列的元素发生了变化，而其他元素与的相同. 特别地，若，由(9.4.5)式中最后的式子知：若取满足关系式 (9.4.6)可使，也就是说，用对进行变换，可将的2个非主对角线元素和化为零. Jacobi方法的一般过程是：记，选择的一对最大的非零非主对角线元素和，使用Givens矩阵对作正交相似变换得，可使的这对非零非主对角线元素. 再选择的一对最大的非零非主
29、对角线元素作上述正交相似变换得，可使的这对非零非主对角线元素化为零. 如此不断地做下去，可产生一个矩阵序列. 注 虽然至多只有对非零非主对角线元素，但是不能期望通过次变换使对角化. 因为每次变换能使一对非零非主对角线元素化为零，例如，和化为零. 但在下一次变换时，它们又可能由零变为非零. 不过可以证明，如此产生的矩阵序列将趋向于对角矩阵，即Jacobi方法是收敛的. 而这个对角矩阵的主对角线元素就是矩阵的特征值.用Jacobi方法求矩阵的特征值的步骤为：(1) 记，在矩阵中找出按模最大的非主对角线元素，取相应的Givens矩阵，记为；(2) 由条件，定出, . 为避免使用三角函数，令 (9.4
30、.7)(3) 按公式(9.4.5)计算的元素；(4) 以代替，重复步骤(1)， (2) ，(3)，求出；依此类推，得. (9.4.8)令，记，则是正交矩阵，且. (9.4.9)若经过步旋转变换，的所有非主对角线元素都小于允许误差时，停止计算. 此时的主对角线元素就是的特征值的近似值. 的列元素就是的对应于上述特征值的全部特征向量.9.4.2 Jacobi方法的收敛性由矩阵理论知：定理9.4.2 设，是正交矩阵，则. (9.4.10)定理9.4.3(收敛性) 设是由Jacobi方法产生的矩阵序列，其中由式(9.4.9)定义，则, (9.4.11)其中就是矩阵的特征值，而正交矩阵的第列就是对应于的
31、特征向量，. 分析：要证明Jacobi方法的收敛性，只要能证明每次变换，总是使主对角线元素的平方和增大，而非主对角线元素的平方和减小即可. 证明 由定理9.4.2知：在正交相似变换下，矩阵元素的平方和不变. 设矩阵的非主对角线元素平方和为， (9.4.12)因此，要证明收敛于对角矩阵，只要证明收敛于零即可. 假设变为时，把的绝对值最大的非主对角线元素和化为零，由式(9.4.5)直接计算，得,且有，所以 (9.4.13)这表明经过变换后，的非主对角线元素的平方和减小了，由定理9.4.2知：与的元素的平方和不变，故的主对角线元素的平方和应增加. 因为是的绝对值最大的非主对角线元素化为零，所以，即
32、(9.4.14)于是因为，所以当时，. 故. (9.4.15)因为相似矩阵的特征值相同，所以就是矩阵的特征值. 由式(9.4.9)和式(9.4.15)得.因为是正交矩阵，所以，若记，其中是的第列，则有，即，得，亦即.这说明正交矩阵的第列就是对应于的特征向量，. 证毕！例9.4.1 利用Jacobi方法求矩阵的全部特征值和特征向量，要求的所有非主对角线元素的绝对值.解 记，因为是中所有非主对角线元素中绝对值最大的元素，取相应的Givens矩阵. 因为，所以， ，于是,.记则.用代替，重复上述过程： 因为是中所有非主对角线元素中绝对值最大的元素，取相应的Givens矩阵.因为，所以；，于是,记则.
33、用代替，重复上述过程： 因为是中所有非主对角线元素中绝对值最大的元素，取相应的Givens矩阵.因为，所以；，于是.记则. 因为的所有非主对角线元素的绝对值，所以迭代停止. 此时的特征值的近似值分别为的对角元素：.相应的特征向量的近似值分别为其中.的特征值的精确值为.相应的特征向量为从例9.4.1可以看出，即使迭代的次数不是很多，迭代矩阵的所有非对角元素的绝对值并不是很小时，用Jacobi方法求得的结果精度都比较高，因此它是求实对称矩阵的全部特征值和特征向量的一个较好的方法.9.4.3 改进的Jacobi方法由于每次旋转变换前选非零非主对角线元素的最大值很费时间，为此介绍如下改进方法.第一种方
34、法：把非主对角线元素按照行的次序依次化为零，称为一次扫描. 一次扫描后，前面已化为零的元素可能成为非零元素，需要再次扫描. 这一方法称为循环Jacobi方法，这种方法的缺点是: 一些已经足够小的元素也作化零处理，浪费了时间.第二种方法：首先对实对称矩阵计算， (9.4.16)设置阀值，按的顺序进行扫描. 若，则选取旋转矩阵作旋转相似变换将和化为零；否则让过关，即不进行旋转相似变换将其化为零. 因为某些绝对值小于的元素的绝对值可能在后面的旋转变换中增长，所以应进行多次扫描，直到的所有非零非主对角线元素的绝对值都小于. 再设置阀值，重复上述过程，直到达到精度要求，即为止(其中是指定的精度). 这种
35、方法称为限值Jacobi方法.9.5 算法程序9.5.1 幂法%幂法% A是方阵，e是精度，输出向量V是A的最大特征值MaxEig对应的特征向量function PowerMethod(A, e)V = ones(size(A,1), 1);for i=1:2Y = A * V;M(i) = norm(Y, Inf);V = Y / M(i);endwhile abs(M(i)-M(i-1) e i = i + 1; Y = A * V; M(i) = norm(Y, Inf); V = Y / M(i);endMaxEig = M(i);disp(sprintf(最大特征值MaxEig %.
36、6fn相应的特征向量, MaxEig)V=V,end例9.5.1 用幂法求方阵按模最大特征值及相应特征向量，要求.解 在MATLAB命令窗口中输入：A=1 2 3;2 1 3; 3 3 5; e=10-3; PowerMethod(A,e) 回车，输出结果：最大特征值MaxEig =8.相应的特征向量V = 0.5598 0.5598 1.00009.5.2 幂法的加速%幂法的加速%A是方阵，e是精度，输出向量V是A的最大特征值MaxEig对应的特征向量function Acc_ PowerMethod (A, e)V = ones(size(A,1), 1);for i=1:2Y = A *
37、 V;M(i) = norm(Y, Inf);V = Y / M(i);endfor i=3:4 Y = A * V; M(i) = norm(Y, Inf); M2(i) = M(i-2) - (M(i-1) - M(i-2)2 / ( M(i)-2*M(i-1)+M(i-2) ); V=Y/M2(i);endwhile abs(M2(i)-M2(i-1) e i = i + 1; Y = A * V; M(i) = norm(Y, Inf); M2(i) = M(i-2) - (M(i-1) - M(i-2)2 / ( M(i)-2*M(i-1)+M(i-2) ); V = Y / M2(
38、i);endMaxEig = M2(i);disp(sprintf(最大特征值MaxEig %.6fn相应的特征向量, MaxEig)V, end9.5.3 反幂法%反幂法% A是方阵，e是精度，输出向量V是A的最大特征值MaxEig对应的特征向量function InversePower(A, e)V = ones(size(A,1), 1);L U = lu(A);for i=1:2Y = U (LV);M(i) = norm(Y, Inf);Vector = Y / M(i);endwhile abs(M(i)-M(i-1) e i = i + 1; Y = U (LV); M(i) =
39、 norm(Y, Inf); V = Y / M(i);endMinEig = 1/M(i);disp(sprintf(该矩阵按模最小的特征值是：%.6fn相应的特征向量是：, MinEig)V, end9.5.4 原点平移下的反幂法%原点平移下的反幂法% A是方阵，e是精度，%输出向量V是A的特征值Eig对应的特征向量，Appr是特征值Eig的近似值function TranInversePower(A, Appr, e)V = ones(size(A,1), 1);B = A - Appr * eye( size(A,1) );L U = lu(B);for i=1:2Y = U (LV)
40、;M(i) = norm(Y, Inf);V = Y / M(i);endwhile abs(M(i)-M(i-1) e i = i + 1; Y = U (LV); M(i) = norm(Y, Inf); V = Y / M(i);endEig = Appr + 1/M(i);disp(sprintf(该矩阵的特征值是：%.6fn相应的特征向量是：, Eig)V=V,end例9.5.2 已知特征值的近似值，用原点平移下的反幂法求方阵的对应特征值的特征向量.解 在MATLAB命令窗口中输入：A=1 2 3;2 1 3; 3 3 5; Appr=-0.3589; e=10-3; TranInversePower(A,Appr,e)输出的结果是：该矩阵的特征值是：-0.相应的特征向量是：V = 0.8931 0.8931 -1.00009.5.5 Householder变换%Householder变换% 参数X是列向量function Result = Householder ( X )d = sign(X(1) * norm(X, 2);e1 = eye( size(X