第三章解线性方程组的直接法ppt课件.ppt

第三章 解线性方程组的直接法,第三章 解线性方程组的直接法,引言Gauss消元法列主元素消元法矩阵三角分解法向量和矩阵的范数误差分析,3.1 引言,小行星轨道问题：,天文学家要确定一小行星的轨道，在轨道平面建立以太阳为原点的直角坐标系。在坐标轴上取天文测量单位(一天文单位为地球到太阳的平均距离：9300万英里，约1.5亿千米)，对小行星作5次观察，测得轨道上5个点的坐标数据如下：,x 5.7640 6.2860 6.7590 7.1680 7.4800 y 0.6480 1.2020 1.8230 2.5260 3.3600,椭圆的一般方程: a1x2 + a2xy + a3y2 + a4x + a5y + 1 = 0,将数据逐个代入，可得五个方程的方程组，求解该线性方程组即可得行星轨道方程。,对一般线性方程组: A x = b, 其中,由以前所学内容知，当且仅当矩阵A行列式不为0时，即A非奇异时，方程组存在唯一解，可根据Cramer法则求解。,其算法设计如下：,（1） 输入系数矩阵A和右端向量b；,（2）计算系数矩阵A的行列式值D，如果D=0，则输出错误信息，结束，否则进行第(3)步；,（3） 对k=1,2,n，用b替换A的第k列数据，并计算替换后矩阵的行列式值Dk；,（4） 计算并输出x1 = D1 / D，x2 = D2 / D，xn=Dn/D， 结束。,但Cramer法则只适用于低阶方程组，高阶方程组工作量太大，故一般用数值方法求解。数值方法分两类：1. 直接法2. 迭代法,3.2 Gauss消元法,第二步: 回代过程， 解上三角形方程组，得原方程组的解。,基本思想：逐步消去未知元，将方程组化为与其等价的上三角方程组求解。,分两步：,第一步: 消元过程，将方程组消元化为等价的上三角形方程组;,Gauss消元的目的：,原始方程组,约化方程组,消元过程(化一般方程组为上三角方程组),以四阶为例：,其系数增广矩阵为：,第一轮消元：,计算3个数: m21 m31 m41T = a21 a31 a41T / a11,用-m21乘矩阵第一行后加到矩阵第二行;,用-m31乘矩阵第一行后加到矩阵第三行;,用-m41乘矩阵第一行后加到矩阵第四行;,其系数增广矩阵变为：,第二轮消元：,计算2个数：m32 m42T = a32(1) a42(1)T / a22(1),用-m32乘矩阵第二行后加到矩阵第三行;,用-m42乘矩阵第二行后加到矩阵第四行;,其系数增广矩阵变为：,第三轮消元：,计算: m43=a43(2)/a33(2),用-m43乘矩阵第三行后加到矩阵第四行;,其系数增广矩阵变为：,其对应的上三角方程组为,若对于一般的线性方程组Ax=b，其消元过程的计算公式为： （k=1,2,n-1）,回代过程(解上三角方程组),上三角方程组的一般形式为： 其中a11ann0,回代过程的计算公式：,工作量计算：,消去过程：,“”：第k步，n-k次，共(n-1)+(n-2)+1=n(n-1)/2,“”：第k步，(n-k)(n-k+1)次，共 (n-1)n+(n-2)(n-1)+12= (n3-n)/3,总工作量：S1=n(n-1)/2+ (n3-n)/3,回代过程：,“”：n,“”：1+2+(n-1)=n(n-1)/2,总工作量：S2=n+ n(n-1)/2= n(n+1)/2,（k=1,2,n-1）,故Gauss消元法的总工作量为：,S=S1+S2 =n(n-1)/2+ (n3-n)/3+ n(n+1)/2= n2+(n3-n)/3,若用克莱姆法则求解，则工作量为：,“”：（n+1个n阶行列式的值）(n+1)(n-1)n!,“”：n,故总工作量为： (n+1)(n-1) n!+n,如当n=6时, Gauss消元法工作量为106 ；而克莱姆法则求解工作量为25206。,1. Newton迭代法及其收敛性,前一次课内容回顾,2. Newton迭代法的变形（简化Newton迭代法、弦截法、Newton下山法）,3. Gauss消元法求解线性方程组,定理： 约化的主元素ak+1,k+1(k) 0 (k=0,1,n-1)的充分必要条件是矩阵A的各阶顺序主子式不为零。即,注：对角线上的元素ak+1,k+1(k)在Gauss消元法中作用突出，称约化的主元素。,推论: 如果A的顺序主子式Dk 0 (k=1,n-1)，则Gauss消元法中的约化主元可以表示为,x1= -13, x2 = 8, x3 = 2,m21=3/2m31=4/2,m32=-3/0.5,例 用高斯消元法求解方程组,矩阵的三角分解:,对线性方程组Ax=b的系数矩阵A施行初等行变换相当于用初等矩阵左乘A，故第一次消元后方程组化为A(1)x=b(1)，即L1Ax=A(1)x，L1b=b(1)，其中,同理LkA(k-1)=A(k) Lkb(k-1)=b (k),其中,将A 分解为单位下三角矩阵 L 和上三角矩阵 U的乘积的算法称为矩阵A的三角分解算法。,重复该过程，最后得,记U=A(n-1)，则,其中,定理：设A为n阶矩阵，若A的顺序主子式Di 0（i=1,2,n-1），则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且这种分解是唯一的。,由Gauss消元过程可推得,U即为Gauss消元后所得的上三角方程的系数矩阵。,例,解,由Gauss消元法可得，,m21=0， m31=2， m32= -1,故,如果已经有 A = L U 则 AX = b = L U X = b，,（1）求解方程组 LY = b 得向量 Y 的值； （L 是下三角矩阵，用顺代算法）,（2）求解方程组 UX = Y 得向量 X 的值。（U 是上三角矩阵，用回代算法）,记UX = Y ， LY = b ，则求解方程组分两步进行：,基本思想：Gauss消元法中，若主元 akk(k) 太小会使误差增大，故应避免采用绝对值小的元素作主元。最好每一步选取系数矩阵中（或消元后的低阶矩阵中）绝对值最大的元素作主元，以具较好的数值稳定性。,3.3 Gauss列主元素消元法,例：求解方程组,（用四位浮点数计算，精确解舍入到4位有效数字为：x1*=-0.4904,x2*=-0.05104,x3*=0.3675）,解：方法一Gauss消元法,（A/b）=,其中， m21=-1.000/0.001=-1000 m31=-2.000/0.001=-2000 m32=4001/2004=1.997 解为x1=-0.4000，x2=-0.09980，x3=0.4000(x1*=-0.4904，x2*=-0.05104，x3*=0.3675)显然，此解并不准确。,方法二交换行，避免绝对值小的主元作除数。,（A/b）=,其中， m21=0.5000 m31=-0.0005 m32=0.6300,解为x1=-0.4900，x2=-0.05113，x3=0.3678,(x1*=-0.4904，x2*=-0.05104，x3*=0.3675),与方法一相比，此解显然要精确得多。,设Axb的增广矩阵为,在A的第一列中选绝对值最大的元素作主元，设该元素所在行为第i1行，交换第一行与第i1行，进行第一次消元；再在第2n行的第二列中选绝对值最大的元素作主元，设该元素所在行为第i2行，交换第二行与第i2行，进行第二次消元，直到消元过程完成为止。,Gauss列主元素消元法的基本思想：,例：用列主元法解,解：第一列中绝对值最大是10，取10为主元,第二列的后两个数中选出主元 2.5,列主元矩阵的三角分解：,解：交换行变换,例：对矩阵A做列主元三角分解：,则列主元的Gauss变换可记为,A(2)=F2P2F1P1A,记 U=A(2)=F2(P2F1P2)P2P1A （因P2P2=I）,P = P2P1,则有,对于一般的n阶矩阵的列主元三角分解，通常令,定理：（列主元素的三角分解定理）若A非奇异，则存在排列阵P使PALU，其中L为下三角阵，U为上三角阵。,矩阵分解关系为,全主元素消元法：,定义：,则称,为全主元素。,经过行列互换，使得 位于经交换行和列后的等价方程组中的 位置，然后再实施消元。,全主元素消元法的基本思想：,若,注：全主元素消元法有可能改变未知数的顺序。,直接三角分解法：若将A分解为LU的积，则求Axb等价于求解两个三角形方程组：（1）Lyb，求y； （2）Uxy，求x。,3.4 矩阵三角分解法,（一） Doolittle分解法,（二 ） Crout分解法,（三） 对称正定矩阵的Cholesky分解法,（四） 三对角方程组的数值解法,设A非奇异，且ALU，L为单位下三角阵，U为上三角阵，即,可得直接三角分解法解Axb的计算公式：,一、Doolittle分解法：,u11=a11, , u1n=a1n,m21u11=a21, , mn1u11=an1,对A的元素aij ，当 jk 和 ik时,m21u12+ u22=a22, , m21u1n+ u2n=a2n,m21=a21/u11, , mn1=an1/u11,u22=a22 - m21u12, , u2n=a2n- m21u1n,m31u12+ m32u22=a32, , mn1u12+ mn2u22=an2,m32=(a32- m31u12)/u22, , mn2=(an2- mn1u12)/u22,矩阵L和矩阵U的元素计算公式：,计算秩序如图所示：,用Doolittle分解法求解线性方程组Axb的具体计算公式如下：,例：用Doolittle法解方程组,解:,由Doolittle分解,解Lyb，得,解Uxy，得,直接三角分解的Doolittle方法可用以下过程表示：,此方法称紧凑格式的Doolittle法。,从上式最后一个矩阵中可知L，U，y，然后解线性方程组Uxy。,例：用紧凑格式的Doolittle分解法解上例方程组,解：,所以,计算公式：,二、Crout分解法：,将矩阵A分解为如下形式：,例：求矩阵A的Crout分解：,所以,1. Gauss消元法中约化主元的作用与矩阵的三角分解,前一次课内容回顾,2. Gauss列主元素消元法,3. 矩阵的直接三角分解（Doolittle分解法、Crout分解法）,若n阶实矩阵A为对称正定矩阵，则：（1）ATA；（2）对任意的X0，有XTAX0；（3） A的各阶顺序主子式大于零。,三、对称正定矩阵的Cholesky分解法（平方根法）,故A可进行LU分解：,记,则 A=LDU1,即 ALDLT,称A的LDLT分解,若记,则,称A的LLT分解,= AT = U1TDLT,= U1=LT,LLT分解的计算公式：,对A的第i行元素, 有 ( j = 1,2,i ),( j =1,2,i-1 ),对于 i = 2,3, n,LDLT分解的计算公式：,d1 = a11, l21= a21/d1, , ln1 = an1/d1,( i = 2,3, n； j = 1,2,i-1 ),无需开方，故称改进的平方根法,应用Cholesky分解可将方程组Axb分解为两个三角形方程组：,而应用改进的Cholesky分解可将方程组Axb分解为下面的方程组：,例：用改进的平方根法解方程组Axb，其中,解：由,d1 = a11, l21= a21/d1, , ln1 = an1/d1,d1=1，l21=2，l31=1，l41=-3,得,又由,得,d2=1，l32=-2，l42=1，d3=9，l43=2/3，d4=1,因此，得到,解方程组Lyb，得,解方程组LTxD-1y，得,最终求得方程组Axb的解为 x=(1,1,1,1)T,四、 三对角方程组的数值解法,且 ai ci 0， (i = 2，3，n - 1),其中,三对角方程组形式如下：,计算公式：,Step 1：系数矩阵的三角分解：ALU,方法一：直接三角分解法,Step 3：求解上三角方程组Uxy。,计算公式：,计算公式：,Step 2：求解下三角方程组Lyf,例： 求解三对角方程组,解：,y=3, 8/5, 4/3, 2T,x = 5, 4, 3, 2T,基本思想：每一轮消元只对增广矩阵中某两行进行将主对角元化为1，主对角元下方元素化为零。,即将其系数增广矩阵化为如下形式：,方法二：（追赶法）,解三对角方程组Axb的追赶法分“追”和“赶”两个环节：,“追”的过程：（消元过程）,按顺序计算系数1- 2- n-1和y1-y2-yn;,“赶”的过程：（回代过程）,按逆序求出xn-xn-1-x2-x1。,其系数增广矩阵为：,i = ci / ( bi ai i 1) , ( i = 2, 3, , n - 1) yi = ( fi ai yi 1) / ( bi aii 1) , ( i = 2, 3, , n),例： 用追赶法求上例三对角方程组,解：,其系数增广矩阵为,追赶法求解,回代可得,3.5 向量和矩阵的范数,向量的范数矩阵的范数谱半径、谱范数与方阵的F范数,定义 设向量XRn ,若X 的某个非负实值函数N(X)=|X|满足条件：（1）非负性： |X|0且|X|=0的充要条件是X=0；（2）齐次性：|kX|=|k| |X|（k R）；（3）三角不等式：对X,YRn , 有 |X+Y|X|+|Y|；则称N(X)=|X|为Rn 上的向量X的范数。,3.5.1 向量的范数,定义 设X= (x1 , x2 , xn)T Rn ,则定义：,（1）向量的“2范数”：,（2）向量的“范数”：,（3）向量的“1范数”：,（4）向量的“p范数”：,例：计算向量X=(1,0,-1,2)T的三种常用范数。,解：,定理 （范数的等价性）对任何向量xRn ,有,证明（2）:,所以,（1）,（2）,（3）,向量序列的收敛性：,定义：在Rn中的一个向量序列X(1)，X(2)，X(k)，称为收敛于一个向量X，当且仅当,或,设X(k) 为Rn中的一向量序列， X*Rn ，记X (k)=(x1(k),x2(k),xn(k)T，X*=（x1*,x2*,xn*)T，若,i=1,2,n,则称X(k)收敛于向量X*，记为,定义 设矩阵ARnn ,若A的某个非负实值函数N(A)=|A|满足条件：（1）非负性： |A|0且|A|=0的充要条件是A=0；（2）齐次性：|kA|=|k| |A|（k R）；（3）三角不等式：对 A,BRnn , 有|A+B|A|+|B|；（4）相容性：|AB|A| |B|；则称N(A)=|A|为Rnn 上的矩阵A的范数。,3.5.2 矩阵的范数,定义 设向量XRn，矩阵ARnn ，且给定一种向量范数|X|p ，则定义(p=1,2,)为矩阵A的范数，并称为A的算子范数。,定义 对于给定的向量范数和矩阵范数，如果对任一个向量XRn和任一个矩阵ARnn 都有不等式|AX|A| |X|成立，则称所给矩阵范数与向量范数是相容的。,向量范数与矩阵范数的关系：,定理 设矩阵ARnn ,A=(aij )nn ,则,（1）,列范数,（2）,（3）,行范数,3.5.3 谱半径、谱范数与方阵的F范数,定义 设n阶方阵A的特征值为1 ,2 , n ，则 称,定理 （特征值上界）设ARnn ,则,为A的谱半径。（| i |是i的模）,即A的谱半径不超过A的任何一种算子范数。,定义 设A=(aij )nn Rnn ，则称,为A的Frobenius范数，简称F范数。,定理 若AR nn 为对称阵，则,故|A|2又称为谱范数。,3.6 误差分析,结论：（1）的常数项b的第二个分量只有1/1000的微小变化，方程组的解变化却很大。,例 记方程组（1）,为Axb，其精确解为：x1*=2，x2*=0,现考察方程组（2）,可将其表示为：A(x+x)=b+b,其中 b= (0,0.0001)T,设x为(1)的解，显然(2)的解为：x+x= (1, 1)T,设Axb的扰动方程组为(A+ A)(x+x)=b+b，其中A叫A的扰动矩阵， x和b叫x和b的扰动向量。,定义 若矩阵A或常数项b的微小变化引起方程组Axb的解的巨大变化，则称此方程组为病态方程组，A为病态矩阵（相对方程组而言）；否则称方程组为良态方程组，A为良态矩阵。,研究方程组中A或b的微小误差对解的影响的分析称“扰动分析”。,(1) A=0，则A (x+x)= b+b，减去Ax=b，得,设Ax=b的扰动方程组为(A+A) (x+x)= b+b，下面进行扰动分析：,(2) b=0，则(A+A) (x+x)= b，同理可得,A x= b，,故 x= A-1 b，,即|x| |A-1 | | b | ，,又由Ax=b，有|b| |A | | x | ，,所以,定义 设A非奇异，称数 cond(A)= |A-1 | |A|为矩阵A的条件数。,说明：,（1）条件数小，扰动引起的解的相对误差一定小；条件数大，扰动引起的解的误差可能很大。（条件数与所取的范数有关，最常用的是|和|2）,（2）由于cond(A)= |A| |A-1 |A A-1 |=|I|=1，故条件数是一放大的倍数，且总以1为下界。,（3）当A为正交矩阵（A-1=AT）时，有cond(A)2=1，所以正交矩阵的方程组是良态的。,可得x= (27.0000 -192.0000 210.0000)T x+x= (30.0000 -210.0000 228.0000)T,例：已知Hilbert矩阵,cond(H3)=748cond(H6)=2.9107cond(H7)=9.85108,计算可得,n越大，Hn的条件数越大，故Hilbert矩阵是一个典型的病态矩阵。,如,定理 （事后误差估计）设A非奇异，x*和x分别是Ax=b的精确解和近似解，r=b-Ax为残差向量，则,残差向量：r=b-Ax (x为Ax=b的近似解),当r很小时，x是否为Ax=b的一个较好的近似解呢？,证明：,由,又,得,P73-75习题三：1，2，3，4，6，8，9，13。,本章作业,