Home > K25104 > LinearSystems > Gauss_Seidel.m

Gauss_Seidel

PURPOSE ^

Solves the linear system Ax=b iteratively using the Gauss-Seidel method.

SYNOPSIS ^

function [ x,k,H,E ] = Gauss_Seidel( A,b,x0,w,imax )

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005