高斯消去法matlab解线性代数,线性代数方程组的数值解法 (1) Gauss 消去法

(1)Gauss

(Demos in Matlab: airfoil in 2D)

Gauss SuperLUJacobi, GS, SOR, SSOR Krylov CG, MINRES , GMRES, QMR, BiCGStabThe Landscape of Ax=b Solvers

Pivoting LU

GMRES,QMR,

Cholesky

Conjugategradient

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 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 of Mathematics, http://www-history.mcs.st-andrews.ac.uk/history/index.html)

Gauss1777-1855

-------

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-1Gauss

Gauss (2)

Gauss A x=b

--- LU A = L U (cost = 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 Algorithmfor k = 1 to n-1 k, (k) for i = k+1 to n ki for j = k to n 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) = A(i,j) - m * A(k,j): k0,

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

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 loopGauss Elimination Algorithm (4)

: 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 , k+1:n) BLAS 2 (rank-1 update)

What we havent told you (A(k)(k,k)0) Sparse LU, Band LU (F.Gustavson & S.Toledo, Recursive Algorithm) ,, : A(k)(k,k)0 :

Matlab Linpack invlu\

C, Fortran, Matlab sgea.fsgefa.f

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 several such columns).% Gaussian elimination with partial pivoting[n, 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+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 Soft/c2-1.c Gauss Elimination */#include #include #include #define TRUE 1/* a[i][j] : 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_init[10][11] = {{1, 2, 3, 6}, {2, 2, 3, 7}, {3, 3, 3, 9}};static double a[10][11];void gauss();/*static int _aini = 1; */

printf( "\nComputer Soft/C2-1 Gauss Elimination \n\n" ); printf( "Augmented matrix\n" ); for( i = 1; i

void gauss(n, a)int n; double a[][11];{int i, j, jc, 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; i

CC PAGE 220-223: NUMERICAL MATHEMATICS AND COMPUTING, CHENEY/KINCAID, 1985CC FILE: GAUSS.FORCC 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.,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.,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(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(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 CONTINUE S(I) = SMAX 3 CONTINUE DO 7 K = 1,N-1 RMAX = 0.0 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),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