矩阵特征与特征向量的计算.doc

《矩阵特征与特征向量的计算.doc》由会员分享，可在线阅读，更多相关《矩阵特征与特征向量的计算.doc（26页珍藏版）》请在三一办公上搜索。

1、第三章 矩阵特征与特征向量的计算3.1 引言在科学技术的应用领域中，许多问题都归为求解一个特征系统。如动力学系统和结构系统中的振动问题，求系统的频率与振型；物理学中的某些临界值的确定等等。设A为n阶方阵，若，有数l使Ax= lx （5.1）则称l为A的特征值，x为相应于l的特征向量。因此，特征问题的求解包括两方面：1求特征值l，满足（5.2）2求特征向量，满足齐方程组（5.3）称j（l）为A的特征多项式，它是关于l的n次代数方程。关于矩阵的特征值，有下列代数理论，定义1 设矩阵A, BR nn，若有可逆阵P，使则称A与B相似。定理1 若矩阵A, BR nn且相似，则（1）A与B的特征值完全相同

2、；（2）若x是B的特征向量，则Px便为A的特征向量。定理2 设AR nn具有完全的特征向量系，即存在n个线性无关的特征向量构成Rn的一组基底，则经相似变换可化A为对角阵，即有可逆阵P，使其中li为A的特征值，P的各列为相应于li的特征向量。定理3 AR nn，l1, , ln为A的特征值，则（1）A的迹数等于特征值之积，即（2）A的行列式值等于全体特征值之积，即定理4 设AR nn为对称矩阵，其特征值l1l2ln，则（1）对任AR n，x0，（2）（3）定理5 （Gerschgorin圆盘定理） 设AR nn，则（1）A的每一个特征值必属于下述某个圆盘之中，（5.4）（5.4）式表示以aii为

3、中心，以半径为的复平面上的n个圆盘。（2）如果矩阵A的m个圆盘组成的并集S（连通的）与其余n m个圆盘不连接，则S内恰包含m个A的特征值。定理4及定理5给出了矩阵特征值的估计方法及界。例1 设有估计A的特征值的范围。解 由圆盘定理，A的3个圆盘为图5.1D1: D2: D3: 见图5.1。D2为弧立圆盘且包含A的一个实特征值l1（因为虚根成对出现的原理），则3l15。而l2，l3D1D2，则，即3.2 乘幂法与反幂法在实际工程应用中，如大型结构的振动系统中，往往要计算振动系统的最低频率（或前几个最低频率）及相应的振型，相应的数学问题便为求解矩阵的按模最大或前几个按模最大特征值及相应的特征向量问

4、题，或称为求主特征值问题。3.2.1 乘幂法乘幂法是用于求大型稀疏矩阵的主特征值的迭代方法，其特点是公式简单，易于上机实现。乘幂法的计算公式为：设AR nn，取初始向量x(0)R n，令x(1) = Ax(0)，x(2) = Ax(1)，一般有（5.5）形成迭代向量序列x(k)。由递推公式（5.5），有（5.6）这表明x(k)是用A的k次幂左乘x(0)得到的，因此称此方法为乘幂法，（5.5）或（5.6）式称为乘幂公式，x(k)称为迭代序列。下面分析乘幂过程，即讨论当k时，x(k)与矩阵A的主特征值及相应特征向量的关系。设A = (aij)nn有完全的特征向量系，且l1, l2, ln为A的n个

5、特征值，满足v1, v2, vn为相应的特征向量且线性无关，从而构成Rn上的一组基底。对任取初始向量x(0) Rn，可由这组基底展开表示为 （5.7）其中a1, a2, an为展开系数。将x(0)的展开式（5.7）代入乘幂公式（5.6）中，得 （5.8）利用（5.8）式为 （5.9）（1）如果A有唯一的主特征值，即，设l1 0，且由（5.9）式，有其中，由于，故当k充分大时，e k 0，此时（5.10）对i = 1, 2, , n，若(a1v1)i 0，考虑相邻迭代向量的对应分量比值，（5.11）即对i = 1, , n（5.12）这表明主特征值l1可由（5.11）或（5.12）式得到。由于迭

6、代序列x(k)，当k充分大时，（5.10）式成立，x(k)与v1只相差一个常数因子，故可取x(k)作为相应于主特征值l1的特征向量的近似值。迭代序列x(k)的收敛速度取决于的大小。（2）如果A的主特征值不唯一，且可分三种情况讨论：a）l1 = l2；b）l1 = - l2；c）情况a）当l1 = l2时，A的主特征值为二重根，根据（5.9）式当k充分小时，由于，j = 3, n，e k 0，则对i = 1, 2, , n，如果，则（主特征值）且x(k)收敛到相应于l1 (=l2)的特征向量的近似值。这种重主特征值的情况，可推广到A的r重主特征值的情况，即当 且时，上述讨论的结论仍然成立。情况b

7、）当l1 = - l2时，A的主特征值为相反数，（5.9）式为当k充分大时，j = 3, 4, , n, e k 0，则（5.13）由于（5.13）式中出现因子(-1) k，则当k变化时，x(k)出现振荡、摆动现象，不收敛，利用(-1) k的特点，连续迭代两步，得从而，对i = 1, 2, , n，若，则（5.14）开方之后，便得到A的以上主特征值l1, l2 = - l1。为计算相应于l1, l2的特征向量，采取组合方式，（5.15）（5.16）可见分别为相应于l1与l2的特征向量。情况c）当时，A的主特征值为共轭复根。因A为实矩阵，于是由有即（v1与v2为互为共轭向量）。设，对任取x(0)

8、 R n，展开式（5.7）可为（5.17）将（5.17）式代入（5.9）式，同理，当k充分大时（5.18）对j = 1, 2, n，设复数表示则（5.18）式的复数表示可为连续迭代，得（5.19）利用三角函数运算性质及l1、l2的复数表示，不难验证。令（5.20）解方程（j = 1, 2, , n）（5.21）求出p, q后，再解出主特征值l1、l2，得 （5.22）同样，采取组合方式求相应于l1、l2的特征向量。由于（5.23）（5.24）则可分别取（5.23）、(5.24）左端的组合表达式作为相应于l1、l2的特征向量的近似值。通过上述分析，有定理6 设A Rnn有完全特征向量系，若l1,

9、 l2, ln为A的n个特征值且满足对任取初始向量x(0) Rn，对乘幂公式确定的迭代序列xk，有下述结论：（1）当时，对i = 1, 2, , n收敛速度取决于的程度，r 1收敛快，r 1收敛慢，且x(k)（当k充分大时）为相应于l1的特征向量的近似值。（2）当时a）若l1 = l2，则主特征值l1及相应特征向量的求法同（1）；b）若l1 = -l2，对i = 1, 2, , n收敛速度取决于的程度。向量、分别为主特征值l1、l2相应的特征向量的近似值。c）若，则连续迭代两次，计算出x(k+1)，x(k+2)，然后对j = 1, 2, , n解方程求出、后，由公式解出主特征值l1、l2。此时

10、收敛速度取决于的程度。向量、分别为相应于l1，l2的特征向量的近似值。从分析乘幂过程可见，乘幂法可用于求矩阵按模最大的一个（或几个）特征值及相应的特征向量，当比值时，收敛速度快，r 1时，收敛速度慢，且计算公式简便，便于上机实现。分析中的假设、，在计算时可不用考虑，如果此条件不满足，则可通过迭代误差自行调整。在用乘幂法求矩阵的主特征值l1及对应的特征向量时，迭代向量的分量可能会出现绝对值非常大的现象，从而造成计算中溢出的可能。为此，需对迭代向量x(k)进行规范化。令max(x)表示向量x分量中绝对值最大者。即如果有某i0，使则max (x) = xi对任取初始向量x(0)，记则一般地，若已知x

11、(k)，称公式（5.25）为规范化的乘幂法公式或改进乘幂法公式，这里，乘幂迭代序列y(k)的分量绝对值最大者1。类似前面的分析乘幂过程，有定理7 设ARnn具有完全特征向量系，l1, l2, , ln为A的n个特征值，且满足则对任初始向量x(0)，由规范化的乘幂法公式（5.25）确定的向量序列y(k)，x(k)满足（1）（5.26）（2）y(k)为相应于主特征值l1的特征向量近似值（5.27）例2 用规范化乘幂法计算矩阵A的主特征值及相应特征向量解 A的特征值l1 = 6，l2 = 3，l3 = 2取初始值x(0) = (1, 1, 1)T，用规范化乘幂法公式（5.25）计算其它结果见表5.1

12、（表中的向量均为转置向量）。表5.1kmax(y(k)x(k) = y(k)/max(x(k)x(k+1) = Ay(k)01(1, 1, 1)(10, 8, 1)110(1, 0.8, 0.1)(7.2, 5.4, -0.8)27.2(1, 0.75, -0.111111)(6.5, 4.75, -1.222222)36.57(1, 0.730769, -0.203704)(6.230766, 4.499997, -1.407408)46.230766(1, 0.722222, -0.225880)(6.111108, 4.388886, -1.1451767)56.111108(1, 0.

13、718182, -0.237561)(6.054548, 4.336336, -1.475122)66.054548(1, 0.716216, -0.243639)(6.027024, 4.310808, -1.487278)76.027024(1, 0.715247, -0.246768)(6.013458, 4.298211, -1.483536)86.013458(1, 0.714765, -0.248366)(6.00671, 4.291945, -1.496732)96.00671(1, 0.714525, -0.249177)(6.00335, 4.28825, -1.496354

14、)106.00335(1, 0.714405, -0.249586)(6.00167, 4.287265, -1.499172)116.00167(1, 0.714345, -0.239792)(6.00083, 4.286485, -1.499584)126.00083(1, 0.714315, -0.249896)取max(x(12) = 6.00083作为主特征值l1的近似值，与真值l1 = 6相比，有较好的近似程度，相应于l1的特征向量的近似值取为y(2) = (1, 0.714315, -0.249896)T。3.2.2 乘幂法的加速及降价当li ( i = 1, 2, , n)为矩

15、阵A Rnn的n个特征值，且 时，乘幂法的收敛速度由决定，r l2 ln时，取则为w(l0)的极小值点。这时原点移位法是一个矩阵变换过程，变换简单且不破坏原矩阵的稀疏性。但由于预先不知道特征值的分布，所以应用起来有一定困难，通常对特征值的分布有一个大略估计，设定一个参数l0进行试算，当所取l0对迭代有明显加速效应以后再进行确定计算。例3 计算A的主特征值解 先用规范化乘幂法计算，得表5.2表5.2ky(k) = x(k)/max(x(k)l1 max (x(k)0(1, 1, 1)1(0.9091, 0.8182, 1)T2.7519(0.7482, 0.6497, 1)T2.53653742

16、0(0.7482, 0.6497, 1)T2.5365323主特征值l1及特征向量v1为（8位有效数字）l1 = 2.5362258v1 = (0.74822116, 0.64966116, 1)T而用规范化乘幂法计算的相应近似值为：l1 max (x(20) = 2.5365323v1 y(20) = (0.7482, 0.6497, 1)T如果采用原点位移的加速法求解，取l0 = 0.75，矩阵b = A - l0I对矩阵B应用规范化乘幂法公式见表5.3。表5.3ky(k) = x(k)/max(x(k)l1 max (x(k)0(1, 1, 1)T9(0.7483, 0.6497, 1)

17、T1.786658710(0.7483, 0.6497, 1)T1.7865914可见此结果与未加速的规范化乘幂法公式计算结果相比，收敛速度要快得多。在已经求出主特征值l1及特征向量v1以后，可将原矩阵进行修改，使修改后的矩阵按模最大特征值是原矩阵的按模次大特征值，再用乘幂法去求按模次大特征值及特征向量，此方法称为降阶过程。为使问题简单，设A Rnn为对称矩阵。假定已求出主特征值l1及特征向量v1，记A(1) =A，构造矩阵（5.29）因A对称，则具完全特征向量系，且特征向量vi可两两相互正交，即满足（i = 2,n），于是 （i = 2, 3, , n）这表明矩阵A(2)除一个特征值l1 =

18、 0与矩阵A(1)不一样之外，其余与A(1)具相同特征值，且A(2)的按模最大特征值是l2，用A(2)代替A(1)进行乘幂迭代即可求得l2及相应的特征向量v2，而l2又为A(1)的按模次大特征值。这种做法称为降阶法。注意，降阶法实际上只可用少数几次，求矩阵的前几个按模最大特征值及特征向量。因为，每降阶一次，计算精度就会损失或降低一些。3.2.3 反幂法设ARnn可逆，则无零特征值，由有即若l为矩阵A的特征值，则必为矩阵A-1的特征值，且特征向量相同。如果A的n个特征值li (i = 1, 2, , n)为则A-1的n个特征值更为因此，若乘幂法可求A的主特征值l1，则用A-1做乘幂矩阵，由乘幂迭

19、代格式（5.30）便可求出A-1的按模最大特征值，取倒数，即为矩阵A的按模最小特征值。因此，对任取初始向量x(0)Rn，称公式（5.30）为求矩阵A按模最小特征值的反幂法。在应用公式（5.30）计算时，由于要计算A的逆矩阵A-1，一方面计算复杂、麻烦，另一方面，有时会破坏A的稀疏性，故改写（5.30）式为： （k = 0, 1, ）（5.31）类似于公式(5.25)的规范化乘幂法公式为（5.32）如果考虑到利用原点移位加速的反幂法，则记B = A - l0I，对任取初始向量x(0)Rn，（5.33）由于反幂法的主要工作量是每迭代一步都要解一个线性方程组(5.31)，且系数矩阵A (或B)是不变

20、的，故可利用矩阵的三角分解A = LU (或B = LU)，则每次迭代只需解二个三角形方程组（5.34）且当时（5.35）同时x(k+1)便为所求的特征向量，收敛速度为。反幂法的主要应用是已知矩阵的近似特征值后，求矩阵的特征向量，且收敛快，精度高，是目前求特征向量最有效的方法之一。3.4 对称矩阵的雅克比 (Jacobi) 旋转法雅克比 (Jacobi) 方法是求实对称矩阵全部特征值及对应的特征向量的方法.它也是一种迭代法，其基本思想是把对称矩阵A经一系列正交相似变换约化为一个近似对角阵，从而该对角阵的对角元就是A的近似特征值，由各个正交变换阵的乘积可得对应的特征向量。1预备知识雅克比方法涉及

21、较多的代数知识，要承认如下一些主要结论：1）若B是上（或下）三角阵或对角阵，则B的主对角元素即是B的特征值。2）若矩阵P满足PTP = I，则称P为正交矩阵。显然PT = P-1，且P1, P2, 是正交阵时，其乘积P = P1P2Pk仍为正交矩阵。3）称矩阵 (5.37）为旋转矩阵，它是在单位阵I的i行，且j行和i列、j列的四个交叉位置上分别置上cosq ，sinq，-sinq 和cosq 而成的。容易验证旋转矩阵是正交矩阵，即RT (i, j) = R-1 (i, j)，所以用它作相似变换阵时十分方便。雅克比方法就是用这种旋转矩阵对实对称阵A作一系列的旋转相似变换，从而将A约化为对角阵的。

22、用Pij(i, j)作旋转变换的几何意义是：在维空间中，以i, j轴形成的平面上，把i, j轴旋转一个角度q 。2雅克比方法先以二阶矩阵为例：旋转矩阵为其中为使A的相似矩阵B成为对角阵，只须适当选取q，使即 ，其中，a11 = a22时，取。由此q 可以确定，从而旋转矩阵P确定。A的特征值为：l1 = b11, l2 = b22关于特征向量的计算设，其中x1, x2为P的行向量，因为 PART = B APT = PTB或 （Ax1, Ax2）= （l1x1, l2x2），即Axi = liZi (i = 1, 2, 3)，所以对应于l1，l2的特征向量是考虑n阶矩阵的情况：设矩阵ARnn是对

23、称矩阵，记A0 = A，对A作一系列旋转相似变换，即（5.38）其中Ak (k = 1, 2,)仍是对称矩阵，Pk的形式也就是对任何角q，可以验证：Pk是一个正交阵，我们称它是（i, j）平面上的旋转矩阵，相应地把变换（5.38）称为旋转变换；Pk和仅在(ii)、(jj)、(ij)和(ji)上不同，PkAk-1只改变Ak-1的第p行，第q行的元素，PkAk-1Pk只改变A的第p行、q行、p列、q列的元素；Ak和Ak-1的元素仅在第P行（列）和第q行（列）不同，它们之间有如下的关系： (5.39) (5.40)我们选取Pk，使得，因此需使q 满足 (5.41)常将q 限制在下列范围内如果，当时，

24、取；当时，取。实际上不需要计算q，而直接从三角函数关系式计算sinq 和cosq，记 (5.42)则当时，有下面三角恒等式：于是cosq 始终取正值。关于sin2q 的计算有几种方法，最简单的一种是利用公式sin2q = 1 cos2q，这个方程有一个缺点，当cos2q 接近于1时，1 cos2q 的有效位数就不多了，为避免这个缺点，采用下面公式计算sinq。由于Ak的对称性，实际上只要计算Ak的上三角元素，而下三角元素由对称性获得，这样即节省了计算量，又能保证Ak是严格对称的。一般地，不能指望通过有限次旋转变换把原矩阵A化为对角阵，因为Ak-1中的零元素（在前面变换中得到的）可能在Ak中成为

25、非零元素，尽管如此，仍可以证明： 当k 时其中li是矩阵A的特征值，但没有一定的大小排列顺序。雅克比方法的优点是可以容易地计算特征向量，如果经过k次旋转变换后，迭代就停止了，即：记则因为Ak可以被看作对角阵（非对角元相当小），所以矩阵Pk的第j列就是特征值所对应的近似特征向量，并且所有特征向量都是正交规范化的。在旋转变换中可以逐步形成Pk，记P0 = I则即这就不需要保存每一次的变换矩阵Pk，若不需要计算特征向量，则（ ）式所示的那一步可以省略。算法：1从A(k-1)中找出绝对值最大元2若，则为对角阵，停若（1）令（2） 当 y = 0时， （3） （4） （5）计算特征向量 （P0 = I）转实习：用Jacobi方法求矩阵A的全部特征值及对应的特征向量。