《第3章 解线性方程组的直接法精选PPT.ppt》由会员分享，可在线阅读，更多相关《第3章 解线性方程组的直接法精选PPT.ppt（52页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第3章 解线性方程组的直接法1第1页，本讲稿共52页解线性方程组有Gram法则，但Gram法则不能用于计算方程组的解，如n100，1033次/秒的计算机要算10120年本章讲解直接法解线性方程组的两类方法:直接法:经过有限次运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法(一般有限步内得不到精确解)2第2页，本讲稿共52页三角形方程组的解三角形方程组的解 (1)上三角方程的一般形式（n1）n/2次运算3第3页，本讲稿共52页 (2)下三角方程的一般形式 （n1）n/2次运算4第4页，本讲稿共52页3.1 顺序顺序Gauss消元
2、法消元法 基本思想基本思想：通过对：通过对(3.1)消元消元,逐步将逐步将(3.1)简化为同解的简化为同解的上三角方程上三角方程组组(称为消元过程称为消元过程),再由下而上地解此方程组，求得解再由下而上地解此方程组，求得解x(称为回代过程称为回代过程)。若若det（A）0，记其增广矩阵记其增广矩阵5第5页，本讲稿共52页步骤如下：第一步：运算量：(n-1)*(1+n)第二步：运算量：(n-2)*(1+n-1)=(n-2)n6第6页，本讲稿共52页第k步：类似的做下去，我们有：运算量：(nk)*(1nk1)=(nk)(nk2)n1步以后，我们可以得到变换后的矩阵为：7第7页，本讲稿共52页因此，
3、总的运算量为：加上 解上述上三角阵的运算量(n+1)n/2，总共为：注意到，计算过程中处在被除的位置，因此整个计算过程要保证它不为0所以，Gauss消元法的可行条件为：就是要求A的所有顺序主子式均不为0，即有以下定理8第8页，本讲稿共52页定理定理约化的主元素 的充要条件是矩阵 的顺序主子式即 证明证明 显然，当 时，定理成立.现设定理充分性对 是成立的，求证定理充分性对 亦成立.首先利用归纳法证明定理的充分性.设于是由归纳法假设有可用高斯消去法将 约化到 ，即 9第9页，本讲稿共52页且有（*）由设利用(*)式，则有定理充分性对 亦成立.显然，由假设利用(*)式亦可推出 10第10页，本讲稿
4、共52页顺序顺序Gauss消去法算法组织消去法算法组织 将增广矩阵将增广矩阵(A,b)放入一个二维数组放入一个二维数组.消去每一步算出的消去每一步算出的aij(k)放放入入aij的位置的位置,bi(k)放入放入bi,mik放在放在aik上上.消去完毕后消去完毕后,数组各位置上的数组各位置上的数据如下示数据如下示:回代过程的计算没有数据组织问题.Ab11第11页，本讲稿共52页 算法算法3.1 顺序顺序Guass消去法消去法 步步1 输入系数矩阵A，右端项b，置k:=1；步步2.2.消元：消元：对 k=1,2,n-1，计算步步3.回代回代对 k=n-1,1，计算程序见程序见P3812第12页，本
5、讲稿共52页 3.2.Gauss列主元消元法列主元消元法 解解:方法方法1 1 小主元小主元a11回代求解回代求解 x2=1.0000,x1=0.0000.计算结果相当糟糕计算结果相当糟糕.原因原因原因原因:求行乘数时用了求行乘数时用了”小主元小主元小主元小主元”0.0001作除数作除数.例例例例3.43.4 在在F(10,4,L,U)上用上用Gauss消去法求解线性方程组消去法求解线性方程组方程组的准确解为方程组的准确解为 x*=(10000/9999,9998/9999)T.13第13页，本讲稿共52页解解:方法方法2若上题在求解时将若上题在求解时将1,2行交换行交换,即即回代后得到回代后
6、得到 x=(1.0000,1.0000)T与准确解与准确解 x*=(9998/9999,10000/9999)T相差不大相差不大.主元主元a11 方法方法2所用的是在所用的是在GaussGauss消去法中加入选主元过程消去法中加入选主元过程,利用换行利用换行,避免避免”小主元小主元”作除数作除数,该方法称为该方法称为GaussGauss列主元消去法。列主元消去法。14第14页，本讲稿共52页 主元素是指主元素是指Gauss消元过程中位于矩阵消元过程中位于矩阵A(k-1)-1)的主对角线的主对角线(k,k)上位置上的元素上位置上的元素akk(k-1)-1)由于计算机和方程组本身的限制，在由于计算
7、机和方程组本身的限制，在Gauss消元过程中会出现零主元和小主元，而使消元过程中会出现零主元和小主元，而使Gauss消去法无法进行或消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素值较大的主元素.为了提高计算的数值稳定性为了提高计算的数值稳定性,在消元过程中采用选择主元在消元过程中采用选择主元的方法的方法.常采用的是常采用的是列主元列主元消去法和消去法和全主元全主元消去法消去法.给定线性方程组Ax=b,记A(1)=A,b(1)=b,列主元列主元GaussGauss消去消去法法的具体过程如下:首先在
8、增广矩阵B(1)=(A(1),b(1)的第一列元素中,取 15第15页，本讲稿共52页Gauss列主元消去法列主元消去法 因为因为|mi,k|1,列主元消去法能避免列主元消去法能避免”小主元小主元”作除数作除数,误差分误差分析与数值试验表明：析与数值试验表明：（Gauss）列主元消去法）列主元消去法”基本稳定基本稳定”.再在矩阵B(2)=(A(2),b(2)的第二列元素中,取 按此方法继续进行下去,经过n-1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回代求解.然后进行第二步消元得增广矩阵B(3)=(A(3
9、),b(3).易证,只要|A|0,列主元Gauss消去法就可顺利进行.16第16页，本讲稿共52页算法算法3.2（列主元消去法）步骤：（列主元消去法）步骤：1、输入矩阵阶数输入矩阵阶数n，增广矩阵增广矩阵 A(n,n+1);2、对于对于(1)按列选主元：选取按列选主元：选取 l 使使(2)如果如果 ，交换，交换 A(n,n+1)的第的第k行与第行与第l 行元素行元素(3)消元计算消元计算 ：3、回代计算回代计算17第17页，本讲稿共52页程序见程序见P42P42列主元Gauss消去法是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素.而全主元Gauss消去法是在每一步消元前,在所
10、有元素中选取绝对值最大的元素作为主元素.但由于运算量大增,实际应用中并不经常使用.18第18页，本讲稿共52页3.3 3.3 解三对角方程组的追赶法解三对角方程组的追赶法Ax=d其中：其中：(3.7)(3.7)19第19页，本讲稿共52页用用Gauss消去法解三对角线性方程组：消去法解三对角线性方程组：其中其中追追赶：赶：20第20页，本讲稿共52页解三对角方程组的追赶法其计算工作量为6n-5次乘除法。工作量小，有以下定理可得证三对角矩阵求解的充分性条件。程序见P4421第21页，本讲稿共52页3.4 3.4 矩阵矩阵LULU分解法分解法 顺序顺序GaussGauss消去法的矩阵表示消去法的矩
11、阵表示矩阵若令 记 22第22页，本讲稿共52页则有 A(2)=记 令若则有A(3)=23第23页，本讲稿共52页如此进行下去,第n-1步得到:A(n)=其中也就是:A(n)=Ln-1A(n-1)=Ln-1Ln-2A(n-2)=Ln-1Ln-2L2L1A(1)24第24页，本讲稿共52页其中所以有:A=A(1)=L1-1L2-1Ln-1-1A(n)=LU其中L=L1-1L2-1Ln-1-1,U=A(n).而且有25第25页，本讲稿共52页 式 A=LU称为矩阵A的三角分解三角分解.设n阶方阵A的各阶顺序主子式不为零,则存在唯一单位下三角矩阵L和上三角矩阵U使A=LU.证明证明 只证唯一性,设有
12、两种分解 A=LU则有=I所以得于是 Ax=b LUx=b 令 Ux=y 得定理定理26第26页，本讲稿共52页LU分解方式直接利用矩阵乘法来计算直接利用矩阵乘法来计算 LU分解分解比较等式两边的比较等式两边的第一行第一行得：得：u1j=a1j比较等式两边的比较等式两边的第一列第一列得：得：比较等式两边的比较等式两边的第二行第二行得：得：比较等式两边的比较等式两边的第二列第二列得：得：(j=1,n)(i=2,n)(j=2,n)(i=3,n)U 的第一行的第一行 L 的第一列的第一列 U 的第二行的第二行 L 的第二列的第二列 27第27页，本讲稿共52页LU分解算法（续）第第 k 步：步：此时
13、此时 U 的前的前 k-1 行和行和 L 的前的前 k-1 列已经求出列已经求出比较等式两边的比较等式两边的第第 k 行行得：得：比较等式两边的比较等式两边的第第 k 列列得：得：直到第直到第 n 步，便可求出矩阵步，便可求出矩阵 L 和和 U 的所有元素。的所有元素。(j=k,n)(i=k+1,n)28第28页，本讲稿共52页LU分解算法（续）算法算法 3.3：(LU 分解分解)For k=1,2,.,nEnd Forj=k,ni=k+1,nMatlab程序参见：程序参见：malu.m为了节省存储空间，通常用为了节省存储空间，通常用 A 的绝对下三角部分来存放的绝对下三角部分来存放 L(对角
14、线元素无需存储对角线元素无需存储)，用，用 A 的上三角部分来存放的上三角部分来存放 U。运算量：运算量：(n3-n)/329第29页，本讲稿共52页由可得这就是求解方程组Ax=b b的Doolittle三角分解方法。30第30页，本讲稿共52页 利用三角分解方法解线性方程组 解解 因为所以例例31第31页，本讲稿共52页先解 ,得再解 ,得 解线性方程组Ax=b的Doolittle三角分解法的计算量约为n3/3,与Gauss消去法基本相同.其优点在于求一系列同系数的线性方程组Ax=bk,(k=1,2,m)时,可大大节省运算量.此外，还有另一种LU分解，称为Crout三角分解.32第32页，本
15、讲稿共52页3.5 对称正定矩阵的Cholesky分解法 A对称：AT=A A正定：A的各阶顺序主子式均大于零。即定理定理：（对称正定矩阵的三角分解）（对称正定矩阵的三角分解）如果如果A为对称正定矩阵为对称正定矩阵，则存在一个实的非奇异则存在一个实的非奇异下三角矩阵下三角矩阵L，使使A=LLT ，且当限定的对角元素为正时且当限定的对角元素为正时，这种分解是唯一的这种分解是唯一的。33第33页，本讲稿共52页U=uij=u11uij/uii111u22unn记为记为 A 对称对称即即记记 D1/2=则则 仍是下三角阵仍是下三角阵 设A为对称正定矩阵,则有唯一分解A=LU,且ukk0.分解A=LL
16、T称为对称正定矩阵的Cholesky分解分解.如何求L34第34页，本讲稿共52页对矩阵对矩阵A进行进行Cholesky分解分解，即即A=LLT，由矩阵乘法由矩阵乘法：对于对于 i=1,2,n 计算（一列一列算）计算（一列一列算）求解下三角形方程组求解下三角形方程组 Ly=b 求解求解LTx=y35第35页，本讲稿共52页Cholesky分解法缺点及优点 2.2.对于对称正定阵对于对称正定阵 A，从，从 可知对任意可知对任意k i 有有 。即即 L 的元素不会增大，误差可控，不需选主元。的元素不会增大，误差可控，不需选主元。缺点：存在开方运算，可能会出现根号下负数。优点：1.可以减少存储单元。
17、为一般矩阵LU分解计算量的一半。36第36页，本讲稿共52页改进Cholesky分解法改进的cholesky分解A=LDLT37第37页，本讲稿共52页改进的cholesky分解逐行相乘，并注意到ij，有由此可得38第38页，本讲稿共52页改进的cholesky分解所以将上述公式改写为解方程，令解方程，令LTX=y，先解下三角形方程组先解下三角形方程组LDY=b得得再解上三角形方程组再解上三角形方程组LTX=Y得得 39第39页，本讲稿共52页算法算法 3.4：(Cholesky 分解法分解法)步步1 输入对称正定矩阵A，右端项b；步步2.Cholesky分解：分解：对 k=2,3,n，计算：
18、步步3.3.解下三角形方程组LDy=bLDy=b步步4.4.解上三角形方程组L LT Tx=yx=y程序见程序见P5340第40页，本讲稿共52页3.6 线性方程组的性态和舍入误差对解的影响线性方程组的性态和舍入误差对解的影响先看一个例子，设线性方程组试分析系数矩阵和右端项有微小扰动，解将产生什么样的变化？解 方程组的精确解为x=(1,1)T设系数矩阵有微小扰动即41第41页，本讲稿共52页若右端项有微小变动则若同时对A,b扰动A，b，则42第42页，本讲稿共52页现计算相对误差：43第43页，本讲稿共52页病态方程组病态方程组由于计算机字长限制，在解Ax=b时，舍入误差是不可避免的。因此我们
19、只能得出方程的近似解x*。x*是方程组(A+A)x=b+b (1)的准确解。在没有舍入误差的解。称方程（1）为方程Ax=b的扰动方程。其中A，b为由舍入误差所产生的扰动矩阵和扰动向量。当A，b的微小扰动，解得（1）的解与Ax=b的解x的相对误差不大称为良态方程，否则为病态方程。44第44页，本讲稿共52页 设线性方程组 Ax=b1)系数矩阵是精确的,常数项有误差b,此时记解为x+x,则 A(x+x)=b+b于是 Ax=b所以 x=A-1b A-1 b又由于 b=Ax A x因此 x b AA-1bx即误差界误差界45第45页，本讲稿共52页 2)再设b是精确的,A有误差A,此时记解为x+x,则
20、 (A+A)(x+x)=b则有 Ax+A(x+x)=0所以 x=-A-1A(x+x)于是 x A-1Ax+x也就是46第46页，本讲稿共52页3)设Ax=b中A和b都有误差(小扰动),此时记仍解为x+x,则 (A+A)(x+x)=b+b由Ax=b，有 Ax+A x=b-Ax即 x+A-1 Ax=A-1 b-A-1 Ax记 Cond(A)=AA-1,称为方程组Ax=b或矩阵A的条件数.47第47页，本讲稿共52页48第48页，本讲稿共52页条件数的性质：条件数的性质：）cond(A)1）cond(kA)=cond(A)k 为非零常数为非零常数）若）若 ，则则经常使用的条件数有 Condp(A)=
21、ApA-1p p=1，2，。当A为对称矩阵时，有 Cond2(A)=|1|/|n|其中1，n分别是A的按绝对值最大和最小的特征值。49第49页，本讲稿共52页 例如，对前面方程组的系数矩阵A有 Cond1(A)=Cond(A)=39601，Cond2(A)39206 由于计算条件数运算量较大,实际计算中若遇到下述情况之一,方程组就有可能是病态的:行列式很大或很小（如某些行、列近似相关）；行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；元素间相差大数量级，且无规则；主元消去过程中出现小主元；主元消去过程中出现小主元；特征值相差大数量级。特征值相差大数量级。50第50页，本讲稿共52页病态病态方程组的处理方程组的处理病态方程组的计算1)用双精度或更高精度计算2)采用数值稳定性较好的算法，如全主元法3)用迭代改善法 设已求得方程组Ax=b的近似解x1,计算剩余向量 r1=b-Ax1则x1的迭代改善解为:再求解余量方程组Ax=r1,得到解x*,Step 1:近似解近似解Step 2:Step 3:Step 4:若若 可被精确解出，则有可被精确解出，则有 就是精确解了。就是精确解了。51第51页，本讲稿共52页练习题练习题第第57-8页页 习题习题3.3-3.73.1552第52页，本讲稿共52页