(内容提要)-6--特征值数值解.doc

第九章 特征值与特征向量的计算 273 矩阵特征值与特征向量的数值计算矩阵的特征值问题有着广泛的应用背景，例如微分方程的刚性比和数值方法的稳定性；动力系统和结构系统中的振动问题；电力系统的静态稳定分析；因子分析模型中的因子载荷、共同度和特殊方差的估计，实质上是矩阵的特征值问题。设阶实矩阵的特征值问题是求和非零向量，使得成立。其中为特征值(Characteristic Value)，为的对应于的特征向量(Characteristic Vector)。一般地，代数方法只适用阶数较低的矩阵，当矩阵阶数很高时，用代数方法求矩阵的特征值和特征向量是极为困难的，本章介绍几种比较有效的数值方法，幂法和反幂法、方法及雅可比方法，我们仅对实矩阵进行讨论。1 幂法和反幂法一、 幂法(The Power Method)幂法主要用于求矩阵的按模最大的特征值与相应的特征向量的数值方法。它是通过迭代产生向量序列，由此计算特征值和特征向量。设阶实矩阵的特征值满足且与对应的特征向量线性无关。幂法的基本思想是给定初始向量，由迭代公式 （9-1）产生向量序列：上述向量称为迭代向量。为简便起见，不妨设。因为线性无关，故必存在个不全为零的数，使得。于是由式（9-1）得 （9-2）设，由得于是故只要k充分大，就有 （9-3）因此，可以近似作为与相应的特征向量。下面我们通过特征向量来计算特征值。用表示的第个分量，由于所以 （9-4）上式这种由已知非零向量及矩阵的乘幂构造向量序列用来计算矩阵按模最大的特征值（按（9-4）式）与对应的特征向量（按（9-3）式）的方法称为幂法。值得注意的是，因为，当时，的各分量要趋于无穷，当时，各分量要趋于零，这在计算机里称为“上溢”和“下溢”。为了克服这个缺点，就需要将迭代向量加以规范化。设有以向量，将其规范化得到向量，其中表示向量的绝对值最大的分量。任取一初始向量，构造向量序列 （9-5）由可知 这说明规范法向量序列收敛到所对应的特征向量。同理可知，综上可知：按上述方法构造的向量序列，。三、反幂法(Inverse Power Method)反幂法是计算矩阵按模最小的特征值及特征向量的方法，也可以用来计算对应于一个给定近似特征值的特征向量。设为阶非奇异矩阵，为的特征值与相应的特征向量，即由于所以的特征值是的特征值的倒数，而相应的特征向量不变。如果的特征值的次序为，则的特征值为，因此，若对矩阵用幂法，即可计算出的按模最大的特征值，从而求得的按模最小的特征值。这就是反幂法的基本思想。因为的计算比较麻烦，而且往往不能保持矩阵的一些好性质（如稀疏性），因此，反幂法在实际运算时以求解方程组 （9-7）代替幂法迭代求得，每迭代一次要解一个线性方程组（9-7）。由于矩阵在迭代过程中不变，故可对先进行三角分解，这样，每次迭代只要解两个三角形方程组。反幂法计算的主要步骤如下：1对进行三角分解2求整数，使得,计算 （9-8）3 解方程组 （9-9）反幂法与原点移位法结合还有更好的应用。设已求得的一个特征值的近似值为，因为接近，并且于是的特征值为，故是矩阵的按模最小的特征值，且由上式可知，比值较小。因此，对用反幂法求一般收敛很快，通常只要经过二、三次迭代就能达到较高的精度。例3 用反幂法求矩阵的最接近的特征值，并求相应的特征向量。解 对矩阵 按公式（9-5）计算的数值结果见下表由表可知 13.22018，相应的特征向量为（1 -0.23510548 -0.17162117）T 。0 1 2 3 4 51 -2.4545450 -4.5908214 -4.54094172 -4.54175138 -4.541738511 0.66666669 1.07818937 1.06764054 1.06779003 1.067787651 0.48484850 0.78750467 0.77934009 0.7746037 0.779458521 1 1 1 1 11 -0.27169496 0.23453777 -0.23511435 -0.23510535 -0.235105481 -0.19753087 -0.17130533 -0.17162521 -0.17162110 -0.171621171 -2.4545450 -4.5908214 -4.5409417 -4.54175138 -4.54173851 -13.40741 -13.21753 -13.22022 -13.22018 -13.22018 2 Jacobi方法Jacobi方法是计算实对称矩阵的全部特征值及相应的特征向量的一种矩阵变换方法，实对称矩阵有如下的性质： （1）的特征值均为实数；（2）任意实对称矩阵可以通过正交相似变换化成对角型，即存在正交矩阵，使得其中是的特征值，中各列即为相应的特征向量。（3）在正交相似变换下，矩阵元素的平方和不变。设,为正交矩阵，记,则Jacobi方法的基本思想是：构造特殊的正交矩阵序列，通过正交变换，使的非零的非对角元素逐次化成零，并且使得非对角元素的平方和减小。从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。一、 矩阵的旋转变换(Rotative Fransformation)Jacobi方法的关键是如何构造正交矩阵？先看一个简单的例子，设 是二阶实对称矩阵，即，其特征值为。我们可以通过对坐标平面作旋转相似变换，使 其中是平面旋转矩阵，为旋转角。显然，是正交矩阵，问题是如何确定中旋转角的大小呢？记则为使对角化，要求非对角元素，即因此，当时， 从而 （9-10）当时，可选。为使旋转角度尽可能小，通常限制。由此可看出，对二阶矩阵，只做依次正交变换，选择适当角，就可将对角化。将这种思想推广到维情况，设为阶实对称矩阵，考虑中平面旋转变换矩阵 （9-11）它是将阶单位阵中对角线上第和第个元换成，非对角元和分别换或和而得到的，容易验证，是正交矩阵，若记则有 （9-12）如果，取使得则有，这样，就得到一个使中非零的非对角元素变成零的正交相似变换。对重复上述过程，可得，如此继续下去，得到一个矩阵序列。可以证明由上述方法构造的旋转矩阵对变换后，就会使非对角元的平方和严格单调递减，而对角元的平方和单调递增。二、 Jacobi方法通过一系列旋转相似变换将变成，求得的全部特征值与特征向量的方法称为Jacobi方法。如果在对作相似变换的过程中，每一步都选绝对值最大的非对角元素，以此确定旋转矩阵，这种方法称为经典Jacobi方法。其计算过程如下：（1）令（2）求整数，使得（3）计算旋转矩阵（4）计算置（5）计算，（令）（6）若，则为特征值，的各列为相应的特征向量；否则，返回2，重复上述过程。一般地，Jacobi方法不能在有限步内将化成对角阵，但有以下收敛性结果。定理9.1 设为阶实对称阵，对用经典Jacobi方法得到序列，其中，则即经典Jacobi法收敛。定理9.1表明，经典Jacobi方法是收敛的，并且进一步分析还可以得出Jacobi方法收敛较快。另外，（1）Jacobi方法对舍入误差有较强的数值稳定性，因而求得的解精度较高。（2）Jacobi方法的不足之处是运算量大，且当是稀疏矩阵时，不能保持稀疏性，因此Jacobi方法是求中小型稠密实对称矩阵的全部特征值与特征向量的较好方法。例4 用Jacobi方法求对称矩阵的全部特征值。解 记(1) 首先取， (2) 取,与(1)类似得到： 重复上面过程，如此继续下去，可得 故得到的特征值： ，3 方法方法是计算一般中小型矩阵全部特征值与特征向量的最有效的方法之一。它是以矩阵正交三角分解为基础的一种变换方法。这里仅讨论实矩阵，并假定矩阵非奇异。一、矩阵的分解定理9.2 任一非奇异实矩阵都可以分解成一个正交矩阵和一个上三角矩阵乘积，而且当的对角元素为正时，分解是惟一的。对矩阵作分解的方法有多种，下面以Schmit正交化方法为例证明。证明 设为阶非奇异实矩阵，将矩阵写成分块形式 其中，因为A非奇异，所以线性无关。取， ，显然，取，则。一般地，取 （9-13）则向量组正交，且。式（9-13）可改写成于是这就是用Schmit正交化方法对矩阵进行分解的过程。二、基本方法基本方法的思想是利用矩阵的分解，通过迭代格式 （9-14）将化成相似的上三角阵（或分块上三角阵），从而求出矩阵的全部特征值与特征向量。具体计算步骤为：令，对作分解令，作分解 重复上述过程，得迭代公式 （9-15）这样可以得到矩阵序列，并且矩阵是相似矩阵。事实上：由可得，于是即与相似。同理可得，故它们有相同的特征值。有上述方法构造正交相似矩阵的方法称为基本方法。在一定条件下，矩阵序列收敛于一个上三角阵或分块上三角阵，且对角块为11矩阵或22矩阵，其中11矩阵对角元素为的实特征值，22矩阵含有的一对复特征值。特别地，如果是实对称阵，则收敛于对角矩阵。另外，当充分大时，的主对角元（或主对角子块的特征值）就可以作为的特征值的近似。例5 用基本方法求矩阵的特征值。解 令，对作分解 然后，求得，作分解，一直下去，得到所以，的一个特征值为4，一个特征值为，还有一对复特征值是方程的根，即，得。事实上，矩阵特征方程为，其特征值为，。基本方法每次迭代都需作一次分解与矩阵乘法，计算量大，而且收敛速度慢。因此实际计算时，总是先用一系列相似变换将化成拟上三角矩阵（称为上Hessenberg矩阵），然后对此矩阵用基本方法。这样可减少运算量。下面介绍一种化为相似的拟上三角阵的方法：Householder变换。4 Mathematica应用实例例1 方阵由前25个相邻整数构成，求的特征多项式、特征方程、特征值与特征向量。解 Mathematica语句：A=PartitionRange25, 5;MatrixForm%CharacteristicPolynomialA,CharacteristicPolynomialA,0EigenvaluesNA/ChopEigenvectorsA/N;MatrixForm%运行结果： 特征多项式 特征方程 特征根 特征向量例2 给出阶Hilbert矩阵特征值的近似值。解 运行结果：附：Mathematica语句n=6;hilbert=Table1/(i+j-1), i,1,n, j,1,n;MatrixForm% NEigenvalueshibert,4例3 用幂法的规范运算求矩阵的按模最大的特征值及对应的特征向量。解 计算公式：，取Mathematica语句：A = 1, -1, 2, -6;MatrixForm%x = -0.5, 1;Doy = A.x; Printk, , y, , x; x = y/MaxAbsy, k, 1, 10运行结果： 特征值的近似值 特征向量的近似值1 -1.5, -7. -0.5, 12 0.785714, 5.57143 -0.214286, -1.3 -0.858974, -5.71795 0.141026, 1.4 0.849776, 5.69955 -0.150224, -1.5 -0.850905, -5.70181 0.149095, 1.6 0.850766, 5.70153 -0.149234, -1.7 -0.850783, -5.70157 0.149217, 1.8 0.850781, 5.70156 -0.149219, -1.9 -0.850781, -5.70156 0.149219, 1.10 -0.850781, -5.70156 0.149219, 1.，相应的特征向量为 0.149219，1 。例 4 用反幂法的规范运算求矩阵的按模最小的特征值及对应的特征向量。解 计算公式：，取Mathematica语句： A = 3,2,1, -1,8,2, 1,4,16;xb = -1,1,0.5;Doxa=LinearSolveA,xb; Printk, , xb, ,xa; xb=xa/MaxAbsxa, k,1,20运行结果： 1 -1, 1, 0.5 -0.390625, 0.0664062, 0.03906252 -1., 0.17, 0.1 -0.325937, -0.0278906, 0.03359373 -1., -0.0855705, 0.103068 -0.307334, -0.0592273, 0.0404574 -1., -0.192713, 0.131638 -0.299819, -0.0728619, 0.04518165 -1., -0.24302, 0.150696 -0.29635, -0.0793666, 0.0477826 -1., -0.267814, 0.161235 -0.294651, -0.0825935, 0.04914137 -1., -0.280309, 0.166778 -0.293798, -0.0842239, 0.0498428 -1., -0.286673, 0.169647 -0.293364, -0.0850551, 0.0502029 -1., -0.28993, 0.171125 -0.293142, -0.0854807, 0.050386910 -1., -0.291602, 0.171886 -0.293028, -0.0856992, 0.050481911 -1., -0.292461, 0.172277 -0.292969, -0.0858115, 0.050530712 -1., -0.292903, 0.172478 -0.292939, -0.0858692, 0.0505559 18 -1., -0.293362, 0.172687 -0.292908, -0.0859292, 0.05058219 -1., -0.293366, 0.172689 -0.