第三章线性方程组解法

1、计算方法 (力学系本科生)3.1 问题的提出 第三章 线性方程组解法 n阶线性方程组3.1 问题的提出11 11221331121 1222233221 12233.nnnnnnnnnnna xa xa xa xba xa xa xa xba xa xa xa xb3.1 问题的提出线性方程组Ax=b，其中A是n维维方阵，x是n维未知数向量，b是n维常数向量。1112111212222212.,.nnnnnnnnaaabxaaabxaaabxAbx3.1 问题的提出如果A是非奇异阵时，方程组有唯一解，且可以用克莱姆(Grammer)法则表示:, (1,2,., )iiDxinD其中xi是解向量

2、x*的第i个分量，D=detA, Di是用b代替A的第i列后得到矩阵的行列式。3.1 问题的提出克莱姆方法求解计算量太大，需要计算(n+1)个n阶行列式，共需要(n+1)!次乘法运算。3.1 问题的提出 求解线性方程组的数值方法有两大类：1)直接法(direct methods)。 经过有限次算术运算可求方程组精确解的方法(实际上，由于舍入误差不可避免，一般得不到精确解)。适合于求解低阶稠密阵方程组。3.1 问题的提出2) 迭代法(iterative methods)。采用极限过程去逐步逼近线性方程组精确解的方法。迭代法需要计算机存储单元较少，对计算机要求不高，程序设计简单，但有收敛性和收敛速

3、度方面的问题。迭代法是求解大型稀疏矩阵方程的重要方法。 3.1 问题的提出 我们在本章将要学习迭代法有：雅可比(Jacobi)迭代法高斯塞德尔(Gauss-Seidel)迭代法超松弛迭代法(Successive overrelaxation method, SOR)。3.1 问题的提出追赶法(Forward elimination and backward substitution)。 我们在本章将要学习直接解法有：高斯消去法(Gauss Elimination)，高斯主元素消去法(Gauss Elemination with pivoting),三角分解法(LU decomposition)

4、，3.1 问题的提出JacobiGauss Seidel迭代法迭代法迭代法超松弛迭代法不选主元高斯消去法选主元（列选，全选）直接法三角分解法追赶法【历史注记】线性代数方程组数值解法有着悠久的历史。我国古代数学著作九章算术(公元1世纪)的“方程”章中就有了较好的线性方程组数值解法相当于现代对方程组的增广矩阵进行初等变换、消去未知数的方法。中世纪的印度数学家也可以求解线性方程组。例如12世纪的婆什迦罗的著作中，也有求解线性方程组的内容。3.1 问题的提出在欧洲，16世纪的比特奥在其算术（1559）中采用了与九章算术类似的消元法。日本数学家关孝和在其解伏题之法一书(1683)中首先采用了类似于现代行

5、列式法求解了三元线性方程组。稍后，莱布尼茨提出关于行列式解线性方程组的思想(1693)。1721年马可劳林用行列式展开式的方法给出了二元、三元、四元线性方程组的解法，3.1 问题的提出但他的符号记法不完善。1750年，克莱姆给出了现在比较通用的线性方程组行列式解法，即克莱姆法则。1764年，贝祖用行列式建立了线性方程组的一般理论。但由于当时计算的效率很低，这一理论几乎只有理论的意义，实际上只能求出未知数很少的线性代数方程组的解。只是在20世纪中叶电子计算机问世并投入应用之后，大型线性代数方程组的数值求解才成为可能。3.1 问题的提出如何利用计算机更精确、更有效地求解大型线性方程组，是计算数学中

6、最重要的课题之一。 现代计算实践中，常用的线性代数方程组数值解法有直接法和迭代法两大类。直接法是在没有舍入误差的假设下，经过有限次运算就可得出方程组的精确解的方法，如各种消元法。迭代法则采用逐次逼近的方法,即从一个初始值出发，按照一定的计算格式(迭代公式)，构造一个向量的无穷序列，其极限才3.1 问题的提出是方程组的精确解，用有限次运算得不到精确解。迭代法是牛顿最先提出来的，1940年经司威尔提出的松弛法也是一种迭代法，共轭梯度法则是另一种迭代法，是弗莱彻等人于20世纪60年代提出来的。3.1 问题的提出例3.13.1 问题的提出52832026xyxy精确解为*2,1xy 0.41.60.1

7、51.3xyyx 将方程写为取(0)(0)0 xy(1)(0)(1)(0)0.41.61.60.151.31.3xyyx 3.1 问题的提出重复以上过程得(2)(1)(2)(1)0.41.62.120.151.31.06xyyx k0123456x(k)01.62.122.0241.99281.998562.000432y(k)0-1.3-1.06-0.982-0.9964-1.00108-1.0002163.1 问题的提出如果把原方程写为6.6678.6672.54.0 xyyx 构造(1)( )(1)( )6.6678.6672.54.0kkkkxyyx k0123x(k)08.66735

8、.335-109.126y(k)04.0-17.668-84.35852832026xyxy3.1 问题的提出例3.2 12312312382124102132516xxxxxxxxx构造(1)( )( )123(1)( )( )213(1)( )( )3120.1250.251.50.40.12.10.60.43.2kkkkkkkkkxxxxxxxxx (0)(0)(0)1230 xxx3.1 问题的提出得 3.1 问题的提出构造(1)( )( )123(1)( )( )213(1)( )( )3122.50.255.25(2)1.52.58.0(3)40.56.0(1)kkkkkkkkkx

9、xxxxxxxx (0)(0)(0)1230 xxx12312312382124102132516xxxxxxxxx由原方程3.1 问题的提出得3.1 问题的提出 如果A非奇异，则线性方程组Ax=b有唯一解x*，将方程化为x=Bx+f，给出初始向量x0，则有： xk+1=Bxk+f, k=0,1,2可以构成一向量序列xk，若向量序列xk收敛于x*，则x*=Bx*+f, 即即x*是方程组的解 。 这种方法称为迭代法, B称为迭代矩阵。 3.1 问题的提出 构造迭代法的中心问题是建立一个由本次近似值计算下一次近似值的规则。用迭代法求解线性方程组时要解决的问题有： 1) 构造一种迭代格式，由xk计算

10、xk+1 3) 证明向量序列xk的收敛性 2) 给出初始向量x04) 如果序列收敛，证明是原方程组的解 5) 给出估计误差和迭代停止判据。 3.1 问题的提出v 定义：在n维空间中给定一个向量序列 ， ，如果对每一个分量 ,当 时都有极限xi, 即 , 则称向量序列 有极限 ,或称 收敛于x。 kx 12(,.)kkkkTnxxxxkixk limkiikxx kx122( ,.)Tx xxx kx3.2 雅可比迭代(Jacobi iteration) 第五章 线性方程组解法3.2 雅可比迭代 最简单的迭代方法是从第i个方程解出未知数xi，i=1,2,n 1,1()niiijjjj iiixb

11、a xa于是雅可比迭代式为 111(),1,2,.,nkkiiijjjj iiixba xina 把系数矩阵分解为A=U+D+L，其中U为由A上三角部分构成的上三角阵，L为由A下三角元素构成的下三角阵，D为由A对角线元素构成的对角阵。 3.2 雅可比迭代 显然，所有aii, i=1,2,n不为零上式才有意义，从线性代数知，对于任何系数方阵非奇异的方程组，通过适当交换方程的顺序总可以使所有方程的0iia 3.2 雅可比迭代121312320.0.0.0nnaaaaaU21313212000.0nnaaaaaL112233.nnaaaaD于是原方程组为 (U+D+L)x=b3.2 雅可比迭代上式两

12、边左乘D-1得 x= D-1 b- D-1(U+L)x=Bx+f其中 B=-D-1(U+L), f=D-1b于是有迭代格式 xk+1=Bxk+f例3.3 用Jacobi迭代格式解下面方程组。3.2 雅可比迭代解：Jacobi迭代格式为 1238111210141153xxx 3.2 雅可比迭代取初始向量x(0)=(0,0,0)T，各次迭代结果(1)( )( )123(1)( )( )213(1)( )( )312(1)/8(42)/10( 3)/( 5)kkkkkkkkkxxxxxxxxx k0123456x1(k)0.00000.12500.25000.22630.22350.22510.2

13、250 x2(k)0.00000.40000.31500.30050.30600.30580.3056x3(k)0.0000-0.6000 -0.4950 -0.4870 -0.4946 -0.4941 -0.49483.3 高斯塞德尔迭代(Gauss-Seidel iteration ) 第五章 线性方程组解法 在雅可比迭代中, 计算第k+1次迭代近似值时用的是上一次即第k次的近似值，从式 3.3 高斯塞德尔迭代11,1(),1,2,.,nkkiiijjjj iiixba xina可以看出，在计算第i个xik+1分量时，前面i-1个分量x1k+1, x2k+1 xi-1k+1已经从上式中计算

14、出来了，于是很自然会想到如果把它们代入用来计算xik+1可能会改进迭代,于是就得到Gauss-Seidel迭代格式： 写成矩阵形式为： 3.3 高斯塞德尔迭代111111(),1,2.,inkkkiiijjijjjj iiixba xa xina 111()kkkxDbLxUx或1()kkDL xbUx如果(D+L)-1存在，则 111()()kkk xDLUxDLbBxf其中3.3 高斯塞德尔迭代11(),() BDLU fDLb【注记】 通常高斯塞得尔方法比雅可比方法有更快的收敛速度，但不是总这样，对于某些方程组，雅可比迭代收敛，而高斯塞得尔方法发散。即，并不是任何时候高斯塞得尔方法都比雅

15、可比方法好。 例3.4 用Gauss-Seidel迭代格式解下面方程组，精确到3位有效数。3.3 高斯塞德尔迭代1238111210141153xxx 解：Gauss-Seidel迭代格式如下3.3 高斯塞德尔迭代取初始近似值x0=(0,0,0)T，各次迭代结果112311213111312(1)/8(42)/10( 3)/( 5)kkkkkkkkkxxxxxxxxx k01234xk10.00000.12500.23440.22450.2250 xk20.00000.37500.30310.30590.3056xk30.0000-0.5000-0.4925-0.4939-0.49363.4

16、逐次超松弛迭代法(SOR) (Successive Overrelaxation Method) 第五章 线性方程组解法 逐次超松弛迭代简称SOR方法，是高斯塞得尔法的一种加速方法。 3.4 逐次超松弛迭代法高斯塞得尔法迭代格式得到 111111(),1,2,.,inkkkiiijjijjjj iiixba xa xina 把xik+1改进为xik与 的加权平均，即 1kix11111(1)(),1,2,.,kkkiiiinkkkiiijjijjjj iiixxxxba xa xina3.4 逐次超松弛迭代法 上式中 时, 就是高斯塞得尔方法，为保证迭代过程收敛,要求 102当 时叫低阶松弛法

17、；当 时叫超松弛法。 121 SOR方法收敛时，希望选择一个最佳的使收敛速度最快。目前还没有最佳松弛因子 的一般算法理论，实际大都由计算经验通过试算确定 的近似值 optopt 超松弛迭代式的矩阵形式 3.4 逐次超松弛迭代法1)直接由分量形式公式写 由1111()inkkkkiiiijjijjjj iiixxba xa xa有111()kkkkkxxDbLxUxDx所以111() ()()kkxDLDDU xDLb(证明)3.4 逐次超松弛迭代法111()kkkkkxxDbLxUxDx111() ()()kkxDLDDU xDLb11()kkkkkDxDxbLxUxDx11kkkkkDxLx

18、DxUxDxb1()()kkDL xDUD xb2)由高斯塞德尔公式推导。 3.4 逐次超松弛迭代法高斯塞德尔迭代公式的矩阵形式是 11kkkDxbLxUx加权平均11(1)kkkxxx11(1)()kkkkDxDxbLxUx111() (1)()kkxDLDU xDLb(证明)3.4 逐次超松弛迭代法11(1)()kkkkDxDxbLxUx111() (1)()kkxDLDU xDLb11(1)kkkkDxLxDxUxb1() (1)kkDL xDU xb3.5 迭代法的收敛性 (convergence) 第五章 线性方程组解法v 矩阵 的特征值 的绝对值最大值称为矩阵A的谱半径，即 3.5

19、 迭代法的收敛性n nRA,1,2,.,iin1( )maxii n A 定理 (迭代法基本定理)：设有方程组x=Bx+f，对于任意初始向量x0及任意f，迭代公式xk+1=Bxk+f收敛的充要条件是 ( )1B3.5 迭代法的收敛性 定理(迭代收敛的充分条件)：设有迭代式xk+1=Bxk+f，如果 ，则对于任意初始向量x0，这个迭代过程收敛于方程组x=Bx+f的唯一解x*，并且有事后估计 1qB*111kkkqxxxx以及事前估计 *101kkqqxxxx3.5 迭代法的收敛性v定义：如果对于方阵 ，有 n nA1,1,2,.,niiijjj iaain则称方阵对角占优。 定理：如果方程组Ax

20、=b的系数阵对角占优，则方程组有唯一解且对任意初始向量x0雅可比迭代和高斯塞德尔迭代都收敛于真解。 3.5 迭代法的收敛性【思考题3.1】如何对方程组进行调整，使用Gauss-Seidel迭代格式求解时收敛？1213123879897xxxxxxx 定理：如果方程组Ax=b的系数阵对称正定，则方程组有唯一解且对任意初始向量x0高斯塞德尔迭代收敛于真解。 3.5 迭代法的收敛性 Jacobi迭代格式的收敛性3.5 迭代法的收敛性 Jacobi迭代矩阵1() JDUL 特征方程0IJ1()0IDUL1()0DDUL3.5 迭代法的收敛性由于10D0DUL所以注意到ADUL所以，将A的对角线元素乘以

21、 后取行列式，令其等于零，就是Jacobi迭代矩阵的特征方程。 例3.5 讨论用Jacobi迭代格式解方程组的收敛性。3.5 迭代法的收敛性1238111210141153xxx解：Jacobi迭代矩阵的特征方程为81121010115展开得3.5 迭代法的收敛性Jacobi迭代格式收敛。34001230解得12,30.146084,0.073041 0.214487i 123( )max,0.226581B Gauss-Seidel迭代格式的收敛性3.5 迭代法的收敛性 Gauss-Seidel迭代矩阵1() GDLU 特征方程0IG1()0IDLU1()()0DLDLU3.5 迭代法的收敛

22、性由于1()0DL()0DLU所以注意到ADUL所以，将A的对角线以下元素乘以 后取行列式，令其等于零，就是Gauss-Seidel迭代矩阵的特征方程。3.5 迭代法的收敛性【思考题3.2】Gauss-Seidel迭代矩阵1() GDLU至少有1个特征值为零，为什么？ 例3.6 讨论用Gauss-Seidel迭代格式解方程组的收敛性。3.5 迭代法的收敛性1238111210141153xxx解： Gauss-Seidel迭代矩阵特征方程为811210105展开得3.5 迭代法的收敛性Gauss-Seidel迭代格式收敛。2(400101)0解得12,30,( 117)/80 123()max

23、,0.0640381G3.5 迭代法的收敛性!作业：(1) 写出下面方程组的Jacobi迭代格式和Gauss-Seidel迭代格式的算法 (2) 讨论它们的收敛性。123211111111121xxx 3.6高斯消去法(Gauss Elimination)第三章 线性方程组解法 高斯消去法是最古老的数值方法之一，现在仍然是一个很有用的方法，它在计算机上容易实现。 3.6 高斯消去法 其基本思想是在各个方程之间进行乘法和加减运算，逐步消去方程中的未知数。它分为消去和回代两个过程。 Axb 给定线性方程组3.6 高斯消去法1112111212222212.nnnnnnnnaaaxbaaaxbaaa

24、xb第一步消元，令 ，用 乘第一个方程再加到第i个方程上作为第i个新方程，消去x1的项，变为1111/,2,3,.,iimaain1 im逐次进行同样过程，最后，经过n-1次消元，得到： 3.6 高斯消去法1112111(1)(1)(1)22222(1)(1)(1)2.0.0.nnnnnnnaaaxbaaxbaaxb以上这些步骤叫消元过程。 3.6 高斯消去法1112111(1)(1)(1)22222(1)(1).0.0.00.nnnnnnnnaaaxbaaxbaxb然后，从第n个方程开始，依次解出 xn, xn-1, x1 3.6 高斯消去法(1)(1)/nnnnnnxba(1)(1)(1)

25、1()/niiiiiijjiij ixbaxa 1,.,1in高斯消去法的计算量3.6 高斯消去法 消去过程的计算量。第一步计算乘数mi1,(i=2,3,n)需要i-1次除法运算，计算aij(2)(i,j=2,3,n)需要(i-1)2次乘法运算以及(i-1)2次加减法运算。 利用求和公式211(1)/2,(1)(21)/6,1nniiin nin nnn 得到3.6 高斯消去法第k步加减法次数乘法次数除法次数1(n-1)2(n-1)2n-12(n-2)2(n-2)2n-2N-1111合计n(n-1)(2n-1)/6 n(n-1)(2n-1)/6n(n-1)/2 于是消去过程乘除法次数为2(1)

26、/3n n 消去过程加减法次数为(1)(21)/6n nn3.6 高斯消去法计算b(n-1)的计算量(1)(2).2 1(1)/2nnn n 乘除法次数为 加减法次数为(1)/2n n3.6 高斯消去法回代计算量 乘除法次数为 加减法次数为(1)/2n n12.(2)(1)(1)/2nnnn n3.6 高斯消去法总的计算量332333nnnn 乘除法次数为 加减法次数为3353263nnnn 高斯消去法还有没有办法进行，为什么？ 3.6 高斯消去法(1)0,1,2,.,(1)iiiain【思考题3.3】如果遇到! 作业：写出高斯消去法Fortran程序。3.7高斯主元消去法(Gauss Eli

27、mination with pivoting )第三章 线性方程组解法3.7 高斯主元消去法 从高斯消去法我们已经看出，为使高斯消去法能顺利进行，必须在每一步消去步都满足条件 ，但若(1)0iiia(1)(1),iiiikiaaki相应的系数 (1)(1)/,iikikiiimaaki计算可能引起很大的舍入误差。 为此，需要改进高斯消去法。有两种改进方法：列选主元，完全选主元。 3.7 高斯主元消去法u列选主元通过交换方程而使得 aki(i-1), k=i,i+1,n,中绝对值最大的一个换到(i, i)位置而成为第i步的消去主元。u完全选主元法就是在系数矩阵的子块 1,11,1,.kkknn

28、kn naaaa中找出绝对值最大的元素，作为第k+1次消去过程的主元。 3.7 高斯主元消去法假设 11maxpqijki nkj naa 则k+1行与p行交换，第k+1列与第q列交换。右端也同时交换，在做列交换时，要注意未知量也作交换，即把xk+1与xp交换。 3.7 高斯主元消去法 列选主元算法： 1. for k=1,2,n-1 2. 找出满足 元素的行位置p ,maxp kikk i naa 3. if , error, stop4. if ( )换行 ,p kapk5. for j=k,k+1,n 6. ak,j与ap,j交换 3.7 高斯主元消去法7. bk与bp交换 8. End

29、for9. Endif 消去计算： /,1,.,ikikkkmaaikn,1,.,ijijikkjaam ai jkn,1,.,iiikkbbm bikn10. Endfor3.7 高斯主元消去法回代计算：/nnnnxba11. for i=n-1,n-2,2,1 1()/niiijjiij ixba xa 12. Endfor 3.7 高斯主元消去法3.8 三角分解法(LU decomposition )第三章 线性方程组解法高斯消去法 3.7 三角分解法第一次消元过程为 111111/,2,3,., ,0iimaaina用-mi1乘上面第一个方程再加到第i个方程，即可消去第i个方程中的未知

30、量x1。这个过程实际上是给系数矩阵A左乘这样一个下三角阵L1： 作第二次消元相当于给系数矩阵A(1)左乘了一个下三角阵L2：211311100.010.00100.00.1nmmm L3.7 三角分解法作第k次消元相当于系数矩阵A(k-1)左乘一个下三角阵Lk2322100.0010.001.0.00.1nmmL(1)(1)(1)222222/,2,., ,0iimaain a3.7 三角分解法|11,2,100.010.1.0.0.0.kKkkkkkn kmmmL3.7 三角分解法于是第n-1次消元过程可表示为：( )1221.nnnLLL L AA( )1221.nnnLLL L bb于是

31、1111( )1221.nnnAL LLLA下三角阵乘积也是一个下三角阵。 记111121.nLL LL3.7 三角分解法 L是单位下三角阵(即对角元素全为1的下三角阵)。 211210.01.0.1nnlllL3.7 三角分解法( )nALA叫矩阵A的三角分解，或LU分解。如果L为单位下三角阵，则叫Doolittle分解，若U为单位上三角阵，则叫Crout分解。 ALU只要A的各顺序主子式不为零，则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积。3.7 三角分解法 设Ax=b，A=LU，则Ax=LUx=b于是令Ux=y，则Ly=b11121212221210.0.1.00.,.0.1

32、00.nnnnnnuuuluulluLU3.7 三角分解法这样原来方程能化为两个简单方程组由 求y11,1,2,.,iiiijjjybl yin然后由 求x1()/niiijjiij ixyu xu 3.7 三角分解法,1,.,1in n 现在显式地给出lij和uij的表达式 111211112121222212221212.1.10.1nnnnnnnnnnnnaaauuuaaaluuaaallu从第一行可以看出，u1j=a1j, j=1,2,n3.7 三角分解法从第一列可以看出 ai1=li1u11, i=2,3,n, 即li1=ai1/u11, i=2,3,n 。 从第二行可以看出，a2j

33、=l21u1j+u2j , j=2,3,n，则u2j=a2j-l21u1j, j=2,3,n 从第二列可以看出，ai2=li1u12+li2u22 , i=2,3,n，则li2=(ai2-li1u12)/u22, i=2,3,n3.7 三角分解法111211112121222212221212.1.10.1nnnnnnnnnnnnaaauuuaaaluuaaallu 一般地，如果U的前k-1行，L的前k-1行已经求出，则第k (k=2,3,n)步 11kkjkmmjkjmal uu(行)3.7 三角分解法11,1,.kkjkjkmmjmual ujk kn(列)11kikimmkikkkmal

34、 ul u11()/,1,.kikikimmkkkmlal uuikn直到第n步，A全部分解成LU。3.7 三角分解法3.9 追赶法( Forward elimination and backward substitution )第三章 线性方程组解法 在很多实际问题中，方程组Ax=b的系数矩阵A是一个稀疏矩阵。 3.6 追赶法11122221110.0000000.0000000.00.00000000.00000000.000000000iiinnnnnnbcxabcxabcabcabx12.nddd假设非零元素集中分布在对角线及其相邻两个次对角线上，且系数阵为对角占优阵，即有：11,2,

35、3,.1iiinnbcbacinba3.6 追赶法把系数矩阵三角分解有11222110.010000.0100,.0000010000nnnnqrqrpqrpqLU3.6 追赶法利用LU分解公式，写出111111,2,3,.,kkkkkkkkkqb rcrcp qaqp cb kn3.6 追赶法得到 11,kkqbrc11,2,3,.,kkkkkkkapqbp cknq于是 ,LydUxy得到 111,2,3,.,kkkkydydp ykn3.6 追赶法1/1(),1,2,.,1nnnkkkkkxyqxyc xknnq3.6 追赶法3.10 其它应用第三章 线性方程组解法 计算行列式值3.6 其他应用ALU1122det( )det( ).nnu uuAU 求矩阵逆1AAI111212122212.1000.0100.0001nnnnnnxxxxxxxxxA分别由 120.,1,2,.,.10iinixxinix A解出 xij, ij=1,2,n 于是1ijxA3.6 其他应用3.11 误差分析(Error analysis)第三章 线