Seidel method algorithm in Matlab

Seidel method algorithm in Matlab :

function [Sol,E,iter,csuf]=seidel(M,b,x0,error,n_iter) % Syntaxis: [Sol,E,iter,csuf]=seidel(M,b,x0,error,n_iter) % ---------------------------------------------------------------------- % Esta funcion resuelve un sistema de ecuaciones lineales (SEL) de la % forma Ax=b por el metodo de Seidel segun % % x0 in V (arbitrario) % x(k+1)=( I-inv(D-E)*A )*x(k) + inv(D-E)*b % % Entrada: % M --> Matriz del sitema % b --> Vector independiente (debe ser un vector columna) % x0 --> Vector de aproximaciones iniciales % error --> Tolerancia del calculo {default | 0.00001} % n_iter --> Numero maximo de iteraciones {default | 1000} % % Salida: % Sol --> Vector solucion del SEL % E --> Error de la aproximacion % csuf --> Cumplimiento de la condición suficiente de convergencia "alpha<1" % % % Ejemplo: % % M=[9 -1 2; 1 8 2; 1 -1 11]; b=[9; 19; 10]; % [Sol,E,iter,csuf]=seidel(M,b) % % Sol = % 1.0000 % 2.0000 % 1.0000 % % E = 1.2249e-06 % % iter = 6 % % csuf = 1 % % % Ver tambien: jacobi % % ---------------------------------------------------------------------- % % Elaborado por: % % Msc. Alexeis Companioni Guerra % Lic. Vianka Orovio Cobo % % Revision: 1.0 % Fecha: 30/09/2007 %% Validacion de las entradas if nargin<2 error('Verifique los datos de entrada.'); else if det(M) == 0 error('M debe ser no singular.'); end if diff(size(M))~=0 error('La matriz de entrada debe ser cuadrada.'); end if min(size(b))>1 error('b debe ser un vector.'); end if max(size(M)) ~= max(size(b)) error('Verifique la dimension del vector independiente respecto a la matriz del sistema.'); end end b = b(:); n = length(b); if nargin<3 x0 = zeros(n,1); else if min(size(x0))>1 error('x0 debe ser un vector.'); end if max(size(M)) ~= max(size(x0)) error('Verifique la dimension del vector x0 respecto a la matriz del sistema.'); end end if nargin<4 error = 0.00001; else if sum(size(error))~=2 error('La tolerancia debe ser un escalar.'); end end if nargin<5 n_iter = 1000; else if sum(size(n_iter))~=2 error('El numero de iteraciones debe ser un escalar.'); end if round(n_iter)-n_iter~=0 error('El numero de iteraciones debe ser un entero.'); end end %% Determinacion del factor de convergencia "beta" Jac = eye(n)-( diag( diag(M).^(-1) )*M ); Jac_up = triu(Jac,1); Jac_lo = tril(Jac,-1); beta = max( sum(abs(Jac_up),2) ./ ( ones(n,1)-sum(abs(Jac_lo),2) ) ); if beta>=1 csuf = 0; else csuf = 1; end %% Ciclo iterativo iter = 0; Sol = zeros(n,1); E = inf; s_inv = inv(diag(diag(M)) - tril(M,-1).*(-1)); S = eye(n)-s_inv*M; % Matriz de Seidel while (E>=error) && (iter<n_iter) Sol = S*x0 + s_inv*b; E = abs(beta/(1-beta))*norm(Sol-x0,inf); x0 = Sol; iter = iter + 1; end