第六章-方程组的数值解法

计算方法,方程组的数值解法,6,1、高斯消去法2、选主元素的高斯消去法3、矩阵的三角分解4、方程组的迭代解法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,主要知识点,第2页，共42页,引言（一）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,快速、高效地求解线性方程组是数值线性代数研究中的核心问题，也是目前科学计算中的重大研究课题之一。,各种各样的科学和工程问题，往往最终都要归结为求解一个线性方程组。,线性方程组的数值解法有：直接法和迭代法。,直接法：在假定没有舍入误差的情况下，经过有限次运算可以求得方程组的精确解；,迭代法：从一个初始向量出发，按照一定的迭代格式，构造出一个趋向于真解的无穷序列。,第3页，共42页,在没有舍入误差的情况下，经过有限次运算可以得到方程组的精确解的方法。,求解,Cramer法则:,所需乘除法的运算量大约为(n+1)!+n,n=20时，每秒1亿次运算速度的计算机要算30多万年！,直接法,引言（二）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第4页，共42页,举例（一）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,解：,例：直接法解线性方程组,第5页，共42页,高斯消去法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,高斯消去法的主要思路：将系数矩阵A化为上三角矩阵，然后回代求解。,考虑n阶线性方程组：,通过初等变换逐步消去未知元，将原方程组化为同解的三角方程组。,第6页，共42页,第一步：消去第一列,依次将增广矩阵的第i行+mi1第1行，得,设，计算,其中,第二步：消去第二列,依次将上述矩阵的第i行+mi2第2行，得,其中,设，计算,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第7页，共42页,高斯消去法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第k步：消去第k列,依此类推，直到第n-1步，原方程化为,设，计算,回代求解：,(i=n-1,1),第8页，共42页,例：用基本Gauss消元法求解下列方程组,解：,增广矩阵,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第9页，共42页,乘除运算量,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,乘除运算量：,由于计算机中做乘除运算的时间远远超过做加减运算时间，故估计运算量时，往往只估计乘除的次数。,第11页，共42页,基本Gauss消元法的工作量,消元过程：,回代过程：,加减法的次数,乘除法的次数,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第12页，共42页,列主元高斯消去法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,高斯消去法有效的条件是：,如果某个主元为0，则可能导致高斯消去法求解失败。,列主元高斯消去法：在第k步消元时,先选取列主元：,ifikkthen交换第k行和第ik行；,列主元高斯消去法比普通高斯消去法要多一些比较运算，但比普通高斯消去法稳定。,消元,第13页，共42页,全主元高斯消去法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,列交换改变了xi的顺序，须记录交换次序，解完后再换回来。,全主元高斯消去法具有很好的稳定性，但选全主元比较费时，故在实际计算中很少使用。,第14页，共42页,举例（二）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,解：,例：采用十进制八位浮点数，分别用Gauss消去法和列主元Gauss消去法求解线性方程组:,Gauss消去法：,9个,第15页，共42页,举例（二）续,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,列主元Gauss消去法：,例：,但列主元高斯消去法有时也会导致求解失败。,第16页，共42页,矩阵三角分解,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,将一个矩阵分解成结构简单的三角形矩阵的乘积称为矩阵的三角分解。,高斯消去过程其实就是一个矩阵的三角分解过程。,第17页，共42页,则A(k)与A(k+1)之间的关系式可以表示为：,将Gauss消去过程中第k-1步消元后的系数矩阵记为：,(k=1,n-1),2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第18页，共42页,LU分解,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,记：，则,其中：L-单位下三角矩阵，U-上三角矩阵,于是有：,(杜利脱尔Doolittle分解),第19页，共42页,LU分解的存在唯一性,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,LU分解存在,普通高斯消去法不被中断,?,定理,普通高斯消去法求解方程组Ax=b时，主元的充要条件是：A的所有顺序主子式不为零。,定理,若A的所有顺序主子式，则A存在唯一的LU分解。,(LU分解的唯一性),证：存在性由上面的定理可得；唯一性可用反证法证明。,第20页，共42页,类似分解,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,定理,若A的所有顺序主子式，则,第21页，共42页,LU分解,LU分解,LU分解紧凑方式,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,直接利用矩阵乘法来计算LU分解,比较等式两边的第一行得：,u1j=a1j,比较等式两边的第一列得：,比较等式两边的第二行得：,比较等式两边的第二列得：,(j=1,n),(i=2,n),(j=2,n),(i=3,n),第25页，共42页,LU分解紧凑算法（续）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,第k步：此时U的前k-1行和L的前k-1列已经求出,比较等式两边的第k行得：,比较等式两边的第k列得：,直到第n步，便可求出矩阵L和U的所有元素。,(j=k,n),(i=k+1,n),第26页，共42页,LU分解紧凑算法（续）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,算法5.1：(LU分解),Fork=1,2,.,n,EndFor,j=k,n,i=k+1,n,运算量：(n3-n)/3,第27页，共42页,LU分解求解线性方程组,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,两次回代过程求出方程组的解：,运算量：n2,LU分解,总运算量：,第28页，共42页,举例（一）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,解：,例：用LU分解求线性方程组Ax=b的解，其中,令A=LU，由LU分解算法6.1可得,回代：解Ly=b得：y=2,8,18,24T,解Ux=y得：x=-1,1,-1,1T,第29页，共42页,迭代法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,直接法的缺点：,迭代法：从一个初始向量出发，按照一定的迭代格式，构造出一个趋向于真解的无穷序列,运算量大，不适合大规模的线性方程组求解,无法保持系数矩阵的稀疏性,迭代解法是目前求解大规模线性方程组的主要方法。,只需存储系数矩阵中的非零元素,运算量不超过O(kn2)，其中k为迭代步数,(1)迭代格式的建立,(3)误差估计和收敛速度,(2)收敛性判断,第30页，共53页,Hilbert矩阵的病态性向量范数与矩阵范数矩阵的条件数,小插曲,引例.Hilbert矩阵的病态性,方程组Ax=b1的解为x1,方程组Ax=b的解为x,xx1=-2.427.0-64.842.0T,数据计算结果,A=hilb(4);b=1;2;1.41;2;b1=b;b1(3)=b1(3)+.01;x=Ab;x1=Ab1;error=x-x1;formatshortex,x1,error,ans=-1.6560e+002-1.6320e+002-2.4000e+0001.8330e+0031.8060e+0032.7000e+001-4.4232e+003-4.3584e+003-6.4800e+0012.8980e+0032.8560e+0034.2000e+001,数值试验程序,定义1.,向量和矩阵的范数,显然,并且由于,例.求下列向量的各种常用范数,解:,向量和矩阵范数/*NormsofVectorsandMatrices*/为了误差的度量,常用向量范数：,注：,可以理解为,可以理解为对任何向量范数都成立。,矩阵范数/*matrixnorms*/,(4)*|AB|A|B|（相容/*consistent*/当m=n时),例.,不难验证其满足定义2的4个条件,称为Frobenius范数,简称F-范数,而且可以验证,tr为矩阵的迹,类似向量的2-范数,定义.,简称为从属范数或算子范数,显然，由定义不难推出,常用矩阵范数：,Frobenius范数,向量|2的直接推广,对方阵以及有,利用Cauchy不等式可证。,算子范数,/*operatornorm*/,由向量范数|p导出关于矩阵ARnn的p范数:,则,矩阵ATA的最大特征根/*eigenvalue*/,例.,求矩阵A的各种常用范数,解:,由于,特征方程为,容易计算,计算较复杂,对矩阵元素的变化比较敏感,不是从属范数,较少使用,使用最广泛,性质较好,我们只关心有相容性的范数，算子范数总是相容的。,即使A中元素全为实数，其特征根和相应特征向量/*eigenvector*/仍可能是复数。将上述定义中绝对值换成复数模均成立。,若不然，则必存在某个向量范数|v使得对任意A成立。,Counterexample?,谱半径/*spectralradius*/,(A),证明：,由算子范数的相容性，得到,将任意一个特征根所对应的特征向量代入,证明：,A对称,若是A的一个特征根，则2必是A2的特征根。,又：对称矩阵的特征根为实数，即2(A)为非负实数，故得证。,对某个A的特征根成立,所以2-范数亦称为谱范数。,矩阵的条件数概念,方程组Ax=b,右端项b有一扰动引起方程组解x的扰动,设x是方程组Ax=b的解,则有,所以,定义条件数:Cond(A)=|A1|A|或C(A)=|A1|A|,当条件数很大时,方程组Ax=b是病态问题;当条件数较小时,方程组Ax=b是良态问题,类似,设方程组Ax=b,矩阵A有一扰动时,将引起方程组解x的扰动,设x是方程组Ax=b的解,则有,化简,得,取范数,Afamousexampleofabadlyconditionedmatrix,迭代法的基本思想,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,迭代格式的建立,Ax=b,k=0,1,2,给定一个初始向量x(0)，可得迭代格式：,若产生的迭代序列x(k)收敛到一个确定的向量x*，则x*就是原方程组的解。,其中B称为迭代矩阵。,第52页，共53页,Jacobi迭代,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,k=0,1,2,则可得雅可比(Jacobi)迭代格式:,令A=D-L-U，其中,称为雅可比(Jacobi)迭代矩阵,第53页，共53页,Jacobi迭代,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,在计算时，如果用代替，则可能会得到更好的收敛效果。此时的迭代公式为,Jacobi迭代的分量形式：,第54页，共53页,GS迭代,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,写成矩阵形式:,称为GS迭代矩阵,此迭代格式称为高斯-塞德尔(Gauss-Seidel)迭代,第55页，共53页,Jacobi、GS算法,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,Jacobi算法,GS算法,第56页，共53页,举例,解：,Jacobi迭代格式,令则迭代得：,举例（续）,2008年10月6日12时7分,沈阳航空工业学院飞机设计教研室,GS迭代格式,得,第58页，共53页,矩阵分裂法,Jacobi迭代,GS迭代,A=M-N,M=D,N=MA=L+U,M=DL,N=U,Jacobi迭代的矩阵表示为Gauss-Seidel迭代的矩阵表示为均属于如下形式的迭代公式.,迭代法的收敛性,定理（迭代法的基本收敛定理）迭代过程X(k+1)=BX(k)+g对于任意初始向量X(0)及右端向量g均收敛的充要条件是迭代矩阵B的谱半径(B)1,并且(B)愈小，收敛速度愈快.,定理（迭代法收敛的充分条件）若迭代法X(k+1)=BX(k)+g的迭代矩阵B满足，B=q1,则对于任意的初始向量X(0)与右端向量g迭代法收敛.,证明,证设X*为方程组X=BX+g的精确解,则X*=BX*+gX*-X(k+1)=B(X*-X(k)X*-X(k+1)BX*-X(k)反复利用此不等式和已知条件有X*-X(k+1)BkX*-X(0)0从而对于任意的X(0)与g迭代收敛.,对角占优阵,称方阵A=(aij)nn为对角占优的，如果或,定理（迭代法收敛的充分条件）若线性代数方程组AX=b的系数方阵A=(aij)nn是按行(或按列)严格对角占优的，则Jacobi迭代法和Gauss-Seidel迭代法都是收敛的.,例,分别用Jacobi迭代法与Gauss-Seidel迭代法求解方程组精确到小数点后四位，并要求分别写出其迭代法的分量形式和矩阵形式.,解,（1）用Jacobi迭代法，其