《高斯消去法.ppt》由会员分享，可在线阅读，更多相关《高斯消去法.ppt（54页珍藏版）》请在三一文库上搜索。
1、2021/4/30,数值分析,第五章 解线性方程组的直接方法,5.1 高斯消去法,5.2 高斯主元素消去法,5.4 误差分析,5.3 矩阵的三角分解,2021/4/30,数值分析,【本章重点】1Gauss 消去法和列主元消去法及其实现条件。2矩阵的三角分解，含LU分解和LLT 分解及三对角方程组的追 赶法。3向量和矩阵范数的定义及性质。4矩阵条件数及病态矩阵定义和解方程组直接法的误差估计。,2021/4/30,数值分析,在自然科学和工程技术中许多问题的解决转化为解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵，一种是高阶稀疏矩阵。,引言,解线性方程组的数值解也有两种： 直接
2、法，就是经过有限步算术运算，可以求得线性方程组的解，但实际计算时有舍入误差的存在和影响，所以求得的结果也只能是近似解对低阶稠密矩阵和部分大型稀疏矩阵有效。 迭代法，就是用某种极限过程去逐步逼近精确解，是解决大型稀疏矩阵的重要方法。从一个初始向量出发，按照一定的迭代格式，构造出一个趋向于真解的无穷序列。,2021/4/30,数值分析,高斯消去法是一个古老的求解线性方程组的方法，但由于它的改进、变形得到的选主元素消去法、三角分解法仍是在计算机上求解系数矩阵为中、低阶稠密矩阵的线性方程组常用的有效方法，所以本节介绍这一方法。,1 高斯消去法,用高斯消去法求解n阶线性方程组Ax=b的基本思想是在逐步消
3、元的过程中，把方程组的系数矩阵化为上三角矩阵，从而将原方程组约化为容易求解的等价三角方程组，然后进行回代求解。,一、 高斯消去法,2021/4/30,数值分析,例1 用消去法求解方程组,还原为方程组,（消元过程）,（回代过程）,x3=3， x2=2， x1=1,2021/4/30,数值分析,设有方程组,改写成矩阵形式：,简记为：,一般线性方程组的高斯消去法,2021/4/30,数值分析,将原方程组记为,其中,则第一步（k=1)，若a11不等于0，则可以计算乘数,用-mi1 乘方程组的第一个方程加到第i个方程，则原方程组同 解方程组为：,其中,简记为,2021/4/30,数值分析,第二步，设经过
4、k-1次消元后的同解方程组为,设,计算乘数,用-mik乘上面的线性方程组的第k个方程加到第i个方程，可 以消去xk元,得到同解方程组,其中,2021/4/30,数值分析,重复上述过程，可以将方程组化为等价的简单方程组,其中,为上梯形阵。,当m=n时，与原方程组等价的方程组为,即,（消元过程）,由上述方程组很快可以求出：,（回代过程）,2021/4/30,乘除运算量,乘除运算量：,由于计算机中做乘除运算的时间远远超过做加减运算时间，故估计运算量时，往往只估计乘除的次数。,高斯消去法总的乘除运算量为：,2021/4/30,数值分析,定理1 设Ax=b,其中AR nn,(1) 如果,则可通过高斯消去
5、法，将方程组化为等价三角形方程组， 计算公式为： (a)消元计算(k=1,2,n),(b)回代计算,(2)如果系数矩阵A为非奇异矩阵则可通过高斯消去法与初等行 变换方法，将方程组约化为三角形方程组。,2021/4/30,数值分析,证明略,2021/4/30,数值分析,function x=gauss(A,b) n=length(b); x=zeros(n,1); for i=1:n-1 %消元过程 for k=i+1:n for j=i+1:n A(k,j)=A(k,j)+A(i,j)*(-A(k,i)/A(i,i); end b(k)=b(k)+b(i)*(-A(k,i)/A(i,i); A
6、(k,i)=0; end end disp(A) disp(b) pause,x(n)=b(n)/A(n,n); for i=n-1:-1:1 %回代过程 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j); end x(i)=(b(i)-sum)/A(i,i); end,程序设计,2021/4/30,数值分析,A=1 1 1;0 4 -1;2 -2 1; b=6 5 1; x=gauss(A,b),程序执行,消元后的A 1 1 1 0 4 -1 0 0 -2 消元后的b 6 5 -6 x = 1 2 3,2021/4/30,数值分析,二、矩阵的三角分解,由上面消去法
7、的分析知，若方程组的系数矩阵A的顺序主子式 不等于0，高斯消元的过程是作一系列初等行变换，等价于矩 阵A左乘一个初等矩阵。,一般第k步消元，A(k)化为A(k+1),b(k)化为b(k+1),相当于,所以第一次消元可表示为：,其中：,2021/4/30,数值分析,重复以上过程，最后得到,令,为单位下三角矩阵。,2021/4/30,数值分析,定理2（矩阵的LU分解）设A为一个n 阶矩阵，若其顺序主子式 Di0(i=1,2,n)。则A可分解为一个单位下三角形矩阵L和一个 上三角形矩阵U的乘积，且这种分解是唯一的。,证明：存在性由定理1 可得； 唯一性：设有两种分解，即A=LU=L1U1,所以由L可
8、逆，U1可逆，知L-1 L1=UU1-1,左边为单位下三角形矩阵，右边为上三角形矩阵，所以 均为单位矩阵，即有 L=L1 ,U=U1 分解唯一。,Ax=b LUx=b 记Ux=y,Ly=b 由Ly=b求出y,再由Ux=y求出x,注意：LU分解的作用是将一个较复杂的矩阵化为两个简单的矩阵。,2021/4/30,数值分析,A=1 1 1;0 4 -1;2 -2 1; L,U=lu(A) L = 0.5000 0.5000 1.0000 0 1.0000 0 1.0000 0 0 U = 2 -2 1 0 4 -1 0 0 1,Matlab函数： L,U=lu(A) U为上三角矩阵 （与理论略有不同
9、）,A=1 1 1;0 4 -1;2 -2 1; L = 1 0 0 0 1 0 2 -1 1 U = 1 1 1 0 4 -1 0 0 -2,2021/4/30,数值分析,高斯消去法在消元的过程中,可能出现 的情况,这时消去法将无法进行;即使主元素 很小,用其作除数,会导致其他元素数量级的严重增长和舍入误差的扩散。 因此，为了减少误差，每次消元选取系数矩阵（或某列）中绝对值最大的元素作为主元素，这就是全（列）主元素消去法。 全主元素消去法在选取主元时要花费较多的机器时间，目前主要使用列主元素消去法。,2 高斯主元素消去法,列主元素消去法的基本思想：依次按列选主元素，然后换行使之变到主元素位置
10、上，再进行消元计算。,2021/4/30,列主元高斯消去法,高斯消去法有效的条件是：,如果某个主元为 0，则可能导致高斯消去法求解失败。,列主元高斯消去法：在第 k 步消元时, 先选取列主元：, if ik k then 交换第 k 行和第 ik 行；,列主元高斯消去法比普通高斯消去法要多一些比较运算，但比普通高斯消去法稳定。, 消元,2021/4/30,全主元高斯消去法,列交换改变了 xi 的顺序，须记录交换次序，解完后再换回来。,全主元高斯消去法具有很好的稳定性，但选全主元比较费时，故在实际计算中很少使用。,2021/4/30,举例（二）,解：,例：采用十进制八位浮点数，分别用Gauss消
11、去法和列主元Gauss消去法求解线性方程组:,Gauss消去法：,9个,2021/4/30,举例（二）续,列主元Gauss消去法：,例：,但列主元高斯消去法有时也会导致求解失败。,2021/4/30,数值分析,function x=gauss_lie(A,b) n=length(b); x=zeros(n,1); c=zeros(1,n); for i=1:n-1 max=abs(A(i,i); m=i; for j=i+1:n if maxabs(A(j,i) max=abs(A(j,i); m=j; end end if(m=i) for k=i:n c(k)=A(i,k); A(i,k)
12、=A(m,k); A(m,k)=c(k); end d1=b(i); b(i)=b(m); b(m)=d1; end,for k=i+1:n for j=i+1:n A(k,j)=A(k,j)+A(i,j)*(-A(k,i)/A(i,i); end b(k)=b(k)+b(i)*(-A(k,i)/A(i,i); A(k,i)=0; end end disp(消元后的A) disp(A) disp(消元后的b) disp(b) pause x(n)=b(n)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j); end x(
13、i)=(b(i)-sum)/A(i,i); end,程序设计,2021/4/30,数值分析,A=1 1 1;0 4 -1;2 -2 1; b=6 5 1; x=gauss_lie(A,b) 消元后的A 2 -2 1 0 4 -1 0 0 1 消元后的b 1 5 3 x = 1 2 3,P14 a=0.00001 1;2 1 b=1 2 gauss(a,b) gauss_lie(a,b),2021/4/30,数值分析,26,很多问题都需要计算方阵的逆阵.线性代数中有多种求逆方法,本节讨论求A-1的数值方法.A非奇异,求A-1,高斯-若当(Gauss-Jordan)消元法,高斯-若当消元法是高斯消
14、元法的一种变形.高斯消元法是消去对角元下方的元素.若同时消去对角元上方和下方的元素,而且将对角元化为1,就是高斯-若当消元法.,设高斯-若当消元法已完成k-1步,于是Ax=b化为等价方程组A(k)x=b(k),增广阵,2021/4/30,数值分析,在第k步计算时,考虑将第k行上下的第k列元素都化为零,且akk化为1.对k=1,2,n 1按列选主元,确定ik使,2021/4/30,数值分析,2换行,交换增广阵第k行和第ik行,(j=k,k+1,n),3计算乘数 mik= - aikakk(i=1,2,n, ik) mkk=1akk. (5.14),4消元 aij=aij+mikakj(i=1,2
15、,n,且ik,j=k+1,n) bi=bi+mikbk (i=1,2,n,且ik) (5.15),2021/4/30,数值分析,5主行计算 akj=akjmkk (j=k,k+1,n) bk=bkmkk,当k=n时,显然xi=bi,i=1,2,n 就是 Ax=b的解.,2021/4/30,数值分析,高斯-若当消元法的消元过程比高斯消元法略复杂,但省去了回代过程,它的计算量约为n3/2,大于高斯消元法.也称为无回代的高斯消元法.,高斯-若当消元法解方程组并不比高斯消元法优越,但用于矩阵求逆是适宜的,实际上它是初等变换方法求逆的一种规范化算法。,定理（高斯-约当法求逆矩阵）,2021/4/30,数
16、值分析,例3,2021/4/30,数值分析,所以,2021/4/30,数值分析,3 矩阵的三角分解,一、直接三角分解,Ax=b,LUx=b,二、列主元的三角分解,令Ux=y,先求y=b,再求Ux=y,PAx=Pb,LUx=Pb,令Ux=y,先求y=Pb,再求Ux=y,2021/4/30,数值分析,三、Cholesky(平方根)法,就是利用对称正定矩阵的三角分解而的到的求解对称正定方程组的一种有效方法。,对称正定矩阵A:AT=A,XTAX0,A=LU,A为对称矩阵，且A的所有各阶顺序主子式均不为零。 A=LU=LDU0, AT=(LDU0)T=U0TDLT=A 由分解的唯一性得：U0T=L 故：
17、 A=LDLT,A为对称正定矩阵，对角矩阵D的元素均大于零。 A=LDLT =LD1/2D1/2LT =(LD1/2)(LD1/2)T =L1L1T 其中：L1为下三角阵,2021/4/30,数值分析,算法1 对正定矩阵A做LLT分解(平方根法),Ax=b LLTx=b 记LTx=y,Ly=b 由Ly=b求出y, 再由LTx=y求出x,公式：故对于j=1,2,n 1）求L,2）由Ly=b求出y,3）由LTx=y求出x,2021/4/30,数值分析,算法2 对A做LDLT分解(改进的平方根法),Ax=b LDLTx=b 记DLTx=y,Ly=b 由Ly=b求出y, 再由DLTx=y求出x,公式：
18、故对于i=1,2,n 1）求L和D,2）由Ly=b求出y,3）由DLTx=y求出x,2021/4/30,数值分析,解:设,由矩阵乘法得,，,，,2021/4/30,数值分析,四、追赶法,三对角占优矩阵,2021/4/30,数值分析,对比得：,2021/4/30,数值分析,计算流程：,Ax=f LUx=f 记Ux=y,Ly=f 先由Ly=f求出y,再由Ux=y求出x,追赶法公式： 1.计算的递推公式：,2.解Ly=f：,3.解Ux=y：,2021/4/30,数值分析,由,例P230.8,由公式,解:设,2021/4/30,数值分析,由,得,2021/4/30,数值分析,function x=th
19、reedia(a,b,c,f) n=length(f); L=eye(n); x=zeros(1,n); y=zeros(1,n); d=zeros(1,n); u=zeros(1,n); d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i+1)*u(i); end y(1)=f(1)/d(1); for i=2:n y(i)=(f(i)-a(i)*y(i-1)/d(i); end x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end,解方程组 a=0 -1 -1 -3; b=2
20、3 2 5; c=-1 -2 -1 0; f=6 1 0 1 threedia(a,b,c,f),X= 5 4 3 2,2021/4/30,数值分析,4 误差分析,一、矩阵范数,设 ，给出向量范数 ，则定义： 为矩阵A的(算子）范数。,定义2 矩阵的算子范数,2021/4/30,数值分析,四种常用的矩阵范数：,（1）列范数,（2）2-范数,（3）行范数,（4）F-范数,2021/4/30,数值分析,例： ，计算A的各种矩阵范数。,解：,2021/4/30,数值分析,二、谱半径：,定义：设 的特征值为（i=1,2,n),称 为A的谱半径。,定理：矩阵范数与谱半径之间的关系为: (A) |A|.,
21、定理：如果 为对称矩阵,则,例： ，计算A的谱半径。,2021/4/30,数值分析,三、矩阵的条件数,定义 设是非奇异阵，称数 为矩阵A的条件数。,常用的条件数有:,当A为对称矩阵时, 为A的绝对值最大最小的特征值.,2021/4/30,数值分析,条件数的性质： 1. A为非奇异阵，有cond(A)v1 2. A为非奇异阵且c0(常数),则cond(cA)v=cond(A)v 3. A为正交阵，有cond(A)2=1,2021/4/30,数值分析,四、误差分析,考虑线性方程组Ax=b中A或b的微小扰动（误差）对解的影响。,精确解为：x=(2,0)T,解为：x=(1,1)T,现考虑Ax=b中b2
22、的微小扰动对解的影响。,定义 如果矩阵A或常数项b的微小的变化，引起方程组Ax=b解的巨大变化，则称此方程组为“病态”方程组,矩阵A称为“病态”矩阵，否则称此方程组为“良态”方程组,矩阵A称为“良态”矩阵。cond(A)1时，方程组是“病态”的(A是“病态”矩阵),2021/4/30,数值分析,2021/4/30,数值分析,例 已知n阶希尔伯特（Hilbert）矩阵,(1)计算H3条件数cond (H3),(2)考虑方程组H3x=b=(11/6,13/12,47/60)T,设H3及b有微小误差(取3为有效数字),分析对解的影响.,2021/4/30,数值分析,解 (1),2021/4/30,数值分析,(2) 方程组H3x=b=(11/6,13/12,47/60)T的解为x=(1,1,1)T,若H3及b有微小误差(取3为有效数字),则:,