Metodo Gauss-Seidel

La iteración de Gauss-Seidel se define al tomar Q como la parte triangular inferior de A incluyendo los elementos de la diagonal:

Si, como en el caso anterior, definimos la matriz R=A-Q

Una ecuacion de la forma: Qx^(k) = -Rx^(k-1) + b

Cualquier elemento i de un vector Qx^(k) está dado por la ecuacion

Si consideramos la forma de las matrices Q y R, resulta que todos los sumandos para los cuales j>i en el lado izquierdo son nulos, mientras que en la parte derevha todos los sumandos los cuales j<= i. Podemos decir que:

De donde x_i^(k) claramente obtenemos que:

En este metodo los valores de xi reemplazan inmediatamente los valores anteriores mientras que en el metodo de Jacobi todoss los componentes nuevos del vector se calculan antes de llevar a cabo la sustitución. Por el contrario, en el método de Gauss-Seidel, los cálculos deben llevarse a cabo en orden, ya que el nuevo valor xi depende de los valores actualizados de x1, x2, ..., xi-1.

Pseudocodigo

CODIGO

clc

clear

fprintf(' METODO ITERATIVO DE GAUSS SEIDEL\n\n\n')

a=input('Ingrese la matriz de coeficientes:\n ');

b=input('\nIngrese los términos independientes:\n ');

x=input('\nIngrese el vector con las aproximacimaciones Iniciales:\n ');

iter=input('\nIngrese el número máximo de iteraciones:\n ');

tol=input('\nIngrese la tolerancia:\n ');

k=norm(a)*norm(a^(-1));

disp('condicional=')

disp(k)

determinante=det(a)

if determinante==0

disp('El determinante es cero, el problema no tiene solución única')

return

end

n=length(b)

d=diag(diag(a));

fprintf('\nLa matriz l:\n')

l=d-tril(a);

disp(l)

fprintf('\nLa matriz u:\n')

u=d-triu(a);

disp(u)

fprintf('\n SOLUCION:\n')

fprintf('\nLa matriz de transicion de gauss seidel:\n')

T=((d-l)^-1)*u;

disp(T)

re=max(abs(eig(T)))

if re>1

disp('Radio Espectral mayor que 1')

disp('el método no converge')

return

end

fprintf('\nEl vector constante es::\n')

C=((d-l)^-1)*b;

disp(C)

i=0;

err=tol+1;

z=[i,x(1),x(2),x(3),x(4),err];

while err>tol & i<iter

xi=T*x+C;

%disp(xi)

err=norm(x-xi);

%err=max(abs(xi-x));

%err=norm(xi-x)/norm(xi);

x=xi;

i=i+1;

z(i,1)=i;

z(i,2)=x(1);

z(i,3)=x(2);

z(i,4)=x(3);

z(i,5)=x(4);

z(i,6)=err;

end

fprintf('\nTABLA:\n\n n x1 x2 x3 x4 Error\n\n ');

disp(z).