第七章矩阵特征值与特征向量 (1).doc

第七章 矩阵特征值与特征向量的计算 

物理、力学中的振动问题，弹性结构中的稳定性问题和许多工程实际问题，最终归结为求矩阵的特征值与特征向量，即求数及非零向量，满足

其中数叫做矩阵特征值，叫做对应特征值的特征向量。
由线性代数的理论知：为特征值的充要条件是为特征方程

的根。特征方程是关于特征值的次代数方程，在复数域内共有个根。并且有
(1) 
(2)
当较大时，如果按行列式展开的方法求的零点，其工作量是巨大的，在实用中往往难以实现。因此，有必要研究有效的计算矩阵特征值的数值方法，本章将介绍迭代方法和相似变换方法。

第一节 幂法与反幂法

一、幂法

矩阵的按模最大的特征值称为主特征值，许多实际问题，往往只要求出矩阵的主特征值。幂法就是求矩阵的主特征值与对应特征向量的一种迭代方法。
设矩阵有一完全的特征向量组，即特征向量线性无关。且所对应的特征值满足主特征值占优
 (7.1)
任取一个非零向量，构造向量序列
 (7.2)
以下利用此序列求主特征值及所对应的特征向量。由于迭代序列实质上是由矩阵各次幂作用于初始向量而形成的，故称此迭代方法为幂法。
因为构成的一个基底，所以对，存在一组不全为零的数，使得

假定，按迭代公式(7.2)得

记，则上式简化为
 (7.3)
由式(7.1)得，则当时，，于是当充分大时
 (7.4)
且

因此，就是的属于主特征值的特征向量的近似值。用表示向量的第个分量，由式(7.3)得

故
 (7.5)
这说明两相邻迭代分量的比值收敛于主特征值。
观察式(7.4)，当时，的不等于零的分量，将随而无限增大，造成“溢出”；而当时，的各分量又将随而趋于零。所以，必须对上述幂法进行修正。

二、实用幂法

从非零向量出发，对迭代序列(7.2)作如下修正
 (7.6)
式中表示向量的绝对值最大的分量。则
 (7.7)
 (7.8)
事实上

由于，所以。于是

同理

例7.1 试用幂法计算实对称矩阵

的主特征值及其所对应的特征向量，要求。
解 取初始向量，应用式(7.6)进行迭代，部分迭代结果见表7.1。

表7.1

0
(1,1,1)

5
(-0.179803366,-0.,1)
-6.26409186
0.43×10-1
10
(-0.081673107,-0.30697771,1)
-6.37858254
0.13×10-1
15
(-0.055174692,-0.35765357,1)
-6.41021590
0.34×10-2
20
(-0.048413347,-0.37058404,1)
-6.41833779
0.87×10-3
25
(-0.046713405,-0.37383503,1)
-6.42038304
0.22×10-3
30
(-0.046287596,-0.37464935,1)
-6.42089555
0.55×10-4
35
(-0.046181038,-0.37485314,1)
-6.42105591
0.14×10-4
40
(-0.046154378,-0.37490412,1)
-6.42105591
0.34×10-5
45
(-0.046147708,-0.37491688,1)
-6.42106394
0.86×10-6
50
(-0.046146040,-0.37492007,1)
-6.42106594
0.21×10-6
55
(-0.046145622,-0.37492086,1)
-6.42106645
0.54×10-7
60
(-0.046145518,-0.37492106,1)
-6.42106657
0.13×10-7
61
(-0.046145509,-0.37492108,1)
-6.42106658
0.10×10-7
62
(-0.046145503,-0.37492109,1)
-6.42106659
0.77×10-8

从以上幂法的论证过程可以看出：其收敛速度依赖于，这个比值愈小，收敛速度愈快；这个比值愈接近1，收敛速度愈慢。例7.1的收敛速度就较慢，经验证。
另一方面，初始向量的选取对迭代次数也有影响。若选取不当，使得中，尽管由于迭代运算中的舍入误差导致中的，并随的增大而逐渐取得优势，使迭代收敛，但迭代次数会大增。实际计算中，可设迭代次数控制量，当仍不成立时，停机检查，根据的变化趋势重置，或运用加速技术。

三、幂法的Aitken(埃特金)加速法

由式(7.6)～(7.8)及推导过程可得


从而可知线性收敛于，即

于是

解得
 (7.9)
同样可得所对应特征向量的如下加速迭代公式
  (7.10)
可见，应用Aitken加速公式(7.9)、(7.10)之前，需应用式(7.6)计算三次。对例7.1进行计算，取，计算结果列于表7.2。与表7.1比较，要达到，幂法需
迭代次，而Aitken加速幂法只需迭代次。

表7.2

0
(1,1,1)

5
(-0.040218749,-0.35237239,1)
-6.43203285
0.34×10-1
10
(-0.046217941,-0.3752841,1)
-6.42093776
0.49×10-3
15
(-0.046144330,-0.37491683,1)
-6.42106876
0.63×10-5
20
(-0.046145497,-0.37492120,1)
-6.42106659
0.95×10-7
25
(-0.046145476,-0.37492111,1)
-6.42106663
0.37×10-7
30
(-0.046145486,-0.37492114,1)
-6.42106661
0.17×10-7
35
(-0.046145482,-0.37492113,1)
-6.42106662
0.67×10-8

若以幂法迭代至时的作为初始值，再应用Aitken加速法，计算结果列于表7.3，其收敛速度增加很快。实际上，要达到的结果，取开始迭代，直接应用幂法需迭代，应用Aitken加速幂法需迭代次。

表7.3

20
(-0.048413347,-0.37058404,1)

25
(-0.046144484,-0.37492304,1)
-6.42106820
-6.42106820
30
(-0.046145483,-0.37492113,1)
-6.42106661
0.000001591
35
(-0.046145483,-0.37492113,1)
-6.42106661
0.000000001

四、反幂法

设矩阵非奇异，则由可推得

若，则的特征值满足

对应用幂法求它的主特征值及所对应的特征向量，这一方法被称为反幂法。换言之，反幂法可用于求的按模最小的特征值与所对应的特征向量的近似值。在式(7.6)中以代替迭代，需要求，为了避免求逆阵，可通过解线性方程组进行迭代：
 (7.11)
由于需要反复求解方程组，所以不宜使用Gauss消去法。若对作分解，则每次迭代只需求解上、下两个简单的三角形方程组：
 (7.12)
当时，有以下极限


于是，当充分大时，。
例7.2 试用反幂法求矩阵

按模最小的特征值及所对应的特征向量，要求。
解 因为的各阶主子式不为零，应用式(3.10)对的分解，得

取初始向量，应用式(7.12)进行迭代，计算结果列于表7.4。只需迭代10次，就可得满足精度要求的及所对应特征向量。

表7.4

0
(1,1,1)

1
(0.43478261,1,-0.47826087)
1.21904762
0.22
2
(0.19018405,1,-0.88343558)
1.01242236
0.79
3
(0.18427518,1,-0.91241509)
1.21279579
0.20
4
(0.18314756,1,-0.91293297)
1.22933399
0.17×10-1
5
(0.18319627,1,-0.91304707)
1.22943487
0.10×10-3
6
(0.18318705,1,-0.91304168)
1.22951369
0.79×10-4
7
(0.18318795,1,-0.91304265)
1.22950865
0.50×10-5
8
(0.18318784,1,-0.91304255)
1.22950941
0.76×10-6
9
(0.18318785,1,-0.91304256)
1.22950933
0.79×10-7
10
(0.18318785,1,-0.91304256)
1.22950934
0.95×10-8

五、用反幂法对近似特征值精确化

若已知的某特征值的近似值，则
 (7.13)
于是就是矩阵的按模最小的特征值，应用反幂法

对进行分解，代入式(7.12)。则有

且当充分大时

即

例7.3 对于例7.1中的实对称矩阵，已知为的一个近似特征值，试用反幂法对精确化，并求对应的特征向量。
解 对作分解得

取代入式(7.12)，迭代结果见表7.5。

表7.5

0
(1,1,1)

1
(-0.04616784,-0.37593814,1)
-474.83801
2
(-0.04614569,-0.37492057,1)
-937.85987
3
(-0.04614548,-0.37492113,1)
-937.54585

只需迭代3次，即可得精度较高的近似特征值及特征向量


这一结果应用幂法需迭代87次。
由式(7.13)知，按模远远小于其它特征值，从而迭代过程收敛极快。另一方面，如果初步求得的全部特征值，则可应用反幂法对每个近似特征值精确化，并求得全部特征向量。

第二节Jacobi法

一、 旋转变换

Jacobi(雅可比)法是用于求实对称矩阵特征值的一种旋转变换方法。由线性代数理论知：对于阶实对称矩阵，一定存在正交矩阵，使

其中的第列向量为对应特征值的特征向量。如何构造数值方法计算正交矩阵及所有的特征值? 回顾平面解析几何中将平面二次曲线通过坐标旋转化为标准型的过程，将二次曲线用矩阵表示，恰是实对称矩阵，而旋转变换对应一个正交矩阵。一般地，记

容易验证:是一正交矩阵，由所构成的线性变换保持向量长度不变。事实上，

显然。相当于把在坐标轴平面中旋转一个角，故称旋转变换，称为旋转矩阵。旋转矩阵有以下性质：
(1) 对任意实对称矩阵，仍是实对称矩阵，且与有相同的特征值。
(2) 、，即无论用左乘或右乘矩阵，其范数保持不变。
(3) 有限个旋转矩阵的乘积仍为正交矩阵。

二、Jacobi 算法

Jacobi 算法是通过一系列旋转变换，把实对称矩阵化为对角阵，从而求得特征值及所对应的特征向量，记，对作一系列旋转变换
 (7.14)
由变到，除第行及列的交叉点元素变化外，其它值不变，计算式为
 
由到，构造旋转矩阵，所需考虑的主要问题有：选取使得
 (7.16)
选取使得。
为了保证，按式(7.15)，应满足

规定，若，则当时，取；当时，。为了式(7.15)计算方便，令

则，又据，，得
 (7.17)

三、Jacobi迭代的收敛性

据式(7.15)容易验证：

若记，，则

即经过一次旋转变换后，所得中对角线元素的平方和增加了，而非对角线元素的平方和减少了。
由式(7.16)，大于或等于的其它非对角线元素的绝对值，则

于是

由于，所以当时，的非对角线元素平方和趋于零，即(为对角阵)。
实际计算时，可给定一个误差控制小量，由控制迭代次数。若

就可视为的全部特征值。
若记，则

于是
其中为的第列向量。说明就是对于的正交特征向量的近似值。
若令(单位矩阵)，则，其迭代计算公式为
 (7.18)

四、实用Jacobi算法

以上所讨论的Jacobi算法中，每次旋转都是为了把非对角元绝对值最大者化为零作出的。对阶数较大的矩阵，为了避免每次变换都寻找最大值，可采用如下改进的Jacobi方法，称为Jacobi过关法：
令

对的非对角元素逐行扫描，如果找到某元素满足

则取，并应用式(7.15)作旋转变换，使中的。再从头开始扫描，若遇到绝对值不小于的元素，就作相应的旋转变换将其化为零，直到所有非对角元素的绝对值全部都小于。
再令，重复上述步骤，直至所有非对角元素的绝对值全部都小于。这样经过一系列“关”(其中)，就可得到与相似的近似对角阵。
例7.4 对例7.1中的实对称矩阵，运用Jacobi过关法求其全部特征值及所对应的标准正交的特征向量(取)。
解 记

令，按Jacobi过关法计算，计算过程如表7.6。

表7.6

1

2

3

4

5

6

7

由表7.6 可得的特征值，实际上满足精度要求的近似特征值为

所对应的特征向量的近似值分别为中的列向量。即

Jacobi过关法算法流程见图7.1。

第三节 QR方法

QR方法是目前计算中小规模矩阵的全部特征值和特征向量的最有效方法之一。QR方法是利用正交相似变换将一个给定矩阵逐步约化为上三角阵或拟上三角阵的一种迭代方法。

一、Householder变换

Householder变换(或矩阵)是指具有如下形式的初等正交矩阵

其中是满足的维实向量。
Householder矩阵具有如下一些简单而又十分重要的性质：
(1) 对称性：；
(2) 正交性：；
(3) 对合性：；
(4) 反射性：对任意的，是关于的垂直超平面的镜像反射。
设是任意给定的维非零实向量，可以通过求得Householder矩阵，使得的若干个指定分量变为零。令

则容易验证，在选取适当的条件下，其所对应的Householder矩阵满足
(7.19)
其中，且有。
记，根据Householder矩阵的正交性有

即，从而由可得。取正或取负，对结果的影响只是符号变了。
令，则，从而
。
例7.5 设，求Householder矩阵满足。
解 取，则，从而有

并且。

二、QR分解

利用Householder变换，对于矩阵可以得到正交阵，使得，其中为上三角阵。下面给出构造的方法。
记，它的第一列记为，则有Householder矩阵，使得
，令，可知其第一列对角线以下方元素全为。一般地，设

其中为阶方阵，为阶方阵，设其第一列为，则可求得Householder矩阵，有

根据构造阶的Householder矩阵

就有

其中为阶方阵，为阶方阵。这样经步运算可得

则知为上三角阵，为正交阵，令即可，并且可以验证矩阵的这种分解是唯一的。
例7.6 求矩阵的QR分解。
解 记，，取，有
则
，
记，，取，有，则
，
正交阵

从而矩阵的QR分解为

三、QR算法

设，对进行QR分解，令矩阵，重复上述过程，对于已获得的进行QR分解，并且有，继续以上递推法则，就会获得矩阵序列，这种方法称为QR算法。只要为非奇异的，则由QR算法可以完全确定矩阵序列。
定义7.1 当时矩阵序列的对角线元素均收敛，且严格下三角元素收敛到零，则称基本收敛到上三角阵。
基本收敛的概念并未指出严格上三角的元素是否收敛，但是对求的特征值而言，基本收敛已经足够了。
定理7.1 设，其特征值满足

为所对应的特征向量，以为列的方阵记为。如果有，其中为单位下三角阵，为上三角阵，则由QR算法产生的序列基本收敛到上三角阵，且
，
例7.7 试用QR算法求解矩阵特征值的近似值。
解 方程组的特征值的精确解为3，-1和1。由上例，则可得

不断重复上述步骤，将矩阵的对角线元素在下表7.7中列出：
表7.7

2
3
4
5
6
7
8
9
10

3.6669
3.2328
3.0752
3.0245
3.0081
3.0026
3.0008
3.0002
3.0000

-1.1961
-1.0798
-1.0256
-1.0088
-1.0029
-1.0015
-1.0003
-1.0007
-1.0001

0.5295
0.8470
0.9500
0.9836
0.9942
0.9983
0.9990
1.0000
0.9996

迭代10次后，所得近似值与精确解已经非常接近。

第四节 应用举例

一、弹簧——质量体系(系统)的基振频率和振型

如图7.2，振动系统由三段弹簧和与其相连的三个物体构成。设表示第段弹的基振频率和振型。

图7.2 弹簧——质量振动系统
设为第个物体的位移，它是以自身的平衡位置为零点。据Newton运动定律：

由于位移与振幅、频率及相角有如下关系：

同时考虑到加速度等于位移的二阶导数，即

代入方程组(7.19)，得

写成矩阵的形式为

简记为

式中，于是将求基振频率转化为求如下特征值问题：

其中

取初始向量，应用反幂法迭代4次，得最小特征值及对应的特征向量的近似值：
，
于是该振动系统的基振频率为

将向量关于矩阵规范化，即满足

则有振型


二、 结构的振动

工程中的振动问题十分普遍、复杂。例如船舶的螺旋浆转轴常被视为一个带有一桨叶盘的悬臂梁的分布参数系统，这个系统有无穷多个自由度，需求解偏微分方程。为简化问题起见，悬臂梁可简化为有限个质点，由无质量的轴段相联接，此时系统为离散参数系统，见图7.3(a)。试计算悬臂梁的自振频率和振型。

若悬臂梁全长15米，弯曲刚度为，质量、、，相互之间的间距为。利用材料力学知识，不同荷载对弯矩的影响是独立的，所以可分别求出单个荷载所产生的挠曲线，然后叠加得静态弯曲方程

式中系数的单位是。简记为

式中称为柔度矩阵。在没有受到其它外力的条件下，作用于一个质量上使它加速的力和作用于梁上同一点的相应力之和必为零，于是有

其中为加速度列向量，质量矩阵为

整理得到

设悬臂梁的振动满足，则有

式中是质点最大位移的列向量，显然是由矩阵的特征值所决定。为了计算方便，令，，，于是计算梁的自振频率和振型问题归结为特征值问题。其中
，
因为，所以的元素的单位是。
应用Jacobi方法迭代7次，得特征值及所对应的特征向量的近似值：

利用关系式，得悬臂梁的自振频率为

它们所对应的特征向量即为所求振型：

小 结

矩阵的特征值与特征向量的计算，是数值计算方法的重要内容之一，内容丰富，振动问题、稳定性问题等许多工程技术问题，往往归结为求矩阵的特征值与特征向量问题。本章重点介绍了三类常用的方法——幂法、Jacobi方法与QR方法。
幂法与反幂法属于迭代方法。幂法主要用来求按模最大的特征值和对应的特征向量，它的主要特点是算法简单，便于在计算机上实现，但收敛速度较慢；反幂法可用来求绝对值最小的特征值及特征向量，如果已初步估计出矩阵A的某特征值的近似值，可应用反幂法对此近似值进行精确化，且迭代速度较快。
 Jacobi方法是把实对称矩阵经一系列正交相似变换化成对角阵，从而求得矩阵的全部特征值及所对应的特征向量，是求矩阵特征值的常用方法，具有收敛快、计算精度较稳定等优点，且所求得的特征向量具有很好的正交性。但Jacobi方法的运算量较大，在进行正交相似变换的过程中，不能保持稀疏矩阵的稀疏性质，因而适宜于求低阶矩阵的全部特征值。
QR方法是求全部特征值的方法，它利用正交相似变换，如Householder变换，对矩阵做化简和QR分解，来求矩阵的特征值。QR方法与Householder方法相结合，具有收敛快，精度高，在中小型稠密矩阵的特征值问题计算中，目前它还是十分有效的方法。

思 考 题

1. 幂法是用来求矩阵的哪些特征值与特征向量?
2.反幂法的思想是什么?可用反幂法求矩阵的哪些特征值?
3. Jacobi法的基本思想是什么?它针对什么样的矩阵?
4. Jacobi过关法的主要优点是什么?
5. 利用QR方法求矩阵特征值时，对矩阵有无具体要求？

习 题 七

1. 用幂法求下列矩阵按模最大的特征值及对应的特征向量。
(1) (2)
2. 用反幂法求第1题中矩阵的按模最小的特征值及对应的特征向量。
3. 已知矩阵

求最接近12的特征值与相应的特征向量。
4. 设是实对称矩阵，为的特征值。
(1) 试证：，与向量的Rayleigh商满足

(2) 对于非零向量，令

此迭代公式称为Rayleigh商加速幂法。试利用该加速公式计算矩阵

的主特征值及所对应的特征向量。
5. 用Jacobi方法求下列对称矩阵
(1) (2)
的全部特征值及所对应的特征向量。
6. 求下列矩阵的QR分解。
(1) (2)
7. 试用QR算法求解第6题中矩阵的全部特征值。

数值实验七

试验目的与要求

加深对求矩阵特征值及特征向量的幂法、Jacobi法理论的理解，掌握幂法、反幂法、Jacobi法的主要步骤和应用范围。

试验内容

1. 应用幂法求下列矩阵的主特征值及所对应的特征向量

(1) 讨论初值对迭代次数的影响，即选取多组进行计算，分别打印迭代次数。
(2) 讨论特征值扰动问题，当的某一元素发生微小变化，例如当、、
时，讨论特征向量的变化。
(3) 比较幂法、Aitken加速幂法、Rayleigh商加速法的迭代次数。
2. 用幂法计算矩阵

的主特征值及对应的特征向量。
(1) 打印、，观察它们的变化规律。
(2) 当矩阵的主特征值是一对相反实数时，随着的增大，及作规律性的摆动，此时

所对应的特征向量为

3. 分别应用Jacobi法、Jacobi过关法及QR方法求矩阵

的全部特征值及所对应的特征向量。并对不同的，比较两种方法所需的变换次数。

实验结果与分析

通过对同一问题运用不同的方法求解，对此不同方法计算机编程的难易程度及所需的运算次数，总结不同方法的优缺点，进一步总结不同方法的特点、主要结论、应用范围，写出实验报告。