《实际问题中解线性方程组的经典解法.docx》由会员分享，可在线阅读，更多相关《实际问题中解线性方程组的经典解法.docx（25页珍藏版）》请在三一办公上搜索。
1、实际问题中解线性方程组的经典解法第二章 在实际问题中解线性方程组的经典解法 直接法与三角形方程组的求解分析 线性方程组求解问题在许多科学计算问题中都会遇到，如应力分析、电学网络、自由振动问题等。在计算机数值方法的课程中，2线性方程组求解在样条插值、数据拟合的最小二乘法以及常微分方程边值问题中都要用到.产生的线性方程组的类型有很多，如按系数矩阵含零元素多少分类，有稠密和稀疏线性方程组之别；如按阶数的高低分类，有高阶和低阶之别；如按系数矩阵的形状和性质分类，又有对称正定、三对角线对角占优等之别，因为在电子计算机上求解，必须要考虑算法的计算复杂性以及算法的数值稳定性问题。所以针对不同类型的线性方程组
2、有不同的解法，但是，基本的方法可归结为两大类，即为直接法和迭代法。 本章介绍的经典解法，都把原方程组化为一个或者两个三角形方程组来求解，主要包括Gauss4消去法和它的变形-直接三角分解法 设有线性方程组 Ax=b 其中 La1nb1x1a22La2nx2b2 x= b= MMMMMbxan2Lannnn 根据线性代数知识，当detA0时，方程组的解存在且唯一，对a11a21 A=Man1a12增广矩阵(A,b)=(A(1)b(1)施行行初等变换，化A(1)为上三角形矩阵A(n)，同时b(1)化为b(n)，这时与增广矩阵(A(n),b(n)相应的线性方程组为上三角形方程组 A(n)x=b(n)
3、 (1)a110其中 (A(n),b(n)=M0(1)a12La1(1n(2)(2)a22La2nM0MM(n)Lannb1(1)(2)b2 M(n)bna(!)11 (A(n)=0M0b1(!)(2)b (b(n)=2 Mb(n)n(!)a12)La1(!n(2)(2)a22La2nMMM(n)0Lann (i)设aii0(i=1.2.3Ln) 则的解为 它便是原方程组的解，实现上述求解过程的方法称为Gauss消去法 如果方程组的系数矩阵可分解为两个形式简单的三角形矩阵L和U的乘积，即 A=LU (1.4) 若L为下三角形矩阵，则U为上三角形矩阵，反之亦然 从而求解Ax=b的问题转化为解三角
4、形方程组 Ly=b (1.5) 和 Ux=y l110L0uu12L其中 L=l21lM02211u22LMMMM U=0MMMln1ln2Llnn00L则Ly=b为下三角形方程组，它的第i个方程为 li1y1+li2y2+L+liiyi=bi(i=1.2Ln) 假定lii0按y1,y2Lyn的顺序解得 yb11= l11i-1 -lijyj+bij=1yi=lij 上三角形方程组Ux=y的第i个方程为 uiixi+uii+1xi+1+L+uinxn=yi假定uii0按xn,xn-1,Lx1的顺序求解得 u1nu2nMunn (i=2,3Ln) (i=1,2Ln) ynx=nunnn -uij
5、xj+yij=i+1yi=uii 直接法的求解过程，在计算过程中无舍入误差的前提下，都可以经有限步算数运算而得到精确解。由于计算机的字长有限，初始数据取浮点数以及运算过程都不断地产生舍入误差，这些误差的传播和积累，会影响计算精度，所以，如何避免舍入误差的增长是设计算法时必须考虑的问题。 直接法通常需要存储系数矩阵的全部元素，当方程组的阶数很高时，需要相当大的内存空间，因此，在算法设计上应当注意节省内存，比如，对称矩阵可以只存其下三角形部分于一维数组中。 综上所述，如果矩阵A非奇异，总可以通过带有行交换或不带行交换的消元过程，将A化为非奇异上三角形矩阵A(n)因此，回代求解过程也可以进行到底，但
6、是，在实际应用中，常常难于事先判断系数矩阵A的奇异性，因此，需要进一步考虑当A为奇异矩阵时计算过程可能发生的情况。一是消元过程的某一(n)(i)步找不到非零的aii，于是计算中断；二是虽然消元过程能进行到底，但ann=0，使回代求解过程无法进行下去，因此，在计算设计中必须考虑到上述两种可能发生的情形，此时应在算法设计中给出计算中断的信息。 1 Gauss列主元素消去法 在Gauss消元过程中位于矩阵A(k)(k=1,2Ln)的主对角线(k,k)位置上的元素(k)称为主元素，因为在计算解的分量时，都做除数，但应当避免用小的数做除akk(k)(k)(k+1in)在数量级上相对小的主元akk数，即避
7、免用较aik做除数，以防止舍入误差的扩大，降低解得精度，下面的例子说明“小主元”对解得精度的影响。 例1用Gauss消去法解线性方程组 10-823x11 -13.7124.623x2=2 -21.0725.643x33用8位十进制尾数的浮点数计算 10-8(1)(1)解 (A,b)=(A,b)=-1-213.7124.6232 1.0725.643323(1)(1)a31a21128m21=(1)=-8=-10 m31=(1)=-8=-2108 a1110a1110(2)a22=fl(a22-m21a12)=fl(3.712+2108) =0.00000000109+0.2109=0.210
8、9 (2)b2=fl(b2-m21b1)=fl(2+108)=0.1109 (2)(2)(2)类似的算出a32我们看到，由于小主元10-8做除数，使得行乘数m21,m31,a33,b3(2)的数量级变得很大，这样在计算a22时，即fl(3.712)与fl(0.2109)在计算机上做代数和时，fl(3.712)由于阶码升为9使尾数右移变成了机器零，这里揭示了数量级相对于分子数做除数，会使Gauss消去法在计算机有限字长的环境下造成较大舍入误差的原因。 经第一步消元得 10-8(2)(2) (A,b)=00经第二步消元得 20.21090.410930.31090.610990.110 0.210
9、9110-8231(3)(3)999(A,b)=00.2100.3100.110 0000这一结果表明原方程组有无穷多个解，在计算机上，将会给出“无唯一解”的信息，停止计算，但实际上，由于detA0，所以原方程组有唯一解，是计算过程的舍入误差使解面目全非了，从前面的分析我们看到，导致这一错误的主要原因是用“小主元”做除数。 为避免小主元做除数，在Gauss消去法中加入选主元过程，即在第k步消元时，首先在A(k)的第k列主对角线元以下元素中挑选出绝对值最大者，并通过交换的第与行对应元素，使位于对角线上，仍记为，并称之为第步消元的列主元素，然后再进行消元计算， 第一步都按列选主元的消去法称为列主元
10、素消去法，由于列主元素满足条件 所以行乘数皆满足，因此，列主元素消去法能避免例1中所出现的问题，误差分析与数值试验表明，用列主元消去法解“良态”线性方程组，在绝大多数情况下，都可以得到令人满意的结果。 例2用列主元素消去法解例1的线性方程组Ax=b 10-8231解 (A,b)=-13.7124.6232 -21.0725.6433-2-110-81.0725.64333.7124.6232 231 m21=0.5,m31=-0.510-8,经第一步消元，得 1.072-2(A(2),b(2)=00.317610100.21015.6430.180151010.31010.5 0.11013，
11、经第二步消元得 m32=0.629722921.072-2(A(3),b(3)=00.3176101005.6430.180151010.186555411010.5 0.685138543回代求解得 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)
12、,b(k+1) 其中Iikk1O10L11 =MOM11L01O11O1Lk=-mk+1,k-mk+2,kM-mnk1 1O1(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
13、=U (2.4) 2令L2=Ii33L2Ii33 L1=Ii33Ii22L1Ii22Ii33 P=Ii33Ii2Ii11 2%=ILI与L的形状相同，也是一个初等下三角容易证明，当i,jk时，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
14、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以及单位下
15、三角形矩阵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)nn的各阶顺序主
16、子式均不为零，A有分解式 La1n1u11u12Lu1na22La2nl211uLu222nLU =MMMMOOMan2Lannln1ln2L1unn 对比等式左边和右边乘积矩阵LU的第r行主对角元右边的对应元素，得 a11a21Man1a12aij=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)假
17、设已求出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,
18、L,n)(3.8)(3.9) uri=ari-lrkukik=1r-1lrr(i=r+1,L,n;r=1,2,L,n-1)0规定lu=0ikkrk=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为上三角形矩阵，并且uii0(i=1,L,n)，根据
19、公式(1.8)和(1.10)得与Doolittle分解相应的求解方程组的公式为 y1=b1r-1yr=br-lriyi(i=2,L,n)i=1及 (3.10) ynx=,nunnn(3.11) yr-urixii=r+1xr=(r=n-1,L,1)urr类似地。可以得出Crout分解相应的求解出公式 b1y=1l11r-1br-lriyii=1yr=(i=2,L,n)l11及 (3.12) xn=yn,n(r=n-1,L,2,1)xr=yr-urixii=r+1(3.13) 综上所述，直接三角分解法的基本方法有Doolittle分解和Crout分解法，他们都包括分解系数矩阵和求解两个三角形方程
20、组的过程。 例1 用Doolittle法解方程组 0-3x110210x-3-4-12132=5 123-4x3-24149-13x47解 由公式(3.4)和(3.5)计算得 (u11,u12,u13,u14)=(a11,a12,a13,a14)=(2,10,0,-3) (l21,l31,l41)T31=-,2 22T对于r=2,3,4，用公式(3.6)和(3.7)计算得 (u22,u23,u24)=11,-12,T17,236-,-(l32,l42)=,111132-,-(u33,u34)=,1111l43=-9,u44=-4.T解Ly=b，由公式(3.11)计算得 7y=10,20,-,-
21、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
22、。在数组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的
23、解y。也就是说，系数矩阵的Doolittle分解与解方程组Ly=b可以同时进行，归结为对数组A=(A,b)的分解计算。然后再解Ux=y,称这种求解方式为紧凑格式。 设N=3,r表示分解的步数，分解过程中数组A的变化情况如下 a11a12a13b1a11a12a13a14 A=a21aaba222323a24a31a32a33b=A=2a21a322a31a32a 33a34u11u12u13u14u13u14u13u14u11r=1laau11u12r=224l21uuuu11u12r=321a22232324u23ul31a32a33a3422l31l32aal=24l21333421u22
24、l31l32u33u34l31分解过程可以缩写为 a11aa13b1u11u12u13u14r=1A=12a21aab22232l21u22u23u24a31a32a33br=2 3l31l32u33u34r=3例2 用紧凑格式的Doolittle法解例1中的方程组 解：分解 100-3102100-3102-311-121720r=12 A=-3-4-121352123-4-21-33217r=2 4149-137211-11-11-11r=32-6r=411-9-4-16u12u13u22u23l32u331-3L=21222100-31711-1212 ,U=332-1-1111116-91-4111020 y= 17-11-16解Ux=y，用公式计算得 x=(1234) T
