function x=gseidel(A,b,x,M,tau); % GSEIDEL: Solves Ax=b using Gauss-Seidel iteration. % USAGE: % x=gseidel(A,b,x,M,tau); % % INPUTS: % A : nxn matrix % b : nxm matrix % [x] : nxm matrix of starting values. Default: b % [M] : maximum number of iterations % [tau] : convergence tolerence % % OUTPUT: % x : nxm matrix of approximate solutions to Ax=b % % Written by: % -- % John L. Weatherwax 2006-12-19 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- if nargin<3 | isempty(x), x=b; end if nargin<4 | isempty(M), M=1000; end if nargin<5 | isempty(tau), tau=sqrt(eps); end % obtain the factors S and T in the splitting A = S - T: S = tril(A); T = S - A; divtol=100*max(abs(b-A*x)); for i=1:M x_new = S \ (T*x + b); % the new value r=b-A*x_new; % the new residual dx=x_new - x; % convergence check if( all(abs(r)divtol ) disp('WARNING: Iterations diverging'); return end x = x_new; end if i==M, disp('WARNING: Maximum iterations exceeded'); end