线性方程组的数值方法

1、线性代数方程组的数值解法(1) Gauss 消去法,(Demos in Matlab: airfoil in 2D),线性代数方程组的数值解法 直接法：Gauss 消去法，SuperLU 迭代法：定常迭代（Jacobi, GS, SOR, SSOR） Krylov 子空间方法（CG, MINRES , GMRES, QMR, BiCGStab）,The Landscape of Ax=b Solvers,刘徽 （约220-280）,Gaussian elimination, which first appeared in the text Nine Chapters on the Mathem

2、atical Art written in 200 BC, was used by Gauss in his work which studied 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 pr

3、ecisely Gaussian elimination on the coefficient matrix. (The MacTutor History of Mathematics, http:/www-history.mcs.st-andrews.ac.uk/history/index.html),Gauss（1777-1855）,今有上禾三秉，中禾二秉，下禾一秉，实三十九斗；上禾二秉，中禾三秉，下禾一秉，实三十四斗；上禾一秉，中禾二秉，下禾三秉，实二十六斗。问上、中、下禾实一秉各几何？答曰：上禾一秉九斗四分斗之一。中禾一秉四斗四分斗之一。下禾一秉二斗四分斗之三。-九章算术,一个两千年前

4、的例子,Basic idea: Add multiples of each row to later rows to make A upper triangular,一个两千年前的例子(2),Solving linear equations is not trivial. Forsythe (1952),After k=1 After k=2 After k=3 After k=n-1,Gauss 消去过程图示,用矩阵变换表达消去过程,利用Gauss 变换矩阵的性质：,用矩阵变换表达消去过程(2),单位下三角形,用Gauss消去法求解 A x=b - LU分解 A = L U (cost =

5、2/3 n3 flops) - 求解 L y = b (cost = n2 flops) - 求解 U x = y (cost = n2 flops),版本一,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 Algorithm,for k = 1 to n-1 对第k列,消去对角线以下元素 (通过每行加上第k行的倍数) for i = k+1 to n 对第k行以下的每一行i for j = k to n

6、 第k行的倍数加到第 i 行 A(i,j) = A(i,j) - (A(i,k)/A(k,k) * A(k,j),版本二: 在内循环中去掉常量 A(i,k)/A(k,k) 的计算,上一版本,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)

7、= A(i,j) - m * A(k,j),版本三: 第k列对角线以下为0,无需计算,上一版本,Gauss Elimination Algorithm (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),版本四: 将乘子 m 存

8、储在对角线以下备用,上一版本,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),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),版本五: Split loop,Gauss Elimination Algorithm (4),上一版本

9、 版本六: 用矩阵运算,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 (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 ,

10、 k+1:n) BLAS 2 (rank-1 update),What we havent told you,选主元策略(当主元A(k)(k,k)为0或很小时) 向后误差分析 并行技术 块算法 Sparse LU, Band LU 最新进展(F.Gustavson for k = 1 : n-1 % find index of largest element below diagonal in column k max = k; for i = k+1 : n if abs(A(i, k) abs(A(max, k) max = i; end end % swap with row k A(k

11、 max, :) = A(max k, :); b(k max) = b(max k); % zero out entries of A and b using pivot A(k, k) A(k+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(

12、i, :) - alpha * A(k, :); %end end % back substitution x = zeros(size(b); for i = n : -1 : 1 j = i+1 : n; x(i) = (b(i) - A(i, j) * x(j) / A(i, i); end,/* Computer Soft/c2-1.c Gauss Elimination */ #include #include #include #define TRUE 1 /* aij : matrix element, a(i,j) n : order of matrix eps : machi

13、ne 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( nComputer Soft/C2-1 Gauss Elimination nn ); printf( Augmented matrixn ); for( i = 1; i = n;

14、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 = n; i+ ) printf( %5d %16.6en, i, ain+1 ); printf( -nn ); exit(0); ,void gauss(n, a) int n; double a11; int i, j, jc,

15、jr, k, kc, nv, pv; double det, eps, ep1, eps2, r, temp, tm, va; eps = 1.0; ep1 = 1.0 ; /* eps = Machine epsilon */ while( ep1 0 ) eps = eps/2.0; ep1 = eps*0.98 + 1; ep1 = ep1 - 1; eps = eps*2; eps2 = eps*2; printf( Machine epsilon=%g n, eps ); det = 1; /* Initialization of determinant */ for( i = 1;

16、 i = 1; nv- ) va = anvn+1; for( k = nv + 1; k = n; k+ ) va = va - anvk*akn+1; anvn+1 = va/anvnv; printf( Determinant = %g n, det ); return; ,C C PAGE 220-223: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985 C C FILE: GAUSS.FOR C C GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING (GAUSS,SO

17、LVE,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.,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 (B

18、2(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.,12./ 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

19、(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) RETURN 10 FORMAT(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 SMAX = 0.0 DO 2 J = 1,N SMAX = AMAX1(SMAX,ABS(A(I,J) 2 CON