Gauss-Seidel method

The Gauss-Seidel iteration is defined by taking Q as the lower triangular part of A including the elements of the diagonal:

If, as in the previous case, we define the matrix R = A-Q

And the equation can be written in the form:

Qx^(k) = -Rx^(k-1) + b

Any element, i, of the vector Qx^(k) will be given by the equation:

If we consider the peculiar form of the matrices Q and R, it turns out that all the summands for which j> i on the left side are null, whereas in the right part all the summands for which . We can then write:

From where x_i^(k) clearing, we obtain

Note that in the Gauss-Seidel method the updated values ​​of xi immediately replace the previous values, whereas in the Jacobi method all new components of the vector are calculated before carrying out the substitution. In contrast, in the Gauss-Seidel method the calculations must be carried out in order, since the new value xi depends on the updated values ​​of x1, x2, ..., xi-1.

Pseudocode

CODE

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).