计算方法 第3章 线性方程组数值解法

第三章 线性方程组的数值解法 线性方程组是应用最为广泛的数学模型，很多复杂问题中都含有线性方程组子问题，因此讨论线性方程组问题的求解很有必要，本章将讨论线性方程组的数值解法。 线性方程组的一般形式： (3.1)如果记： 这里A称为系数矩阵，x称为解向量，b称为右端项，则得线性方程组的矩阵形式： (3.2) 求解线性方程组问题(3.1)或(3.2)的数值方法可分为两类：直接解法和迭代解法，其直接解法是通过有限次初等运算，求得其解，虽然直接解法的推导过程都是无误差的，但是由于计算机的运算都是有舍入误差的，所求解其实是一个有误差的近似解；迭代解法则是从某个初始近似解出发，按照一个确定的迭代公式得到一个更好的近似解，反复迭代，直到求得一个满足精度要求的近似解。 本章首先讨论线性方程组的直接解法然后再介绍迭代解法。3.1 消去法1.顺序Gauss消去法首先回顾一下线性代数中所讲的线性方程组消去法过程，然后归纳出消去法的数值算法，请看如下的例子：求解线性方程组解：求解线性方程组的第一阶段称为消元过程，其方法是：第2个方程减去第1个方程的倍，第3个方程减去第1个方程的2倍，得第3个方程减去第2个方程的倍，得这一过程就是消元过程，即把方程化为等价的上三角方程(对角线下变为0)。第一个阶段完成后，进入第二个阶段，称为回代过程，其方法是：先由第3个方程解出，将代入第2个方程解出，再将和代入第1个方程解出也就解出所有的未知量。如下就是所求的解： 归纳以上求解方法，求解线性方程组包括两个过程，消元过程和回代过程。首先给出回代过程的算法，回代过程其实是一个特殊形式的方程组的求解方法，就是一个上三角线性方程组的求解方法，如： (3.3)第1步：根据第n个方程解出，得 第2步：根据已求出的和第n-1个方程求，得 假如已经求出，代入第k个方程可以求出，得： (3.4)回代过程就是对实施这一公式，注意必须从后向前计算方可，所以此过程叫做回代过程，(3.4)式也是求解下三角方程的一般公式，可看作一个独立的算法。 算法3.1(上三角线性方程组的回代算法)[初始化]设置上三角方程系数矩阵A，右端项向量b[回代过程]对于循环对于循环[算法结束]再看Gauss消去法的消元过程，对于线性方程组的消去法本质上可看作将其增广矩阵用初等行变换化为梯形矩阵的过程。为清楚的表示每次消元前后系数矩阵和右端项的状态，通常以和表示系数矩阵和右端项，其元素分别记作，其上标表示在第k次消元前的状态，其初始增广矩阵为： (3.5)首先对第1列消元，也就是对增广矩阵做初等行变换将元素约化为0，希望消元之后形如： (3.6*)为此，对，令消元因子，将第i行减去第1行的倍，则可将(3.5)式变为(3.6*)式，之后的论述中，在不引起混淆的情况下，将(3.6*)式记为： (3.6)第2列消元，对，令，将第i行减去第2行的倍，则有对第3、4、…、n-1列消元，最后得到上三角方程组的增广矩阵为以上的消元过程可用一个统一公式表示为： (3.7)通常在编程求解线性方程组时，消元因子存储于矩阵对角线以下相应的位置上，因为此时这些位置的元素将被消元为0。算法3.2(顺序Gauss消去法)[初始化]置系数矩阵A，右端项向量b；[消元过程]对于循环对于循环[回代过程]执行算法3.1，即[算法结束]定理3.1顺序Gauss消去法能够顺利进行的充要条件是系数矩阵A的顺序主子式，且，特别有。证：从前述的算法可知，消元过程能够顺利进行的充要条件是步骤1.1)中对角线元素，以下只要证明顺序主子式为即可。事实上，顺序消元过程所做行变换并不改变顺序主子式的值，而在完成第k次消元过程后，增广矩阵被约化成：显然第k个顺序主子式为：2．顺序Gauss算法的计算量通常我们在度量一个有限步终止的算法的计算量时，主要考虑整个计算过程中使用乘除法运算的总次数，在以后统计算法的计算量时，我们仅统计算法中乘除运算的次数。消元过程的计算量 回代过程的计算量 总的乘除法运算的计算量为总计算量=可以估计出总的加减法运算的计算量也约为列主元Gauss消去法前节所讲顺序Gauss消去法，有以下两个问题：在消元过程中当某个元素对角线为零时，算法无法进行；如果某个对角线元素的绝对值非常小，将会引起很大的误差。因此顺序Gauss消去法不具实用性，但由此方法改进的列主元消去法是最为广泛使用的方法，以下我们通过一个实例体会一下顺序Gauss消去法的问题和稍作改进之后的效果：用顺序Gauss消去法求解线性方程组(运算过程中保留4位有效数字)解：写出方程组的增广矩阵消元过程：，第2行减去第1行的倍，得回代过程得到计算解为：，此方程的准确解为，显然此计算解已经没有任何意义。如果将两个方程交换一下位置再求解，增广矩阵：消元过程：，第2行减去第1行的倍，得回代过程得到计算解为：，与精确解完全相同！ 第1种解法产生非常大的误差的原因是的绝对值太小，做为分母将产生太大的误差而导致最后的解面目全非；用第2种方法是通过交换方程的位置使做分母的的绝对值尽可能的大，正如第一章误差分析所述，避免绝对值太小的数做分母，可以提高计算结果的精度。这一过程称为选主元(所选的对角线元素称为主元)，一般的做法如下。 设进行第k步计算时增广矩阵如下：为了在第k次消元过程中使消元因子中的分母绝对值尽可能的大，在第k列中取第k到第n个元素中绝对值最大者： (3.8)在选定主元之后，考虑如下情况如果，则线性方程组的行列式的值为0，其解不定(无解或有无穷多组解，停止计算)；如果，则的绝对值最大，以作为主元；如果则交换矩阵的第行，交换之后第k行对角线元素绝对值变为最大。而交换增广矩阵的两行相当于交换两个方程的位置，显然不改变原方程的解；另一个需要指出的是，由于误差的原因判断的方法采取：对于指定的充分小的，如果，则认为。增加选主元机制的Gauss消去法即为列主元Gauss消去法，列主元Gauss消去法的计算过程包括：选主元，交换主元行，消元。以下我们写出适合编程的Gauss算法。算法3.3(列主元Gauss消去法)[初始化]置系数矩阵A，右端项向量b，取充分小的数；[消元过程]对于循环[选主元]如果，停止计算，算法失败如果交换方程，[消元]对于循环1.4.1)1.4.2)1.4.3)[回代过程]执行算法3.1，即计算[算法结束]关于列主元Gauss消去法的计算量，从列主元Gauss消去法的算法可以看出，仅比顺序Gauss消去法多找主元部分，这一部分的操作主要是比较运算及必要时交换两个方程，其他的部分算法完全一样，因此影响不大。全主元Gauss消去法全主元Gauss消去法与列主元Gauss消去法的区别在于取主元的方法是如果交换增广矩阵的第行，如果交换增广矩阵的第列(本质上是交换未知量与，计算结束时应该交换回来)。显然全主元消去法取主元的方法范围更大，主元的绝对值会更大，计算过程也就更可靠更稳定，但计算量稍大。实际计算效果表明，列主元Gauss消去法已经很稳定，全主元法的改进甚微，是没必要的，所以通常使用列主元Gauss消去法。Gauss-Jordan消去法列主元Gauss消去法在消元过程中是将对角线以下的元素约化为0，将系数矩阵化为上三角矩阵，通过回代得到方程的解；Gauss-Jordan消去法是将系数矩阵约化成单位矩阵，无需回代就直接得到了方程组的解，也称为无回代过程的消去法。设线性方程组的增广矩阵为： 第1步计算，设，第1行除以，第行减去新第1行的倍，则增广矩阵变为如下形式 第2步计算，设，第2行除以，第行减去新第2行的倍，则增广矩阵变为如下形式设第k-1步计算结束时增广矩阵为如下形式 第k步计算，设，第k行除以，第k)行减去新第k行的倍，则增广矩阵变为完成n-1步计算之后，增广矩阵变为显然最后一列就是方程组的解。 注意到在求解过程中，总是假设对角线上的主元，若算法将无法进行下去，如同列主元消去法，在消元过程中也可用选列主元的策略来避免这种情况发生，除非系数矩阵奇异。 算法3.4(求解线性方程组的Gauss-Jordan消去法)[初始化]置系数矩阵A，右端项向量b，取充分小的数；[消元过程]对于循环[选主元]如果，停止计算，算法失败！如果交换方程 [规格化主元行][消元]即对于循环1.5.1)1.5.2)[算法结束]6.Gauss-Jordan消去法求逆矩阵及算法设计 首先通过一个求逆过程来写出相应算法。 例3.3用Gauss-Jordan消去法求下矩阵的逆阵： 解：在矩阵A的右边添加单位矩阵进行初等行变换此时对应初始状态的单位矩阵就是所求的逆矩阵：通常编程实现求逆矩阵时并不添加右边的单位矩阵，而是对原始矩阵直接进行变换，在完成计算时原矩阵就变成了逆矩阵。其实我们注意一下求逆的过程，有如下两个特征：首先、从初始状态到最终状态，增广矩阵的每个状态都有n个列是单位列向量，这也就意味着每次变换，都会将原矩阵的一列变为单位列向量，而将原单位矩阵的一个列向量变为逆矩阵的一列，每次变换时可让这一对列向量占用共同的列，也就是可将右边新产生的一列直接放在左边将要变为单位列向量的位置上。其次、在求逆过程中如果没有行交换，则增广矩阵中单位列向量的消失与新增，在左右两部分都是按自然顺序递增的。从例3.3可以看出，因交换主元行，右边的单位列向量顺序随主元行而变，因此在完成计算后需要换回原位置，这正说明原矩阵的行交换对应逆矩阵的列交换这一性质。以下过程是通常编程时的实际计算过程：整个计算有个消元过程，但是最后一次消元不须选主元，所以主元行标分别是{2,3}，按逆序分别交换列，则得到正确的逆矩阵，此过程直观性太差，仅作为调试程序时核对之用。 算法3.5(矩阵求逆)[初始化]输入矩阵A，并以向量p记录主元行下标(共n-1个主元行下标)，取充分小的数；[行变换]对循环[选主元]，如果，矩阵不可逆，算法失败；否则计算[交换主元行]如果交换行[规格化主元行]，[消元]对循环1.5.1)1.5.2)[交换列]对循环如果，交换列[算法结束]以下是计算机求逆验证矩阵 (3.9) ， (3.10)3.2 LU分解 1.消去法的矩阵运算表示形式 从矩阵计算的角度讲Gauss消去法是通过一系列初等行变换将矩阵约化成上三角矩阵；Guass-Jordan消去法是通过一系列初等行变换将矩阵约化成单位矩阵。矩阵的初等行变换本质上是矩阵左乘一个初等矩阵，因此消去法也可以用矩阵乘积形式进行描述。 Gauss消去法的矩阵形式： Gauss-Jordan消去法的矩阵形式： 其中(3.9)式和(3.10)式中的就是我们前述消去法的消元因子，如果记： (3.11)则Gauss消去法的整个过程可以表示为：注意到形如(3.11)式的单位下三角矩阵具有以下两个性质：逆阵仍然是单位下三角阵；任意两个单位下三角阵的乘积仍然是单位下三角阵，所以上式也可以表示为：因此Gauss消去法的本质上是将系数矩阵A进行如下的矩阵分解： (3.12)其中是单位下三角矩阵，是上三角矩阵，此种分解也称为杜里特尔(Doolittle)分解，如果是下三角矩阵，是单位上三角矩阵，则称为克洛特(Crout)分解，以下内容讨论的是杜里特尔分解方法。LU分解上述(3.12)式也可以写成： 实现矩阵的LU分解算法既可以仿照消去法设计算法也可以直接从该式推导出计算公式，以下我们讨论直接LU分解方法。首先用L的第1行乘以U的各列，则有：，由此可以得到U的第1行：再让L的各行乘以U的第1列得到，由此又可以得到L的第1列假如我们已经得到了L的前k-1列和U的前k-1行：按如下方法确定U的第k行和L的第k列，将L的第k行与U的第各列相乘，有 从该式可以得到U的第k行 将L的第i(i>k)各行与U的第k列相乘，有 从该式可以得到L的第k列 则有如下结果：按如此方法，从执行到n-1就完成了矩阵的LU分解。综上所述，我们得到直接LU分解的计算公式： (3.13)定理3.2设n阶矩阵A的顺序主子式均不为0，则A存在LU分解，且分解是唯一的。事实上，如果分解成功，A的第k个顺序主子式的值就是，因此结论是显然的。本定理表明LU分解可以正常进行的前提是顺序主子式不为0，表现在分解过程中就是所有。 通常可以将下三角矩阵L和上三角矩阵U放在一个矩阵内，即： (3.14)例3.4求矩阵A的LU分解 解：设对于k=1，利用(3.13)式，先计算U的第1行 再根据的计算公式计算L的第1列 ， ， 当前的分解为 对于k=2，根据的计算公式计算U的第2行 ， ， 再根据的计算公式计算L的第2列 ， 当前的分解为对于k=3，根据的计算公式计算 最终得到分解为：将下三角矩阵L和上三角矩阵U放在一个矩阵内，如：用LU分解求解线性方程组考虑线性方程组如果系数矩阵A存在LU分解A=LU，代入原方程组令，则可先求解然后再求解即可得到原线性方程组的解。此处两个方程组，一个是下三角方程，一个是上三角方程。上三角方程的求解可以使用算法3.1，现在考虑下三角方程的求解。将下三角方程展开求解上三角线性方程组是按的次序求解所以称为回代过程，求解此下三角线性方程组时次序正好相反是所以称为顺代过程，计算公式如下： 算法3.6(求解线性方程组的直接LU分解算法)[初始化]输入系数矩阵A，右端项b[LU分解]对于循环[计算U的第k行][计算L的第k列][解下三角方程][解上三角方程][算法结束]例3.5用LU分解方法求解下列线性方程组 解：先对系数矩阵LU分解得到：解下三角方称：解上三角方程：前述的LU分解相当于顺序Gauss消去法，如同顺序Gauss消去法不具实用性一样，算法3.6也不具实用性，应该使用更稳定的列主元策略。由线性代数知识可知，交换两行相当于系数矩阵左乘一个交换矩阵，将所有交换矩阵的乘积记为P，就得到如下定理：定理3.3如果矩阵A非奇异(矩阵可逆)，则存在交换矩阵P，单位下三角阵L和上三角阵U，使得 (3.15) 用列主元实现矩阵的直接LU分解远不如对算法3.3稍作改进更容易，以下我们给出这一改进的算法。 算法3.7(列主元LU分解算法)[初始化]设置系数矩阵A，向量存放主元行下标，取充分小的数；[分解过程]对于循环[选主元]如果，停止计算，算法失败！如果交换第行对于循环1.4.1)1.4.2)[算法结束]本算法的几点说明：向量p记录主元行的下标，必要时可以生成交换矩阵P，完成矩阵或方程的交换功能，在求解方程时根据主元行下标对右端项作相应的交换；交换主元行时，要对n个元素进行交换，这是因为此时矩阵L也要做相应的交换；每次循环都对k行k列以后的元素进行更新，步骤1.4.1)是对的最后计算结果，对的计算则隐含在步骤1.4.2)中；完成算法之后，下三角矩阵L除主对角上的1之外其他元素置于矩阵A对应的位置上，上三对角矩阵U置于矩阵A的对应位置上。3.3 三对角方程组的追赶法本节考虑如下形式的线性方程组 (3.16)称上述方程组为三对角方程组，其系数矩阵称为三对角矩阵。其矩阵形式为： (3.17)其中系数矩阵A是三对角矩阵。在一些实际问题中，例如解常微分方程边值问题、偏微分方程的差分方程、解热传导方程、三次样条函数、数值微分等，都会遇到此类线性方程组，求解此类方程组，用Gauss消去法显然会浪费内存和计算量。以下我们通过LU分解，导出求解三对角方程更为简便有效的算法，称其为追赶法。 对三对角矩阵进行求解LU分解：右端两矩阵相乘我们可以得到：从而得到以下三对角阵LU分解的递推公式： (3.18)在得到三对角阵的LU分解之后，先解下三角方程：由此得到y的计算公式： 再求解上三角方程：由此得到x的计算公式： 定理3.4设三对角线性方程组中系数矩阵A满足：则系数矩阵A非奇异，且有：证明：注意到结论1)与系数矩阵A非奇异是等价的。下面我们用归纳法来证明结论1)和结论2)。对于，，显然结论成立；先假设结论对成立，则对有从而可以得到结论1)和结论2)均成立。考虑结论3)，根据结论2)和的计算公式有：算法3.8(求解三对角方程的追赶法)[初始化]输入数据，令。[矩阵分解][解下三角方程][解上三角方程][算法结束]例3.6求解三对角方程(计算过程保留小数点以后3位数字)解：首先对三对角矩阵进行LU分解，结果为：先解下三角矩阵再解上三角矩阵用追赶法求解三对角方程组的计算量，也就是乘除法运算次数仅为，所以对于三对角方程组通常不使用Guass消去法而是使用追赶法也就很容易理解了。有时还会遇到如下形式称为广义三对角方程组(如第4章的三次样条中第三类边界条件形成的三弯矩方程)： (3.19)对(3.19)式也可以设计出计算量的算法，仿照Gauss消去法，可用如下方法求解。 首先消元第一列：显然仅需要消元，在完成对应的两个方程的消元之后会增加两个位置的非零元素，我们仍然用表示，成为如下形式：一般情况，对于，每列仅两个位置的元素应该消元，同时处在的位置上会新增两个非零元素，我们仍记为本次被消元的元素符号，当完成了的消元过程时是如下形式 (3.20)显然只要再完成(3.20)式对的消元就可以将其变为如下形式 (3.21)剩下的问题就是一个简单的回代算法。算法3.8A求解广义三对角方程算法[初始化]输入数据对于执行，，，，对于做，，，，对于，[回代过程][算法结束]显然容易算出该算法的计算量约为次乘除法运算。3.4 平方根法平方根法本节我们讨论对称正定矩阵的一个更为有效的分解方法—平方根法或Cholesky分解，首先回顾线性代数中对称正定矩阵的有关结论。 定义3.1设矩阵满足，则称矩阵A是对称矩阵，如果对任意非零向量均有 则称A是正定矩阵，如果有则称A是半正定矩阵。 定理3.5设A是对称正定矩阵，则有以下结论成立：所有特征值为正；所有顺序主子式为正；如果对称正定矩阵可逆，则逆矩阵也是对称正定矩阵。定理3.6设，且的所有顺序主子式均大于0，则存在唯一的单位下三角阵和对角阵，使得 (3.22)证明：因为的顺序主子式不为零，所以存在唯一的LU分解A=LU，注意到上三角矩阵U可以分解为：将其带入则有，由对称性得再根据LU分解的唯一性，则有得证。 注意到如果矩阵A对称正定的，则有，令：将(3.22)式可以改写成如下形式：则有如下结论。定理3.7(Cholesky分解)设是正定对称矩阵，则存在唯一的非奇异下三角阵，使得 (3.23)其中的对角元素均为正数，称为Cholesky(乔莱斯基)分解或平方根分解。 至此我们仅给出了Cholesky分解的存在性条件，以下讨论分解算法。如同LU分解，也考虑直接的分解方法，从下式推导出计算公式： (3.24)我们写出下三角部分的表达式为假设已得到L的第1列到第k-1列，现要计算第k列，则有 (3.25) 例3.7求以下对称矩阵的Cholesky分解，并求其行列式的值 解： 2．改进平方根法 平方根算法是稳定的，缺点是算法中有开平方运算。一次开平方运算的CPU时间远远大于乘除运算(开平方运算一般是迭代法实现的)，因此通常我们应该尽量避免使用平方根运算，改进平方根法就是无开方运算的平方根法。 其实，所谓无开方运算的平方根法就是基于上节所述的(3.22)式进行求解，将(3.22)展开以作为中间矩阵，则有，直接对按矩阵乘法规则推导出计算公式为： (3.26)如果求解线性方程组Ax=b，以下只要分别求解下三角方程和上三角方程，即用如下计算公式： (3.27)在编程实现改进平方根算法时，对角阵可以存放在矩阵的对角线上，中间矩阵可以借用矩阵的上三角部分，改进平方根算法如下： 算法3.9改进平方根矩阵分解算法[初始化]输入系数矩阵和右端项[分解]对于循环[计算L]对[计算D][算法结束]注：(1)在编程时，如果希望将下三角矩阵L和中间矩阵T存放在原矩阵A内，只需要将算法3.9中的所有改为即可，其它不须作任何修改；(2)如果利用改进平方根算法求解线性方程组，可以在算法3.9的基础上按照(3.27)公式继续计算即可以得到方程组的解。例3.8用改进平方根算法分解以下矩阵解：按算法3.9步骤1)可求得如下结果：即：迭代法以上所述Guass消去法、LU分解求解线性方程组的方法均属于直接方法，还有一类求解线性方程组的方法属于迭代法，迭代法在求解大型稀疏问题、病态方程组等特殊情况时是一个非常有效的方法。求解线性方程组的直接方法和迭代法的主要区别是：直接法可以用有限次初等运算求出其解，比如Gauss消去法可以用大约次乘除法运算和次加减法运算求出方程组的解；而迭代法思想如同非线性方程的迭代法一样，先给出一个初始近似解，根据一个迭代公式得到一个更好的近似解，通过反复执行迭代过程逐步逼近到精确解，当近似解达到指定的精度时结束计算。如同求解非线性方程的迭代法，求解线性方程组迭代法的也是先把所给的线性方程组，化成同解的线性方程组 (3.28)这里称B为迭代矩阵，给定初始向量并按迭代格式 (3.29)进行迭代计算。如果按上述迭代格式所得到的向量序列收敛于某一向量，则就是方程组的解，并称此迭代法收敛；否则称为不收敛或发散。本课程仅介绍Jacobi迭代法和Gauss-Seidel迭代法，并简单讨论其收敛条件。Jacobi迭代法设有线性方程组写成其紧凑形式为 假设系数矩阵非奇异且主对角线元素，将方程组的非对角线项移至等式的右端并两端除以对角线元素，则方程组的等价形式为： (3.30)按此式建立迭代格式： (3.31)则称为Jacobi迭代格式，其相应的迭代方法称为Jacobi迭代法。为便于今后的收敛性分析，现将迭代公式改写成矩阵形式。记则，方程组改写成：则有：相应的矩阵形式的迭代公式为：简记为： (3.32)其中称为Jacobi迭代矩阵，。在非线性方程的迭代法中，迭代过程的终止准则是两相邻迭代近似解之差的绝对值小于给定精度，而对于线性方程组的迭代法，其终止准则是两次相邻迭代近似解之间的某种范数小于给定精度，即 例3.9用Jacobi迭代法求解方程组 解：等价形式的方程组为：对应的迭代格式为：取初始点，终止精度为，迭代过程如下表：序号00.0000.0000.00010.444-0.1251.00020.417-0.8610.90930.918-0.7971.12740.851-1.0590.96651.043-0.9411.05960.954-1.0480.97171.035-0.9701.02780.977-1.0260.98191.019-0.9831.014100.987-1.0130.990111.010-0.9901.008120.993-1.0070.994131.005-0.9951.004140.996-1.0040.997151.003-0.9971.002160.998-1.0020.998171.002-0.9981.001180.999-1.0010.999191.001-0.9991.001200.999-1.0010.999211.000-1.0001.000经21次迭代，收敛到线性方程组的精确解，收敛速度较慢。 算法3.10求解线性方程组的Jacobi迭代法[初始化]置初始近似解，最大迭代次数，精度；[迭代]当时循环执行如下步骤对于循环如果，则停止计算，取为方程组的近似解；令，继续执行下一次循环。[算法结束]Guass-Seidel迭代法前述的Jacobi迭代的优点是算法简单，缺点是收敛速度较慢。我们稍加观察就会注意到，如果迭代是收敛的，则有理由认为一个新的近似解应该优于前一个近似解，但在迭代过程中并不是充分利用了最新的计算结果。比如在计算时仍然在使用，其实此时已经计算出来，而近似值应该比更好，如果用参与新的计算结果会更好些。Guass-Seidel迭代法就是基于此思想而设计的，因此我们将(3.31)式改写为： (3.33)称其为Guass-Seidel迭代格式，相应的迭代法称为Guass-Seidel迭代法。以下我们给出Gauss-Seidel迭代格式的矩阵形式，沿用前述的矩阵记号，(3.33)式可以写成等价的形式：迭代格式的矩阵形式为： (3.34)这里称为Gauss-Seidel迭代矩阵，。 例3.10请用Gauss-Seidel迭代法求解例3.9中的线性方程组，初始点和终止精度仍按原例所取。 解：首先写出Gauss-Seidel迭代格式为：如同例3.9，取初始点，终止精度为，迭代过程如下表：序号00.0000.0000.00010.444-0.2360.94020.497-0.8371.09730.881-1.0311.04341.016-1.0311.00451.020-1.0080.99661.006-0.9990.99871.000-0.9991.00080.999-1.0001.000显然仅经过8次迭代，已经达到了终止精度，迭代速度远胜过Jacobi迭代法。 算法3.11(求解线性方程组的Gauss-Seidel迭代法)[初始化]置初始近似解，最大迭代次数，精度；[迭代]当时循环执行如下步骤对于循环如果，则停止计算，取为方程组的近似解；令，继续执行下一次循环；[算法结束]读者如果编程实现Jacobi迭代法和Gauss-Seidel迭代法会发现，Gauss-Seidel迭代法更易编程且使用更少的存储。一般认为，Gauss-Seidel迭代法的收敛速度要优于Jacobi迭代法，数值实验表明对绝大多数线性方程组，确实如此，但并无法保证对每一个例子均是如此，读者可以把下列两个矩阵作为系数矩阵构成线性方程组，随便取非0的右端项，前述的两种迭代法，各有一个收敛而另一个不收敛。， 超松弛(SOR)迭代解线性方程超松驰迭代法是目前解大型方程组的一种最常用方法，该算法是对Gauss-Seidel迭代法的一种加速方法，其思想是根据两次相邻迭代的近似解构造一个新的更好的近似解。在Gauss-Seidel迭代过程中，对Gauss-Seidel迭代法进行如下修正： (3.35)该迭代法称为逐次超松弛(SOR)迭代法，其中系数称为松弛因子，当时称为低松弛迭代，取为就是Gauss-Seidel迭代，当称为超松弛迭代。 例3.11用超松弛迭代法求解方程组 解：取初始点，精度均取。松弛因子的迭代结果序号01.0001.0001.00015.2503.813-5.04723.1413.883-5.02933.0883.927-5.01843.0553.954-5.01153.0343.971-5.00763.0213.982-5.00473.0133.989-5.00383.0083.993-5.00293.0053.996-5.001103.0033.997-5.001113.0023.998-5.000123.0013.999-5.000松弛因子的迭代结果序号01.0001.0001.00016.3133.520-6.65022.6223.959-4.60033.1334.010-5.09742.9574.007-4.97353.0044.003-5.00662.9964.001-4.99873.0004.000-5.000 算法3.12求解线性方程组的逐次超松弛(SOR)迭代法[初始化]置初始近似解、松弛因子、迭代终止，最大迭代次数N，精度，；[迭代]当时循环执行如下步骤对于循环如果，则停止计算，取为方程组的近似解；令，继续执行下一次循环[算法结束]迭代法的收敛性我们使用任何迭代算法，都要清楚该算法必须具有收敛性、了解收敛速度如何？使用发散的迭代算法或者是收敛速度太慢的迭代算法是没有实际意义的，因此我们要研究一下前述的Jacobi迭代法和Gauss-Seidel迭代法的收敛性、收敛速度问题。 定义3.2设是上的任一向量范数。若存在使则称序列收敛于。无论在任何一种范数下，序列收敛于，根据向量范数的等价性，必可推出在其它范数下也收敛于，即收敛性与范数的选择无关。 定理3.8(迭代法收敛的充分条件)对于迭代形式的线性方程组，如果则对于任意的初始向量及任意的右端向量，解此线性程组的迭代格式是收敛的，且 (3.36) (3.37) 证明：由于，根据定理1.5矩阵非奇异，线性方程组有唯一解，即有，与迭代公式相减两边取范数，并利用矩阵范数与矩阵范数的相容性则有反复使用该不等式，则有注意到对上式两端取极限也就是迭代格式收敛。同样根据前边分析还有从而有(3.36)式成立，将(3.36)式再综合下式可推得(3.37)式成立，即 本定理说明迭代法的收敛性与右端项无关，也与所取初始点无关，但初始点的选取与收敛的快慢有关。例3.12设线性方程组的系数矩阵为试判断使用Jacobi迭代法和Gauss-Seidel迭代法求解该方程组的收敛性。 解：分解系数矩阵为对于Jacobi迭代法，有对迭代矩阵取1-范数，则有，根据定理3.8可知Jacobi迭代法收敛。对于Gauss-Seidel迭代法，有对迭代矩阵取-范数，则有，同样是根据定理3.8可得Gauss-Seidel迭代法收敛。 定理3.9(迭代法收敛的充要条件)对于迭代形式的线性方程组，迭代法收敛的充要条件是谱半径小于1，即。 证明：如果方程组的解为，则应有显然的充要条件是，根据定理1.7可知，极限成立的充要条件是，证毕。 与定理3.8相比，定理3.9的条件不容易验证，以下我们再给出一个容易验证的收敛条件，首先给出如下定义。 定义3.3如果矩阵满足则称矩阵A为对角占优矩阵，如果严格不等式成立，则称A是严格对角占优的。 定理3.10在求解线性方程组的迭代法中：如果矩阵A严格对角占优，则Jaocbi迭代法和Gauss-Seidel迭代法均是收敛的，当松弛因子满足时，SOR迭代法是收敛的；如果矩阵A是对称正定的，则Gauss-Seidel迭代法是收敛的，当松弛因子时，SOR迭代法是收敛的。在以上给出的判断迭代法收敛条件中，主要有矩阵对角占优、矩阵范数小于1、谱半径小于1、矩阵对称正定这四个条件，显然前两个判断方法简单，而后两个因为要计算特征值或顺序主子式则不便于使用。3.6 线性方程组的误差分析 在本章的前几节，我们给出了一些线性方程组的求解方法，一个还没有涉及到而且很重要的问题是，这些方法所得到的解可靠吗？如果是近似解其近似程度如何？在何种情况下所求的解是可靠的？何种方法求出的解更为可靠？这就是本节所要研究的问题。病态线性方程组在求解出方程组的计算解后，直观上会认为可用如下方法度量其计算解的准确性，计算其解的残向量：如果的值很小，可能认为比较精确。但有时并非如此，看下述线性方程组： (3.38)此方程组的精确解是，但将带入方程组计算残量的-范数值为此时残向量确实已经很小，但是无论如何也不可以看作是的近似解，其实是下列方程组的精确解而这两个线性方程组的差别仅仅是元素相差，这种差别可以当做有一个微小的扰动。定义3.4如果线性方程组的系数矩阵或右端项有一个微小的扰动将导致扰动前后两方程组的解有一个很大的变化，我们称此线性方程组是病态的，其系数矩阵称为病态矩阵；反之称此方程组是良态的。众所周知，计算机求解线性方程组的每步计算都会产生舍入误差，而当前的舍入误差又会传播到之后的求解过程中去，导致误差不断的积累，而这些计算过程中的误差及误差积累势必影响到最终计算解误差的大小。目前对这类误差的分析采用一种称之为向后误差分析的方法，其思想是在误差分析过程中，把舍入误差导致最终解的误差看作是初始数据有一个扰动，把有舍入误差的计算过程看作是精确的运算。之后的讨论就是基于如此思想的误差分析方法，为了度量方程组的病态或良性程度，首先给出关于矩阵条件数的概念。矩阵的条件数定义3.5设是非奇异矩阵，对于某一矩阵范数，称 (3.39)为矩阵A求解线性方程组的条件数。特别如果取矩阵的2-范数，则可记作：如果矩阵A是对称矩阵，分别是A的最大、最小特征值，则有 考虑(3.38)式的方程组，系数矩阵为，而的两个特征值分别为，，则其系数矩阵A在2-范数下的条件数为，显然该条件数相比较于系数矩阵的元素已经非常大了。 通常我们使用(3.39)式定义的条件数来度量线性方程组的病态程度，其中的矩阵范数可取为2-范数、1-范数、范数中的任一范数，条件数愈大方程组的病态程度就愈严重。 矩阵的条件数具有如下性质：；；；。容易验证单位矩阵的条件数为1，且无论是2-范数、1-范数、范数中的任一范数，单位矩阵的条件数均达到最小值。对于任取的非奇异矩阵A，有 特别需要指出的，2-范数下正交矩阵Q其条件数也为1，即也达到条件数的最小值。线性方程组的误差分析如前所述，可以把求解过程的误差看作是原始数据的扰动问题，首先考虑系数矩阵的扰动对方程组解的影响。设是方程组的精确解，是方程组的精确解，即假设是系数矩阵有一个扰动之后的方程组的精确解：展开并化简得到则，从而有，得最后写成如下形式 (3.40)该式给出了一个系数矩阵的扰动对最终所求解的误差界，显然影响解的精度的主要因素是系数矩阵的条件数，从该式可以看出，系数矩阵的一个微小的扰动，在计算解中被放大了倍，其实条件数就是误差的放大因子。 再看右端项的扰动，设是方程组右端项有一个扰动后方程组的精确解，即则有，从而有，再由，以下两不等式相除可得 (3.41)对于右端项的扰动，其相对误差也被放大了倍。 最后看系数矩阵与右端项均有扰动的情况。设是系数矩阵有一个扰动同时右端项也有一个扰动之后的方程组的解，即则有，从而，当扰动较小时，可以在分析中忽略扰动的二次项，即忽略上式中的项，两边取范数有从而有 (3.42)根据向后误差分析的观点，在实际计算过程中产生的误差及误差的传播、积累可以看作是原始数据的两个扰动导致的误差，从(3.42)式可以看出这些误差被放大了cond(A)倍。因此我们可以用条件数cond(A)来度量求解线性方程组的解的病态程度，条件数愈大愈病态，条件数愈小愈良性。 现在我们看一个严重病态的矩阵—Hilbert矩阵 (3.43)其逆矩阵为 (3.44)此矩阵的条件数见下表：Hilbert矩阵的条件数表n581012此表说明，以Hilbert矩阵为系数矩阵的线性方程组，将是严重病态的，当方程组的阶数大于15时，不论是用列主元Gauss消去法还是用正交分解方法求解所得解已经面目全非，读者可以编程一试。病态线性方程组的求解线性方程组是否为病态是由其系数矩阵的条件数确定的，但是求矩阵的条件数将是十分困难的，通常并不会去计算矩阵的条件数，而是根据计算过程中的某些现象来判断方程组是否为病态，比如消元过程出现小的主元，则有可能是病态的；接近奇异的系数矩阵基本肯定是病态的等等。关于病态方程组的求解方法，有以下几种可选择的方法：最理想的方法是提高计算机的字长，增加基本运算精度，用双精度、四倍精度甚至更高的精度进行计算，这些已经超出本课程讨论的范围;改善系数矩阵的条件数、即寻找预处理矩阵P，对矩阵做预处理： (3.45)这里，使得方程组变成良态的，但预处理矩阵的求解是非常困难的；迭代改善、迭代改善是最为经典的求解病态方程组的算法，其思想是用单精度求得方程组的解，之后用双精度计算其残量，然后再依此残量为右端项用单精度求解新构成的方程