Solves the linear system Ax=b iteratively using the Gauss-Seidel method. Input arguments: A, square matrix b, column vector x0, initial estimate of the solution w, relaxation parameter ( w=1 gives unrelaxed iterations ) imax, maximum number of iterations Output arguments: x, solution k, number of iterations H, iteration matrix E, eigenvalues of the iteration matrix
0001 function [ x,k,H,E ] = Gauss_Seidel( A,b,x0,w,imax ) 0002 % Solves the linear system Ax=b iteratively using the Gauss-Seidel method. 0003 % Input arguments: 0004 % A, square matrix 0005 % b, column vector 0006 % x0, initial estimate of the solution 0007 % w, relaxation parameter ( w=1 gives unrelaxed iterations ) 0008 % imax, maximum number of iterations 0009 % Output arguments: 0010 % x, solution 0011 % k, number of iterations 0012 % H, iteration matrix 0013 % E, eigenvalues of the iteration matrix 0014 0015 [m,n]=size(A); % find the size of A 0016 [p,q]=size(b); % find the size of b 0017 if n~= m 0018 error('input is not a square matrix'); 0019 elseif p~=m || q~=1; 0020 error('input dimensions do not agree'); 0021 end 0022 0023 B=A; % initialize B to be A 0024 for i=1:n; 0025 for j=1:i; 0026 B(i,j)=0; % set B to be the strictly upper triangular part of A 0027 end 0028 end 0029 0030 J=A-B; % initialize J to be the lower triangular part of A. 0031 0032 y=-B*x0+b; % calculate matrix vector product on right hand side 0033 temp=Forward(J,y); % solve the lower triangular system by forward 0034 % substitution. 0035 x1=w*temp+(1-w)*x0; % implement relaxation 0036 k=1; % Set iteration counter to 1 0037 0038 while (abs(x1-x0)>1e-12) & (k <= imax); % convergence if the approximations 0039 % are close or the maximum number of 0040 % iterations is reached 0041 x0=x1; 0042 y=-B*x0+b; % calculate matrix vector product on right hand side 0043 temp=Forward(J,y); % solve the lower triangular system by forward 0044 % substitution. 0045 x1=w*temp+(1-w)*x0; % implement relaxation 0046 k=k+1; % increment counter 0047 end 0048 x=x1; 0049 0050 [H,E]=iteration_analysis(J,B,w); % this function calculates the iteration 0051 % matrix and its eigenvalues for 0052 % illustrative purposes. 0053 end 0054