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

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)消元,逐步将(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消元法的可行条件为：,就是要求A的所有顺序主子式均不为0，即有以下定理,9,定理,约化的主元素的充要条件,是矩阵的顺序主子式,即,证明,显然，当时，定理成立.,现设定理充分性对是成立的，求证定理充分性对亦成立.,首先利用归纳法证明定理的充分性.,设,于是由归纳法假设有,可用高斯消去法将约化到，,即,10,且有,（*）,由设,利用(*)式，,则有,定理充分性对亦成立.,显然，由假设,利用(*)式亦可,推出,11,顺序Gauss消去法算法组织将增广矩阵(A,b)放入一个二维数组.消去每一步算出的aij(k)放入aij的位置,bi(k)放入bi,mik放在aik上.消去完毕后,数组各位置上的数据如下示:,回代过程的计算没有数据组织问题.,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,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消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素.,为了提高计算的数值稳定性,在消元过程中采用选择主元的方法.常采用的是列主元消去法和全主元消去法.,给定线性方程组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)是与原方程组等价的上三角形方程组,可进行回代求解.,然后进行第二步消元得增广矩阵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消去法是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素.但由于运算量大增,实际应用中并不经常使用.,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(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=bLUx=b,令Ux=y得,定理,27,LU分解方式,直接利用矩阵乘法来计算LU分解,比较等式两边的第一行得：,u1j=a1j,比较等式两边的第一列得：,比较等式两边的第二行得：,比较等式两边的第二列得：,(j=1,n),(i=2,n),(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分解),Fork=1,2,.,n,EndFor,j=k,n,i=k+1,n,Matlab程序参见：malu.m,为了节省存储空间，通常用A的绝对下三角部分来存放L(对角线元素无需存储)，用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=AA正定：A的各阶顺序主子式均大于零。即,定理：（对称正定矩阵的三角分解）,如果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，从可知对任意ki有。即L的元素不会增大，误差可控，不需选主元。,缺点：存在开方运算，可能会出现根号下负数。,优点：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.解上三角形方程组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为由舍入误差所产生的扰动矩阵和扰动向量。当A，b的微小扰动，解得（1）的解与Ax=b的解x的相对误差不大称为良态方程，否则为病态方程。,45,设线性方程组Ax=b,1)系数矩阵是精确的,常数项有误差b,此时记解为x+x,则A(x+x)=b+b,于是Ax=b,所以x=A-1bA-1b,又由于b=AxAx,因此xbAA-1bx,即,误差界,46,2)再设b是精确的,A有误差A,此时记解为x+x,则(A+A)(x+x)=b,则有Ax+A(x+x)=0,所以x=-A-1A(x+x),于是xA-1Ax+x,也就是,47,3)设Ax=b中A和b都有误差(小扰动),此时记仍解为x+x,则(A+A)(x+x)=b+b,由Ax=b，有(A+A)(x)=b-Ax,即x+A-1Ax=A-1b-A-1Ax,记Cond(A)=AA-1,称为方程组Ax=b或矩阵A的条件数.,48,49,条件数的性质：,）cond(A)1,）cond(kA)=cond(A)k为非零常数,）若，则,经常使用的条件数有Condp(A)=ApA-1pp=1，2，。,当A为对称矩阵时，有Cond2(A)=|1|/|n|其中1，n分别是A的按绝对值最大和最小的特征值。,50,例如，对前面方程组的系数矩阵A有,Cond1(A)=Cond(A)=39601，Cond2(A)39206,由于计算条件数运算量较大,实际计算中若遇到下述情况之一,方程组就有可能是病态的:,行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；主元消去过程中出现小主元；特征值相差大数量级。,51,病态方程组的处理,病态方程组的计算1)用双精度或更高精度计算2)采用数值稳定性较好的算法，如全主元法3)用迭代改