矩阵A紧凑格式的Doolittle分解为PPT幻灯片

第三章线性方程组的解 概述 消元法 三角分解法 平方根法 向量与范数 方程组的性态与误差分析 迭代法3.1概述大量实际计算问题＝>一组线性方程组＝>如何求解？1.术语 非齐次线性方程组： 若记： 则：(1.1)式可 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示为Am×nx=b－－－线性方程组的矩阵表示 分类：据方程的个数与未知数的个数间关系，分为： m<n，即方程的个数小于未知数的个数，称为亚定方程组，有无穷多组解 m>n，即方程的个数大于未知数的个数，称为超定方程组，或矛盾方程组。它没有一般意义下的解，但可以求其广义解。 m=n，一般意义下的方程组，本章主要讨论的重点 面临的问题： 方程组Ax=b有没有解？ 有多少解？ 如何求解？ 2.Crammer法则求解An×nx=b 方程组Ax=b有唯一解的充要条件是：|A|≠0A可逆A是非奇异矩阵 Ax=b在A可逆时，存在唯一解 结论：Crammer法则计算量非常大，需算n+1个n阶行列式。如：n＝100，1000次/秒的计算机要算10120年 3线性方程组的数值解法： 直接法：思想：经有限步运算求得精确解的方法：舍入误差，因此也是近似解高斯消元法及其变形特点：它可靠且效率高，但它只适用于中小型方程组。 迭代法：思想：构造适当的初始近似解向量x，按一定的法则，使之逐步逼近精确解，得到一个满足精度要求的近似解。即是用某种极限过程逐步逼近准确解的方法。主要方法：Jaccobi迭代法，Gauss-Sidel迭代法等3.1.1GAUSS消元法一.几种可以直接求解的线性方程组 对角矩阵： ＝>n次除法 下三角矩阵： 即：l11x1=b1 l21x1+l22x2=b2 l31x1+l32x2+l33x3=b3 …… li1x1+li2x2+……＋liixi=bi …… ln1x1+ln2x2+……+lnnxn=bn 上三角矩阵： 即： u11x1+u12x2+………+u1nxn=b1 u22x2+………+u2nxn=b2 …… uiixi+ui,i+1xi+1+…+ui,nxn=bi un-1,n-1xn-1+un-1,nxn=bn-1 un,nxn=bn二、同解变换 初等方阵： :rirj:ri*k 注 对矩阵A实施初等行变换对A左乘初等方阵 初等方阵是可逆的，多个初等方阵的积仍然可逆 可逆阵A经有限次初等行变换单位矩阵:rj+k*ri 2.对方程组Ax=b作如下的变换，解不变: 交换两个方程的次序 一个方程的两边同时乘以一个非0数 一个方程的两边………………………，再将之加到另一个方程 将方程组Ax=b对应的增广矩阵（A,b)，作如下变换，解不变 交换矩阵的两行 某行乘以一个非0数 某行乘以一个非0数，加到另一行增广矩阵第i行第i个方程三、Gauss消元法消元法基本思想：对Ax=b的增广矩阵（A|b)进行初等变换，变成可直接求解的三种形式之一，再求解 Gauss消元法：将A化为上三角阵回代求解用高斯消去法解下列线性方程组 原方程组对应的增广矩阵为：1.消元步骤方程组AX=b的矩阵表示为：A(1)x=b(1) 初始状态： 第一次消元：消去对角线下第1列元素（a11≠0)方法： 消去对角线下第1列元素为0，即ri-r1×ai1/a11：（i=2,3,… 第二次：若a22≠0，则消去对角线下第二列为0， 即：ri-r2×ai2/a22(i=3,4,…..) 经过k-1次消元后，增广矩阵变为： 第k次消元：若akk≠0，则消去对角线下第k列元素 即：ri-rk×aik/akk(i=k+1,k+2,……) 2.经n次消元后得到同解方程组A(n)x=b(n)，且A(n)为上三角矩阵，可逐步回代求解即：UX=Y3.算法步骤将方程组用增广矩阵（A|b)＝（aij)n×(n+1)表示，注意c语言中数组下标从0开始，故i=0,1,…,n-1，j=0,1,…,n 消元：对k＝0，1，……，n-2（消去第k列对角线下的元素） 计算li,k=ai,k/ak,,k 第i行－第k行×li,k（i=k+1,k+2,…n-1) 回代求解： xi=(ai,n-∑xjai,j)/ai,i.i=n-1,….,0,j=i+1,i+2,…..n-14.时间复杂度消元过程中： 第k次消元每行均需计算第k行要乘以的常数：lik=aik/akk，消去第i行（第k列）时，i行各列元素加上k行各列元素乘以常数lik：aij=aij-akj*lik(j=k,k+1,…,n)，故需做n-k+1次乘法、加减法一共要做n-k行综上：第k次消元，除法：n-k次，乘法：（n-k)*(n-k+1)次，加减法：（n-k)*(n-k+1)次 整个消元过程：除法：乘法、加减法： 综上：算法时间度量为O(n3)5.Gauss顺序消元法的适用条件GAUSS直接消元法的适用条件：akk(k)≠0 A的各阶顺序主子式均不为0时，有akk(k)≠0（证明略） A是严格对角占优矩阵，则：akk(k)≠0（证明略） 举例说明A的每行主对角元的绝对值>同行其余元素的绝对值之和 例2：用Gauss消元法解方程组（4位十进制机）0.0001x1+x2=1x1+x2=2 因要对阶，损失有效数字，求得x1=0,x2=1 显然x1,x2不是方程的解，怎么办？？？ x1+x2=2 0.0001x1+x2=1 解之，得到x1=1,x2=1，接近与精确解 3.1.2列主元Gauss消元法 当采用Gauss消元法求解时，若|akk|《|aik|,会导致误差扩散 第k步消元时，ri-rk*aik/akk=aij-akj*aik/akk 因为|akk|<<|aik|，若rk行的误差为e=(ek+1,….,en+1)，则e将被扩大（aik/akk)倍Gauss消元法在特殊情况下是不稳定的 解决办法： 选择列主元：即选出第k列对角线下绝对值最大者设|atk|=max|aik|k≤i≤nrkrt以新的akk作为主对角元进行消元 故列主元消元法只是在每步消元之前选出主对角元 举实例说明列主元Gauss消元法的计算过程 列主元算法，只需在消元的第1步前增加以下几行代码i=k;//选择第k列的主元 for(j=k+1;j<n;j++) if(fabs(a[j][k])>fabs(a[i][k]))i=j; if(i!=k) for(j=k;j<=n;j++) { tmp=a[k][j]; a[k][j]=a[i][j]; a[i][j]=tmp; }3.1.3Gauss-Jordan消元法 算法思想： 在Gauss消元的第k次消元时，只是利用akk，将第k行以下各行的第k列消元为0，最终将系数矩阵消元为s上三角阵。 将上述算法改进：第k次消元时，用akk将第k列的所有元素（akk除外）均消元为0，则系数矩阵最终成为对角矩阵（或单位矩阵） 第k步消元：若akk≠0，则消去除对角线外的第k列元素 方法： 计算：lik，即 消去第k列元素（akk除外）为0，即ri-rk×lik第k-1次消元结果第k次消元结果经过n次消元，原方程的增广矩阵变为：则，方程的根为： 示例：Gauss-Jordan消元法求解Gauss-Jordan消元法的归一化 消元后对角线元素全为1，即单位矩阵E 方法：第k次消元时，先将akk化为1(第k行除以akk)，再用第k行消去各行k列的元素 经过n-1次消元，则原方程的增广矩阵变为：则，方程的根为： 归一化Gauss-Jordan消元法步骤： 消元：（c/c++实现）对k＝0，1，……，n-1（消去除对角线外第k列） 第k行除以akk，将akk化为1 第i行－第k行×aik（i=0,1,…n-1，且i≠k) 增广矩阵最后一列ai,n即为方程组的解，无需回代 示例：用归一化的Gauss-Jordan消元法求解 归一化的Gauss－Jordan消元法的应用 可逆阵A经有限次初等行变换单位矩阵 故可利用Gauss－Jordan消元法，构造矩阵 （A|E)Gauss-Jordan消元(E|A-1) 示例：用Gauss-Jordan消元法求矩阵的逆矩阵3.2矩阵的三角分解法Gauss消元法的矩阵表示 第一次消元 即：第一步消元，相当于左乘初等变换矩阵L1，相当于：L1A(1)=A(2)，L1b=b(2) 第二次消元：即：L2A(2)=A(3)，L2L1A(1)=A(3) 经过n-1次消元后 Ln-1Ln-2…L2L1A=A(n)，Ln-1Ln-2…L2L1b=b(n) 故：A=L1-1L2-1…Ln-2-1Ln-1-1A(n），令L=L1-1L2-1…Ln-2-1Ln-1-1，则： 故：A=LA(n）=LU Gauss消元法的实质：将系数矩阵分解为两个三角矩阵的乘积 将矩阵A分解成一个上三角、下三角矩阵的乘积，即：A=LU，矩阵的三角分解、LU分解 三角分解是从A的元素直接得到L、U的元素，不用中间步骤，即不用消元，所以称直接三角分解，或直接分解。 A的LU分解只要求L是一个下三角阵、U是一个上三角阵，当A有LU分解时，这种分解不是唯一的。当L是一个单位下三角阵或者U是一个单位上三角阵时，A的这两种LU分解是唯一的。 上述是A的两种不同的三角分解 一般地，若A=LU是一个三角分解，任取与A同阶的 非奇异对角矩阵D，则A=(LD)(D-1U)=L1U1 也是A的三角分解。 常用的两种分解： L为单位下三角阵而U为一般上三角阵的分解称为Doolittle分解一般上三角单位下三角 L为一般下三角阵而U为单位上三角阵的分解称为Crout分解 一般下三角单位上三角定理：n阶矩阵A存在唯一的杜利特尔分解、克洛特分解的充要条件是A的顺序主子矩阵Ak(k=1,2,…,n-1)非奇异（各阶顺序主子式非0）证明略（反证法）若不唯一，则可设A=L1U1=L2U2，推出 因下三角矩阵的积仍然是下三角阵，上三角矩阵的积仍然是上三角阵 下三角阵=上三角阵，则两个必为单位矩阵1.Doolittle分解 据矩阵的乘法，知： A的第一行元素：，故 第一列元素：，故： A的第二行元素： 故： 第二列元素： 故：＝A=LU 一般地： 推而广之：比较A第k行： 求U的第k行 比较A第k列求L的第k列 即： 计算顺序为： 为了便于记忆，将A的LU分解写成紧凑格式：123456.....对矩阵A进行Doolittle分解 解：矩阵A紧凑格式的Doolittle分解为：对矩阵A进行Doolittle分解，并求|A| 解：-12-3三角分解法解线性方程组方程组Ax=b，若A非奇异，且有三角分解A=LU，则原方程组等价于：LUx=b 令：Ux=y,则Ly=b 故可先解出y，再求x 所以等价于求U阵第k行元素ukj的公式，故可将增广矩阵（A|b）同时做LU分解，则最后一列即为y，再用Ux=y求解例：用Doolittle分解法求解线性方程组： 解：对增广矩阵进行Doolittle分解 故：Ux=y为： 所以x3=1，x2=-2，x1=22用三角分解法求矩阵A的逆矩阵A-1若AX=E，则X=A-12.Crout分解类同于Doolittle分解 比较两边可得： 紧凑格式的Crout分解为：214365.....原方程：AX=b变为：LUX=b 令Ux=y，则Ly=b， 可见yk的表达式完全类同于ukj的表达式，故可将常向量放在系数矩阵之后，对整个增广矩阵进行Crout三角分解 则UX=y的解即为线性方程组的解： 即：例3-11：分别用Doolittle、Crout分解法求解线性方程组 解一：Doolittle分解：1/2y向量 则解Ux=y可得线性方程组的解 解之，得：x3=2，x2=-1，x1=1解二：Crout分解法 则Ux=y的解即为原方程组的解 解之，得：x3=2，x2=-1，x1=13/22注意：此处与Doolittle方法的区别3.追赶法解三对角线性方程组概念： 三对角矩阵：除主对角线、次主对角线上的元素外，其余元素均为0对系数矩阵进行Crout分解： 比较矩阵两边可得： 故：|A|=|L|*|U|=l1*l2*……*ln 定理3-6：定理：若A为对角占优的三对角阵，即：满足 则方程组有唯一的LU分解 证明：要证三对角方程组具有唯一的LU分解 即要证：三对角方程组的各阶顺序主子式均≠0 也即要证：该方程组组的Crout分解中lii≠0 由于三对角矩阵的特殊性，无需用二维数组存放系数矩阵，而采用一维数组A,B,C,D存放所有的非0元 求解步骤： 三角分解：A=LU 则方程组Ax=fLUx=f，令：Ux=y，则Ly=f 故先解方程组Ly=f，求出y，再解Ux=y即可 追：由Ly=f，得： 赶：再从Ux=y，回代求解：示例： 23.3平方根法： 系数矩阵为对称正定矩阵的方程组称为对称正定方程组 对称正定方程组可用高斯消去法、LU分解法求解 也可导出计算量更小的平方根法 利用对称正定矩阵的三角分解（Cholesky）求解对称正定方程组的方法称为平方根法 对称矩阵：A=AT 对称正定矩阵A=AT，且对任意非零向量x有：f(x)=xTAx>0 对称正定阵A的几个重要性质 A1亦对称正定，且aii>0 A的顺序主子阵Ak亦对称正定 A的特征值i>0 A的全部顺序主子式det(Ak)>0 对称正定矩阵的乔累斯基分解 对称正定矩阵A=Rn×n，存在唯一的单位下三角矩阵L和对角矩阵D，使A=LDLT证明：前面已经证明了对称正定矩阵A做Doolittle分解的唯一性，设A=LU（其中L为单位下三角矩阵，U为上三角矩阵）令 ，其中di为U阵对角线元素，故di>0 则：A=LD(D-1U)，AT=(D-1U)TDTLT=(D-1U)TDLT，又A=AT，故：LD(D-1U)=(D-1U)DLT 由上式知：L=(D-1U）T，LT=D-1U 所以A=LDLT 对称正定矩阵A=Rn×n，存在唯一的主对角元素均为正数的下三角矩阵L，使A=LLT，此即：对称正定矩阵的乔累斯基分解 证明：因A可唯一分解为A=LDLT，且D是di>0的对角矩阵故可将D进行拆分令：比较的等式两边，可得第j列元素： 计算顺序为：l11li1li2…lin |A|=|L||LT|=l112…lnn2示例：对矩阵做乔累斯基分解，并求|A| 乔累斯基分解的缺点：开方改进的平方根法： 利用已经证明的结论：对称正定矩阵A，可唯一分解为：A=LDLT 故：d1=u11，d1l12=u12，…，dilik=uik则实对称正定矩阵A=LDLT，可视为A=LU分解，由于A的对称性，致使L、U之间具有以下关系： 则写成紧凑格式为： 此即为改进平方根法：在LU分解过程中，利用系数矩阵的对称性，作Doolittle分解，只需计算U矩阵，L矩阵第k列的值可由U矩阵k行的值除以ukk算得，使得LU分解的计算量减少了一半 其他可完全按Doolittle分解求解线性方程组的解的方法进行示例3-17：用改进平方根法解线性方程组 解：因系数矩阵是对称正定矩阵（需判断），用改进平方根法对方程组的增广矩阵做Doolittle分解2则A=LU，Ax=bLUx=b，令y=Ux，则Ly=b，A矩阵改进平方根法LU分解的最后一列即为y向量 所以解Ux=y 得：x3=0，x2=-1，x1=1