实际问题中解线性方程组的经典解法.docx

实际问题中解线性方程组的经典解法第二章 在实际问题中解线性方程组的经典解法 直接法与三角形方程组的求解分析 线性方程组求解问题在许多科学计算问题中都会遇到，如应力分析、电学网络、自由振动问题等。在计算机数值方法的课程中，2线性方程组求解在样条插值、数据拟合的最小二乘法以及常微分方程边值问题中都要用到.产生的线性方程组的类型有很多，如按系数矩阵含零元素多少分类，有稠密和稀疏线性方程组之别；如按阶数的高低分类，有高阶和低阶之别；如按系数矩阵的形状和性质分类，又有对称正定、三对角线对角占优等之别，因为在电子计算机上求解，必须要考虑算法的计算复杂性以及算法的数值稳定性问题。所以针对不同类型的线性方程组有不同的解法，但是，基本的方法可归结为两大类，即为直接法和迭代法。 本章介绍的经典解法，都把原方程组化为一个或者两个三角形方程组来求解，主要包括Gauss4消去法和它的变形-直接三角分解法 设有线性方程组 Ax=b 其中 La1nöæb1öæx1ö÷ç÷ç÷a22La2n÷xç2÷çb2÷ x=ç÷ b=ç÷ MMM÷MM÷ç÷ç÷çb÷çx÷an2Lann÷øènøènø 根据线性代数知识，当detA¹0时，方程组的解存在且唯一，对æa11çça21 A=çMççaèn1a12增广矩阵(A,b)=(A(1)b(1)施行行初等变换，化A(1)为上三角形矩阵A(n)，同时b(1)化为b(n)，这时与增广矩阵(A(n),b(n)相应的线性方程组为上三角形方程组 A(n)x=b(n) (1)æa11çç0其中 (A(n),b(n)=ççMç0è(1)a12La1(1n(2)(2)a22La2nM0MM(n)Lannb1(1)ö÷(2)b2÷ ÷M÷(n)÷bnøæa(!)ç11ç (A(n)=ç0çMç0èæb1(!)öç(2)÷çb÷ (b(n)=ç2÷ çM÷çb(n)÷ènø(!)a12)La1(!n(2)(2)a22La2nMMM(n)0Lannö÷÷÷ ÷÷ø(i)设aii¹0(i=1.2.3Ln) 则的解为 它便是原方程组的解，实现上述求解过程的方法称为Gauss消去法 如果方程组的系数矩阵可分解为两个形式简单的三角形矩阵L和U的乘积，即 A=LU (1.4) 若L为下三角形矩阵，则U为上三角形矩阵，反之亦然 从而求解Ax=b的问题转化为解三角形方程组 Ly=b (1.5) 和 Ux=y æçl110L0öæuu12L其中 L=çl21lM0÷ç22÷ç11u22LçMMMM÷ U=ç0çMMMçèln1ln2Ll÷çnn÷øçè00L则Ly=b为下三角形方程组，它的第i个方程为 li1y1+li2y2+L+liiyi=bi(i=1.2Ln) 假定lii¹0按y1,y2Lyn的顺序解得 ìïyb11= ïïl11íi-1 ï-ålijyj+biïj=1ïyîi=lij 上三角形方程组Ux=y的第i个方程为 uiixi+uii+1xi+1+L+uinxn=yi假定uii¹0按xn,xn-1,Lx1的顺序求解得 u1nöu÷2n÷M÷u÷nn÷ø (i=2,3Ln) (i=1,2Ln) ynìx=nïunnïïn í-åuijxj+yiïj=i+1ïyi=ïuiiî 直接法的求解过程，在计算过程中无舍入误差的前提下，都可以经有限步算数运算而得到精确解。由于计算机的字长有限，初始数据取浮点数以及运算过程都不断地产生舍入误差，这些误差的传播和积累，会影响计算精度，所以，如何避免舍入误差的增长是设计算法时必须考虑的问题。 直接法通常需要存储系数矩阵的全部元素，当方程组的阶数很高时，需要相当大的内存空间，因此，在算法设计上应当注意节省内存，比如，对称矩阵可以只存其下三角形部分于一维数组中。 综上所述，如果矩阵A非奇异，总可以通过带有行交换或不带行交换的消元过程，将A化为非奇异上三角形矩阵A(n)因此，回代求解过程也可以进行到底，但是，在实际应用中，常常难于事先判断系数矩阵A的奇异性，因此，需要进一步考虑当A为奇异矩阵时计算过程可能发生的情况。一是消元过程的某一(n)(i)步找不到非零的aii，于是计算中断；二是虽然消元过程能进行到底，但ann=0，使回代求解过程无法进行下去，因此，在计算设计中必须考虑到上述两种可能发生的情形，此时应在算法设计中给出计算中断的信息。 1 Gauss列主元素消去法 在Gauss消元过程中位于矩阵A(k)(k=1,2Ln)的主对角线(k,k)位置上的元素(k)称为主元素，因为在计算解的分量时，都做除数，但应当避免用小的数做除akk(k)(k)(k+1£i£n)在数量级上相对小的主元akk数，即避免用较aik做除数，以防止舍入误差的扩大，降低解得精度，下面的例子说明“小主元”对解得精度的影响。 例1用Gauss消去法解线性方程组 æ10-823öæx1öæ1öç÷ç÷ç÷ ç-13.7124.623÷çx2÷=ç2÷ ç-21.0725.643÷çx÷ç3÷èøè3øèø用8位十进制尾数的浮点数计算 æ10-8ç(1)(1)解 (A,b)=(A,b)=ç-1ç-2è1ö÷3.7124.6232÷ 1.0725.6433÷ø23(1)(1)a31a21128m21=(1)=-8=-10 m31=(1)=-8=-2´108 a1110a1110(2)a22=fl(a22-m21a12)=fl(3.712+2´108) =0.00000000´109+0.2´109=0.2´109 (2)b2=fl(b2-m21b1)=fl(2+108)=0.1´109 (2)(2)(2)类似的算出a32我们看到，由于小主元10-8做除数，使得行乘数m21,m31,a33,b3(2)的数量级变得很大，这样在计算a22时，即fl(3.712)与fl(0.2´109)在计算机上做代数和时，fl(3.712)由于阶码升为9使尾数右移变成了机器零，这里揭示了数量级相对于分子数做除数，会使Gauss消去法在计算机有限字长的环境下造成较大舍入误差的原因。 经第一步消元得 æ10-8ç(2)(2) (A,b)=ç0ç0è经第二步消元得 20.2´1090.4´10930.3´1090.6´109ö÷90.1´10÷ 0.2´109÷ø1æ10-8ö231ç÷(3)(3)999(A,b)=ç00.2´100.3´100.1´10÷ ç0÷000èø这一结果表明原方程组有无穷多个解，在计算机上，将会给出“无唯一解”的信息，停止计算，但实际上，由于detA¹0，所以原方程组有唯一解，是计算过程的舍入误差使解面目全非了，从前面的分析我们看到，导致这一错误的主要原因是用“小主元”做除数。 为避免小主元做除数，在Gauss消去法中加入选主元过程，即在第k步消元时，首先在A(k)的第k列主对角线元以下元素中挑选出绝对值最大者，并通过交换的第与行对应元素，使位于对角线上，仍记为，并称之为第步消元的列主元素，然后再进行消元计算， 第一步都按列选主元的消去法称为列主元素消去法，由于列主元素满足条件 所以行乘数皆满足，因此，列主元素消去法能避免例1中所出现的问题，误差分析与数值试验表明，用列主元消去法解“良态”线性方程组，在绝大多数情况下，都可以得到令人满意的结果。 例2用列主元素消去法解例1的线性方程组Ax=b æ10-8231öç÷解 (A,b)=ç-13.7124.6232÷ ç-21.0725.6433÷èøæ-2ç®ç-1ç10-8è1.0725.6433ö÷3.7124.6232÷ 231÷ø m21=0.5,m31=-0.5´10-8,经第一步消元，得 1.072æ-2ç(A(2),b(2)=ç00.3176´101ç00.2´101è5.6430.18015´1010.3´101ö÷0.5÷ 0.1´101÷ø3，经第二步消元得 m32=0.629722921.072æ-2ç(A(3),b(3)=ç00.3176´101ç00è5.6430.18015´1010.18655541´101ö÷0.5÷ ÷0.68513854ø3回代求解得 T x=(-0.49105,8 9-02.005088,60.0376725)73方程组具有9位有效数字的解 x=(-0.491058227,-0.050886075,0.367257384)T 可见用Gauss列主元素消去法计算，得到了一个高精度的近似解。 2-2 带有行交换的矩阵分解 设第k步消元时，先交换的第行与行后再消元，用矩阵运算表示为 LkIikk(A(k),b(k)=(A(k+1),b(k+1) 其中Iikkæ1öç÷çO÷ç÷1ç÷0L1ç÷ç÷1ç÷÷ =çMOMç÷1ç÷ç÷1L0ç÷1ç÷çO÷ç÷ç1÷èøæ1ççOç1çLk=ç-mk+1,kç-mk+2,kççMç-mnkèö÷÷÷÷1÷ ÷1÷O÷÷1ø(k)aik mik=(k),(i=k+1,Ln) akk当ik=k时Iikk=I,Iikk称为置换矩阵。 于是带行交换的消元过程用矩阵运算表示为 Ln-1I1n-1n-1LL2Ii22L1Ii11(A,b)=(A(n),b(n) 由此得Ln-1I1n-1n-1LL2Ii22L1Ii11A=A(n)=U 比如，n=4则有 L3Ii33L2Ii22L1Ii11A=U 可改写为 L3(Ii33L2Ii33)(Ii33Ii22L1Ii22Ii33)(Ii33Ii2Ii11)A=U (2.4) 2令L2=Ii33L2Ii33 L1=Ii33Ii22L1Ii22Ii33 P=Ii33Ii2Ii11 2%=ILI与L的形状相同，也是一个初等下三角容易证明，当i,j>k时，Lkkijkij%和形矩阵，只是交换了Lk的第k列对角元以下的第i行与第j行的元素，因此，L1%仍然是初等下三角形矩阵，从而 L2%=L% L3L2L1为单位下三角形矩阵，于是可以写成 %PA=U L 将对n=4的处理方法推广到一般情形，令 %=I Lki-1-n1LIk+i2n+2kIk+1i+1kLIkk+1+i1kk+2+i2ILk n In-1-1i P=Iin-1n-1Iin-nLIi2I2i 2-211称P为排列矩阵，于是(2.3)式可以改写为 %LL%L% Ln-1Ln-22PA1=U 从而 %LL%L%)-1U=LU PA=(Ln-1Ln-221%LL%L%)-1 其中L=(Ln-1Ln-221为单位下三角形矩阵，它与(1.4)中的L区别仅在于第k列的主对角元一下元素的排序可能不同，取决于消元过程所做的行交换，若设行交换的次数为m，则还有detA=(-1)mdetU 式表明，带行交换的消元过程产生了矩阵PA的LU分解，相当于先对系数矩阵A施行消元过程中所需的行交换后再分解，只要A非奇异，带行交换的消元过程就可以进行到底 定理2.1 设A为非奇异矩阵，则必存在排列矩阵P以及单位下三角形矩阵L和非奇异上三角形矩阵U，使 PA=LU (2.8) 3 直接三角分解法 本节以矩阵分解定理和矩阵乘法规则为基础，研究矩阵分解的计算公式以及消去法的变形直接三角分解法 当实现了系数矩阵的三角分解之后，解线性方程组原则上化为解三角形方程组(1.5)和(1.6)其解法在本章§1中已经基本解决了。定理2.1表明，A在满足一定条件下必有A=LU或PA=LU成立，同时借助于Gauss消去法得到了L和U的表达式,这一节着重研究在A=LU或PA=LU时，L，U的元素与A的元素之间的直接关系，并在此基础上给出求解矩阵方程组AX=B的方法。 3-1 基本的三角分解法 设矩阵A=(aij)n´n的各阶顺序主子式均不为零，A有分解式 La1nöæ1öæu11u12Lu1nö÷ç÷ç÷a22La2n÷çl211uLu222n÷ç÷ºLU =MM÷çMMO÷çOM÷÷ç÷ç÷an2Lannøèln1ln2L1øèunnø 对比等式左边和右边乘积矩阵LU的第r行主对角元右边的对应元素，得 æa11çça21çMçèan1a12aij=ålrkukik=1r(i=r,L,n;r=1,L,n-1)(3.2) 再对比等式两边第r列主对角元以下的对应元素，得 air=ålikukrk=1r(i=r+1,L,n;r=1,L,n-1)(3.3) 当r=1时，由(3.2)和(3.3)式分别解出 u1i=a1ilir=ai1u11(i=1,2,L,n)(i=1,2,L,n)(3.4)(3.5)假设已求出U的第1至r-1行，L的第1至r-1列，由(3.2)和(3.3)式分别得出计算U的第r行、L的第r列元素计算公式： uri=ari-ålrkukik=1r-1(i=r,L,n;r=2,L,n)(3.6)lir=air-ålikukrk=1r-1urr(i=r+1,L,n;r=2,L,n-1)(3.7)称由(3.4)(3.7)式所表示的矩阵分解为Doolittle分解，也称LU分解。 如果在分解式(3.1)中，设U为单位上三角矩阵，L为下三角矩阵，用类似于推导Doolittle分解的方法，可以得到递推地计算L和U的公式为 lir=air-ålikukrk=1r-1(i=r,L,n;r=1,2,L,n)(3.8)(3.9) uri=ari-ålrkukik=1r-1lrr(i=r+1,L,n;r=1,2,L,n-1)0æö规定lu=0åikkrç÷k=1èø称上述分解为Crout分解。注意第r步分解是先算L的第r列，后算U的第r行，与Crout分解的计算次序恰好相反。可以证明，当矩阵A的顺序主子式全不为零，A可以分解为一个下三角矩阵L和一个单位上三角矩阵U的乘积A=LU,且分解是唯一的。 实现了系数矩阵A的Doolittle分解或Crout分解之后，方程组Ax=b的求解可以分解为Ly=b和Ux=y这两个三角形方程组求解。 设L为单位下三角形矩阵，U为上三角形矩阵，并且uii¹0(i=1,L,n)，根据公式(1.8)和(1.10)得与Doolittle分解相应的求解方程组的公式为 y1=b1ìïr-1íïyr=br-ålriyi(i=2,L,n)i=1î及 (3.10) ynìx=,nïunnïïn(3.11) íyr-åurixiïi=r+1ïxr=(r=n-1,L,1)urrïî类似地。可以得出Crout分解相应的求解出公式 b1ìy=1ïl11ïïr-1íbr-ålriyiïi=1ïyr=(i=2,L,n)lïî11及 (3.12) xn=yn,ìïní(r=n-1,L,2,1)ïxr=yr-åurixii=r+1î(3.13) 综上所述，直接三角分解法的基本方法有Doolittle分解和Crout分解法，他们都包括分解系数矩阵和求解两个三角形方程组的过程。 例1 用Doolittle法解方程组 0-3öæx1öæ10öæ210ç÷çx÷ç÷-3-4-1213ç÷ç2÷=ç5÷ ç123-4÷çx3÷ç-2÷ç÷ç÷ç÷4149-13èøèx4øè7ø解 由公式(3.4)和(3.5)计算得 (u11,u12,u13,u14)=(a11,a12,a13,a14)=(2,10,0,-3) (l21,l31,l41)Tæ31ö=ç-,2÷ è22øT对于r=2,3,4，用公式(3.6)和(3.7)计算得 (u22,u23,u24)=æç11,-12,èT17ö÷,2ø36ö-,-(l32,l42)=æç÷,1111èø32ö-,-(u33,u34)=æç÷,è1111øl43=-9,u44=-4.T解Ly=b，由公式(3.11)计算得 7æöy=ç10,20,-,-16÷ 11èø解Ux=y x=(1,2,3,4) TT用电子计算机算时，可以用二维数组A(1:n;1;n+1)按行存放增广矩阵(A,b)，数组元素记为A(i,j)(i=1,2,L,n;j=1,2,L,n+1).分解计算所得L及U的元素冲掉数组A的相应位置的元素。如用Doolittle法，则 A(r,i)¬uriA(i,r)¬lir(r=1,2,L,n;i=r,r+1,L,n+1),， r=1,2,L,n-1;i=r,r+1,L,n()值得注意的是在计算U和解y的公式中，从参与计算的量在数组A。中所处的相对位置看，它们遵循着相同的规则，如果用yr冲掉br(r=1,L,n)，那么，可以视yr为ui,n+1。在数组A中，计算uri(i=r,L,n+1)和lir(i=r+1,L,n)的公式可以按如下格式计算，即 A(r,i)¬uri=A(i,r)-第r行左方数组元素A(r,k)(k=1,L,r-1)与第i列上方数组元素A(k,i)(k=1,L,r-1)对应元乘积之和； (3.14) Ai,r)¬lir=éûA(r,r)ëA(i,r)-第行左方数组元素与第列上方数组元素对应元乘积之和ù(3.15) 利用上述格式计算，按第i步先算行数组元素、后算列数组元素的次序，在完成第n步计算时，在数组A上三角部分得到矩阵U的上三角形部分，在数组A的严格下三角部分得到L的严格下三角部分，而在数组A的最后一部列得到的是Ly=b的解y。也就是说，系数矩阵的Doolittle分解与解方程组Ly=b可以同时进行，归结为对数组A=(A,b)的分解计算。然后再解Ux=y,称这种求解方式为紧凑格式。 设N=3,r表示分解的步数，分解过程中数组A的变化情况如下 æa11a12a13b1öæa11a12a13a14ö A=çça21aab÷aç222323a÷24èa31a32a33b÷=A=ç2ça21a3÷øç22èa31a32a÷ 33a34÷øæu11u12u13u14öu13u14öu13u14öæu11¾¾r=1®ççlaa÷æu11u12r=2ç24÷¾¾®çl21uuu÷æu11u12r=3®çç21a22232324u23u÷èl31a32a33a34÷øç22èl31l32aa÷¾¾çl÷=ç24çl213334÷øç21u22èl31l32u33u34÷øçèl31分解过程可以缩写为 æa11aa13b1öæu11u12u13u14ör=1A=ç12ça21aab÷çç22232÷®çl21u22u23u÷24èa31a32a33b÷r=2 3÷øçèl31l32u33u34÷ør=3例2 用紧凑格式的Doolittle法解例1中的方程组 解：分解 æ100-310öæ2100-310öç2ç-3ç11-1217÷20÷r=12÷ A=ç-3-4-12135÷ç2çç123-4-2÷÷®ç1-33217÷r=2 è4149-137÷çøç211-11-11-11÷ç÷r=3çè2-6÷r=411-9-4-16÷øu12u13u22u23l32u33æç1çç-3çL=ç2ç1ç2çç2èöæö2100-3÷ç÷÷ç÷17÷ç11-12÷1÷ç2÷ ,U=÷ç÷332÷ç-1-÷÷ç111111÷÷ç÷6-91÷ç-4÷11øèøæö10ç÷ç÷ç20÷ç÷ y=ç÷ 17ç-÷ç11÷ç÷ç-16÷èø解Ux=y，用公式计算得 x=(1234) T