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

北京科技大学数理学院 卫宏儒 W,科学与工程计算,第4章:线性方程组的直接解法,关键词： 高斯消元法 主元消元法,高斯消元法与主元消元法 高斯消元法是一个古老的直接法,由它改进得到的选主元的消元法,是目前计算机上常用于求低阶稠密矩阵方程组的有效方法,其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题,(一)引言 为便于以下讨论,把涉及到的有关名词及问题的引出暂介绍如下: 如果未知量的个数为 n ,而且关于这些未知量x1,x2, ,xn 的幂次都是一次的(线性的),那么, n 个方程 a11x1+a12x2+ +a1nxn=b1 (1) an1x1+an2x2+ +annxn=bn 构成一个含 n 个未知量的线性方程组,称为 n 阶线性方程组。其中,系数a11,a1n,a21, ,a2n, ,an1, ,ann 是给定的常数;b1, ,bn 也是给定的常数,通常称为常数项,或称为方程组的右端. 方程组(1)也常用矩阵的形式表示,写为 Ax=b,其中,A是由系数按次序排列构成的一个n阶矩阵，称为方程组的系数矩阵，x和b都是n维向量，称为方程组的右端向量。 . 使方程组(1)中每一个方程都成立的一组数x1*,x2*, ,xn* 称为式(1)的解，把它记为向量的形式，称为解向量。,我们总是希望方程组有解,且有唯一解.由线性代数的克莱姆(cramer)规则可知,如果方程组(1)的系数矩阵A的行列式(一般记为D=A)不等于零,那么,这个方程组有唯一解,而且它们可以表示为 xi=Di/D (i=1,n) 这里,Di是指D中第i列元素用右端b1, bn代替构成的行列式. 如果方程组(1)有唯一解,我们按上面的等式求解,就必须计算n+1个n阶行列式.由行列式的定义,n阶行列式包含有n!项,每一项含有n个因子,计算一个n阶行列式就需要做(n-1)n!次乘法.而我们一共要计算n+1个n阶行列式,共需做(n2-1)n!次乘法.此外,还要做n次除法才能算出xi(i=1, n).也就是说,用这个办法求解就要做 N=(n2-1)n!+n 次乘除法运算,这个计算量是大得惊人的.例如,当n=10(即求解一个含10个未知量的方程组),乘除法的运算次数共为32659210次;,当n=40，乘除法运算次数可达3.181049次。对于上百个未知量的方程组，次数运算量就更大了。因此可莱姆规则在理论上尽管是完善的，但在实际计算中却没有什么实用价值。我们将重点讨论求解线性方程组的一种有效的数值方法。 (二)求解线性方程组的消元法 消(元)去法是求解线性方程组 Ax=b (2) 和满秩矩阵的逆阵A-1的一种直接方法.尽管它比较古老,但它具有演算步骤,推算公式都系统化的特点(对其中选主元消去法,还可以证明是稳定的).因此,它至今仍然是求解方程组的一种有效的方法. 消去法可以引出几种计算方法,下面按三角形方程组和一般线性方程组的顺序来讨论。,1.三角形方程组的解法 三角形方程组包括上三角形方程组和下三角形方程组，是最简单的线性方程组之一，实际上消元法就是通过简化一般线性方程组为三角形方程组后再求解的。上三角方程组的一般形式是：,2.Gauss消元法,(一) 高斯消去法的求解过程：分为两个阶段:首先，把原方程组化为上三角形方程组，称之为“消去”过程；然后，用逆次序逐一求出三角方程组(原方程组的等价方程组)的解，并称之为“回代”过程。下面分别写出“消去”和“回代” 两个过程的计算步骤。 消去过程: 第一步: 设a110,取 做(消去第i个方程组的x1) mi1第一个方程+第i个方程 i=2,3,n 则第i个方程变为：,因为 可得第一步消元后的方程组为：,i,j=2,3,n,i,j=2,3,n,第二步: 设 ,取 做(消 去第i个方程组的x2,i=3,4,n) mi2第二个方程 + 第i个方程 i=3,4,n 类似可得第二步消元后的方程组为：,第k步: 设 ,取 做(消去第 i个方程组的xk,i=k+1,k+2,n) mik第k个方程 + 第i个方程 i= k+1,k+2,n 类似可得第k步消元后的方程组为：,继续下去到第n-1步消元，可将线性方程组化为如下上三角方程组:,3. 主元消元法,例1:考虑如下线性方程组的Gauss消元法 求解性,2x2=1 2x1+3x2=2 解:因为a11= 0 ,故此题不能用Gauss消元法求解,但交换方程组的顺序后,就可用Gauss消元法求解了.,例题2:讨论下面方程组的解法,0.0001x1+x2=1 x1+x2=2,假设求解是在四位浮点十进制数的计算机上进行。,0.100010-3 x1 + 0.1000 101 x2 = 0.1000 101 0.1000 101 x1 +0.1000 101 x2 = 0.2000 101,解:本题的计算机机内形式为,因为a11 =0.00010,故可用Gauss消元法求解,进行第一次消元时有 a22(1)= 0.1000 101 - 104 0.1000 101 (m21=a21/a11=1/0.0001= 104) = 0.00001 105 - 0.1000 105 (对阶计算) = 0.0000 - 0.1000 105 = - 0.1000 105 ,得三角方程组 0.100010-3 x1 + 0.1000 101 x2 = 0.1000 101 -0.1000 105 x2 = -0.1000 105,回代解得x2=1 ,x1=0 严重失真! (因为本题的准确解为x1=10000/9999, x2=9998/9999）,列主元基本思想及程序框图 用高斯消去法求解线性方程组时,应避免小的主元.在实际计算中,进行第k步消去前,应该在第k列元素aik (i=k,n)中找出绝大值最大者,例如 a = max a 再把第p个方程与第k个方程组进行交换,使apk(k-1)成为主元.我们称这个过程为选主元.由于只在第k列元素中选主元,通常也称为按列选主元(或称部分选主元). 如果在第k步消去前,在第k个方程到第n个方程所有的xk到xn的系数a (i=k,n;j=k,n)中,找出绝对值最大者,例如 a =maxa 再交换第k,p两个方程和第k,q两个未知量的次序,使a 成为主元. 称这个过程为完全选主元。 不论是哪种方式选出主元,而后再按上面介绍的计算步骤进行消去的计算,一般都称为选主元的高斯消去法。在实际计算中,常用按列选主元的高斯消去法。,(k-1),(k-1),pk,(k-1),ik,kin,(k-1),pq,(k-1),ij,ki,jn,(k-1),ij,(k-1),pq,用列主元消去法解线性方程组AX=b的程序框图见右图.图中w作控制常数,通常为比较小的正实数,当某个选出的主元或完成消元后的系数ann的绝对值小于w时,就认为detA0,从而终止计算.为了节省工作单元图中还用A(k) 冲掉A,b(k)冲掉b,并注意A到的下三角部分都应变为零,而且在第k次消元计算式算出乘数mik后aik不再起作用,故用mik冲掉aik,回代后所得到的解放在数组b内。,kn-1,是,否,是,是,否,例三 用列主元消去法解方程组,解 第一次消元对,因列主元素为a31(1),故先作行变换r1_r3,然后进行消元计算可得,第二次消元对A(2) |b(2) ,因列主元素为a32(2) ,故先作行变换r2_r3,然后进行消元计算可得,由此回代,得x=(1.9272,-0.69841,0.90038)T与精确解 x=(1.9273,-0.69850,0.90042)T相比较是比较准确的.,4、高斯消去法的计算量分析 高斯消去法的乘除总运算分析如下： 消元次数k 消元乘法次数 消元除法次数 回代乘除法总次数 1 n(n-1) n-1 2 (n-1)(n-2) n-2 . . . k (n-k+1)(n-k) n-k . . n-1 2*1 1 n(n+1)/2 故高斯消去法的计算量为 N=n(n2-1)/3+n(n-1)/2+n(n+1)/2=n3/3+n2-n/3 当 N 充分大时为 n3/3,5、方法的特点 由具体计算结果可知,全主元素法的精度优于主元素法,这是由于全主元素是在全体系数中选主元,故它对控制舍入误差十分有效.但全主元素法在计算过程中,需同时作行与列的互换,因而程序比较复杂,计算时间较长.列主元素法的精度虽稍低于全主元素法,但其计算简单工作量大为减少,且计算经验与理论分析均表明,它与全主元素法同样具有良好的数值稳定性,列主元素法是求解中小型稠密性方程组的最好方法之一。 选主元消去法(包括解线性方程组的所有直接的方法)比较适用于中小型方程组.对高阶方程组,即使奇数矩阵是稀疏的,但在计算中很难保持稀疏性,因而有存储量大，程序复杂等不足,所幸的是这一缺点可用迭代法解决. 另外,高斯选主元消去法还可技巧性的解决一些特殊线性方程组。 由于误差的影响,计算过程中可能出现一些坏现象,要尽可能防止,表现在求解线性方程组的消元法比较上,则应该注意要简化运算,减小运算次数,提高效率;提高数值稳定性等.,（三）线性方程组解对系数的敏感性,在计算过程中，由于计算机字长的有限性，不可避免地产生所谓的舍入误差。同时，由于所求问题的初始数据（例如线形方程组的系数矩阵和右端项系数）往往是带有一定误差的。因此计算结果总是不可避免地带有误差，或者说，如果初始数据有扰动，势必将带来具有一定误差的计算结果。就拿A x = b来说，由于观测或计算等原因，线性方程组两端的系数A和b都带有误差A和b，这样实际建立的方程组是近似方程组(A+A)(x+x)=b+b。对近似方程组求出的解是原问题的真解x加上误差x,即x+x。而 x是由A及b引起的，它的大小将直接影响所求解的可靠性。 这种解依赖于方程组系数的误差A及b的问题，称为线性方程组解对系数的敏感性。,方程组的系数矩阵发生微小扰动，就有可能引起方程组性质上的变化，这是方程组本身的“条件问题”。 相对误差关系式： 设原线形方程组 A x= b 和近似方程组 (A+A)(x+x)=b+b 在 1、A=0,b0 2、A0,b=0,一般情形,由这些关系式可看到，带有扰动的近似方程组中，扰动的大小直接影响着所求解的相对误差，故可作如下定义： 定义：设A非奇异，称|A-1|A| 为矩阵A的条件数，记为Cond (A)，即Cond (A)= |A-1|A|. Cond (A)可反映出方程组解对系数的敏感性。对此，可通过下面的例子加以理解。,方程组 此方程组的准确解为x1=0, x2=-1。现将其右端加以微小的扰动使之变为： 经计算可得准确解为x1=2, x2=-3. 这两个方程组的解相差很大，说明方程组的解对常数项b的扰动很敏感。 同时注意到|A|=4.0001, |A-1|=3.0001104 ,故有Cond (A)=1.23 105 ,可见条件数很大。,一般来说，方程组的条件数越小，求得的解就越可靠；反之，解的可靠性就越差。 通常当从方程组A x = b求出计算解x*后，有时用残向量 r = b- Ax* 的大小来检验x*的精度，认为如果对某种范数|r|很小，就说明是好的。不过要记住这种处理在某些情况下是不准确的。 例如，对上述举例的方程组，当用某种方法求得计算解后，则有，显然|r|是很小的，但与其真实解相差很大。为说明这方面的原因，我们可以把残向量r看成是A x = b的右端有扰动的b，由相对误差关系式，有：,不可轻易用残向量|r|的大小来判断计算解的精度。,（四）病态方程组：,如果方程组AX=b由于A或b的小扰动而导致解严重失真，则此方程组称为病态方程组，否则称为良态方程组。 判定一个病态方程组的简单方法； 病态方程组一般不能用解方程组的常有方法求解，而采用“迭代求精法”来计算。（列主元消元法的应用）,(五)、LU分解法 1、LU分解法的基本思想 假定我们能把矩阵A写成下列两个矩阵相乘的形式：A=LU 其中L为下三角矩阵，U为上三角矩阵。这样我们可以把 线性方程组 Ax=b写成 Ax=（LU)x=L( U x ) = b Ly=b 令 U x=y，则原线性方程组 Ax=b Ux=y 于是可首先求解向量y使 Ly=b 然后求解 Ux=y,从而求解线性方程组 Ax=b的目的。,定义： 1.LU分解 :将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是下三角和上三角矩阵 ，而且要求U的对角元素都是1。 2.紧凑格式:由于可以把L和U两个矩阵压缩到一个数组中，而且还可以存储在原来的系数矩阵A的数组中。这种LU分解常被称为紧凑格式.,由LU=A及对L和U的要求可以得到分解的计算公式。根据下式(Doolittle分解)：,1 l21 1 l31 l32 1, ,ln1 ln2 lnn-1 1,u11 u12 u13 u1n u22 u23 u2n,un-1n u(n-1)n unn,=, ,a,ann,L U,第j个分量,第i个分量,根据矩阵乘法及相等的定义,有,得公式 u1j=a1j j=1,2,n li1=ai1 / u11 i=2,3,n,在计算机程序中常常用这种方法解线性代数方程组。 它的优点是存储量很省。L和U中的三角零元素都不 必存储，就是U的对角元素也因为都是1没有必要再 记录在程序中，这样只用一个n阶方阵就可以把L和 U贮存起来。即：下三角（包括对角元）存储L各元 素 而上三角存储U的元素。 再考察公式S会发现A中任一元素 只在计算 （ji)中用到一次以后就不再出现了，因而完全 可以利用原始数组A的单元，一个个逐次贮存L或U中 的相应元素，即： a11 a12 a13 a1n u11 u12 u13 u1n a21 a22 a23 a2n l21 u22 u23 u2n a31 a32 a33 a3n l31 l32 u33 u3n an1 an2 an3 ann ln1 ln2 ln3 unn,., ,.,.,(1),(3),(5),(2n-1),(2) (4) (6) (2n),采用LU分解有如下特点： （1）LU分解与右端向量无关。先分解，后 回代。一般说来，分解的运算次数正比于n 回代求解正比与n。求解遇到多次回代时，分 解的工作不必重新做。这样节省计算时间。 （2）分解按步进行，前边分解得到的信息 为后边所用。 （3）A阵的存储空间可利用，节省存储。,3,2,2、特殊方程组的解法,（1）追赶法 （2） 分解法-平方根法,（1）追赶法,-追赶法与稀疏线性方程组 追赶法仍然保持LU分解特性，它是一种特殊的LU分解。充分利用了系数矩阵的特点，而且使之分解更简单，得到对三对角线性方程组的快速解法。 因三对角矩阵的