计算方法 第三章 线性方程组解法.ppt

1、计算方法 (力学系本科生),第三章 线性方程组解法 (Solution for Linear Algebraic Equations ),3.1 问题的提出,第三章 线性方程组解法,n阶线性方程组,3.1 问题的提出,3.1 问题的提出,线性方程组Ax=b，其中A是n维方阵，x是n维未知数向量，b是n维常数向量。,3.1 问题的提出,如果A是非奇异阵时，方程组有唯一解，且可以用克莱姆(Grammer)法则表示:,其中xi是解向量x*的第i个分量，D=detA, Di是用b代替A的第i列后得到矩阵的行列式。,3.1 问题的提出,克莱姆方法求解计算量太大，需要计算(n+1)个n阶行列式，共需要(n

2、+1)!次乘法运算。,3.1 问题的提出,求解线性方程组的数值方法有两大类：,直接法(direct methods)。 经过有限次算术运算可求方程组精确解的方法(实际上，由于舍入误差不可避免，一般得不到精确解)。适合于求解低阶稠密阵方程组。,3.1 问题的提出,2) 迭代法(iterative methods)。采用极限过程去逐步逼近线性方程组精确解的方法。迭代法需要计算机存储单元较少，对计算机要求不高，程序设计简单，但有收敛性和收敛速度方面的问题。迭代法是求解大型稀疏矩阵方程的重要方法。,3.1 问题的提出,我们在本章将要学习迭代法有：,雅可比(Jacobi)迭代法,高斯塞德尔(Gauss-

3、Seidel)迭代法,超松弛迭代法(Successive overrelaxation method, SOR)。,3.1 问题的提出,追赶法(Forward elimination and backward substitution)。,我们在本章将要学习直接解法有：,高斯消去法(Gauss Elimination)，,高斯主元素消去法(Gauss Elemination with pivoting),三角分解法(LU decomposition)，,3.1 问题的提出,【历史注记】线性代数方程组数值解法有着悠久的历史。我国古代数学著作九章算术(公元1世纪)的“方程”章中就有了较好的线性方程

4、组数值解法相当于现代对方程组的增广矩阵进行初等变换、消去未知数的方法。中世纪的印度数学家也可以求解线性方程组。例如12世纪的婆什迦罗的著作中，也有求解线性方程组的内容。,3.1 问题的提出,在欧洲，16世纪的比特奥在其算术（1559）中采用了与九章算术类似的消元法。日本数学家关孝和在其解伏题之法一书(1683)中首先采用了类似于现代行列式法求解了三元线性方程组。稍后，莱布尼茨提出关于行列式解线性方程组的思想(1693)。1721年马可劳林用行列式展开式的方法给出了二元、三元、四元线性方程组的解法，,3.1 问题的提出,但他的符号记法不完善。1750年，克莱姆给出了现在比较通用的线性方程组行列式

5、解法，即克莱姆法则。1764年，贝祖用行列式建立了线性方程组的一般理论。但由于当时计算的效率很低，这一理论几乎只有理论的意义，实际上只能求出未知数很少的线性代数方程组的解。只是在20世纪中叶电子计算机问世并投入应用之后，大型线性代数方程组的数值求解才成为可能。,3.1 问题的提出,如何利用计算机更精确、更有效地求解大型线性方程组，是计算数学中最重要的课题之一。 现代计算实践中，常用的线性代数方程组数值解法有直接法和迭代法两大类。直接法是在没有舍入误差的假设下，经过有限次运算就可得出方程组的精确解的方法，如各种消元法。迭代法则采用逐次逼近的方法,即从一个初始值出发，按照一定的计算格式(迭代公式)

6、，构造一个向量的无穷序列，其极限才,3.1 问题的提出,是方程组的精确解，用有限次运算得不到精确解。迭代法是牛顿最先提出来的，1940年经司威尔提出的松弛法也是一种迭代法，共轭梯度法则是另一种迭代法，是弗莱彻等人于20世纪60年代提出来的。,3.1 问题的提出,例3.1,3.1 问题的提出,精确解为,将方程写为,取,3.1 问题的提出,重复以上过程得,3.1 问题的提出,如果把原方程写为,构造,3.1 问题的提出,例3.2,构造,3.1 问题的提出,得,3.1 问题的提出,构造,由原方程,3.1 问题的提出,得,3.1 问题的提出,如果A非奇异，则线性方程组Ax=b有唯一解x*，将上式化为x=

7、Bx+f，给出初始向量x0，则有：,xk+1=Bxk+f, k=0,1,2,可以构成一向量序列xk，若向量序列xk收敛于x*，则x*=Bx*+f, 即x*是方程组的解 。,这种方法称为迭代法, B称为迭代矩阵。,3.1 问题的提出,构造迭代法的中心问题是建立一个由本次近似值计算下一次近似值的规则。用迭代法求解线性方程组时要解决的问题有：,构造一种迭代格式，由xk计算xk+1,证明向量序列xk的收敛性,给出初始向量x0,如果序列收敛，证明是原方程组的解,给出估计误差和迭代停止判据。,3.1 问题的提出,定义：在n维空间中给定一个向量序列 ， ，如果对每一个分量 ,当 时都有极限xi, 即 , 则

8、称向量序列 有极限 ,或称 收敛于x。,3.2 雅可比迭代 (Jacobi iteration),第五章 线性方程组解法,3.2 雅可比迭代,最简单的迭代方法是从第i个方程解出未知数xi，i=1,2,n,于是雅可比迭代式为,把系数矩阵分解为A=U+D+L，其中U为由A上三角部分构成的上三角阵，L为由A下三角元素构成的下三角阵，D为由A对角线元素构成的对角阵。,3.2 雅可比迭代,显然，所有aii, i=1,2,n不为零上式才有意义，从线性代数知，对于任何系数方阵非奇异的方程组，通过适当交换方程的顺序总可以使所有方程的,3.2 雅可比迭代,于是原方程组为 (U+D+L)x=b,3.2 雅可比迭代

9、,上式两边左乘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迭代格式为,3.2 雅可比迭代,取初始向量x(0)=(0,0,0)T，各次迭代结果,3.3 高斯塞德尔迭代 (Gauss-Seidel iteration ),第五章 线性方程组解法,在雅可比迭代中, 计算第k+1次迭代近似值时用的是上一次即第k次的近似值，从式,3.3 高斯塞德尔迭代,可以看出，在计算第i个xik+1分量时，前面i-1个分量x1k+1,

10、x2k+1 xi-1k+1已经从上式中计算出来了，于是很自然会想到如果把它们代入用来计算xik+1可能会改进迭代,于是就得到Gauss-Seidel迭代格式：,写成矩阵形式为：,3.3 高斯塞德尔迭代,或,如果(D+L)-1存在，则,其中,3.3 高斯塞德尔迭代,【注记】 通常高斯塞得尔方法比雅可比方法有更快的收敛速度，但不是总这样，对于某些方程组，雅可比迭代收敛，而高斯塞得尔方法发散。即，并不是任何时候高斯塞得尔方法都比雅可比方法好。,例3.4 用Gauss-Seidel迭代格式解下面方程组，精确到3位有效数。,3.3 高斯塞德尔迭代,解：Gauss-Seidel迭代格式如下,3.3 高斯塞

11、德尔迭代,取初始近似值x0=(0,0,0)T，各次迭代结果,3.4 逐次超松弛迭代法(SOR) (Successive Overrelaxation Method),第五章 线性方程组解法,逐次超松弛迭代简称SOR方法，是高斯塞得尔法的一种加速方法。,3.4 逐次超松弛迭代法,高斯塞得尔法迭代格式得到,把xik+1改进为xik与 的加权平均，即,3.4 逐次超松弛迭代法,上式中 时, 就是高斯塞得尔方法，为保证迭代过程收敛,要求,当 时叫低阶松弛法； 当 时叫超松弛法。,SOR方法收敛时，希望选择一个最佳的使收敛速度最快。目前还没有最佳松弛因子 的一般算法理论，实际大都由计算经验通过试算确定

12、的近似值,超松弛迭代式的矩阵形式,3.4 逐次超松弛迭代法,直接由分量形式公式写,由,有,所以,(证明),3.4 逐次超松弛迭代法,由高斯塞德尔公式推导。,3.4 逐次超松弛迭代法,高斯塞德尔迭代公式的矩阵形式是,加权平均,(证明),3.4 逐次超松弛迭代法,3.5 迭代法的收敛性 (convergence),第五章 线性方程组解法,矩阵 的特征值 的绝对值最大值称为矩阵A的谱半径，即,3.5 迭代法的收敛性,定理 (迭代法基本定理)：设有方程组x=Bx+f，对于任意初始向量x0及任意f，迭代公式xk+1=Bxk+f收敛的充要条件是,3.5 迭代法的收敛性,定理(迭代收敛的充分条件)：设有迭代

13、式xk+1=Bxk+f，如果 ，则对于任意初始向量x0，这个迭代过程收敛于方程组x=Bx+f的唯一解x*，并且有事后估计,以及事前估计,3.5 迭代法的收敛性,定义：如果对于方阵 ，有,则称方阵对角占优。,定理：如果方程组Ax=b的系数阵对角占优，则方程组有唯一解且对任意初始向量x0雅可比迭代和高斯塞德尔迭代都收敛于真解。,3.5 迭代法的收敛性,【思考题3.1】如何对方程组进行调整，使用Gauss-Seidel迭代格式求解时收敛？,定理：如果方程组Ax=b的系数阵对称正定，则方程组有唯一解且对任意初始向量x0高斯塞德尔迭代收敛于真解。,3.5 迭代法的收敛性,Jacobi迭代格式的收敛性,3

14、.5 迭代法的收敛性,Jacobi迭代矩阵,特征方程,3.5 迭代法的收敛性,由于,所以,注意到,所以，将A的对角线元素乘以 后取行列式，令其等于零，就是Jacobi迭代矩阵的特征方程。,例3.5 讨论用Jacobi迭代格式解方程组的收敛性。,3.5 迭代法的收敛性,解：Jacobi迭代矩阵的特征方程为,展开得,3.5 迭代法的收敛性,Jacobi迭代格式收敛。,解得,Gauss-Seidel迭代格式的收敛性,3.5 迭代法的收敛性,Gauss-Seidel迭代矩阵,特征方程,3.5 迭代法的收敛性,由于,所以,注意到,所以，将A的对角线以下元素乘以 后取行列式，令其等于零，就是Gauss-S

15、eidel迭代矩阵的特征方程。,3.5 迭代法的收敛性,【思考题3.2】Gauss-Seidel迭代矩阵,至少有1个特征值为零，为什么？,例3.6 讨论用Gauss-Seidel迭代格式解方程组的收敛性。,3.5 迭代法的收敛性,解： Gauss-Seidel迭代矩阵特征方程为,展开得,3.5 迭代法的收敛性,Gauss-Seidel迭代格式收敛。,解得,3.5 迭代法的收敛性,作业：(1) 写出下面方程组的Jacobi迭代格式和Gauss-Seidel迭代格式的算法 (2) 讨论它们的收敛性。,3.6高斯消去法 (Gauss Elimination),第三章 线性方程组解法,高斯消去法是最古

16、老的数值方法之一，现在仍然是一个很有用的方法，它在计算机上容易实现。,3.6 高斯消去法,其基本思想是在各个方程之间进行乘法和加减运算，逐步消去方程中的未知数。它分为消去和回代两个过程。,给定线性方程组,3.6 高斯消去法,第一步消元，令 ，用 乘第一个方程再加到第i个方程上作为第i个新方程，消去x1的项，变为,逐次进行同样过程，最后，经过n-1次消元，得到：,3.6 高斯消去法,以上这些步骤叫消元过程。,3.6 高斯消去法,然后，从第n个方程开始，依次解出 xn, xn-1, x1,3.6 高斯消去法,高斯消去法的计算量,3.6 高斯消去法,消去过程的计算量。第一步计算乘数mi1,(i=2,

17、3,n)需要i-1次除法运算，计算aij(2)(i,j=2,3,n)需要(i-1)2次乘法运算以及(i-1)2次加减法运算。,利用求和公式,得到,3.6 高斯消去法,于是消去过程乘除法次数为,消去过程加减法次数为,3.6 高斯消去法,计算b(n-1)的计算量,乘除法次数为,加减法次数为,3.6 高斯消去法,回代计算量,乘除法次数为,加减法次数为,3.6 高斯消去法,总的计算量,乘除法次数为,加减法次数为,高斯消去法还有没有办法进行，为什么？,3.6 高斯消去法,【思考题3.3】如果遇到,作业：写出高斯消去法的Fortran程序。,3.7高斯主元消去法 (Gauss Elimination wi

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

19、端也同时交换，在做列交换时，要注意未知量也作交换，即把xk+1与xp交换。,3.7 高斯主元消去法,列选主元算法：,for k=1,2,n-1,找出满足 元素的行位置p,if , error,无唯一解，stop,if ( )换行,for j=k,k+1,n,ak,j与ap,j交换,3.7 高斯主元消去法,bk与bp交换,Endfor,Endif,消去计算：,Endfor,3.7 高斯主元消去法,回代计算：,for i=n-1,n-2,2,1,Endfor,3.7 高斯主元消去法,3.8 三角分解法 (LU decomposition ),第三章 线性方程组解法,高斯消去法,3.7 三角分解法,

20、第一次消元过程为,用-mi1乘上面第一个方程再加到第i个方程，即可消去第i个方程中的未知量x1。,这个过程实际上是给系数矩阵A左乘这样一个下三角阵L1：,作第二次消元相当于给系数矩阵A(1)左乘了一个下三角阵L2：,3.7 三角分解法,作第k次消元相当于系数矩阵A(k-1)左乘一个下三角阵Lk,3.7 三角分解法,3.7 三角分解法,于是第n-1次消元过程可表示为：,于是,下三角阵乘积也是一个下三角阵。,记,3.7 三角分解法,L是单位下三角阵(即对角元素全为1的下三角阵)。,3.7 三角分解法,叫矩阵A的三角分解，或LU分解。如果L为单位下三角阵，则叫Doolittle分解，若U为单位上三角

21、阵，则叫Crout分解。,只要A的各顺序主子式不为零，则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积。,3.7 三角分解法,设Ax=b，A=LU，则Ax=LUx=b,于是令Ux=y，则Ly=b,3.7 三角分解法,这样原来方程能化为两个简单方程组,由 求y,然后由 求x,3.7 三角分解法,现在显式地给出lij和uij的表达式,从第一行可以看出，u1j=a1j, j=1,2,n,3.7 三角分解法,从第一列可以看出 ai1=li1u11, i=2,3,n, 即li1=ai1/u11, i=2,3,n 。,从第二行可以看出，a2j=l21u1j+u2j , j=2,3,n，则u2j=a

22、2j-l21u1j, j=2,3,n,从第二列可以看出，ai2=li1u12+li2u22 , i=2,3,n，则li2=(ai2-li1u12)/u22, i=2,3,n,3.7 三角分解法,一般地，如果U的前k-1行，L的前k-1行已经求出，则第k (k=2,3,n)步,(行),3.7 三角分解法,(列),直到第n步，A全部分解成LU。,3.7 三角分解法,3.9 追赶法 ( Forward elimination and backward substitution ),第三章 线性方程组解法,在很多实际问题中，方程组Ax=b的系数矩阵A是一个稀疏矩阵。,3.6 追赶法,假设非零元素集中分

23、布在对角线及其相邻两个次对角线上，且系数阵为对角占优阵，即有：,3.6 追赶法,把系数矩阵三角分解有,3.6 追赶法,利用LU分解公式，写出,3.6 追赶法,得到,于是,得到,3.6 追赶法,3.6 追赶法,3.10 其它应用,第三章 线性方程组解法,计算行列式值,3.6 其他应用,求矩阵逆,分别由,解出 xij, ij=1,2,n,于是,3.6 其他应用,3.11 误差分析 (Error analysis),第三章 线性方程组解法,3.6 误差分析,定义： 设 为方程组Ax=b的精确解，近似解与精确解之间的误差向量定义为,一种衡量近似解精确度的方法是把近似解 代入方程组，看残差向量的大小。,于是,3.6 误差分析,由范数性质有,而,上面两个式子相乘得,3.6 误差分析