第3章解线性方程组的直接法

1、1,Ax = b (3.2),（3.1）,第3章 解线性方程组的直接法,2,解线性方程组有Gram法则，但Gram法则不能用于计算方程组的解， 如n100，1033次/秒的计算机要算10120年,本章讲解直接法,解线性方程组的两类方法: 直接法: 经过有限次运算后可求得方程组精确解的方法(不计舍入误差!) 迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法(一般有限步内得不到精确解),3,三角形方程组的解 (1) 上三角方程的一般形式,（n1）n/2次运算,4,(2) 下三角方程的一般形式,（n1）n/2次运算,5,3.1 顺序Gauss消元法 基本思想：通过对(3.1)消元

2、,逐步将(3.1)简化为同解的上三角方程组(称为消元过程),再由下而上地解此方程组，求得解x ( 称为回代过程)。,若det（A）0，记其增广矩阵,6,步骤如下：,第一步：,运算量： (n-1)*(1+n),第二步：,运算量： (n-2)*(1+n-1)=(n-2)n,7,第k步：,类似的做下去，我们有：,运算量： (nk)*(1nk1)=(nk)(nk2),n1步以后，我们可以得到变换后的矩阵为：,8,因此，总的运算量为：,加上 解上述上三角阵的运算量(n+1)n/2，总共为：,注意到，计算过程中,处在被除的位置，因此整个计算,过程要保证它不为0,所以，Gauss消元法的可行条件为：,就是要

3、求A的所有顺序主子式均不为0，即有以下定理,9,定理,约化的主元素 的充要条件,是矩阵 的顺序主子式,即,证明,显然，当 时，定理成立.,现设定理充分性对 是成立的，求证定理充分性对 亦成立.,首先利用归纳法证明定理的充分性.,设,于是由归纳法假设有,可用高斯消去法将 约化到 ，,即,10,且有,（*）,由设,利用(*)式，,则有,定理充分性对 亦成立.,显然，由假设,利用(*)式亦可,推出,11,顺序Gauss消去法算法组织 将增广矩阵(A,b)放入一个二维数组.消去每一步算出的aij(k)放入aij的位置,bi(k) 放入bi ,mik放在aik上.消去完毕后,数组各位置上的数据如下示:,

4、回代过程的计算没有数据组织问题.,A,b,12,算法3.1 顺序Guass消去法 步1 输入系数矩阵A，右端项b，置k:=1； 步2. 消元：对 k=1,2,n-1，计算,步3. 回代,对 k=n-1,1，计算,程序见P38,13,3.2. Gauss列主元消元法,解: 方法1,小主元a11,回代求解 x2=1.0000, x1=0.0000. 计算结果相当糟糕. 原因:求行乘数时用了”小主元”0.0001作除数.,例3.4 在F(10,4,L,U)上用Gauss消去法求解线性方程组,方程组的准确解为 x*=(10000/9999,9998/9999)T.,14,解: 方法2若上题在求解时将1

5、,2行交换,即,回代后得到 x =(1.0000 ,1.0000 )T,与准确解 x*=(9998/9999,10000/9999)T相差不大.,主元a11,方法2所用的是在Gauss消去法中加入选主元过程,利用换行,避免”小主元”作除数, 该方法称为Gauss列主元消去法。,15,主元素是指Gauss消元过程中位于矩阵A(k-1)的主对角线(k,k)上位置上的元素akk(k-1)由于计算机和方程组本身的限制，在Gauss消元过程中会出现零主元和小主元，而使Gauss消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素.,为了提高计算的数值稳定性,在消元过程

6、中采用选择主元的方法.常采用的是列主元消去法和全主元消去法.,给定线性方程组Ax=b, 记A(1)=A,b(1)=b,列主元Gauss消去法的具体过程如下:,首先在增广矩阵B(1)=(A(1),b(1)的第一列元素中,取,16,Gauss列主元消去法,因为|mi,k|1,列主元消去法能避免”小主元”作除数,误差分析与数值试验表明：（Gauss）列主元消去法”基本稳定”.,再在矩阵B(2)=(A(2),b(2)的第二列元素中,取,按此方法继续进行下去, 经过n-1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回

7、代求解.,然后进行第二步消元得增广矩阵B(3)=(A(3),b(3).,易证,只要|A|0,列主元Gauss消去法就可顺利进行.,17,算法3.2（列主元消去法）步骤：,1、输入矩阵阶数n，增广矩阵 A(n,n+1);,2、对于,(1) 按列选主元：选取 l 使,(2) 如果 ，交换 A(n,n+1) 的第k行与第l 行元素,(3) 消元计算 ：,3、回代计算,18,程序见P42,列主元Gauss消去法是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素.而全主元Gauss消去法是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素.但由于运算量大增,实际应用中并不经常使用.

8、,19,3.3 解三对角方程组的追赶法,Ax = d,其中：,(3.7),20,用Gauss消去法解三对角线性方程组：,其中,追,赶：,21,解三对角方程组的追赶法,其计算工作量为6n-5次乘除法。工作量小，有以下定理可得证三对角矩阵求解的充分性条件。,程序见P44,22,3.4 矩阵LU分解法,顺序Gauss消去法的矩阵表示,矩阵,若,令,记,23,则有,A(2)=,记,令,若,则有,A(3)=,24,如此进行下去, 第n-1步得到:,A(n)=,其中,也就是:,A(n)=Ln-1A(n-1),=Ln-1Ln-2A(n-2),= Ln-1Ln-2L2L1A(1),25,其中,所以有:,A=A

9、(1)= L1-1L2-1Ln-1-1A(n)=LU,其中L=L1-1L2-1Ln-1-1, U=A(n) .,而且有,26,式 A=LU称为矩阵A的三角分解.,设n阶方阵A的各阶顺序主子式不为零,则存在唯一单位下三角矩阵L和上三角矩阵U使A=LU .,证明 只证唯一性, 设有两种分解 A=LU,则有,=I,所以得,于是 Ax=b LUx=b,令 Ux=y 得,定理,27,LU分解方式,直接利用矩阵乘法来计算 LU分解,比较等式两边的第一行得：,u1j = a1j,比较等式两边的第一列得：,比较等式两边的第二行得：,比较等式两边的第二列得：,( j = 1, n ),( i = 2, n ),

10、( j = 2, n ),( i = 3, n ),28,LU分解算法（续）,第 k 步：此时 U 的前 k-1 行和 L 的前 k-1 列已经求出,比较等式两边的第 k 行得：,比较等式两边的第 k 列得：,直到第 n 步，便可求出矩阵 L 和 U 的所有元素。,( j = k, , n ),( i = k+1, , n ),29,LU分解算法（续）,算法 3.3：( LU 分解 ),For k =1,2,.,n,End For,j = k, , n,i = k+1, , n,Matlab程序参见：malu.m,为了节省存储空间，通常用 A 的绝对下三角部分来存放 L (对角线元素无需存储)

11、，用 A 的上三角部分来存放 U。,运算量：(n3 - n)/3,30,由,可得,这就是求解方程组Ax=b的Doolittle三角分解方法。,31,利用三角分解方法解线性方程组,解 因为,所以,例,32,先解,得,再解,得,解线性方程组Ax=b的Doolittle三角分解法的计算量约为n3/3,与Gauss消去法基本相同.其优点在于求一系列同系数的线性方程组Ax=bk ,(k=1,2,m)时,可大大节省运算量.,此外，还有另一种LU分解，称为Crout三角分解.,33,3.5 对称正定矩阵的Cholesky分解法,A对称：AT=A A正定：A的各阶顺序主子式均大于零。即,定理：（对称正定矩阵的

12、三角分解）,如果A为对称正定矩阵，则存在一个实的非奇异 下三角矩阵L，使A=LLT ，且当限定的对角元素为正时， 这种分解是唯一的。,34,即,则 仍是下三角阵,设A为对称正定矩阵, 则有唯一分解A=LU, 且ukk0.,分解A=LLT称为对称正定矩阵的Cholesky分解.,如何求L,35,对矩阵A进行Cholesky分解，即A=LLT，由矩阵乘法：,对于 i = 1, 2, n 计算（一列一列算）,求解下三角形方程组 Ly=b,求解LTx = y,36,Cholesky分解法缺点及优点,2. 对于对称正定阵 A ，从 可知对任意k i 有 。即 L 的元素不会增大，误差可控，不需选主元。,

13、缺点：存在开方运算，可能会出现根号下负数。,优点：1.可以减少存储单元。为一般矩阵LU分解计算量的一半。,37,改进Cholesky分解法,改进的cholesky分解A=LDLT,38,改进的cholesky分解,逐行相乘，并注意到ij，有,由此可得,39,改进的cholesky分解,所以将上述公式改写为,解方程，令LTX = y，先解下三角形方程组LDY = b得,再解上三角形方程组LTX = Y得,40,算法 3.4：( Cholesky 分解法 ),步1 输入对称正定矩阵A，右端项b；,步2. Cholesky分解：,对 k=2,3,n，计算：,步3.解下三角形方程组LDy = b,步4

14、.解上三角形方程组LTx = y,程序见P53,41,3.6 线性方程组的性态和舍入误差对解的影响,先看一个例子，设线性方程组,试分析系数矩阵和右端项有微小扰动，解将产生什么样的变化？,解 方程组的精确解为x=(1,1)T,设系数矩阵有微小扰动,即,42,若右端项有微小变动,则,若同时对A,b扰动A ，b，则,43,现计算相对误差：,44,病态方程组,由于计算机字长限制，在解Ax=b时，舍入误差是不可避免的。因此我们只能得出方程的近似解x*。,x*是方程组(A+A) x=b+b (1)的准确解。,在没有舍入误差的解。称方程（1）为方程Ax=b的扰动方程。其中A， b为由舍入误差所产生的扰动矩阵

15、和扰动向量。当A， b的微小扰动，解得（1）的解与Ax=b的解x的相对误差不大称为良态方程，否则为病态方程。,45,设线性方程组 Ax=b,1) 系数矩阵是精确的, 常数项有误差b, 此时记解为x+x ,则 A(x+x) =b+b,于是 Ax =b,所以 x = A-1b A-1 b,又由于 b= Ax A x,因此 x b AA-1bx,即,误差界,46,2) 再设b是精确的, A有误差A, 此时记解为x+x ,则 (A+A)(x+x) =b,则有 Ax +A(x+x) =0,所以 x =-A-1A(x+x),于是 x A-1Ax+x,也就是,47,3) 设Ax=b中A和b都有误差(小扰动)

16、 ,此时记仍解为x+x ,则 (A+ A)(x+x) =b+b,由Ax=b，有 Ax + A x =b- Ax,即 x+A-1 Ax =A-1 b- A-1 Ax,记 Cond(A)=AA-1, 称为方程组Ax=b或矩阵A的条件数.,48,49,条件数的性质：,）cond ( A )1,）cond ( kA )= cond ( A ) k 为非零常数,）若 ， 则,经常使用的条件数有 Condp(A)=ApA-1p p=1，2，。,当A为对称矩阵时，有 Cond2(A)=|1|/|n| 其中1，n分别是A的按绝对值最大和最小的特征值。,50,例如，对前面方程组的系数矩阵A有,Cond1(A)= Cond(A)= 39601， Cond2(A)39206,由于计算条件数运算量较大, 实际计算中若遇到下述情况之一,方程组就有可能是病态的:, 行列式很大或很小（如某些行、列近似相关）； 元素间相差