线性方程组直接解法.pptx

《线性方程组直接解法.pptx》由会员分享，可在线阅读，更多相关《线性方程组直接解法.pptx（79页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、西安电子科技大学理学院 主讲:王卫卫通过某种迭代系统（公式）求得近似解，优点：编程简单 缺点：存在收敛性和收敛速度问题How to get the solution？Coefficient matrix A低阶稠密阵 高阶稀疏阵small dense matrix large sparse matrixDirect methodsIteration methods Gaussian elimination列列/行行/完全主元素完全主元素(pivoting)消去法消去法Gauss-Jordan eliminationSquare root/improved square root methods

2、追赶法追赶法Jaccobi iterationGauss-Sidel iterationSORExistence and uniqueness of the solution？Cramer rule:Computation cost:(n+1)!上万阶，零元素很多，非零元素很少非零元素较多，零元素较少经过有限步算术运算直接求得精确解（在没有舍入误差的情况下），但实际上机器总存在舍入误差，因此求得的是近似解第1页/共79页西安电子科技大学理学院 主讲:王卫卫Gaussian elimination-通过初等变换将原方程组化成三角方程租来求解-2*(1)+(3)(2)+(3)回代求解第2页/共79

3、页西安电子科技大学理学院 主讲:王卫卫Step 1.Denote Ax=b as Suppose let-li1*(1)+(i),i=2,nGeneral procedure of Gaussian elimimation第3页/共79页西安电子科技大学理学院 主讲:王卫卫初等行变换,相当于左乘初等行变换矩阵-li1*row 1+row i,i=2,nformulaeDirectly replace with 0Leave aloneNeed to be updatedConvenient in using Matlab,Mathmatica,Maple第i行第1列第4页/共79页西安电子科技

4、大学理学院 主讲:王卫卫Step k.After k-1 eliminations,we have letSuppose-lik*(k)+(i),i=k+1,n第5页/共79页西安电子科技大学理学院 主讲:王卫卫-lik*row k+row i,i=k+1,n初等行变换,相当于左乘初等行变换矩阵Leave aloneNeed to be updatedDirectly replace with 0第i行第k列第6页/共79页西安电子科技大学理学院 主讲:王卫卫formulaek=1,n-1Convenient in using Matlab,Mathmatica,Maplek=1,n-1消元过

5、程第7页/共79页西安电子科技大学理学院 主讲:王卫卫对角线元素全为1，非对角线上除第i行第k列元素为lik外其余元素全为0第8页/共79页西安电子科技大学理学院 主讲:王卫卫要使要使Gauss消去法能够进行下去，必消去法能够进行下去，必须有约化后的主对角元素非零，即须有约化后的主对角元素非零，即矩阵矩阵A在什么条件下能够保证此条件在什么条件下能够保证此条件成立？成立？Lemma 1A的顺序主子式第9页/共79页西安电子科技大学理学院 主讲:王卫卫(n-k)kkk,det=Dikk,det=1kk,det=0k(n-k)证明要点证明要点第10页/共79页西安电子科技大学理学院 主讲:王卫卫In

6、ferenceTheorem 1 可通过Gauss消去法将线性方程组化为三角方程组第11页/共79页西安电子科技大学理学院 主讲:王卫卫 In Gaussian elimination,if a certain diagonal element vanishes,then you cant continue eliminating.Choose column/row/overall main element!Right.Even if it is very small,you cant do that because it will cause error Expand and the sca

7、le increase dramatically then what should we do?第12页/共79页西安电子科技大学理学院 主讲:王卫卫列主元Gauss column/row/complete pivoting 行主元完全主元第13页/共79页西安电子科技大学理学院 主讲:王卫卫-1/22行+3行2,3行互换列主元列主元例1用Gauss列主元消去法解方程组Solution:消元回代求解：1,2行互换-2/51行+2行-1/51行+3行第14页/共79页西安电子科技大学理学院 主讲:王卫卫方程的根x对角线下方元素对角线上方元素假设已是列主元，否则在对角线下方选列主元，做行置换Gau

8、ss-Jordan消去法无回代过程的列主元消去法：Gauss列主元消去法只消去对角线下方的元素，Gauss-Jordan消去法同时消去对角线下方和上方的元素设已经过k-1步G-J消元，得到第k步消元formulaeG-J消元过程k=1,n第15页/共79页西安电子科技大学理学院 主讲:王卫卫第16页/共79页西安电子科技大学理学院 主讲:王卫卫Gauss-Jordan消去法能够进消去法能够进行下去的条件和前面一样！行下去的条件和前面一样！计算量高于计算量高于Gauss列主元列主元 消消去法，通常用去法，通常用G-J消去法求消去法求矩阵的逆矩阵的逆(初等行变换初等行变换)第17页/共79页西安电

9、子科技大学理学院 主讲:王卫卫列主元Example 2用Gauss-Jordan消去法解方程组Solution:消元列主元1,2行互换-2/51行+2行-1/51行+3行2,3行互换-5/142行+1行-1/22行+3行53行+1行-8.43行+2行各行/对角线元素第18页/共79页西安电子科技大学理学院 主讲:王卫卫example3用Gauss-Jordan消去法求矩阵的逆 solution1,2行互换-1/21行+2行-1/21行+3行2,3行互换-0.42行+1行-0.62行+3行-43行+1行-153行+2行各行/对角线元素第19页/共79页西安电子科技大学理学院 主讲:王卫卫Cram

10、er ruleGauss elimination(LU factorization)Gauss-Jordan elimination(n+1)!n3/3n3/239916800n=10430500Computation cost第20页/共79页西安电子科技大学理学院 主讲:王卫卫detLk=1,Lk invertible,k=1,2,n-1矩阵的直接三角分解(LU分解，不选主元)k=1,2,n-1上三角阵第21页/共79页西安电子科技大学理学院 主讲:王卫卫单位下三角阵单位下三角阵上三角阵第22页/共79页西安电子科技大学理学院 主讲:王卫卫LU分解定理A的各阶顺序主子式单位下三角阵上三角阵

11、!L,U证明要点存在性不必证,见前面的分析,仅证唯一性(采用反正法)L,U,L1U1都可逆单位下三角上三角阵第23页/共79页西安电子科技大学理学院 主讲:王卫卫矩阵的直接三角分解(LU分解)对解线性方程组有什么帮助？三角方程组，易于求解第24页/共79页西安电子科技大学理学院 主讲:王卫卫L矩阵直接三角分解(LU分解)中如何求L和U的元素？矩阵乘法A的第一行U的第一行A的第一列L的第一列U第25页/共79页西安电子科技大学理学院 主讲:王卫卫Lrk=0,k=r+1,n;Lrr=1设已知U的前r-1行，L的前r-1列，求U的第r行，L的第r列 A的第r行U的第r行Ukr=0,k=r+1,n;A

12、的第r列L的第r列LU factorization r=2,n -Doolittle factorization第26页/共79页西安电子科技大学理学院 主讲:王卫卫A24Example 4用紧凑格式解线性方程组solution1 2 3 41112 6 12376 246UL第27页/共79页西安电子科技大学理学院 主讲:王卫卫平方根法和改进的平方根法A对称正定DU0A symmetric and positive definite单位下三角diagonal单位上三角第28页/共79页西安电子科技大学理学院 主讲:王卫卫单位下三角上三角单位下三角上三角LU分解的唯一性Cholesky fac

13、torization第29页/共79页西安电子科技大学理学院 主讲:王卫卫用直接分解法确定 的元素的递推公式Obviously,if j i,then lij=0Else jj,then ljk=0 第30页/共79页西安电子科技大学理学院 主讲:王卫卫Given Ax=b,where A is symmetric and positive definite,what does Cholesky factorization to do with solution of Ax=b三角方程组，易于求解对系数矩阵为对称正定的线性方程组，先对A作Chloesky分解，然后将方程组化成两个三角方程组求解

14、的方法称为解线性方程组的平方根法，计算量是n3/6,是Gauss法的一半，只要存储下半个三角即可Given Ax=b,where A is symmetric and positive definite,what does Cholesky factorization to do with solution of Ax=b第31页/共79页西安电子科技大学理学院 主讲:王卫卫Example 5用平方根法解方程组solutionA 对称正定第32页/共79页西安电子科技大学理学院 主讲:王卫卫diagonal改进的平方根法避免平方根法中的开方运算单位下三角Obviously,if j i,the

15、n lij=0 and ljj=1;else j0 on P0 xRn,For x=0,the equation holds obviously,so we only need to prove the case x0Take M=b,m=a第46页/共79页西安电子科技大学理学院 主讲:王卫卫Proof of property4第47页/共79页西安电子科技大学理学院 主讲:王卫卫Def 2 矩阵范数ARnn,若有非负实值函数N(A)满足1.A0,A0A=0;(正定性)2.cA=cA,cR;(齐次)3.A+BA+B;(三角不等式)4.ABAB 则称N(A)是Rnn上的矩阵范数。example

16、Frobenius norm第48页/共79页西安电子科技大学理学院 主讲:王卫卫Def 4 矩阵范数与向量范数相容xRn,ARnn,若矩阵范数和向量范数满足 Ax Ax -相容性条件则称矩阵范数与向量范数相容。Def 5 矩阵的算子范数xRn,ARnn,给定一种向量范数x p，例如p=1,2,相应的定义一个矩阵的非负函数满足矩阵范数条件和相容性条件，称为A的算子范数。checkWhats the relation between Ax and A,x?Rn ,x xAxA第49页/共79页西安电子科技大学理学院 主讲:王卫卫check第50页/共79页西安电子科技大学理学院 主讲:王卫卫第5

17、1页/共79页西安电子科技大学理学院 主讲:王卫卫第52页/共79页西安电子科技大学理学院 主讲:王卫卫第53页/共79页西安电子科技大学理学院 主讲:王卫卫第54页/共79页西安电子科技大学理学院 主讲:王卫卫常用的矩阵算子范数行范数,p=1列范数,p=22范数,p=3checkxRn,ARnn第55页/共79页西安电子科技大学理学院 主讲:王卫卫check1.1 In case A=0,the equation holds obviously,so we only need to check for the case A0.1.2 y,s.t.第56页/共79页西安电子科技大学理学院 主讲

18、:王卫卫1.1 第57页/共79页西安电子科技大学理学院 主讲:王卫卫1.2 y,s.t.suppose第58页/共79页西安电子科技大学理学院 主讲:王卫卫Similarly,we can prove the result in 21.1 1.2 y,s.t.第59页/共79页西安电子科技大学理学院 主讲:王卫卫1.1 In case A=0,the equation holds obviously,so we only need to check for the case A0.1.2 y,s.t.第60页/共79页西安电子科技大学理学院 主讲:王卫卫Symmetric and posit

19、ive definite,So characteristic values of AtA are all nonnegative real numbers:Let the corresponding characteristic vectors are 1.1 第61页/共79页西安电子科技大学理学院 主讲:王卫卫第62页/共79页西安电子科技大学理学院 主讲:王卫卫1.2 y,s.t.第63页/共79页西安电子科技大学理学院 主讲:王卫卫例第64页/共79页西安电子科技大学理学院 主讲:王卫卫Def 6 矩阵的谱半径谱半径的性质1.矩阵A的谱半径不超过它的任何一种算子范数，包括F-范数Let

20、 be any eigen-alue of A,x the corresponding characteristic eigen-vector,第65页/共79页西安电子科技大学理学院 主讲:王卫卫Let i be any eigen-value of A,ui the corresponding eigen-vectorSo is an eigen-value of AtA,ui the corresponding eigen-vector第66页/共79页西安电子科技大学理学院 主讲:王卫卫Th1Suppose det(IB)=0,then Contradict to 第67页/共79页西

21、安电子科技大学理学院 主讲:王卫卫perturbationperturbationIts funny that such small perturbations in the coefficients lead to so big change in the solution!Thats right!its due to the nature of A and Such a system is said to be ill-posed.Error analysis第68页/共79页西安电子科技大学理学院 主讲:王卫卫由实际问题建立起来的线性方程组Ax=b本身存在模型误差和观测误差，或者是由计算

22、得到的，存在舍入误差等。总之，A,b都会有一定扰动A,b,因此实际处理的是A+A或b+b，我们需要分析A或b的扰动对解的影响。第69页/共79页西安电子科技大学理学院 主讲:王卫卫 x exact solutionx+x exact solution当方程组的系数矩阵A和齐次项b受到扰动A,b后，其解x会受到怎样的扰动？分三种情况讨论1.b有扰动b,A无扰动(A=0)解的相对误差齐次项b的相对误差当b受到扰动b时，引起的解的相对误差不超过b的相对误差的AA-1倍第70页/共79页西安电子科技大学理学院 主讲:王卫卫 2.A有扰动A,b无扰动(b=0)解的相对误差A的相对误差设A-1 A1Th1

23、A+A=A(I+A-1A)可逆A-1A1(不易计算)2.系数矩阵A 的三角约化中出现小主元3.|max(A)|/|min(A)|14.Det(A)很小，或A的某些行近似线性相关5.A的元素间数量级相差很大，且无一定规则。第77页/共79页西安电子科技大学理学院 主讲:王卫卫对病态方程组可采取以下措施：对病态方程组可采取以下措施：1.采用高精度运算，减轻病态影响，例如采用双倍字长运采用高精度运算，减轻病态影响，例如采用双倍字长运算。算。2.用预处理方法改善用预处理方法改善A的条件数，即选择非奇异矩阵的条件数，即选择非奇异矩阵P,Q,s.t.PAQ(Q-1)x=Pb 与与 Ax=b 等价，而等价，而 B=PAQ 的条的条件数比件数比 A 改善，令改善，令y=Q-1x,f=Pb,则求出则求出By=f的解的解y,原方原方程组的解为程组的解为x=Qy.一般选一般选P,Q为对角阵或三角阵。为对角阵或三角阵。3.平衡方法，当平衡方法，当A中元素数量级相差很大时，可采用行或中元素数量级相差很大时，可采用行或列平衡法改善列平衡法改善A的条件数。设的条件数。设A非奇异，计算非奇异，计算 求求Ax=b等价于求等价于求DAx=Db,s.t.DA的条件数得到改善。的条件数得到改善。第78页/共79页西安电子科技大学理学院 主讲:王卫卫感谢您的观看！第79页/共79页