数值分析 线性方程组直接解法

1、第三 章,线性方程组 直接解法,第三章目录,1. Gauus 消元法 2. 主元素法 2.1 引入主元素法的必要性 2.2 列主元素法 2.3 全主元素法 2.4 解三对角方程组的追赶法 3. 矩阵分解法 3.1 Gauss消去法的矩阵形式 3.2 矩阵的三角分解 3.3 直接三角分解法 4. 平方根法与改进的平方根法 5. 矩阵求逆 6.方程组的性态和条件数,设n阶线性方程组：,其矩阵形式为：,Ax=b (2-2),其中：,在科学研究和工程技术中所提出的计算问题中，线性 方程组的求解问题是基本的，常见的，很多问题如插值函 数，最小二乘数据拟合，构造求解微分方程的差分格式等 ，都包含了解线性方

2、程组问题，因此，线性方程组的解法 在数值计算中占有较重要的地位。,求解Ax = b，曾经学过高斯（Gauss）消元法，克莱姆（Cramer）法则，矩阵变换法等，但已远 远满足不了实际运算的需要，主要体现两个方面：一是运算的快速和准确，其次是方程组的个数增大时的计算问题。如何建立能在计算机上可以实现的有效而实用的解法，具有极其重要的意义，我们也曾指出过，Cramer法则在理论上是绝对正确的，但当n较大时，在实际计算中却不能用。,如果线性方程组Ax = b的系数行列式不为零，即det(A) 0，则该方程组有唯一解。,线性方程组的数值解法,解线性方程组的数值方法大致分为两类：,请注意：由于在计算中某

3、些数据实际上只能用有限位小 数，即不可避免地存在着舍入误差的影响，因 而即使是准确解法，也只能求到近似解。 直接法在求解中小型线性方程组（100个），特别是系数矩阵为稠密型时，是常用的、非常好的方法。,直接法：指假设计算过程中不产生含入误差，经过有 限步四则运算可求得方程组准确解的方法。,2. 迭代法：从给定的方程组的一个近似值出发，构造某种算法逐步将其准确化，一般不能在有限步内得到准确解。,这一章介绍计算机上常用的直接法，它们都是以Gauss消元法为基本方法，即先将线性方程组化为等价的三角形方程组，然后求解。,1 Gauss消元法,Gauss消元法是最基本的一种方法，下例说明其基本思想:,例

4、1,解线性方程组：,解：消去x1，进行第一次消元：首先找乘数，以 -12乘第一个方程加到第二个方程，以18乘第一个方程加到第三个方程上可得同解方程组：,例1（续）,上述Gauss消元法的基本思想是：先逐次消去变量，将方程组化成同解的上三角形方程组，此过程称为消元过 程。然后按方程相反顺序求解上三角形方程组，得到原 方程组的解，此过程称为回代过程。,再消一次元得：,二次消元后将方程化为 倒三角形式，然后进行 回代容易解出： x3 = 3, x2 = 2, x1 = 1。,我们的目的，是要总结归纳出一般情况下的n阶线性方程 组的消元公式和回代求解公式，从而得到求解n阶线性方 程组的能顺利在计算机上

5、实现的行之有效的算法。,为能更清楚地得到算法，下面以4阶线性方程组为例总结求解步骤，并且很容易地可推广至一般的n阶线性方程组。,可以检查，分别以li1乘第一个方程加到第i个方程上可以完成第一次消元，得同解方程组：,变化以后的方程组系数及右边的常数项可总结出如下的计算公式：,Gauss消元法的基本步骤3（4阶）,以方程组中第i个方程减去第二个方程乘li2 (i = 3,4)，完 成第二次消元。,上标为3的系数 和右端项可由 下面公式计算：,第三步：消元（4阶方程组需进行3次消元） 将上述 A (3) X = b(3) 中最后一个方程中的x3消为零：,然后可回代求解：由于A(4)为上三角形，所以可

6、按变量的逆序逐步回代求原方程组的解：,上述 消元、回代求解过程 很容易推广到一般的n阶线 性方程组。,经过上述消元步骤， 得到同解的上三角形方程组：A (4) x = b(4),Gauss消元法的消元过程1、2（n阶）,一般地，设 n阶方程组：,消元过程为：,第k步消元后同 解方程组中上标 为k+1的元素的 计算公式见下屏,照此消元下去，完成n1次 消元后，可将原方程组化成 同解的上三角形方程组如下：,Gauss消元法的回代过程（n阶）,回代过程：逐步回代求得原方程组的解,Gauss消元法的计算量,由于在计算机中作乘除运算量所需时间远大于作加减运算所需时间，故只考虑作乘除运算量。 由消元法步骤

7、知，第k次消元需作nk次除法，作 (n k)(n k + 1)次乘法，故消元过程中乘除法运算量为：,所以Gauss 消去法的 乘除法总运算量为：,Gauss法与Cramer法则的计算量比较,Gauss 消元法的乘 除法总运算量为：,与我们曾经介绍的Cramer法则的乘除法总运算量 (n21)n!+n 相比，由下表可知：当阶数越高时，Gauss消 元法所需乘除法次数比Cramer法则要少得多：,Gauss 消元法的优缺点：,但其计算过程中，要求akk(k)（称为主元素）均不为零，因而适用范围小，只适用于从1到n 1阶顺序主子式均不为零的矩阵A，计算实践还表明，Gauss消元法的数值稳定性差，当出

8、现小主元素时，会严重影响计算结果的精度，甚至导出错误的结果。,Gauss消元法简单易行。,2 主元素法,2.1 引入主元素的必要性 对线性方程组AX = b，若其系数行列式 det(A) 0，则该方程组有唯一 解，但是这一条件 不能保证所有主元素都不等于零，只要某一主元 素等于零，就不能用Gauss消元法求解该方程组， 即使所有主元素不等于零，但 某一主元素的绝对 值很小时，Gauss消元法也是不适用的。如下例：,例2,例2（续1）,解：为减小误差，计算过程中保留3位有效数字。 按Gauss消元法步骤：,第一次消元后得同解方程组：,第二次消元后得同解方程组,回代得解，x3 = 2.02，x2

9、= 2.40，x1 = 5.80。 容易验证，方程组（3-8）的准确解为： x1 = 2.60， x2 = 1.00，x3 = 2.0。,显然两种结果相差很大。,若在解方程组前，先交换方程的次序， 如将（2-8）交换一行与二行改写成如下所示：,再用Gauss消元法，顺序消元后得同解方程组：,回代得解：x3 = 2.00，x2 = 1.00，x1 = 2.60 与准确解相同。,例2（续2）,例2两种解法的误差分析,在例2中，对(2-8)的方程进行顺序消元时， 主元a(1)11=0.50，a(2)22=0.100都比较小，以它们作 除数就增长了舍入误差,而导致计算结果不准确。,产生上述现象的原因在

10、于舍入误差，当|a(k)kk| 很小时，进行第k次消元，要用|a(k)kk|作除数，这 样就可能增大舍入误差造成溢出停机，或者导致 错误的结果。,为了在计算过程中，抑制舍入误差的增长， 应尽量避免小主元的出现。如例2的第二种解法, 通过交换方程次序，选取绝对值大的元素作主元 基于这种思想而导出主元素法。,2.2 列主元素法,为简便起见，对方 程组（2-1），用 其增 广矩阵：,表示，并直接在增广 矩阵上进行运算。,列主元素法的具体步骤如下：,列主元素法,如此经过n1步，增广矩阵（2-9）被化成上三角形，然后由回代过程求解。 在上述过程中，主元是按列选取的，故称为列主元法。例2中的第二种解法就是

11、按列主元法进行的。,2.3 全主元素法,经过nk次消元后，得到与方程组（2-1）同解的上三角形方程组，再由回代过程求解。,例6,计算过程保留三位小数。,2.4 解三对角方程组的追赶法,在很多问题中，需要解如下形式的三对角方程组：,三对角方程组的系数矩阵为三对角阵，对于这种特殊而又简单的方程组，用前面介绍的方法求解由于有大量的零元素既占内存又浪费计算时间，显然很不经济。充分注意到三对角方程组的特点，根据顺序消元的思想导出一个简便的算法追赶法。,首先进行顺序消元， 且每步将主元系数 化为1，将方程组化为：,其中系数按下式计算 ：,然后回代 求解，得：,上述追赶法能进行到底 。,追赶法举例,用追赶法

12、解下 列三对角方程组：,例4,解： 首先将方程组化为（先追）：,然后回代（赶）求解： x5 = 0，x4 = 30/7， x3 = 6/7，x2 = 12/7， x1 = 0,可以看出，追赶法本质上还 是顺序消元法，但由于计算过程 中只涉及系数矩阵的非零元，因此大大节约了计算机内存与计算量，按乘除法次数进行比较，Gauss消元法约为n3/3，全主元法为n3/2，而追赶法仅为5n-3次，可见追赶法是求解三对角方程组的非常好的方法。,3 矩阵分解法,如果用矩阵形式表示，Gauss消元法的消元过程是对方程组（2-1）的增广矩阵（A、b）进行一系列的初等行变换，将系数矩阵A化成上三角矩阵的过程，也等价

13、于用一串初等变换阵去左乘增广矩阵，因此，消元过程可以通过矩阵运算来实现。,紧接下屏：,3.1 Gauss消元法的矩阵形式,事实上，Gauss消元法的 第一次消元相当于用初等矩阵：,第二次消元相当于用初等矩阵：,第k次消元相当于用初等矩阵：,Gauss消元法的矩阵形式,经过n1步消元后得到：,因为Lk (k = 1,2, n1) 均为非奇异阵，故它们 的逆矩阵存在。容易求出：,这说明：在 的条件下，消元过程实际上是把系数矩阵A分解成单位下三角阵与上三角矩阵的乘积的过程。,Gauss消元法的矩阵形式（续）,杜利特尔（Doolittle）分解LU分解,事实上，只要A满足一定条件，由上述结论，它一定可

14、以分解成两个三角形矩阵的乘积，即：A=LU 。,上述分解称为杜利特尔（Doolittle）分解，也称为LU分解，当系数矩阵完成三角分解后，对于求解方程组：Ax = b 。,消元过程相当于分解A = LU及求解三角形方程组Ly = b，回代过程则是求解另一个三角形方程组Ux = y，因此，解线性方程组问题可转化为矩阵的三角分解问题。,其中：L为单位下三角矩阵，,U为上三角矩阵：,3.2 矩阵的三角分解,正如Gauss消元法要在一定条件下才能进行到底一样，矩阵A也必须满足一定条件才能进行三角分解。,设A为n阶方阵，若A的顺序主子式Ai (i = 1,2,n )均不为零,则矩阵A存在唯一的Dooli

15、ttle分解。,定理2.1,下面讨论如何对A进行LU分解： （1行1列）,由于两个矩阵相等就是它们的对应元素都相等，因此通过比较A与LU的对应元素，即可得到直接计算L、U的元素的公式：A=LU 即：,紧接下屏：,对A进行LU分解,由矩阵乘法规则及比较（2-11）两端的元素，得：,即可由（2-12）求出U的第一行，L的第一列。,下面讨论一般情况，即：U的第i 行，L的第j列：,一般情况下，可由：,对A进行LU分解 （r 行r 列）,计算过程应按U第1行、L第1列（第1框），U第2行、L第2列（第2框）， 的顺序,具体分解步骤见下屏：,得到计算uij 和 lij 的公式：,对A进行LU分解的具体步

16、骤,1. 计算U的第1行，L的第1列，亦称为计算 第1框；,2. 计算U的第r 行，L的第r 列（r =2,n）， 即第r 框 ：,矩阵A的LU分解举例,解：按分解公式（2-13）， 一框一框分解，每框计算时先行后列 ：,所以：,例5,3.3 直接三角分解法,若线性方程组Ax = b的系数矩阵 A完成三角分解，A = LU，那么解方程组Ax = b等价于求解两个三角形方程组Ly = b，Ux = y，即由：,再由：,可解得：,直接三角分解法（续）,容易看出，式（2-14）与式（2-13）的运算规律相同， 或者由：,解线性方程组举例,例6,解 ：按表2-3计算：,三角分解法的几点说明,1、用三角

17、分解法求线性方程的乘除法运算量也是n3/3 数量级。由于在求出uij, lij和yi后，aij和bi无需保留， 故上机计算时，可把L，U和y存在A，b所占的单元， 回代时x取代y，整个计算过程中不需要增加新的存 贮单元。,3、完成A=LU分解后可以较容易地求出行列式|A|的值：,2、从三角分解法的推导及例中可以看出，系数矩阵的 三 角分解与右端项无关。因 而在计算多个系数矩阵 为A而右端不同的线性方程组系时，用三角分解法更 为 简便（如可用于求逆矩阵）。,三角分解法的几点说明（续）,6、分解法的优点除上述2、3外，还有： a. 可求解A2z = b,因为算A2计算量大，可用,b. 可根据A的形

18、状设计算法，当A为大型稀疏，且非零 元素有规律如带状，三对角等，作分解时能充分利用A的 特点，L，U能保持A的原状，即L与A的下三角相同，U与 A的上三角形状相同。,4、三角分解法一般也可采用选主元的技术，以使算法更 具数值稳定性。,5、也可以把矩阵A分解成一个下三角矩阵与一个单位上三 角矩阵的乘积，矩阵的 这种分解称为克劳特（Crout） 分解。,解线性方程组举例,例：下述矩阵能否分解为LU（其中L为单位下三角阵，U为 上三角阵）？若能分解，那么分解是否惟一？,4 平方根法与改进的平方根法,对称正定矩阵概念：,工程实际问题的计算中，线性方程组的系数矩阵常 常具有对称正定性，即其各阶顺序主子式

19、及全部特征值均 大于零。矩阵的这一特性使它的三角分解也有更简单的形 式，从而导出一些特殊的解法。,5. A的所有顺序主子式均为正数， 即det(Ak)0, k = 1,2,n) 。,4. A的所有特征值0；,3. A的对角线元素aii0(i = 1,2,n)；,2. A是非奇异阵，且A1也是对称正定阵；,1. A的所有顺序主子矩阵Ak(k=1,2,n)也是对称正定阵；,若A为对称正定矩阵，则：,4.1 平方根法,定理 2.2,证明： 首先A可直接作LU分解，记为A =LU1 ，,其中：,定理 2.2（续）,其中U0是单位上三角阵，则由A的对称性可得：,再由唯一性得U0 = LT，所以A = L

20、DLT，并且这种分解是唯一的，又由于U1的对角线元素Ukk就是Gauss消元法的主元素：,由于LD1/2是下三角阵，因此有推论：,Choleskg分解1,推论,若A是对称正定矩阵，则存在唯一的主对角线元素都为正的下三角阵L，使得： A=LLT,矩阵的这种分解称为Choleskg分解。 用比较法可以导出L的计算公式。设：,比较A与LLT的相应元素， 即由A = LLT可得：,计算顺序 按列进行。,因此：,Choleskg分解2,当完成矩阵A的Cholesky分解后，求解方 程组AX = b就转化求：,求解对称 正定线性方 程组的方法 称为平方根 法，也称为 Cholesky分 解法，这种 方法无需选 主元，计算过程也是稳定的。由于A