第4章线性方程组和矩阵特征值迭代解法

1、 第4章 线性方程组和矩阵特征值的迭代解法 线性代数计算方法中的迭代解法（即迭代法）是一类重要方法。其基本思想是构造适当的矩阵序列或向量序列，使其逐步逼近所求问题的精确解，故又称矩阵迭代方法。在求解阶数较高且零系数较多的大型稀疏线性代数方程组时，迭代法是很有效的。矩阵特征值问题的求解通常也要用迭代法。本章着重介绍求解线性代数方程组常用的简单迭代法及其收敛条件，并对计算矩阵特征值问题的雅可比方法和QR方法作一些介绍。 4.1 线性代数方程组的迭代解法线性方程组(3.1)的迭代解法其基本思想与一元非线性方程的迭代解法类似，即构造适当的迭代公式，任选一个初始向量进行迭代计算，使生成的向量序列，,收敛

2、于方程组的精确解。4.1.1 简单迭代法的一般形式设方程组(3.1)的系数矩阵非奇异，把它化为等价的方程组 (4.1)其中 按(4.1)构造迭代公式 (4.2)其中。任取初始向量，用(4.2)逐次计算近似解向量这种方法称为简单迭代法，称(4.2)为简单迭代公式，为迭代矩阵。公式的分量形式是即 (4.3)如果的各分量存在极限 (4.4)则称向量序列收敛于向量，并记为 (4.5)这时，称简单迭代法(4.2)是收敛的，否则就是发散的。 当(4.2) 收敛时，容易看出满足方程组(4.1)，从而必定是原方程组(3.1)的解。实际计算时，一般用 (4.6)来控制迭代结束，取为满足要求的近似解，其中e是指定

3、的精度要求。这样做的理论根据将在后面(4.23)中指出。例如下面左边的方程组可改写为右边的同解方程组 取，可构造迭代公式 ， 这是分量形式。若写成矩阵形式就是(4.2),其中就是迭代矩阵。容易计算出，设精度要求为，计算结果如表4.1所示。表4.1 k01234567045.55.14.954.995.0055.001032.21.91.982.012.0021.999由于,故求得方程组的解为上例中的迭代公式是收敛的，当时，和的值越来越靠近准确解5和2。对给定的方程组，可以构造各种迭代公式，下面介绍两种常用的简单迭代公式。 4.1.2 雅可比迭代法 设方程组(3.1)的系数矩阵非奇异，其主对角元

4、素，并且绝对值相对来说比较大，我们从方程组的第个方程 中解出，得到等价的方程组 (4.7)即 (4.8)其中 由(4.7)构造迭代公式 (4.9)其矩阵形式为 (4.10)我们称(4.9)或(4.10)为雅可比(Jacobi)迭代公式，为雅可比迭代矩阵。任取初始向量，按上述迭代公式逐次计算，这就是雅可比迭代法，这是一种简单迭代法。 例4.1 用雅可比迭代法求解方程组 精度要求为e=0.005，用5位有效数字计算。 解 将方程组写成等价的方程组 构造雅可比迭代公式 取初始向量进行迭代，计算结果如表4.2表4.2 k 0 12 3 4 5 0 0.90.98 0.994 0.999 2 0.999

5、 80 0 0.70.96 0.99 0.998 0.999 64 0 0.80.94 0.992 0.998 0.999 60 由于 所以满足精度要求，即得 图4.1 雅可比迭代法的计算框图 与精确解相比较，上述计算解都具有3位有效数字。简单迭代法在计算机上实现时，注意到计算只用到，无须保存及以前的计算结果，故一般只用两个一维数组和分别存放相继两次的迭代向量和，当进行下一步迭代时，将的值存入取代旧的结果，将新的计算值存入，依次循环进行计算。为了防止迭代发散时陷入死循环，可以事先设一个最大迭代次数N0，若迭代N0次仍达不到精度要求，则输出表示迭代失败的标志，并停止计算。 雅可比迭代法的计算框图

6、见图4.1。对于一般形式的简单迭代法(4.3)，计算框图类似，只需将图4.1中“输入”改为“输入”，并将计算的部分改为 (4.11)即可。当是大型稀疏矩阵时，还可以只存储迭代矩阵的非零元素，以便节约存贮单元及提高运算速度。 4.1.3 高斯-赛德尔迭代法 在雅可比迭代公式(4.9)中，计算的第i(>1)个分量时，所用的全是的各个分量，对新算出的分量并没有利用。若迭代收敛，我们设想把算出后立即代替用于后面分量的计算，当算出后立即代替用于后面分量的计算，期望这样会收敛得快些。根据这种思想得到高斯-赛德尔(Seidel)迭代公式 (4.12)若令 ，(4.12)即为 (4.12)其矩阵形式是

7、(4.13)其中 比如例4.1中根据迭代公式容易看出B以及L、U为， (4.14)容易看出雅可比迭代矩阵B与L、U的关系为 B=L+U (4.15) 例 4.2 用高斯-赛德尔迭代法求解例4.1的方程组。 解 按例4.1的方式化方程组为等价的方程组，构造高斯-赛德尔迭代公式为 表4.3 k 0 12 3 0 0.9 0.997 6 0.999 94 0 0.88 0.997 12 0.999 93 0 0.9760.999 42 0.999 99 取初始向量=0=0,0,0T，迭代结果如表4.3。由于 故即为满足精度的近似解，得 =0.999 94 , =0.999 93, =0.999 99

8、 图4.2 高斯-赛德尔迭代法的计算框图高斯-赛德尔迭代法(简称GS迭代法)中，由于算出后立即代替用于后面的计算，所以在计算机上只需用一个一维数组存放迭代结果，每次用新值取代旧值，即 (4.16)但为了判断是否达到精度，需计算，为此可设它为，并用一个变量临时寄存的旧值，用如下步骤计算和：(1)(2) 对作 。按(4.16)计算新的。高斯-赛德尔迭代法的计算框图见图4.2。 需要注意的是，尽管在例4.2中用GS迭代法比用雅可比迭代法收敛得快一些，但是一般情况并非总如此理想。对任一个给定的方程组，在两种方法都收敛的情况下，可能GS迭代法收敛得快，也可能雅可比迭代法收敛得快。在某些情况下，可能一种方

9、法收敛而另一种方法不收敛，或者两者都不收敛。例如对方程组直接构造雅可比迭代公式或GS迭代公式，都是发散的；但如果交换两个方程的次序后再构造相应的迭代公式，则两种方法都收敛(请读者自行验算之)。为此有必要了解迭代法的收敛条件。雅可比迭代法是一种简单迭代法，它的公式与简单迭代法的一般形式(4.2)、(4.3)一致。GS迭代公式看起来却与一般形式不同，它是否属于简单迭代法呢?在(4.13)中移项可得 两端左乘得 与(4.2)对照可知，高斯-赛德尔迭代法也是一种简单迭代法，迭代矩阵为 (4.17)称为高斯-赛德尔迭代矩阵(简称GS迭代矩阵)。 如由(4.14)知，例4.2的GS迭代矩阵为 (4.18)

10、下一节将简要地介绍简单迭代法的收敛条件，理论证明从略。 4.2 迭代法的收敛性 定理 4.1 简单迭代法(4.2)收敛的充分必要条件是迭代矩阵M的谱半径。其中M的谱半径是指 (4.19)是矩阵M的特征值。(证明略,见5)例如在例4.1中，雅可比迭代矩阵B（见(4.14) ）的特征值为0.2,0.1+0.1i,0.10.1i, ,所以雅可比迭代是收敛的。易见例4.2的高斯-赛德尔迭代矩阵（见(4.18)）具有谱半径，所以GS迭代法也是收敛的。 在具体问题中，矩阵的谱半径一般不易计算，所以通常用其他较强的充分条件来判断迭代收敛，为此先介绍范数的概念。 4.2.1 向量和矩阵的范数 定义1 如果对任

11、意实n维向量(或n阶矩阵) ，都有确定的实数与之对应，记为，并且满足以下三条性质： (1)非负性，即，等号当且仅当=0时成立； (2)齐次性，即为任意实数； (3)三角不等式，即和是任意n维向量(或n阶矩阵)。则称实数为向量(或矩阵) 的范数。 向量范数是向量长度的推广。设为n维向量，常用的向量范数有 分别称为向量的1范数、2范数、无穷范数，可以验证它们都具有定义1的三条性质。 例如设 x=1,2,3T，则容易算出 设为n阶矩阵，常用的矩阵范数有 它们分别称为矩阵A的1范数、无穷范数。例如设 则有 矩阵的1范数和无穷范数都满足相容性，即 (1) 为任意n阶矩阵。 (2) 为任意n维向量，a为1

12、或¥。还要指出，对于n维向量(或n阶矩阵)所定义的任意两种范数是等价的。即存在数m、M，使对一切n维向量(或n阶矩阵)x有。例如对于n维向量x有等等。显然，向量序列收敛于向量的概念等价于的各分量收敛于0，即这相当于收敛于0，由于向量范数有等价性，因此可以用 (4.20)来表示收敛于，这里可以是任意向量范数。在迭代法中也就可以用任一种范数来表示近似解向量的误差。当然较方便的还是使用无穷范数。 4.2.2 迭代法收敛的充分条件 矩阵的谱半径和矩阵范数有如下关系: 定理 4.2 设为n阶矩阵，则。 证 设为的特征值，是对应于的特征向量，则。上式两端取范数，利用范数的有关性质可得 ，由于，故

13、，从而结论成立。 根据定理4.1和4.2，我们就可以用来作为判断迭代法(4.2)收敛的充分条件。 定理4.3 若迭代矩阵的范数，则迭代公式(4.2)收敛，且有误差估计式 (4.21)及 (4.22)对于例4.1的雅可比迭代矩阵，容易计算，或，故可根据定理4.3判定相应的迭代公式收敛。显然这比用谱半径来判断方便得多。那么，能否根据来判定迭代公式(4.2)发散呢？请自己作出结论（见习题）。定理4.3提供了控制迭代结束的方法。在(4.21)中用无穷范数便是 (4.23)所以，为使迭代解与真解足够接近，可事先指定适当小的正数e，用(4.6)式控制迭代结束。下面不加证明地列出常用的几个判断简单迭代法收敛

14、的充分条件：设对于方程组构造简单迭代公式，则 (1)若迭代矩阵满足或，则简单迭代公式收敛。 (2)若系数矩阵或其转置是严格对角占优矩阵，则雅可比迭代公式和高斯-赛德尔迭代公式都收敛。 (3)若系数矩阵对称正定，则高斯-赛德尔迭代公式收敛。（证略） 例4.1的系数矩阵是严格对角占优的，所以用雅可比迭代法和GS迭代法都收敛。 需要注意，当上述的充分条件不满足时，迭代公式也是有可能收敛的，这时可用定理4.1来判断，当时迭代公式(4.2)收敛，当时迭代发散。实际计算中，往往事先调整方程组，如改变方程或未知量的次序，对若干方程加以组合等，使得尽可能满足收敛条件。比如要构造雅可比迭代公式或GS迭代公式，可

15、以事先调整方程组使系数矩阵的主对角元绝对值尽可能相对地大，这样构造出的迭代公式就很可能收敛。 4.3 矩阵特征值问题的计算方法 许多工程技术实际问题中需要计算矩阵的特征值和特征向量。当矩阵A的阶数较高时，若首先构造特征多项式再求其根，不仅实现起来比较困难，而且往往舍入误差的影响很大。因此通常采用其它行之有效的数值方法。 4.3.1 雅可比方法 设是n阶实对称矩阵，则一定存在正交矩阵，使 (4.24)其中是对角矩阵，其对角元是的全部特征值，正交矩阵的第j列就是对应于lj的特征向量。 雅可比(Jacobi)方法的基本思想是通过一系列特殊的正交相似变换，逐步构造上述正交矩阵，使化为对角矩阵，从而求得

16、的特征值和特征向量。 设的一对非对角元素，令 (4.25)只有四个元素与单位矩阵不同，容易验证是正交矩阵，称它为平面旋转矩阵。用对作一次正交相似变换后得 (4.26)通过计算知，的第i、j两行和第i、j两列发生了变化，其它元素与相同。如以三阶矩阵为例，设，则i=1,j=3,可得= 一般n阶情形有 (4.27)如果取使满足 (4.28)即若，则当时取,当时取；否则取满足 (4.29)于是有，即把的一对非对角元素和化成了0。 上述变换可以逐次进行，每次选一对绝对值最大(不为零)的非对角元素，用这种变换把它们化为零，这样得到一系列实对称矩阵。一般情况下不可能只用有限次变换就把化成对角形，因为每次变换

17、都有可能使先前已经化为零的非对角元变为非零。但是可以证明，非对角元的平方和会逐次减小而趋于零，对角元的平方和则逐次增加，矩阵序列收敛于一个对角矩阵(证明略，可参阅8)。 在计算机上计算时，一般给定精度要求e，当矩阵的非对角元素的最大绝对值小于e时，把最终矩阵的对角元作为的近似特征值。假设一共经过m次正交相似变换，已满足精度要求，各次所用的平面旋转矩阵设为，则 其中已近似等于对角矩阵，其对角元就是的近似特征值，令 就有 所以的第j列就是的特征值的近似特征向量。 综上所述得雅可比方法的计算步骤为: (1)输入，精度，令 (单位矩阵)； (2)确定中绝对值最大的非对角元素；若 ，则转(6)； (3)

18、计算使满足(4.28）、(4.29）； (4)按式(4.27)、(4.25)计算、的元素； ； (5)；返(2) (6)输出的对角元，输出。 例 4.3 用雅可比方法求矩阵的特征值和特征向量。取四位有效数字计算。 解 绝对值最大的非对角元为，即i=1，j=2，因，故取，,按(4.25)、(4.26)得 ，故由(4.29)得 经计算得 ， 如此继续计算，7次旋转相似变换后得 ，故求得的近似特征值为2.415，-0.414 3，1.000；V的各列即为相应的特征向量。 事实上的精确特征值为，可见上述结果是令人满意的。 雅可比方法的优点是可以同时求出实对称矩阵的特征值和特征向量，并且算法是稳定的，结

19、果的精度一般都比较高；缺点是当为稀疏矩阵时，计算过程中不能保持原有的零元素分布特征，因此这种方法一般用于阶数不很高的“稠密”对称矩阵情形。此外，每次挑选绝对值最大的非对角元素比较费时，这方面有一些改进的措施，具体方法可参看有关书籍。 4.3.2 QR方法简介 QR方法是求解一般矩阵全部特征值的最有效的方法之一，它基于下述定理： 定理 4.4 任意n阶矩阵总可以分解成一个正交矩阵Q和一个上三角矩阵R的乘积，即 (4.30) 称(4.30)为n阶矩阵的QR分解。例4.4 设矩阵 ，判断是否具有QR分解式。解 显然是上三角阵，是方阵，且故为正交阵。经验算得，所以这就是的一个QR分解式。一般，矩阵的Q

20、R分解式是不唯一的，但总可通过一定的算法找出一个QR分解式来。计算矩阵特征值的QR方法，就是利用矩阵的QR分解，构造一个矩阵迭代序列，使与保持相似且基本收敛于上三角形或分块上三角形，从而容易求得全部特征值。QR方法的计算公式是: (4.31)即先对作QR分解，再令 ；然后对作分解，再令；。显然这是一种迭代法，由于 故与相似，依此可知每一个 (k=2,3,)都与相似。 可以证明，若的特征值是两两不同的实数，则在一定条件下，当时，的主对角元下方各元素收敛于零，而主对角元素收敛于的特征值。实际计算中，当的主对角元下方各元素的绝对值足够小时停止迭代，将的对角元素作为近似特征值。当实矩阵有复的特征值时，

21、在一定条件下收敛于分块上三角矩阵，对角线上的子块是一阶或二阶的，由每个二阶子块可以求出一对共轭的复特征值。例如可看成一个分块上三角阵，它的特征值就是3，5，以及，后者是2阶子矩阵的特征值。 QR算法计算量较大，具体的计算方法这里不作介绍。 4.4 上机实验参考程序 例4.5 据雅可比迭代法计算框图4.1编制程序求解方程组 , , 。其中的计算流程为：(1) 。(2) 对作:;若max<d则。程序的结构：在程序的说明部分将方程组的系数矩阵及常数项分别存入数组aa和bb。并指定迭代的精度。在主函数main中提示人工键入最大迭代次数N0,根据框图4.1进行迭代计算。最后输出方程组的解或者迭代失

22、败信息。程序中的主要常量,变量：n,eps全局常量。分别为方程组的阶数及指定的精度。aann,bbn分别存放方程组的系数矩阵及常数项,并保持不变。an+1n+1,bn+1程序中用a11ann,b1bn存放系数矩阵及常数项。int N0最大迭代次数,人工键入.当迭代达到N0次以上时程序输出失败信息,并结束。xn+1,yn+1分别用x1xn,y1yn存放相继两次的迭代向量。sum,max用sum计算连加和。max计算。程序清单如下：#include <math.h>#include <stdio.h>#define n 3 /* 方程组的阶数 */#define eps 0

23、.5e-4 /* 给定精度要求 */static double aann=10,-1,-2,-1,10,-2,-1,-1,5;static double bbn=7.2,8.3,4.2;main() int i,j,k,N0; double an+1n+1,bn+1,xn+1,yn+1; double d,sum,max; for(i=1;i<=n;i+) for(j=1;j<=n;j+) aij=aai-1j-1; bi=bbi-1; printf("n Please enter N0:"); scanf(" %d",&N0); p

24、rintf("n");/* 键入最大迭代次数N0 */ for(i=1;i<=n;i+) xi=0; k=0; do for(i=1;i<=n;i+) sum=0.0;for(j=1;j<=n;j+) if(j!=i) sum=sum+aij*xj;yi=(-sum+bi)/aii; max=0.0; for(i=1;i<=n;i+) d=fabs(yi-xi);if(max<d) max=d;xi=yi; k+; while(max>=eps)&&(k<N0); /* 当max<eps 或 k>=N0时

25、结束迭代 */ printf("n k=%dn",k); if(k>=N0) printf("n fail!n"); else for(i=1;i<=n;i+) printf("n x%d=%f",i,xi); 运行结果： k=11 x1=1.099993 x2=1.199993 x3=1.299991例4.6 根据高斯-赛德尔迭代法计算框图4.2编制程序求解上题的方程组。程序中的主要常量及变量: n,eps,aa,bb,a,b,x,N0,sum均同上例。每次迭代计算时,用变量s临时存放xi的旧值。用d计算xi与其旧值之差

26、的最大绝对值。程序清单：#include <math.h>#include <stdio.h>#define n 3#define eps 0.5e-4static double aann=10,-1,-2,-1,10,-2,-1,-1,5;static double bbn=7.2,8.3,4.2;main() int i,j,k,N0; double an+1n+1,bn+1,xn+1,yn+1; double d,s,sum,max; for(i=1;i<=n;i+) for(j=1;j<=n;j+) aij=aai-1j-1; bi=bbi-1; printf("n Please enter N0:"); scanf(" %d",&N0); printf("n"); for(i=1;i<=n;i+) xi=0; k=0; do d=0.0; for(i=1;i<=n;i+) s=xi; /* s临时存放x