线性方程组的数值方法

《线性方程组的数值方法》由会员分享，可在线阅读，更多相关《线性方程组的数值方法（28页珍藏版）》请在装配图网上搜索。

1、线性代数方程组的数值解法线性代数方程组的数值解法(1)Gauss 消去法消去法bxAnnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111(Demos in Matlab:airfoil in 2D)线性代数方程组的数值解法线性代数方程组的数值解法直接法直接法：Gauss 消去法，SuperLU迭代法迭代法：定常迭代（Jacobi,GS,SOR,SSOR）Krylov 子空间方法（CG,MINRES,GMRES,QMR,BiCGStab）The Landscape of Ax=b Solvers Pivoting LU GMRES,QMR,Cho

2、lesky Conjugate gradient DirectA=LUIterativey=AyNon-symmetricSymmetricpositivedefiniteMore RobustLess StorageMore RobustMore General刘徽刘徽（约220-280）Gaussian elimination,which first appeared in the text Nine Chapters on the Mathematical Art written in 200 BC,was used by Gauss in his work which studied

3、the orbit of the asteroid Pallas.Using observations of Pallas taken between 1803 and 1809,Gauss obtained a system of six linear equations in six unknowns.Gauss gave a systematic method for solving such equations which is precisely Gaussian elimination on the coefficient matrix.(The MacTutor History

4、of Mathematics,http:/www-)Gauss（1777-1855）263234323923zyxzyxzyx今有上禾三秉，中禾二秉，下禾今有上禾三秉，中禾二秉，下禾一秉，实三十九斗；上禾二秉，一秉，实三十九斗；上禾二秉，中禾三秉，下禾一秉，实三十四中禾三秉，下禾一秉，实三十四斗；上禾一秉，中禾二秉，下禾斗；上禾一秉，中禾二秉，下禾三秉，实二十六斗。问上、中、三秉，实二十六斗。问上、中、下禾实一秉各几何？答曰：上禾下禾实一秉各几何？答曰：上禾一秉九斗四分斗之一。中禾一秉一秉九斗四分斗之一。中禾一秉四斗四分斗之一。下禾一秉二斗四斗四分斗之一。下禾一秉二斗四分斗之三。四分斗之三。-

5、九章算术九章算术一个两千年前的例子263213413239123130803912338343135437,417,411xyz263234323923zyxzyxzyx138392338343135zyzyzyx533512313583923zzyzyx5335123135008039123vBasic idea:Add multiples of each row to later rows to make A upper triangular一个两千年前的例子(2)nnnnnnaaaaaaaaa212222111211bxAnnbbbxxx2121Solving linear equati

6、ons is not trivial.Forsythe(1952)AaaaaaaaaaAnnnnnn)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1()2()2(2)2(2)2(22)1(1)1(12)1(11)2(00nnnnnaaaaaaaA)3()3(3)3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11)3(00000nnnnnnaaaaaaaaaaaA),3,2,(),3,2()1(11)1()2()1(11)1(11njialaaniaaljiijijii),4,3,(),4,3()2(22)2()3()2(22)2

7、(22njialaaniaaljiijijii),1,(),1()()1()1()()(nkjialaankiaalkkjikkijkijkkkkikik)()(1,)()(,1)(1,1)(,1)()(1,)()1(,1)1(1,1)1(,1)1(1,1)1(1)1(1,1)1(1)1(1,1)1(11)(000000knnkknknkknkkkkkkkkknkkkkkkknkkkkkkkkkknkkkkaaaaaaaaaaaaaaaaaaA)1()1(1,)1(,1)1(1,1)()(1,)()1(,1)1(1,1)1(,1)1(1,1)1(1)1(1,1)1(1)1(1,1)1(11)1

8、(00000000knnkknknkkkkkknkkkkkkknkkkkkkkkkknkkkkaaaaaaaaaaaaaaaaA)()3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11nnnnnnaaaaaaaaaa)()3(3)2(2)1(1321nnnbbbbxxxx After k=1 After k=2 After k=3 After k=n-1Gauss 消去过程图示)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1(nnnnnnaaaaaaaaaA)2()2(2)2(2)2(22)1(1)1(12)1(11)2(00n

9、nnnnaaaaaaaA)3()3(3)3(3)3(33)2(2)2(23)2(22)1(1)1(13)1(12)1(11)3(00000nnnnnnaaaaaaaaaaaA111)1(11def)1(AcraAAT222)2(221)1(11def11)1(11)1(1)2(00)1(1111AcraraAraALATTarcTT22)2(221)1(11)2(2)3(00)2(2222AraraALAarcTTTTnacelIIL22def22)2(22211Tnace lIIL11def11)1(1111用矩阵变换表达消去过程UaaaaaaALLLAnnnnnnndef)()2(2)2(

10、22)1(1)1(12)1(11)1(121)(LUULLLAndef111211利用Gauss 变换矩阵的性质：TnnTTnelelelILLLL112211111211)()2()1(11jielelIelIelILLelIelILTjjTiiTjjTiijiTiiTiii用矩阵变换表达消去过程(2)单位下三角形1111)3(33)3(3)2(22)2(2)1(11)1(1)2(22)2(32)1(11)1(31)1(11)1(21aaaaaaaaaaaannn用Gauss消去法求解 A x=b -LU分解 A=L U (cost=2/3 n3 flops)-求解 L y=b (cost=

11、n2 flops)-求解 U x=y (cost=n2 flops)v版本一for k=1 to n-1 for i=k+1 to n m=A(i,k)/A(k,k)for j=k to n A(i,j)=A(i,j)-m*A(k,j)算法实现：Gauss Elimination Algorithmfor k=1 to n-1 对第对第k列列,消去对角线以下元素消去对角线以下元素 (通过每行加上第通过每行加上第k行的倍数行的倍数)for i=k+1 to n 对第对第k行以下的每行以下的每一一行行i for j=k to n 第第k行的倍数加到第行的倍数加到第 i 行行 A(i,j)=A(i,

12、j)-(A(i,k)/A(k,k)*A(k,j)v版本二:在内循环中去掉常量 A(i,k)/A(k,k)的计算v上一版本Gauss Elimination Algorithm(2)for k=1 to n-1 for i=k+1 to n m=A(i,k)/A(k,k)for j=k to n A(i,j)=A(i,j)-m*A(k,j)for k=1 to n-1 for i=k+1 to n m=A(i,k)/A(k,k)for j=k+1 to n A(i,j)=A(i,j)-m*A(k,j)v版本三:第k列对角线以下为0,无需计算v上一版本Gauss Elimination Algor

13、ithm(3)for k=1 to n-1 for i=k+1 to n m=A(i,k)/A(k,k)for j=k+1 to n A(i,j)=A(i,j)-m*A(k,j)for k=1 to n-1 for i=k+1 to n A(i,k)=A(i,k)/A(k,k)for j=k+1 to n A(i,j)=A(i,j)-A(i,k)*A(k,j)v版本四:将乘子 m 存储在对角线以下备用v上一版本for k=1 to n-1 for i=k+1 to n A(i,k)=A(i,k)/A(k,k)for j=k+1 to n A(i,j)=A(i,j)-A(i,k)*A(k,j)f

14、or k=1 to n-1 for i=k+1 to n A(i,k)=A(i,k)/A(k,k)for i=k+1 to n for j=k+1 to n A(i,j)=A(i,j)-A(i,k)*A(k,j)v版本五:Split loopGauss Elimination Algorithm(4)v上一版本v版本六:用矩阵运算for k=1 to n-1 for i=k+1 to n A(i,k)=A(i,k)/A(k,k)for i=k+1 to n for j=k+1 to n A(i,j)=A(i,j)-A(i,k)*A(k,j)Gauss Elimination Algorithm

15、(5)for k=1 to n-1 A(k+1:n,k)=A(k+1:n,k)/A(k,k)BLAS 1(scale a vector)A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n)BLAS 2(rank-1 update)What we havent told youl 选主元策略选主元策略(当主元A(k)(k,k)为0或很小时)l 向后误差分析向后误差分析l 并行技术并行技术l 块算法块算法 l Sparse LU,Band LUSparse LU,Band LUl 最新进展最新进展(&,Recursive Algorithm)l 还

16、可用于矩阵求逆还可用于矩阵求逆,求行列式求行列式,秩秩l 定理定理:主元A(k)(k,k)不为0的充要条件是顺序主子矩阵非奇异l 定理定理:分解的存在性和唯一性 Matlab 中的相应函数Linpack 中对应的函数 invluC,Fortran,Matlab 代码function x=lsolve(A,b)%x=lsolve(A,b)returns the solution to the equation Ax=b,%where A is an n-by-b matrix and b is a column vector of%length n(or a matrix with severa

17、l such columns).%Gaussian elimination with partial pivotingn,n=size(A);for k=1:n-1%find index of largest element below diagonal in column kmax=k;for i=k+1:nif abs(A(i,k)abs(A(max,k)max=i;endend%swap with row kA(k max,:)=A(max k,:);b(k max)=b(max k);%zero out entries of A and b using pivot A(k,k)A(k+

18、1:n,k)=A(k+1:n,k)/A(k,k);b(k+1:n)=b(k+1:n)-A(k+1:n,k)*b(k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);%for i=k+1:n%alpha=A(i,k)/A(k,k);%b(i)=b(i)-alpha*b(k);%A(i,:)=A(i,:)-alpha*A(k,:);%end end%back substitutionx=zeros(size(b);for i=n:-1:1j=i+1:n;x(i)=(b(i)-A(i,j)*x(j)/A(i,i);end/*Computer

19、Soft/c2-1.c Gauss Elimination*/#include#include#include#define TRUE 1/*aij:matrix element,a(i,j)n:order of matrix eps:machine epsilon det:determinant */void main()int i,j,_i,_r;static n=3;static float a_init1011=1,2,3,6,2,2,3,7,3,3,3,9;static double a1011;void gauss();/*static int _aini=1;*/printf(n

20、Computer Soft/C2-1 Gauss Elimination nn);printf(Augmented matrixn);for(i=1;i=n;i+)for(j=1;j=n+1;j+)aij=a_initi-1j-1;printf(%13.5e,aij);printf(n);gauss(n,a);printf(Solutionn);printf(-n);printf(i x(i)n);printf(-n);for(i=1;i 0)eps=eps/2.0;ep1=eps*0.98+1;ep1=ep1-1;eps=eps*2;eps2=eps*2;printf(Machine eps

21、ilon=%g n,eps);det=1;/*Initialization of determinant*/for(i=1;i=(n-1);i+)pv=i;for(j=i+1;j=n;j+)if(fabs(apvi)fabs(aji)pv=j;if(pv!=i)for(jc=1;jc=(n+1);jc+)tm=aijc;aijc=apvjc;apvjc=tm;det=-det;if(aii=0)/*Singular matrix*/printf(Matrix is singular.n);exit(0);for(jr=i+1;jr=n;jr+)/*Elimination of below-di

22、agonal.*/if(ajri!=0)r=ajri/aii;for(kc=i+1;kc=(n+1);kc+)temp=ajrkc;ajrkc=ajrkc-r*aikc;if(fabs(ajrkc)eps2*temp)ajrkc=0.0;/*If the result of subtraction is smaller than*2 times machine epsilon times the original*value,it is set to zero.*/for(i=1;i=1;nv-)va=anvn+1;for(k=nv+1;k=n;k+)va=va-anvk*akn+1;anvn

23、+1=va/anvnv;printf(Determinant=%g n,det);return;CC PAGE 220-223:NUMERICAL MATHEMATICS AND COMPUTING,CHENEY/KINCAID,1985CCC GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING(GAUSS,SOLVE,TSTGAUS)C DIMENSION A1(4,4),A2(4,4),A3(4,4),B1(4),B2(4),B3(4)DIMENSION L(4),S(4),X(4)DATA(A1(I,J),I=1,4),J=1,4)/3.,

24、1.,6.,0.,4.,5.,3.,0.,3.,-1.,7.,A 0.,0.,0.,0.,0./DATA(B1(I),I=1,4)/16.,-12.,102.,0./DATA(A2(I,J),I=1,4),J=1,4)/3.,2.,1.,0.,2.,-3.,4.,0.,-5.,1.,-1.,A 0.,0.,0.,0.,0./DATA(B2(I),I=1,4)/4.,8.,3.,0./DATA(A3(I,J),I=1,4),J=1,4)/1.,3.,5.,4.,-1.,2.,8.,2.,2.,1.,6.,A 5.,1.,4.,3.,3./DATA(B3(I),I=1,4)/5.,8.,10.,1

25、2./C CALL TSTGAUS(3,A1,4,L,S,B1,X)CALL TSTGAUS(3,A2,4,L,S,B2,X)CALL TSTGAUS(4,A3,4,L,S,B3,X)END SUBROUTINE TSTGAUS(N,A,IA,L,S,B,X)DIMENSION A(IA,N),B(N),X(N),S(N),L(N)PRINT 10,(A(I,J),J=1,N),I=1,N)PRINT 10,(B(I),I=1,N)CALL GAUSS(N,A,IA,L,S)CALL SOLVE(N,A,IA,L,B,X)PRINT 10,(X(I),I=1,N)RETURN10 FORMAT

26、(5X,3(F10.5,2X)END SUBROUTINE GAUSS(N,A,IA,L,S)DIMENSION A(IA,N),L(N),S(N)DO 3 I=1,N L(I)=I DO 2 J=1,N SMAX=AMAX1(SMAX,ABS(A(I,J)2 CONTINUE S(I)=SMAX 3 CONTINUE DO 7 K=1,N-1 DO 4 I=K,N R=ABS(A(L(I),K)/S(L(I)IF(R.LE.RMAX)GO TO 4 J=I RMAX=R 4 CONTINUE LK=L(J)L(J)=L(K)L(K)=LK DO 6 I=K+1,N XMULT=A(L(I),

27、K)/A(LK,K)DO 5 J=K+1,N A(L(I),J)=A(L(I),J)-XMULT*A(LK,J)5 CONTINUE A(L(I),K)=XMULT 6 CONTINUE 7 CONTINUE RETURN END SUBROUTINE SOLVE(N,A,IA,L,B,X)DIMENSION A(IA,N),L(N),B(N),X(N)DO 3 K=1,N-1 DO 2 I=K+1,N B(L(I)=B(L(I)-A(L(I),K)*B(L(K)2 CONTINUE 3 CONTINUE X(N)=B(L(N)/A(L(N),N)DO 5 I=N-1,1,-1 SUM=B(L(I)DO 4 J=I+1,N SUM=SUM-A(L(I),J)*X(J)4 CONTINUE X(I)=SUM/A(L(I),I)5 CONTINUE RETURN END