线性方程组的直接法.doc

第二章 线性方程组的直接法

　　在近代数学数值计算和工程应用中,求解线性方程组是重要的课题。例如，样条插值中形成的关系式，曲线拟合形成的法方程等，都落实到解一个元线性方程组，尤其是大型方程组的求解，即求线性方程组(2.1)的未知量的数值.
　　　　　　　 （2。1）
　　其中ai j，bi为常数.上式可写成矩阵形式Ax = b,即
　　　　 　　 　　　　 （2。2）
　　其中，为系数矩阵,为解向量，为常数向量。当detA=D0时,由线性代数中的克莱姆法则,方程组的解存在且惟一,且有

　　为系数矩阵的第列元素以代替的矩阵的行列式的值。克莱姆法则在建立线性方程组解的理论基础中功不可没，但是在实际计算中，我们难以承受它的计算量。例如，解一个100阶的线性方程组，乘除法次数约为（101·100！·99）,即使以每秒的运算速度，也需要近年的时间。在石油勘探、天气预报等问题中常常出现成百上千阶的方程组，也就产生了各种形式方程组数值解法的需求。研究大型方程组的解是目前计算数学中的一个重要方向和课题.
　　解方程组的方法可归纳为直接解法和迭代解法。从理论上来说，直接法经过有限次四则运算，假定每一步运算过程中没有舍入误差,那么，最后得到方程组的解就是精确解。但是，这只是理想化的假定，在计算过程中，完全杜绝舍入误差是不可能的，只能控制和约束由有限位算术运算带来的舍入误差的增长和危害，这样直接法得到的解也不一定是绝对精确的.
　　迭代法是将方程组的解看作某种极限过程的向量极限的值，像第2章中非线性方程求解一样,计算极限过程是用迭代过程完成的，只不过将迭代式中单变量换成向量而已。在用迭代算法时,我们不可能将极限过程算到底，只能将迭代进行有限多次，得到满足一定精度要求的方程组的近似解。
　　在数值计算历史上，直接解法和迭代解法交替生辉。一种解法的兴旺与计算机的硬件环境和问题规模是密切相关的.一般说来，对同等规模的线性方程组,直接法对计算机的要求高于迭代法。对于中等规模的线性方程组，由于直接法的准确性和可靠性高一般都用直接法求解。对于高阶方程组和稀疏方程组(非零元素较少）,一般用迭代法求解.
§1  消元法
　　一、 三角形方程组的解
　　形如下面三种形式的线性方程组较容易求解。
　　对角形方程组
　　　　　 　　　　 　　（2。3）
　　设，对每一个方程,。
　　显然,求解n阶对角方程的运算量为 。
　　下三角方程组
　　　 　　　　　（2。4)
　　按照方程组的顺序，从第一个方程至第个方程，逐个解出。
　　由方程，得。将的值代入到第二个方程

　　得　　　　　　　　　
　　将的值代入到第个方程

　　得　　　　　　
　　计算需要次乘法或除法运算,。因此,求解过程中的运算量为

　　上三角方程组
　　 （2。5）
　　与计算下三角方程组的次序相反，从第个方程至第一个方程,逐个解出。
　　由第个方程。将的值代入到第个方程

　　得　　　　　　
　　将的值代入到第个方程

　　得解的通式

　　计算需要次乘法或除法运算。因此求解过程中的运算量为

　　消元法的基本思想就是通过对方程组做初等变换，把一般形式的方程组化为等价的具有上述形式的易解方程组.
　　二、 高斯消元法与列主元消元法
　　高斯消元法
　　高斯消元法是我们熟悉的古老、简单而有效的解方程组的方法。下面是中学阶段解二元方程组（高斯消元法）的步骤：

（2.6）
（2.7）
　　方程（2。6）乘以－3加到第（2。7)个方程中得

　　代入（2。6)得。
　　其方法相当于对方程组的增广矩阵做行的初等变换:

　　已是上三角矩阵,而

　　为原方程组的等价方程组，已化成易解的方程组形式.再用回代方法求解,得到：
　　这就是高斯消元法解方程组的消元和回代过程.
　　一般地，可对线性方程组（2。1)施行以下一系列变换；
　　（1）对换某两个方程的次序；
　　(2）对其中某个方程的两边同乘一个不为零的数；
　　（3）把某一个方程两边同乘一个常数后加到另一个方程的两边。
　　记变换后的方程组为：
　　　　 　 　2。8）
　　显然方程组（2。1）与（2。8）是等价方程组，或者说它们有相同的解.分别记方程组(2.1）与（2。8）的增广矩阵为：

　　可以看出,实际上是由按一系列初等换后得到的
　　（1)对换某两行元素；
　　（2）中的某行乘一个不为零的数；
　　（3）把的某一行乘一个常数后加到另一行。
　　高斯消元法就是通过以上（3)的变换，把化为等价的上三角形式。
　　下面我们以为例演示消元过程。
　　设方程组：
　　　　　　 　　　　　 (2.9)
　　其增广矩阵为：

　　（1）若，则将第一行乘以加到第二行上；将第一乘以
　　加到第三行上；将第一行乘以 加到第四行上得到
　　　　　　　　　　 　　　　（2.10）
　　即  
　　其中：

　　（2）若则将第二行乘以加到第三行上；将第二行乘以加到第四行上，得到
　　　　　　　　　　　　　　（2.11）
　　其中：

　　（3)若则将第三行乘以加到第四行上，得到
　　　　　　　　　　 (2。12)
　　其中:

　　已是上三角矩阵,于是得到了与原方程等价的易解形式的方程组：
　　 　 （2.13）
　　再对方程组（2。13）依次回代解出。
　　从式（2.12）可以得到系数矩阵的行列式的值为的对角元素的乘积。即

　　这也正是在计算机上计算方阵的行列式的常规方法
　　要将上述解方程步骤推广到阶方程组，只需将对控制量“4”的操作改成对控制量的操作即可。
　　元方程组高斯消元法的步骤如下：
　　对于约定有
（2.14)
　　经过以上消元，我们已得到与等价的方程组,其中已是一个上三角矩阵.为简单起见，仍记的元素为
　　　　　　　 　　（2.15）
　　即已得到原方程组的解。
　　高斯消元法算法
　　在算法中略去了变量，矩阵和向量的定义部分.在消元过程中，将仍放在元素的位置上。
　　1．输入：方程组阶数n，方程组系数矩阵A和常数向量b。
　　2．FOR  k:=1 TO n—1                    //消元过程
　　       ｛   FOR  i：=k+1 TO n
　　          {                  //  假定
　　            FOR  j：=k+1 TO n
　　               {  }     //   ENDFORJ

　　           ｝                            //   ENDFORI
　　       }                                //   ENDFORK
　　3．FOR i:=n TO1          //回代求解
　　    {  s:=bi
　　         FOR  j:=i+1  TO  n  DO

　　｝
　　4．输出方程组的解。
　　高斯消元法的运算量
　　整个消元过程即式(2.14）的乘法和除法的运算量为

　　回代过程即式（2.15）的乘除运算量为

　　故高斯消元法的运算量为
　　　　　　　　 　　　　　　　　　（2.16）
　　高斯消元法的可行性
　　在上面的消元法中，未知量是按照在方程组中的自然顺序消去的，也叫顺序消元法。
　　在消元过程中假定对然元素，消元步骤才能顺利进行,由于顺序消元不改变的主子式值，故高斯消元法可行的充分必要条件为的各阶主子式不为零。但是，实际上只要方程组就有解。故高斯消元法本身具有局限性。
　　另一方面，即使高斯消元法可行，如果很小,在运算中用它作为除法的分母,会导致其它元素数量级的严重增长和舍入误差的扩散。这是高斯消元法的另一缺陷。
　　例2。1 方程组

(2。17)
（2。18）
的精确解为：。在高斯消元法计算中取5位有效数字。
　　解:方程（2。17）×(－1）/0.0003+方程（2。18)得：

　　，代入方程（2.17）得。由此得到的解完全失真，如果交换两个方程的顺序，得到等价方程组

　　经高斯消元后有

　　由此可看到，在有些情况下，调换方程组的次序对方程组的解是有影响的，对消元法中抑制舍入误差的增长是有效的。
　　如果不调换方程组的次序，取6位有效数字计算方程组的解，得到

　　取9位有效数字计算方程组的解，得到

　　由此可见有效数字在数值计算中的作用.
　　列主元消元法
　　如果在一列中选取按模最大的元素，将其调到主干方程位置再做消元，则称为列主元消元法.调换方程组的次序是为了使运算中做分母量的绝对值尽量地大，减少舍入误差的影响.
　　用列主元方法可以克服高斯消元法的额外限制，只要方程组有解,列主元消元法就能畅通无阻地顺利求解，同时又提高了解的精确度。
　　更具体地，第一步在第一列元素中选出绝对值最大的元素，交换第一行和第行的所有元素,再做化简为零的操作。
　　对于每个在做消元前，选出中绝对值最大的元素，对行和行交换后，再做消元操作,这就是列主元消元法的操作步骤.由于，可证中至少有一个元素不为零，从理论上保证了列主元消元法的可行性.列主元消元法与高斯消元法相比，只增加了选列主元和交换两个方程组(即两行元素)的过程。
　　列主元消元法算法
　　1．输入：方程组阶数，方程组系数矩阵和常数向量项。
　　2．                //选主元的消元过程
　　     ｛//选择

　　 　　　    //  交换第行和第行

　　         }                                    //  ENDFORI
　　 ｝                                        //  ENDFORK
　　3．FOR i:=n TO 1                  //  回代求解

　　4．输出方程组的解 .
　　如果对于第步，从行至行和从列至列中选取按模最大的，对第行和第行交换，对第列和第v列交换，这就是全主元消元法.在k列和第v列交换时，还要记录下v的序号,以便恢复未知量xk和xv的位置.
　　3.1。3 高斯－若尔当(Gauss－Jordan）消元法
　　高斯消元法将系数矩阵化为上三角矩阵,再进行回代求解；高斯－若尔当消元法是将系数矩阵化为对角矩阵，再进行求解,即对高斯消元法或列主元消元法得到的等价增广矩阵:

　　用第4行乘加到第3行上，用第4行乘加到第2行上，用第4行乘加到第1行上，得到

　　用第3行乘加到第2行上,用第3行乘加到第1行上，再用第2行乘加到第1行上，得到
　　　　　                  (2.19）
　　为方便起见,仍记式（2。19)系数矩阵元素为 ，常数项元素为.那么

　　用初等变换化系数矩阵为对角矩阵的方法称为高斯－若尔当消元法。从形式上看对角矩阵（2。15）比上三角矩阵(2.11)更为简单，易于计算 ，但是将上三角矩阵（2。11）化为对角矩阵（2。15 ）的工作量略大于上三角矩阵回代的工作量.
　　高斯—若尔当消元法求逆矩阵
　　设为非奇异矩阵，方程组的增广矩阵为。如果对应用高斯-若尔当消元法化为，则。
　　例2。2 用高斯—若尔当消元法求的逆矩阵。
　　解：

　　解得:
§2  直接三角分解法
　　仍以为例,在高斯消元法中,从对方程组进行初等变换的角度观察方程组系数矩阵的演变过程.第一次消元步骤将方程组（2.9）化为方程组（2.10）,相当于用了三个初等矩阵左乘 和.记，容易验证，

　　由得到其中

　　将方程组(2.10）化为方程组（2。11），相当于用了两个初等矩阵左乘 和。
　　记有

　　·

　　同理，将方程组（2.11）化为方程组（2。12），相当于用一个初等矩阵左乘 和。
　　记,有

　　完成了消元过程,有

　　亦有所有消元步骤表示为：左乘一系列下三角初等矩阵。容易验证,这些下三角矩阵的乘积仍为下三角矩阵并有

　　于是有

　　或

　　这里仍为下三角矩阵,其对角元素为1，称为单位下三角矩阵,而已是上三角矩阵。
　　记，则有

　　结果表明若消元过程可行，可以将分解为单位下三角矩阵与上三角矩阵的乘积。由此派生出解方程组的直接分解法。
　　由高斯消元法得到启发，对消元的过程相当于将分解为一个上三角矩阵和一个下三角矩阵的过程。如果直接分解得到和,。这时方程化为，令,由解出；再由，解出。这就是直接分解法.
　　将方阵分解为，当是单位下三角矩阵，是上三角矩阵时,称为多利特尔（Doclittle）分解；当是下三角矩阵，是单位上三角矩阵时,称为库朗（Courant）分解。分解的条件是若方阵的各阶主子式不为零,则多利特尔分解或库朗分解是可行和惟一的。
　　一、 多利特尔分解
　　多利特尔分解步骤
　　若的各阶主子式不为零,可分解为，其中为单位下三角矩阵，为上三角矩阵,即
　　 （2.20）
　　矩阵和共有个未知元素,按照的行元素的列元素的顺序，对每个列出式(2.16)两边对应的元素关系式，一个关系式解出一个或的元素。下面是计算和的元素的步骤:
　　（1）计算的第一行元素
　　要计算,则列出式（2。20）等号两边的第1行第1列元素的关系式:

　　故。一般地，由的第一行元素的关系式

　　得到

　　（2）计算的第一列元素

　　要计算，则列出式(2.20)等号两边的第2行第1列元素的关系式：

　　故.一般地，由的第1列元素的关系式

　　得到

　　（3）计算的第2行元素
　　(4）计算的第2列元素
　　……
　　若已算出的前行,的前列的元素，则
　　（5）计算的第行元素
　　由的第行元素的关系式：

　　是上三角矩阵，列标行标.

　　得到
　　　　　　　　　　　（2.21）
　　（6）计算的第列元素
　　由的第列元素的关系式:

　　∵是下三角矩阵，∴行标标行标.

　　得到
　　　  　 （2.22）
　　一直做到的第列,的第行为止。
　　用直接分解方法求解方程组所需要的计算量仍为，和用高斯消元法的计算量基本相同。
　　可以看到在分解中的每个元素只在式（2.21）或（2。22）中做而且仅做一次贡献，如果需要节省空间，可将 以及的元素直接放在矩阵相应元素的位置上。
　　用直接分解法解方程，首先作出分解令，解方程组.由于是单位下三角矩阵，容易得到
　　　　             　   (2。23）
　　再解方程组,其中
　　　 　　　 （2.24)
　　对进行分解时，并不涉及常数项。因此，若需要解具有相同系数的一系列线性方程组时,使用直接分解法可以达到事半功倍的效果。
　　多利特尔直接分解算法
　　1。输入:方程组阶数,系数矩阵和常数项.
　　2。
　　                 //   计算的第行元素

　　                    //计算的第列元素

　　　　　　｝                           //   完成分解
　　3。                    //  解方程组

　　4.                      //  解方程组

　　5.输出方程组的解
　　例3.2 用多利特尔分解求解方程组：

　　解：
　　对，有

　　对，有

　　对,有

　　解,即

　　 解，即

　　二、 追赶法
　　很多科学计算问题中，常常所要求解的方程组为三对角方程组
　　　　　　　　　　　　　　　　　　　　　　　　(3。25)
　　其中
　　　 　　　 　（2.26）
　　并且满足条件：

　　称为对角占优的三对角线矩阵。其特征是非零元素仅分布在对角线及对角线两侧的位置.在样条插值函数的关系式中就出现过这类矩阵。事实上，许多连续问题经离散化后得到的线性方程组，其系数矩阵就是三对角或五对角形式的对角矩阵。
　　利用矩阵直接分解法,求解具有三对角线系数矩阵的线性方程组十分简单而有效。现将分解为下三角矩阵,及单位上三角矩阵的乘积。即，其中
　　　,　　（3。27）
　　利用矩阵乘法公式,比较两边元素,可得到

　　于是有
　　   　 　　　 （3.28)
　　由些可见将分解为及，只需计算及 两组数，然后解 ，计算公式为：
　　　  　　 　 （3.29)
　　再解则得
　　　　　　    　    （3。30）
　　整个求解过程是先由（3。28)及(3.29）求， 及 ，这时 是“追”的过程，再由（3。24)求出,这时 是往回赶，故求解方程组(3.25)的整个过程称为追赶法。
　　三、 平方根法
　　对称正定矩阵也是很多物理问题产生的一类矩阵，正定矩阵的各阶主子式大于零.可以证明，若正定，则存在下三角矩阵，使，直接分解的分解方法，称为平方根法。由于在平方根法中含有计算平方根的步骤，因此很少采用平方根法的分解形式.对于对称矩阵，常用的是分解。
　　对作多利特尔分解，即
（提出矩阵的对角元素）

　　由对称正定,可得，令

　　由的对称性和分解的惟一性可证

　　即。
　　(3。31）
　　是对角元素为1的单位下三角矩阵。
　　对矩阵作多利特尔或库朗分解,共计算个矩阵元素；对称矩阵的分解，只需计算个元素,减少了近一半的工作量。借助于多利特尔或库朗分解计算公式，容易得到分解计算公式。
　　设有多利特尔分解形式：

　　其中
　　在分解可套用多利特尔分解公式，只要计算下三角矩阵和的对角元素。计算中只需保存的元素,的行列的元素用的表示。由于对称正定矩阵的各阶主子式大于零，直接调用多利特尔或库朗分解公式可完分解计算，而不必借助于列主元的分解算法。
　　对于有

　　　　　　　 　　　　　　　　　　　　（2。32）

　　　　　　　     （3。33）
　　由,解方程组可分三步完成：
　　(1）解方程组，其中。
　　　　　                  （2。34）
　　（2）解方程组其中。
　　　　            　　 　　 (3.35)
　　3）解方程
　　　      　    （2。36）
　　对称矩阵的分解算法
　　1．输入：方程组阶数 ，系数矩阵和常数项。
　　2．

　　3．略去解方程组步骤
　　从矩阵分解角度看,直接分解法与消元法本质上没有多大区别，但实际计算时它们各有所长。一般来说,如果仅用单字长进行计算,列主元消元法具有运算量较少、精度高的优点，故是常用的方法。但是,为了提高精度往往采取单字长数双倍内积的办法（即作向量内积计算时，采用双倍加法，最终结果再舍入成单字长数)，这时直接分解法能获得较高的精度.
　　例3.4 用分解求解方程组:

　　解：

　　由,有

小 结
本章介绍了适合解系数矩阵稠密、低阶的线性方程组的直接方法：Gauss消去法及三角分解法,这两种方法本质上是一样的。直接法的优点是：计算量小、精度高，是一种精确地求线性方程组的方法（如果每步计算没有舍入误差）；缺点是:程序较复杂，占用内存大，所以它适用于解中小型（）线性方程组。在实际应用中Gauss列主元消去法是一种稳定的算法。当方程组的系数矩阵是三对角阵时，特别是严格对角占优，追赶法是一种既稳定、有快速的方法.