数值分析 Numerical Analysis.ppt

第六章第六章 解线性方程组的直接方法解线性方程组的直接方法6.1 6.1 高斯消去法高斯消去法6.2 6.2 高斯高斯主元素消去法主元素消去法6.4 6.4 误差分析误差分析6.3 6.3 矩阵的三角分解矩阵的三角分解1/5/2023数计学院数值计算课程建设组QAB【本章重点本章重点本章重点本章重点】1 1 1 1Gauss Gauss Gauss Gauss 消去法和列主元消去法及其实现条件。消去法和列主元消去法及其实现条件。消去法和列主元消去法及其实现条件。消去法和列主元消去法及其实现条件。2 2 2 2矩阵的三角分解，含矩阵的三角分解，含矩阵的三角分解，含矩阵的三角分解，含LULULULU分解和分解和分解和分解和LLLLLLLLT T T T 分解及三对角方程组的追分解及三对角方程组的追分解及三对角方程组的追分解及三对角方程组的追 赶法。赶法。赶法。赶法。3 3 3 3向量和矩阵范数的定义及性质。向量和矩阵范数的定义及性质。向量和矩阵范数的定义及性质。向量和矩阵范数的定义及性质。4 4 4 4矩阵条件数及病态矩阵定义和解方程组直接法的误差估计。矩阵条件数及病态矩阵定义和解方程组直接法的误差估计。矩阵条件数及病态矩阵定义和解方程组直接法的误差估计。矩阵条件数及病态矩阵定义和解方程组直接法的误差估计。1/5/2023数计学院数值计算课程建设组QAB 在自然科学和工程技术中许多问题的解决转化为解线性方在自然科学和工程技术中许多问题的解决转化为解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵，一种是高阶稀疏矩阵。解线性方程组的数值解也有两密矩阵，一种是高阶稀疏矩阵。解线性方程组的数值解也有两种，一种是种，一种是直接法直接法，就是经过有限步算术运算，可以求得线性，就是经过有限步算术运算，可以求得线性方程组的解，但实际计算时有舍入误差的存在和影响，所以求方程组的解，但实际计算时有舍入误差的存在和影响，所以求得的结果也只能是近似解对低阶稠密矩阵和部分大型稀疏矩阵得的结果也只能是近似解对低阶稠密矩阵和部分大型稀疏矩阵有效，另一种是有效，另一种是迭代法迭代法，就是用某种极限过程去逐步逼近精确就是用某种极限过程去逐步逼近精确解，是解决大型稀疏矩阵的重要方法。解，是解决大型稀疏矩阵的重要方法。引言引言1/5/2023数计学院数值计算课程建设组QAB 高斯高斯消去法是消去法是一个古老的求解线性方程组的方法，但由一个古老的求解线性方程组的方法，但由于它的改进、变形得到的选主元素消去法、三角分解法仍是于它的改进、变形得到的选主元素消去法、三角分解法仍是在计算机上求解系数在计算机上求解系数矩阵为中、低阶稠密矩阵的线性矩阵为中、低阶稠密矩阵的线性方程组方程组常用的有效方法，所以本节介绍这一方法。常用的有效方法，所以本节介绍这一方法。1 1 高斯消去法高斯消去法 用高斯消去法求解用高斯消去法求解n阶线性方程组阶线性方程组Ax=b的的基本思想基本思想是在是在逐步消元的过程中，把方程组的系数逐步消元的过程中，把方程组的系数矩阵化为上三角矩阵，矩阵化为上三角矩阵，从而将原从而将原方程组约化为容易求解的等价三角方程组，然后进方程组约化为容易求解的等价三角方程组，然后进行回代求解。行回代求解。一、一、高斯消去法高斯消去法1/5/2023数计学院数值计算课程建设组QAB例例1 用消去法求解方程组用消去法求解方程组1 1 1 6 4 -1 52 -2 1 11 1 1 6 4 -1 50 -4 -1 -111 1 1 6 4 -1 50 0 -2 -6还原为方程组还原为方程组（消元（消元过程）过程）（回代（回代过程）过程）x3=3，x2=2，x1=11/5/2023数计学院数值计算课程建设组QAB设有方程组设有方程组改写成矩阵形式：改写成矩阵形式：简记为：简记为：一般线性方程组的一般线性方程组的高斯消去法高斯消去法1/5/2023数计学院数值计算课程建设组QAB将原将原方程组记为方程组记为其中其中则则第一步（第一步（k=1)，若，若a11不等于不等于0，则可以计算乘数，则可以计算乘数用用-mi1 乘方程组的第一个方程加到第乘方程组的第一个方程加到第i个方程，则原方程组同个方程，则原方程组同解方程组为：解方程组为：其中其中简记为简记为1/5/2023数计学院数值计算课程建设组QAB第二步，设经过第二步，设经过k-1次消元后的同解方程组为次消元后的同解方程组为设设计算乘数计算乘数用用-mik乘上面的线性方程组的第乘上面的线性方程组的第k个方程加到第个方程加到第i个方程，可个方程，可以消去以消去xk元元,得到同解方程组得到同解方程组其中其中1/5/2023数计学院数值计算课程建设组QAB重复上述过程，可以将方程组化为等价的简单方程组重复上述过程，可以将方程组化为等价的简单方程组其中其中为为上梯形阵。上梯形阵。当当m=n时，与原方程组等价的方程组为时，与原方程组等价的方程组为即即（消元（消元过程）过程）由由上述方程组很快可以求出：上述方程组很快可以求出：（回代（回代过程）过程）1/5/2023数计学院数值计算课程建设组QAB定理定理5 设设Ax=b,其中其中AR nn(1)如果如果则可通过高斯消去法，将方程组化为等价三角形方程组，则可通过高斯消去法，将方程组化为等价三角形方程组，计算公式为：计算公式为：(a)消元计算消元计算(k=1,2,n)(b)回代计算回代计算(2)如果系数矩阵如果系数矩阵A为非奇异矩阵则可通过高斯消去法与初等行为非奇异矩阵则可通过高斯消去法与初等行变换方法，将方程组约化为三角形方程组。变换方法，将方程组约化为三角形方程组。1/5/2023数计学院数值计算课程建设组QAB定理定理6 约化的主元素约化的主元素的充要条件的充要条件是是矩阵矩阵A的的顺序主子式顺序主子式即即证明略（证明略（P171）推论推论 如果矩阵如果矩阵A的顺序主子式的顺序主子式则则1/5/2023数计学院数值计算课程建设组QABfunction 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(k,i)=0;endenddisp(A)disp(b)pausex(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程序设计程序设计1/5/2023数计学院数值计算课程建设组QABA=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 -6x=1 2 31/5/2023数计学院数值计算课程建设组QABa=10-100 1;2 1b=1 2gauss(a,b)a=0.0000 1.0000 2.0000 1.0000b=1 2消元后的A 1.0e+100*0.0000 0.0000 0 -2.0000消元后的b 1.0e+100*0.0000 -2.0000ans=0 1a=2 1;10-100 1b=2 1gauss(a,b)a=2.0000 1.0000 0.0000 1.0000b=2 1消元后的A 2 1 0 1消元后的b 2 1ans=0.5000 1.00001/5/2023数计学院数值计算课程建设组QAB二、矩阵的三角分解二、矩阵的三角分解由由上面消去法的分析知，若方程组的系数矩阵上面消去法的分析知，若方程组的系数矩阵A的的顺序主子式顺序主子式不等于不等于0，高斯消元的过程是作一系列初等行变换，等价于矩，高斯消元的过程是作一系列初等行变换，等价于矩阵阵A左乘一个初等矩阵。左乘一个初等矩阵。一般第一般第k步消元，步消元，A(k)化为化为A(k+1),b(k)化为化为b(k+1),相当于相当于所以第一次消元可表示为：所以第一次消元可表示为：其中：其中：1/5/2023数计学院数值计算课程建设组QAB重复以上过程，最后得到重复以上过程，最后得到重复以上过程，最后得到重复以上过程，最后得到令令为单位下三角矩阵。为单位下三角矩阵。1/5/2023数计学院数值计算课程建设组QAB定理定理7（矩阵的（矩阵的LU分解）分解）设设A为一个为一个n 阶矩阵，若其顺序主子式阶矩阵，若其顺序主子式Di0(i=1,2,n)。则。则A可分解为一个单位下三角形矩阵可分解为一个单位下三角形矩阵L和一个和一个上三角形矩阵上三角形矩阵U的乘积，且这种分解是唯一的。的乘积，且这种分解是唯一的。证明：存在性由定理证明：存在性由定理6 6 可得；可得；唯一性：设有两种分解，即唯一性：设有两种分解，即A=LU=L1U1所以由所以由L可逆，可逆，U1可逆，知可逆，知L-1 L1=UU1-1左边为单位下三角形矩阵，左边为单位下三角形矩阵，右边为上三角形矩阵，所以右边为上三角形矩阵，所以均为单位矩阵，即有均为单位矩阵，即有 L=L1 ,U=U1分解唯一。分解唯一。Ax=bAx=bLUxLUx=b=b记记UxUx=y,Lyy,Ly=b=b由由Ly=bLy=b求出求出y,y,再由再由UxUx=y=y求出求出x x注意：注意：LULU分解的作用是将一个较复杂的矩阵化为两个简单的矩阵。分解的作用是将一个较复杂的矩阵化为两个简单的矩阵。1/5/2023数计学院数值计算课程建设组QABA=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 0U=2 -2 1 0 4 -1 0 0 1Matlab函数：函数：L,U=lu(A)U为上三角矩阵为上三角矩阵（与理论略有不同）（与理论略有不同）A=1 1 1;0 4-1;2-2 1;L=1 0 0 0 1 0 2 -1 1U=1 1 1 0 4 -1 0 0 -21/5/2023数计学院数值计算课程建设组QAB 高斯消去法实质上将高斯消去法实质上将A A分解为一个下三角矩阵分解为一个下三角矩阵L L和一个和一个上三角上三角U U矩阵，即矩阵，即A=LUA=LU（1 1）L L为单位下三角矩阵，称为为单位下三角矩阵，称为DoolittleDoolittle分解（默认）；分解（默认）；（2 2）U U为单位下三角矩阵，称为为单位下三角矩阵，称为CroutCrout分解。分解。1/5/2023数计学院数值计算课程建设组QAB 高斯消去法在消元的过程中高斯消去法在消元的过程中,可能出现可能出现 的情况的情况,这时消去法将无法进行这时消去法将无法进行;即使主元素即使主元素 很小很小,用其作除数用其作除数,会导致其他元素数量级的严重增长和色、舍入误差的扩散。会导致其他元素数量级的严重增长和色、舍入误差的扩散。因此，为了减少误差，每次因此，为了减少误差，每次消元选取系数矩阵（或某列）消元选取系数矩阵（或某列）中绝对值最大的元素作为主元素，这就是全（列）主元素消中绝对值最大的元素作为主元素，这就是全（列）主元素消去法。去法。全主元素消去法在选取主元时要花费较多的机器时间，全主元素消去法在选取主元时要花费较多的机器时间，目前主要使用目前主要使用列主元素消去法列主元素消去法。2 2 高斯主元素消去法高斯主元素消去法 列主元素消去法的基本思想：依次按列选主元素，然列主元素消去法的基本思想：依次按列选主元素，然后换行使之变到主元素位置上，再进行消元计算。后换行使之变到主元素位置上，再进行消元计算。1/5/2023数计学院数值计算课程建设组QAB记作记作1/5/2023数计学院数值计算课程建设组QABPA=LU，U为上三角，为上三角，L为单位下三角。为单位下三角。1/5/2023数计学院数值计算课程建设组QAB例例2 用列主元消去法求解方程组用列主元消去法求解方程组1 1 1 6 4 -1 52 2 -2 1 1还原为方程组还原为方程组（消元（消元过程）过程）（回代（回代过程）过程）x3=3，x2=2，x1=12 -2 1 1 4 -1 51 1 1 62 -2 1 1 4 4 -1 5 2 0.5 6.52 -2 1 1 4 -1 5 1 31/5/2023数计学院数值计算课程建设组QABfunction 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 max0A=LUA为对称矩阵，且为对称矩阵，且A的所有各的所有各阶顺序主子式均不为零。阶顺序主子式均不为零。A=LU=LDU0,AT=(LDU0)T=U0TDLT=A由分解的唯一性得：由分解的唯一性得：U0T=L故：故：A=LDLTA为对称正定矩阵，对角矩为对称正定矩阵，对角矩阵阵D的元素均大于零。的元素均大于零。A=LDLT =LD1/2D1/2LT =(LD1/2)(LD1/2)T =L1L1T其中：其中：L1为下三角阵为下三角阵1/5/2023数计学院数值计算课程建设组QAB算法算法1 1 对正定矩阵对正定矩阵A A做做LLLLT T分解分解(平方根法平方根法)Ax=bLLTx=b记记LTx=y,Ly=b由由Ly=b求出求出y,再由再由LTx=y求出求出x公式：故对于公式：故对于j=1,2,j=1,2,n,n1 1）求）求L L2 2）由）由Ly=b求出求出y3 3）由）由LTx=y求出求出x1/5/2023数计学院数值计算课程建设组QAB算法算法2 2 对对A A做做LDLLDLT T分解分解(改进的平方根法改进的平方根法)Ax=bAx=bAx=bAx=bLDLLDLLDLLDLT T T Tx x x x=b=b=b=b记记记记DLDLDLDLT T T Tx x x x=y,Lyy,Lyy,Lyy,Ly=b=b=b=b由由由由Ly=bLy=bLy=bLy=b求出求出求出求出y,y,y,y,再由再由再由再由DLDLDLDLT T T Tx x x x=y=y=y=y求出求出求出求出x x x x公式：故对于公式：故对于i=1,2,i=1,2,n,n1 1）求）求L L和和D D2 2）由）由Ly=b求出求出y3 3）由）由D DLTx=y求出求出x1/5/2023数计学院数值计算课程建设组QAB解解:设设 由矩阵乘法得由矩阵乘法得 ，例例P231.9P231.91/5/2023数计学院数值计算课程建设组QAB四、四、追赶追赶法法三对角占优矩阵三对角占优矩阵1/5/2023数计学院数值计算课程建设组QAB对比得：对比得：1/5/2023数计学院数值计算课程建设组QAB计算流程：计算流程：Ax=fAx=fLUxLUx=f=f记记UxUx=y,Lyy,Ly=f=f先由先由Ly=fLy=f求出求出y,y,再由再由UxUx=y=y求出求出x x追赶法公式：追赶法公式：追赶法公式：追赶法公式：1.1.计算的递推公式：计算的递推公式：2.2.解解Ly=fLy=f：3.3.解解UxUx=y=y：1/5/2023数计学院数值计算课程建设组QAB由 例例P230.8P230.8由公式解:设1/5/2023数计学院数值计算课程建设组QAB由 得1/5/2023数计学院数值计算课程建设组QABfunction x=threedia(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);endy(1)=f(1)/d(1);for i=2:n y(i)=(f(i)-a(i)*y(i-1)/d(i);endx(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 3 2 5;c=-1-2-1 0;f=6 1 0 1 threedia(a,b,c,f)X=5 4 3 2X=5 4 3 21/5/2023数计学院数值计算课程建设组QAB4 4 4 4 误差分析误差分析误差分析误差分析一、矩阵范数一、矩阵范数设设 ，给出向量，给出向量范数范数 ，则，则定义：定义：为矩阵为矩阵A A的的(算子）算子）范数。范数。定义定义2 2 矩阵的矩阵的算子算子范数范数定义定义1 1 矩阵范数矩阵范数如果矩阵如果矩阵 的某个非负的实值函数的某个非负的实值函数N(A)=N(A)=满足条件：满足条件：则称则称N(A)N(A)是是 上的一个矩阵范数（或模）上的一个矩阵范数（或模）1/5/2023数计学院数值计算课程建设组QAB四种常用的矩阵范数：四种常用的矩阵范数：（1 1）列范数）列范数（2 2）2-2-范数范数（3 3）行范数）行范数（4 4）F-F-范数范数1/5/2023数计学院数值计算课程建设组QAB例：例：，计算，计算A A的各种矩阵范数。的各种矩阵范数。解：解：1/5/2023数计学院数值计算课程建设组QAB二、谱半径：二、谱半径：定义：定义：设设 的特征值为（的特征值为（i=1,2,n),i=1,2,n),称称为为A A的谱半径。的谱半径。定理定理：矩阵范数与谱半径之间的关系为：矩阵范数与谱半径之间的关系为:(A)|A|.定理定理：如果：如果 为对称矩阵为对称矩阵,则则例：例：，计算，计算A A的谱半径。的谱半径。1/5/2023数计学院数值计算课程建设组QAB三、矩阵的条件数三、矩阵的条件数定义定义 设是非奇异阵，称数设是非奇异阵，称数 为矩为矩阵阵A A的条件数。的条件数。常用的条件数有常用的条件数有:当当A为对称矩阵时为对称矩阵时,为为A的绝对值最大最小的特征值的绝对值最大最小的特征值.1/5/2023数计学院数值计算课程建设组QAB条件数的性质：条件数的性质：1.A1.A为非奇异阵，有为非奇异阵，有cond(A)v12.A2.A为非奇异阵且为非奇异阵且c0(常数常数),),则则cond(cA)v=cond(A)v3.A3.A为正交阵，有为正交阵，有cond(A)2=11/5/2023数计学院数值计算课程建设组QAB四、误差分析四、误差分析考虑线性方程组考虑线性方程组Ax=bAx=b中中A A或或b b的微小扰动（误差）对解的影响。的微小扰动（误差）对解的影响。精确解为：精确解为：x=(2,0)T解为解为：x=(1,1)T现考虑现考虑Ax=bAx=b中中b b2 2的微小扰动对解的影响。的微小扰动对解的影响。定义定义定义定义 如果矩阵如果矩阵A A或常数项或常数项b b的微小的变化，引起方程组的微小的变化，引起方程组Ax=bAx=b解的巨大变化，则称此方程组为解的巨大变化，则称此方程组为“病态病态”方程组方程组,矩阵矩阵A A称为称为“病态病态”矩阵，否则称此方程组为矩阵，否则称此方程组为“良态良态”方程组方程组,矩阵矩阵A A称称为为“良态良态”矩阵。矩阵。cond(Acond(A)1)1时，方程组是时，方程组是“病态病态”的的(A(A是是“病态病态”矩阵矩阵)1/5/2023数计学院数值计算课程建设组QAB定理定理1 设是非奇异阵设是非奇异阵，Ax=b0,且且A(x+x)=b+b,则则：定理定理2 2 设是非奇异阵设是非奇异阵，Ax=b0,且且(A+A)(x+x)=b,若若 则则：1/5/2023数计学院数值计算课程建设组QAB例例 已知已知n n阶希尔伯特（阶希尔伯特（HilbertHilbert）矩阵）矩阵(1)(1)计算计算H H3 3条件数条件数condcond(H(H3 3)(2)(2)考虑方程组考虑方程组H H3 3x=b=(11/6,13/12,47/60)x=b=(11/6,13/12,47/60)T T,设设H H3 3及及b b有有微小误差微小误差(取取3 3为有效数字为有效数字),),分析对解的影响分析对解的影响.1/5/2023数计学院数值计算课程建设组QAB解解 (1)(1)1/5/2023数计学院数值计算课程建设组QAB(2)(2)方程组H3x=b=(11/6,13/12,47/60)T的解为x=(1,1,1)T若H3及b有微小误差(取3为有效数字),则:1/5/2023数计学院数值计算课程建设组QAB