数值分析第5章解线性方程组的直接方法

1、5 .1 引言与预备知识 5 .2 高斯消去法 5 .3 高斯主元素消去法 5 .4 矩阵三角分解法 5 .5 向量和矩阵的范数 5 .6 误差分析,Ch5 解线性方程组的直接方法,在自然科学和工程技术中有很多问题的解决常常归结为解线性代数方程组如三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程的边值问题等都导致求解线性代数方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵，另一种是大型稀疏矩阵。,关于线性方程组的数值解法一般有两类： 1.直接法 2.迭代法,5.1引言与预备知识,5.1.1 引言,本章讨论n

2、元线性方程组,（5.1）,的直接解法。方程组(5.1)的矩阵形式为,Ax=b,其中,若矩阵A非奇异,即det(A)0,则方程组(2.1)有唯一解。,所谓直接解法是指，若不考虑计算过程中的舍入误差，经过有限次算术运算就能求出线性方程组的精确解的方法。,但由于实际计算中舍入误差的存在，用直接解法一般也只能求出方程组的近似解。,Cramer法则是一种不实用的直接法，本章将介绍几种实用的直接法。,5.1.2 预备知识,M行n列矩阵.,n维列向量.,矩阵的基本运算: (1)矩阵的加法,(7)矩阵的行列式 行列式性质: (a) det(AB)=det(A)det(B),(6)非奇异矩阵,(5)单位矩阵,(

3、4)转置矩阵,(3)矩阵与矩阵的乘法,(2)矩阵与标量的乘法,5.1.3矩阵特征值与谱半径,定义1设,若存在一个数(实数或复数)和非零向量,使,(1.1),则称为A的特征值,x为A对应的特征向量,A的全体特征值称为,A的谱，记作,称为A的谱半径.,(1.2),由式(1.1)知, 可使齐次方程组,有非零解,故系数行列式,记,称为矩阵的特征多项式，方程(1.3)称为的特征方程,(1.3),在复数域中有n个根,故,由行列式(1.3)展开可知：,的n个特征值,故,是它的特征,方程(1.3)的几个根,并有,(1.4),(1.5),A的迹.,A的特征值和特征向量x还有以下性质:,(1) AT与A有相同的特

4、征值 及相同的特征向量x .,(2) 若A非奇异，则A-1的特征值为-1,特征向量为x.,(3) 相似矩阵 B=S-1AS有相同的特征多项式.,例1 求,的特征值及谱半径.,解: A的特征方程为,故A的特征值为,A的谱半径为,5.1.4 特殊矩阵,定理1.,定理2.,定理3.,定理4 (Jordan标准型) 设A为n阶矩阵,则存在一个非奇异矩阵P使得,其中:,(1)当A的若当标准型中所有若当块Ji均为一阶时,此标准型变成 对角矩阵.,返回主页,求解, 高斯消去法(逐次消去法)及消去法和矩阵三角分解之间的关系：,5.2 高斯消去法,例2 用消去法解方程组,解 第1步.将方程(2.2)乘上-2加到

5、方程(2.4)上,消去未知数x1,得到,(2.2),(2.3),(2.4),第2步.将方程(2.3)加到方程(2.5)上去,消去方程(2.5)中的x2, 得到与原方程组等价的三角形方程组,解为:,首先举一个简单的例子来说明消去法的基本思想.,上述过程相当于,消元,记,其中,Step k：设 ，计算因子 且计算,回代,注意1：只要 A 非奇异，即 A1 存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。,共进行 ？ 步,n 1,注意2:设Ax=b,其中A为非奇异矩阵,如果,由于A为非奇异矩阵,所以A的第一列一定有元素不等于零.例如,定理5 设Ax=b,其中,(1)如果,则可通

6、过高斯消去法将,Ax=b约化为等价的三角形线性方程组(2.10),且计算公式为:,消元计算(k=1,2,n-1),回代计算,(2)如果A为非奇异矩阵,则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为方程组(2.10).,定理6 约化的主元素aii(i) 0(i=1,2,k)的充要条件是矩阵A的所有顺序主子式 /* determinant of leading principal submatrices */ Di0(i=1,2,k) .即,(2.12),推论 如果A的顺序主子式Dk0(k=1,2, ,n-1),则,5.2.2 三角分解法 /* Matrix Factorizat

7、ion */, 高斯消元法的矩阵形式 /* Matrix Form of G.E. */：,Step 1:,单位下三角阵 /* unitary lower-triangular matrix */,证明：由1中定理可知，LU 分解存在。下面证明唯一性。,若不唯一，则可设 A = L1U1 = L2U2 ，推出,Upper-triangular,Lower-triangular With diagonal entries 1,注： L 为一般下三角阵而 U 为单位上三角阵的分解称为Crout 分解。 实际上只要考虑 A* 的 LU 分解，即 ，则 即是 A 的 Crout 分解。,例3 对于例2

8、,系数矩阵,由高斯消去法,返回主页,5.3.1 列主元消去法:,在顺序消元过程中,只要,即可进行计算,但如果,很小,则将导致舍入误差增长,使解的误差很大.,例4 用Gauss 消去法求解方程组,5 .3 高斯主元素消去法,解:因,故方程有唯一解,且精确解为,若用Gauss消去法取四位有效数字计算,可得解,比较,误差很大,若将两个方程互换为,用Gauss消去法取四位有效数字计算,可得解,本例表明通过行交换可避免舍入误差增长,这就是列主元消去法的基本思想. 其计算步骤如下:,第1步，在,中的第1列选主元，即,行为主元,若,将,的第i1行与第1行互换，再按消元,公式计算得到,假定上述过程已进行(k-

9、1)步，得到,第k步，在,中的第k列选主元，即,若,则在,中将ik行与第k行互换，再按消元,公式计算得到,对k=1,2, ,n-1,重复以上过程则得,如果某个k出现主元,方程没有唯一解或严重病态，否则可由(3.2.4)求得解.,则表明detA=0，,它也表明当非奇异时，存在排列矩阵(若干初等排列矩阵的乘积),使PA=LU,其中L为单位下三角矩阵,其元素|lij|=1,U为上三角矩阵.,上述每步行交换后再消元相当于,其中,是指标为k的初等下三角矩阵,为初等排列矩阵,时,表示不换行,经过(n-1)步换行与消元,A,化为上三角矩阵.,即:,解:,例5用列主元消去法解x=b,其中,消元,消元,消元结束

10、由回代公式求得解,此例的精确解为,可见结果精度较高若不选列主元auss消去法，求得解,，误差较大,除列主元消去法外，还有一种消去法，是在的所有元素aij中选主元，称为全主元消去法因计算量较大且应用列主元已能满足实际要求，故不再讨论目前很多数学软件库都有列主元消去法，可直接调用,注:,为了减少计算的舍入误差,使用消去法通常都要选主元.目前最常用的是列主元消去法,也就是每步消元之前选主元,当A=(aij) 第一步选A中第1列的主元,即max|ai1|=ai1.然后将i1行与第行互换，再进行消元，以后每步消元做法类似，先选主元，再消元,5.3.2 高斯若当消去法,消去对角线上方和下方的元素.,假设已

11、经完成k-1步,得到与方程Ax=b等价的方程组,返回主页,高斯消去法有很多变形,有的是高斯消去法的改进、改写，有的是用于某一类特殊性质矩阵的高斯消去法的简化。,5.3.1 直接三角分解法.,将高斯消去法改写为紧凑形式,可以直接从A的元素得到计算L,U元素的递推公式,而不需要任何中间步骤,这就是所谓直接三角分解法.,一旦实现A的LU分解,那么求解Ax=b的问题就等价于求解两个三角形方程组 (1) Ly=b,求y; (2) Ux=y,求x.,.4 矩阵三角分解法 /* Matrix Factorization */,.不选主元的三角分解法 设A为非奇异矩阵,且有分解式A=LU,其中L为单位下三角,

12、U为上三角即,L,U元素可以由n步直接计算定出,其中第r步定出U的第r行和L的第r列元素.由上式有:,故,同样有:,设已经定出U的第1行到第r-1行元素与L的第1列到第r-1列元素.利用矩阵乘法(注意当rk时,lrk=0),有,得,固定 r : 对 i = r, r+1, , n 有,lii = 1,a,将 r ，i 对换，对 r = i, i+1, , n 有,b,结论:用直接三角分解法解Ax=b(要求A的所有顺序主子式都不等于0)的计算公式如下.,Step 1: u1i = a1i; li1 = ai1 / u11; ( i = 2, , n ),计算U的第r行,L的第r列元素(r=2,3

13、,，n）,Step 2:,求解Ly=b, Ux=y 的计算公式：,Step 3:,Step 4:,Step 5:,例5 用直接三角分解法解,解 用分解公式计算得,求解,2.选主元的三角分解法 当urr=0时计算中断，或者当urr绝对值很小时，按分解公式计算可能引起舍入误差的积累。 但如果A非奇异，可以通过交换A的行实现矩阵A的LU分解，因此可采用与列主元消去法类似的方法，将直接三角分解法修改为（部分）选主元的三角分解法。,设第r-1步分解已完成，这时有,第r步分解需用到（3.2）及（3.3）式，为避免用小的数urr做除数，引进量,于是有：,取,交换A的r行与ir行元素，将,调到(r,r)位置,

14、(将（i，j）位置的新元素仍记为lij 及 aii),于是有| lir |=1,(i=r+1,，n）.由此再进行第r步分解计算。,5.3.2 平方根法,当A对称正定时，A的顺序主子式,故由定理知，A=LU的分解存在且唯一，其中L为单位下三角,为了A利用对称性,其中D为对角阵,U0为单位上三角阵,于是,又,代入到上式,就得到对称矩阵A的分解式,矩阵，U为上三角矩阵，且,定理9(对称阵的三角分解定理),设A为n阶对称阵,且A的顺序主子式,则A可唯一分解为,，其中L为单位下三角矩阵，,.,D为对角矩阵,定理10(对称正定矩阵的三角分解或Cholesky分解) 如果A为n阶对称正定矩阵,则存在一个实的

15、非奇异下三角阵L使A=LLT,当限定L的对角元素为正时,这种分解是唯一的.,将求方程组的解转化为求方程的解LLTx=b.令,，求得方程的解,由,根据矩阵乘法，由,得,i=j有,当ij,得,例6用平方根法求以下方程组的解.,解先验证系数矩阵A对称正定，对称显然，,故A对称正定，可用Cholesky分解计算，求得,求解,得,再求,得, 5.3.3 追赶法解三对角方程组 /* Crout Reduction for Tridiagonal Linear System */,在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对

16、角线方程组.,简记为Ax=f.其中,当|i-j|1时,aij=0,且：,Step 1: 对 A 作Crout 分解,直接比较等式两边的元素，可得到计算公式。,为下三角，为单位上三角,注意当j=1时有,对j=2，3，n求得L的元素,，这就是A的Cholesky分解，然后再解两个三角方程组，得,这就是对称正定方程组的平方根法,另外，由于,故有,这表明分解过程中矩阵L中元素,因此平方根法计算是数值稳定的.,的数量级不增长，,Step 2: 追即解,Step 3: 赶即解 ：,求解xf等价于,Step :计算,的递推公式,将计算系数,为追的过程,将计算方程组的解,为赶的过程,定理：设有三对角线方程组x

17、f，其中满足对角占优的条件，,则为非奇异矩阵且追赶法计算公式中的,满足：,注： 如果 A 是严格对角占优阵，则不要求三对角线上的所有元素非零。 根据不等式 可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。 运算量为 O(6n)。,返回主页,5.5.1 内积与向量范数,为了研究方程组Ax=b解的误差和迭代法收敛性，需对向量,及矩阵,的大小引进一种度量，就要定义范数，,它是向量长度概念的直接推广，通常用,表示n维实向量,空间，,表示n维复向量空间.,定义2 设,将实数,称为向量x,y的数量积.,非负实数,称为向量x的欧氏范数或2-范数.,5.5 向量和矩阵范数,定理12设,则内积有以下性质：

18、,(1),(2),(3),(4),(5) (柯西-施瓦茨不等式),等式当且仅当x与y线性相关时成立;,(6) 三角不等式,定义3(向量范数)如果向量,的某个实值函数,满足条件:,则称,对于,由内积性质可知它满足定义2的三个条件，,故它是一种向量范数.此外还有以下几种常用的向量范数.,容易验证,均满足定义2的三个条件.更一般的还可定义,但只有p=1,2,时的三种范数是常用的向量范数.,例如给定,则可求出,定理14设,是,则N(x)是向量x的分量,上任一种向量范数，,的连续函数.,定理15(向量范数的等价性)设,是,上任意两种向量范数，则存在常数,使,5.5.2 矩阵范数,矩阵,可看成nn维向量，

19、如果直接将向量,的2-范数用于矩阵A，则可定义,称为矩阵A的Frobenius范数，简称F-范数.它显然满足向量范数的三条性质，但由于矩阵还有乘法运算，因此矩阵范数的定义中应增加新条件.,定义4 如果,的某个非负实函数N(A)，记作A，,满足条件：,则称,显然,满足定义中的四个条件, (3)，(4)两条均可由,Cauchy-Schwarz不等式证明，故,是一种矩阵范数.,除矩阵自身的运算外，在解方程中矩阵乘向量的运算即Ax，也是必不可少的.因此要求所引进的范数应满足条件：,上式称为相容性条件.为使引进的矩阵范数满足条件(4.5)，我们给出以下定义.,(4.5),定义6(矩阵的算子范数) 设,当

20、给定向量范数,时可定义,称为矩阵的算子范数或从属范数.,(4.6),定理17 设,上的一种向量范数，,则由(4.6)定义的,是一种矩阵范数，且满足相容性条件,证明 因,中有界闭集,上的连续函数，故,在D上有最大值，即,使,而对,故,所以,从而当,成立，而x=0时,显然也成立.,定理17 设,则,这里,为矩阵的谱半径.,例7 已知,解,从定理可以看出，计算,较容易，而计算,时因为要求,的特征值，所以较为困难.但当A对称时，有,定理19,定理18 对任何,为任一种从属范数则,反之，对任意0，至少存在一种从属范数,使,证明:设,为A的特征值，则,由,得,非奇异，且,证明用反证法.假定(I+B)奇异，

21、则齐次方程,有非零解,而,与B1的假设矛盾，故(I+B)非奇异.,又,得,取范数,得,定理20设,返回主页,5.6.1矩阵条件数与扰动方程组误差界在解方程组Ax=b时，由于各种原因，A或b往往有误差，从而使得解也产生误差.,例8 方程组,的准确解为,当A与b有微小变化时，如变为方程,则准确解为,它表明A,b的微小扰动引起方程解x的很大,变化，这就是病态方程.,5.6 误差分析与病态方程组,定义7 求解线性方程组Ax=b时，若A或b有微小扰动,解x的误差,很大, ，则称此方程组为,病态方程组，相应的系数矩阵A称为病态矩阵.,反之，若此时,很小, ，则称此方程组为,良态方程组，相应的系数矩阵A称为

22、良态矩阵.,注意方程组是否病态与用什么数值方法无关，它是由方程自身性质决定的.,在例8中因为行列式,因此出现病态.但有时A从表面上看性质很好，也可能是病态的.,那么如何判断A是否病态?先给出如下定义.,例9 方程组Ax=b表示为,它的准确解,A对称正定且,表面看性质较好，但若对右端b作微小变化，如方程改为,则解变为,这里b的相对误差大约只有,但解的相对误差却很大，,故A也是病态矩阵.,定义8 设,非奇异，v为矩阵的任一种从属范数，则,称为矩阵A的条件数.,从定义看到矩阵条件数依赖于范数的选取，如范数为2-范数，,则记为,同理有,等等.,条件数有以下性质：,(1),(2),(3) U为正交矩阵，

23、则,(4) 若,与,为A的按模最大与最小特征值，则,若A对称，则,下面给出扰动方程组解的误差分析.先考察b有扰动,则扰动方程为,由于Ax=b,故得,于是,再由Ax=b,有,即,故得,下面再研究方程Ax=b，当A有扰动,时，其解,的误差分析.,此时扰动方程为,因Ax=b,故有,因,存在,若假定,则由定理20可知,非奇异,并有:,(5.6),由(5.6)可得,因此,(5.7),定理22 设A为非奇异矩阵,Ax=b0,且,如果,则(5.7)式成立.,从(5.7)看到，当A的条件数Cond(A)很大时，解的相对误差,也很大，故方程组为病态.在例9中,而,于是,条件数很大，故方程是严重病态的.,例10

24、Hibert矩阵是一个著名的病态矩阵，记作,它是一个对称正定矩阵，当n3时它是病态矩阵.例如,故,另外还有,等等.,因此Hn是严重的病态矩阵，且n越大 Cond(Hn)越大.,例11 在例10的方程组中可算出A的特征值,故,例中,实际相对误差是,而根据(5.6)的误差估计为,这与实际相差不大，即相对误差放大了将近3 000倍.故方程为病态方程组.,定理23(事后误差估计)设方程组,，则,若实际求得解为,证明记剩余,则,它表明如果方程组病态，即使剩余r很小，解的相对误差仍可能很大.,5.6.2 病态方程组的解法,如果A的条件数Cond(A)1，则Ax=b为病态方程，但计算Cond(A)时需要求

25、A-1,计算量很大，相当于解方程组， 在实际中常可通过求解过程直观地判断方程组的病态性质，如果解方程时出现下述情况之一，则可能是病态方程组.,(1) 在列主元消去法中出现小主元；,(2) 在计算过程中行或列几乎线性相关或三角分解中对角元出现近似零的元素；,(3) 矩阵A的元素数量级相差很大且无规律；,(4) 剩余,很小，而解,很大，又达不到精度要求.,(1) 采用高精度运算，减轻病态影响，例如用双倍字长运算.,对病态方程组求解可采用以下措施：,(2) 用预处理方法改善A的条件数，即选择非奇矩阵,与Ax=b等价,而,(3) 平衡方法，当A中元素的数量级相差很大，可采用行均衡或,列均衡的方法改善A

26、的条件数.设,非奇异，计算,于是求Ax=b等价于求,的条件数可得到改善，这就是行均衡法.,例12 给定方程组Ax=b为,A的条件数,若用行均衡法可取,则平衡后的方程,用三位有效数字的列主元消去法求解得,function x=threedia(a,b,c,f) N=length(f); x=zeros(1,N); y=zeros(1,N); beta=zeros(1,N); gramma=zeros(1,N); beta(1)=b(1); for i=1:N-1 gramma(i)=c(i)/beta(i); beta(i+1)=b(i+1)-a(i+1)*gramma(i); end %追的过

27、程 y(1)=f(1)/beta(1); for i=2:N y(i)=(f(i)-a(i)*y(i-1)/beta(i); end %赶的过程 x(N)=y(N); for i=N-1:-1:1 x(i)=y(i)-gramma(i)*x(i+1); end,a=0,-1,-1,-3; b=2,3,2,5; c=-1,-2,-1,0; f=6,1,0,1; x=threedia(a,b,c,f),追赶法求解三对角方程组,Cholesky方法：,A=4,-2,4;-2,17,10;4,10,9;,b=8.7,13.7,-0.7;,x,L,D=Chol_decompose(A,b),L = 1.

28、0000 0 0 -0.5000 1.0000 0 1.0000 0.7500 1.0000,D = 4 16 -4,x = -5.1457 -3.1727 5.7344,%用Cholesky分解求解 %A是对称矩阵 %L是单位下三角阵 %D是对角阵 %对称阵A进行三角分解：A=LDL,function x,L,D=Chol_decompose(A,b) N=length(A); L=zeros(N,N); D=zeros(1,N); for i=1:N L(i,i)=1; end D(1)=A(1,1); for i=2:N for j=1:i-1 if j=1 L(i,j)=A(i,j)/

29、D(j); else sum1=0; for k=1:j-1 sum1=sum1+L(i,k)*D(k)*L(j,k); end L(i,j)=(A(i,j)-sum1)/D(j); end end sum2=0; for k=1:i-1 sum2=sum2+L(i,k)2*D(k); end D(i)=A(i,i)-sum2; end,%分别求解线性方程组Ly=b;Lx=y/D y=zeros(1,N); y(1)=b(1); for i=2:N sumi=0; for k=1:i-1 sumi=sumi+L(i,k)*y(k); end y(i)=b(i)-sumi; end x=zero

30、s(1,N); x(N)=y(N)/D(N); for i=N-1:-1:1 sumi=0; for k=i+1:N sumi=sumi+L(k,i)*x(k); end x(i)=y(i)/D(i)-sumi; end,用Dollittle三角分解法求解方程组：, A=0.001,2,3;-1,3.712,4.623;-2,1.072,5.643; b=1,2,3; x,L,U=lu_decompose(A,b),x = -0.4904 -0.0510 0.3675 L = 1.0e+003 * 0.0010 0 0 -1.0000 0.0010 0 -2.0000 0.0020 0.001

31、0 U = 1.0e+003 * 0.0000 0.0020 0.0030 0 2.0037 3.0046 0 0 0.0059,function x,L,U=lu_decompose(A,b) %用Dollittle三角分解 %A是系数矩阵 %b表示方程组右边的向量 %L是单位下三角阵 %U是一般上三角阵 n=length(b); L=eye(n); U=zeros(n,n); for i=1:n U(1,i)=A(1,i); if i=1 L(i,1)=1; else L(i,1)=A(i,1)/U(1,1); end end,for i=2:n for j=i:n sum=0; for

32、k=1:i-1 sum=sum+L(i,k)*U(k,j); end U(i,j)=A(i,j)-sum; if j=n sum=0; for k=1:i-1 sum=sum+L(j+1,k)*U(k,i); end L(j+1,i)=(A(j+1,i)-sum)/U(i,i); end end end,%分别求解线性方程组Ly=b; y(1)=b(1); for k=2:n sum=0; for j=1:k-1 sum=sum+L(k,j)*y(j); end y(k)=b(k)-sum; end %解方程组Ux=y x(n)=y(n)/U(n,n); for k=n-1:-1:1 sum=

33、0; for j=k+1:n sum=sum+U(k,j)*x(j); end x(k)=(y(k)-sum)/U(k,k); end,function x,L,U=lr_decompose(A,b) %用Dollittle三角分解 %A是系数矩阵 %b表示方程组右边的向量 %L是单位下三角阵 %U是一般上三角阵 n = length(b); L = zeros(n,n); U = eye(n,n); %U的对角元素为1 L(1:n,1) = A(1:n,1); %L的第一列 U(1,1:n) = A(1,1:n)/L(1,1); %U的第一行 for k=2:n for i=k:n L(i,k) = A