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

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,研究生学位课程,数值分析,Ax,=,b,(3.2),（,3.1,）,第,3,章 解线性方程组的直接法,1,解线性方程组有,Gram,法则，但,Gram,法则不能用于计算方程组的解，,如,n,100,，,10,33,次,/,秒的计算机要算,10,120,年,本章讲解直接法,解线性方程组的两类方法,:,直接法,:,经过有限次运算后可求得方程组精确解的方法,(,不计舍入误差,!),迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法,(,一般有限步内得不到精确解,),2,三角形方程组的解,(1),上三角方程的一般形式,（,n,1,）,n/2,次运算,3,(2),下三角方程的一般形式,（,n,1,）,n/2,次运算,4,3.1,顺序,Gauss,消元法,基本思想,：通过对,(3.1),消元,逐步将,(3.1),简化为同解的,上三角方程组,(,称为消元过程,),再由下而上地解此方程组，求得解,x,(,称为回代过程,),。,若,det,（,A,）,0,，,记其增广矩阵,5,步骤如下：,第一步：,运算量：,(n-,1,)*(,1,+n),第二步：,运算量：,(n-,2,)*(,1,+n-,1,)=,(,n-,2),n,6,第,k,步：,类似的做下去，我们有：,运算量：,(n,k)*(,1,n,k,1,)=,(,n,k,)(,n,k,2),n,1,步以后，我们可以得到变换后的矩阵为：,7,因此，总的运算量为：,加上 解上述上三角阵的运算量,(,n,+1),n,/2,，总共为：,注意到，计算过程中,处在被除的位置，因此整个计算,过程要保证它不为,0,所以，,Gauss,消元法的可行条件为：,就是要求,A,的所有顺序主子式均不为,0,，即有以下定理,8,定理,约化的主元素 的充要条件,是矩阵 的顺序主子式,即,证明,显然，当 时，定理成立,.,现设定理充分性对 是成立的，求证定理充分性对 亦成立,.,首先利用归纳法证明定理的充分性,.,设,于是由归纳法假设有,可用高斯消去法将 约化到 ，,即,9,且有,（,*,）,由设,利用,(*),式，,则有,定理充分性对 亦成立,.,显然，由假设,利用,(*),式亦可,推出,10,顺序,Gauss,消去法算法组织,将增广矩阵,(,A,b,),放入一个二维数组,.,消去每一步算出的,a,ij,(,k,),放入,a,ij,的位置,b,i,(,k,),放入,b,i,m,ik,放在,a,ik,上,.,消去完毕后,数组各位置上的数据如下示,:,回代过程的计算没有数据组织问题,.,A,b,11,算法,3.1,顺序,Guass,消去法,步,1,输入系数矩阵,A,，,右端项,b,，,置,k,:=1,；,步,2.,消元：,对,k,=1,2,n,-1,，计算,步,3.,回代,对,k,=,n,-1,1,，计算,程序见,P38,12,3.2.Gauss,列主元消元法,解,:,方法,1,小主元,a,11,回代求解,x,2,=1.0000,x,1,=0.0000.,计算结果相当糟糕,.,原因,:,求行乘数时用了,”,小主元,”,0.0001,作除数,.,例,3.4,在,F(10,4,L,U),上用,Gauss,消去法求解线性方程组,方程组的准确解为,x,*=(10000/9999,9998/9999),T,.,13,解,:,方法,2,若上题在求解时将,1,2,行交换,即,回代后得到,x=(1.0000,1.0000),T,与,准确解,x*=(9998/9999,10000/9999),T,相差不大,.,主元,a,11,方法,2,所用的是在,Gauss,消去法中加入选主元过程,利用换行,避免,”,小主元,”,作除数,该方法称为,Gauss,列主元消去法。,14,主元素是指,Gauss,消元过程中位于矩阵,A,(,k,-1),的主对角线,(,k,k,),上位置上的元素,a,kk,(,k,-1),由于计算机和方程组本身的限制，在,Gauss,消元过程中会出现零主元和小主元，而使,Gauss,消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素,.,为了提高计算的数值稳定性,在消元过程中采用选择主元的方法,.,常采用的是,列主元,消去法和,全主元,消去法,.,给定线性方程组,Ax,=,b,记,A,(1),=,A,b,(1),=,b,列主元,Gauss,消去法,的具体过程如下,:,首先在增广矩阵,B,(1),=(,A,(1),b,(1),),的第一列元素中,取,15,Gauss,列主元消去法,因为,|,m,i,k,|,1,列主元消去法能避免,”,小主元,”,作除数,误差分析与数值试验表明：,（,Gauss,）,列主元消去法,”,基本稳定,”,.,再在矩阵,B,(2),=(,A,(2),b,(2),),的第二列元素中,取,按此方法继续进行下去,经过,n,-1,步选主元和消元运算,得到增广矩阵,B,(n)A,(n),b,(n),),.,则,方程组,A,(n),x,=,b,(n),是与原方程组等价的上三角形方程组,可进行回代求解,.,然后进行第二步消元得增广矩阵,B,(3),=(,A,(3),b,(3),),.,易证,只要,|,A,|,0,列主元,Gauss,消去法就可顺利进行,.,16,算法,3.2,（列主元消去法）步骤：,1,、,输入矩阵阶数,n,，,增广矩阵,A,(,n,n,+1);,2,、,对于,(1),按列选主元：选取,l,使,(2),如果 ，交换,A,(,n,n,+1),的第,k,行与第,l,行元素,(3),消元计算：,3,、,回代计算,17,程序见,P42,列主元,Gauss,消去法,是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素,.,而,全主元,Gauss,消去法,是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素,.,但由于运算量大增,实际应用中并不经常使用,.,18,3.3,解三对角方程组的追赶法,Ax=d,其中：,(3.7),19,用,Gauss,消去法解三对角线性方程组：,其中,追,赶：,20,解三对角方程组的追赶法,其计算工作量为,6n-5,次乘除法。工作量小，,有以下定理可得证三对角矩阵求解的充分性条件。,程序见,P44,21,3.4,矩阵,LU,分解法,顺序,Gauss,消去法的矩阵表示,矩阵,若,令,记,22,则有,A,(2),=,记,令,若,则有,A,(3),=,23,如此进行下去,第,n,-1,步得到,:,A,(n,),=,其中,也就是,:,A,(n),=,L,n-1,A,(n-1),=,L,n,-1,L,n,-2,A,(,n,-2),=,L,n-1,L,n-2,L,2,L,1,A,(1),24,其中,所以有,:,A,=,A,(1,),=,L,1,-1,L,2,-1,L,n-1,-1,A,(n),=,LU,其中,L,=,L,1,-1,L,2,-1,L,n-1,-1,U,=,A,(n),.,而且有,25,式,A,=,LU,称为矩阵,A,的,三角分解,.,设,n,阶方阵,A,的各阶顺序主子式不为零,则存在唯一单位下三角矩阵,L,和上三角矩阵,U,使,A,=,LU,.,证明,只证唯一性,设有两种分解,A,=,LU,则有,=,I,所以得,于是,Ax,=,b,LUx,=,b,令,Ux,=,y,得,定理,26,LU,分解方式,直接利用矩阵乘法来计算,LU,分解,比较等式两边的,第一行,得：,u,1,j,=a,1,j,比较等式两边的,第一列,得：,比较等式两边的,第二行,得：,比较等式两边的,第二列,得：,(,j=,1,n,),(,i=,2,n,),(,j=,2,n,),(,i=,3,n,),U,的第一行,L,的第一列,U,的第二行,L,的第二列,27,LU,分解算法（续）,第,k,步：,此时,U,的前k,-1,行和,L,的前,k,-1,列已经求出,比较等式两边的,第,k,行,得：,比较等式两边的,第,k,列,得：,直到第,n,步，便可求出矩阵,L,和,U,的所有元素。,(,j=k,n,),(,i=k,+1,n,),28,LU,分解算法（续）,算法,3.3,：,(,LU,分解,),For,k,=1,2,.,n,End For,j=k,n,i=k,+1,n,Matlab,程序参见：,malu.m,为了节省存储空间，通常用,A,的绝对下三角部分来存放,L,(,对角线元素无需存储,),，用,A,的上三角部分来存放,U,。,运算量：,(,n,3,-,n,)/3,29,由,可得,这,就是求解方程组,Ax,=,b,的,Doolittle,三角分解方法。,30,利用三角分解方法解线性方程组,解,因为,所以,例,31,先解,得,再解,得,解线性方程组,Ax,=,b,的,Doolittle,三角分解法的计算量约为,n,3,/3,与,Gauss,消去法基本相同,.,其优点在于求一系列同系数的线性方程组,Ax,=,b,k,(,k,=1,2,m,),时,可大大节省运算量,.,此外，还有另一种,LU,分解，称为,Crout,三角分解,.,32,3.5,对称正定矩阵的,Cholesky,分解法,A,对称：,A,T,=,A,A,正定：,A,的各阶顺序主子式均大于零。即,定理,：,（对称正定矩阵的三角分解）,如果,A,为对称正定矩阵,，,则存在一个实的非奇异,下三角矩阵,L,，,使,A,=,LL,T,，,且当限定的对角元素为正时,，,这种分解是唯一的,。,33,U,=,u,ij,=,u,11,u,ij,/,u,ii,1,1,1,u,22,u,nn,记为,A,对称,即,记,D,1/2,=,则 仍是下三角阵,设,A,为对称正定矩阵,则有,唯一,分解,A,=,LU,且,u,kk,0,.,分解,A,=,LL,T,称为对称正定矩阵的,Cholesky,分解,.,如何求,L,34,对矩阵,A,进行,Cholesky,分解,，,即,A=LL,T,，,由矩阵乘法,：,对于,i,=1,2,n,计算（一列一列算）,求解下三角形方程组,Ly,=,b,求解,L,T,x,=y,35,Cholesky,分解法缺点及优点,2.,对于对称正定阵,A,，,从 可知对任意,k,i,有 。,即,L,的元素不会增大，误差可控，不需选主元。,缺点：存在开方运算，可能会出现根号下负数。,优点：,1.,可以减少存储单元。为一般矩阵,LU,分解计算量的一半。,36,改进,Cholesky,分解法,改进的,cholesky,分解,A,=,LDL,T,37,改进的,cholesky,分解,逐行相乘，并注意到,i,j,，有,由此可得,38,改进的,cholesky,分解,所以将上述公式改写为,解方程，令,L,T,X=y,，,先解下三角形方程组,LDY=b,得,再解上三角形方程组,L,T,X=Y,得,39,算法,3.4,：,(,Cholesky,分解法,),步,1,输入对称正定矩阵,A,，,右端项,b,；,步,2.,Cholesky,分解：,对,k,=2,3,n,，计算：,步,3.,解下三角形方程组,LDy,=b,步,4.,解上三角形方程组,L,T,x,=y,程序见,P53,40,3.6,线性方程组的性态和舍入误差对解的影响,先看一个例子，设线性方程组,试分析系数矩阵和右端项有微小扰动，解将产生什么样的变化？,解 方程组的精确解为,x,=(1,1),T,设系数矩阵有微小扰动,即,41,若右端项有微小变动,则,若同时对,A,b,扰动,A,，,b,，则,42,现计算相对误差：,43,病态方程组,由于计算机字长限制，在解,Ax,=,b,时，,舍入误差,是不可避免的。因此我们只能得出方程的近似解,x,*,。,x,*,是方程组,(,A,+,A,),x,=,b,+,b,(1),的准确解。,在没有舍入误差的解。称方程（,1,）为方程,Ax,=,b,的扰动方程。其中,A,，,b,为由舍入误差所产生的扰动矩阵和扰动向量。当,A,，,b,的微小扰动，解得（1）的解与,Ax,=,b,的解,x,的相对误差不大称为,良态方程,，否则为,病态方程,。,44,设线性方程组,Ax,=,b,1),系数矩阵是精确的,常数项有误差,b,此时记解为,x,+x,则,A,(,x,+,x,)=,b,+,b,于是,A,x,=,b,所以,x,=,A,-1,b,A,-1,b,又,由于,b,=,Ax,A,x,因此,x,b,A,A,-1,b,x,即,误差界,45,2),再设,b,是精确的,A,有误差,A,此时记解为,x,+,x,则,(,A,+,A,)(,x,+,x,)=,b,则有,A,x,+,A,(,x,+,x,)=,0,所以,x,=-,A,-1,A,(,x,+,x,),于是,x,A,-1,A,x,+,x,也就是,46,3),设,Ax,=,b,中,A,和,b,都有误差,(,小扰动,),此时记仍解为,x,+,x,则,(,A,+,A,)(,x,+,x,)=,b,b,由,Ax,=,b,，,有,A,x,+,A,x,=,b-,Ax,即,x+A,-,1,A,x,=,A,-,1,b-A,-,1,Ax,记,Cond(,A,)=,A,A,-1,称为方程组,Ax,=,b,或矩阵,A,的条件数,.,47,48,条件数的性质：,）,cond,(A)1,）,cond,(,k,A)=,cond,(A),k,为非零常数,）若 ，则,经常使用的条件数有,Cond,p,(,A,)=,A,p,A,-1,p,p=1,，,2,，。,当,A,为对称矩阵时，有,Cond,2,(,A,)=,|,1,|/|,n,|,其中,1,，,n,分别是,A,的按绝对值最大和最小的特征值。,49,例如，对前面方程组的系数矩阵,A,有,Cond,1,(,A,)=,Cond,(,A,)=39601,，,Cond,2,(,A,)39206,由于计算条件数运算量较大,实际计算中若遇到下述情况之一,方程组就有可能是病态的,:,行列式很大或很小（如某些行、列近似相关）；,元素间相差大数量级，且无规则；,主元消去过程中出现小主元；,特征值相差大数量级。,50,病态,方程组的处理,病态方程组的计算,1),用双精度或更高精度计算,2),采用数值稳定性较好的算法，如全主元法,3),用迭代改善法,设已求得方程组,Ax,=,b,的近似解,x,1,计算剩余向量,r,1,=,b,-,Ax,1,则,x,1,的迭代改善解为,:,再求解余量方程组,Ax,=,r,1,得到解,x,*,Step 1,:,近似解,Step 2,:,Step 3,:,Step 4,:,若 可被精确解出，则有,就是精确解了。,51,练习题,第,57-8,页 习题,3.3-3.7,3.15,52,