Método de Gauss-Seidel

Los métodos iterativos o aproximados proveen una alternativa aceptable a los métodos de pivoteo como el método de Gauss y Gauss-Jordan. En el caso que les expongo en esta ocasión, el método de Gauss-Seidel emplea un enfoque similar al utilizado por los métodos para obtener las raíces de ecuaciones simples lineales o no lineales.

Dichos planteamientos se fundamentan en suponer un valor inicial de la raíz y luego utilizar un método sistemático para obtener una aproximación más precisa. En este caso extenderemos este concepto a sistemas de ecuaciones lineales.

El método de Gauss-Seidel, es el método iterativo más comúnmente utilizado para obtener la solución de sistemas de ecuaciones  lineales, Supongamos que se nos da en conjunto de $n$ ecuaciones con $n$ incógnitas:

$$[A]\{x\}=\{b\}$$
en donde:

\begin{eqnarray}
A&=&\left [
\begin{array}{ccccccc}
a_{11}    &a_{12}   &a_{13} &+&\cdots&+&a_{1n}\\
a_{21}   &a_{22}  &a_{23} &+&\cdots&+&a_{1n}\\
\vdots      & \vdots   &  \vdots & \vdots  & \vdots  & \vdots \\
a_{n1}   &a_{n2}  &a_{n3}&+&\cdots&+&a_{nn}
\end{array}
\right ]
\end{eqnarray}

\begin{eqnarray}
b&=&\left \{
\begin{array}{c}
b_1\\b_2\\\vdots\\b_n
\end{array}
\right \}
\end{eqnarray}



Si suponemos que los elementos de la diagonal principal son distintos de cero y  $n=3$, el sistemas se pude resolver de la siguiente forma:

\begin{eqnarray}
x_{1}&=&\frac{b_1-a_{12}x_2-a_{13}x_3}{a_{11}} \label{ec2}\\
x_{2}&=&\frac{b_2-a_{21}x_1-a_{23}x_3}{a_{22}} \label{ec1}\\
x_{3}&=&\frac{b_3-a_{31}x_1-a_{32}x_2}{a_{33}}\label{ec3}
\end{eqnarray}

Ahora es posible obtener la solución de los sistemas de ecuaciones presentados en \eqref{ec2}, \eqref{ec1} y \eqref{ec3} a partir de un valor inicial $x$, el cual representa un vector con los valores iniciales para $x_1$, $x_2$ y $x_3$. Lo que se acostumbra a hacer para fijar estos valores iniciales es asumirlos todos iguales a $0$, de esta forma se remplaza $x_2=0$ y $x_3=0$ en la ecuación \eqref{ec2}, obteniéndose de esta manera el primer valor para $x_1$. Luego se reemplaza este nuevo valor  de $x_1$ y el valor  $x_3=0$ en la ecuación \eqref{ec2} obteniéndose así el nuevo valor de $x_2$. Y por último los nuevos valores de $x_1$ y $x_2$ se reemplazan en la ecuación \eqref{ec3} y se obtiene el nuevo valor estimado de $x_3$. El método iterativo se repite hasta alcanzar una aproximación lo suficientemente precisa.

Como criterio para verificar que el método converge a una solución se utiliza el error relativo, esto es:
\begin{equation}
|\epsilon_{r,i}|=\left | \frac{x_i^j-x_i^{j-1}}{x_i^j}\right |
\end{equation}
Para todas las $i$, en donde $j$ y $j-1$ son las iteraciones actuales y previas respectivamente.

Ahora volvamos al caso general en el que se cuenten con $n$ ecuaciones y $n$ incógnitas. ¿Cómo podemos determinar una expresión general para evaluar los valores de $x$, para de esta forma realizar un algoritmo que luego nos permita generar un programa?. Veamos.

Podemos escribir un conjunto de expresiones como las siguientes:
\begin{eqnarray}
x_i&=&\frac{b_i-\sum_{j=i+1}^n a_{ij}x_j}{a_{ii}}& i=1 \label{ec4}\\
x_i&=&\frac{b_i-\sum_{j=i+1}^n a_{ij}x_j-\sum_{j=1}^{i-1} a_{ij}x_j}{a_{ii}}& i=2,3,..,n-1 \label{ec5}\\
x_i&=&\frac{b_i-\sum_{j=1}^{n-1} a_{ij}x_j}{a_{ii}}& i=n \label{ec6}
\end{eqnarray}

De esta forma se nos es fácil implementar el método en un lenguaje de programación. verifiquemos  que sucede con las ecuaciones \eqref{ec4}, \eqref{ec5} y \eqref{ec6} para $n=3$.

De la ecuación \eqref{ec4} tenemos:
\begin{eqnarray}
x_1&=&\frac{b_1-\sum_{j=2}^3 a_{1j}x_j}{a_{11}} \nonumber\\
x_1&=&\frac{b_1-a_{12}x_{2}+a_{13}x_3}{a_{11}} \nonumber
\end{eqnarray}
De la ecuación \eqref{ec5}:
\begin{eqnarray}
x_2&=&\frac{b_2-\sum_{j=3}^3 a_{2j}x_j-\sum_{j=1}^1 a_{2j}x_j}{a_{22}}\nonumber\\
x_2&=&\frac{b_2-a_{23}x_3-a_{21}x_1}{a_{22}}\nonumber
\end{eqnarray}
y de la ecuación \eqref{ec6}:
\begin{eqnarray}
x_3&=&\frac{b_3-\sum_{j=1}^2 a_{3j}x_j}{a_{33}}\nonumber\\
x_3&=&\frac{b_3-a_{31}x_1-a_{32}x_2}{a_{33}}
\end{eqnarray}
Adicionalmente recordemos que  los elementos de la diagonal principal deben ser distintos de cero para evitar una división por cero. De modo que se debe tener esto en cuenta a la hora de especificar la matriz de coeficientes $A$.

Cabe mencionar que al ser un método numérico, pudiera darse el caso que el mismo no convergiera. De tal forma que para garantizar una convergencia la matriz  $A$ debe cumplir lo siguiente:

Lo cual quiere decir que el elemento perteneciente a la diagonal principal en cada ecuación debe ser mayor que la suma del resto de elementos de la misma ecuación. De esta forma la convergencia se garantiza si se cumple la condición, aunque pueden existir algunos casos en que el método converja sin cumplir con ella. Los sistemas que cumplen con la condición antes mencionada son conocidos como diagonalmente dominantes.

Iteración de Jacobi
Como en el caso del método de Gauss-Seidel cada valor de $x_i$ que se calcula es utilizado  inmediatamente para calcular el valor siguiente, de tal forma que si la solución es convergente se utilizará la mejor combinación posible de los valores de $x_i$. Pero también existe otra forma de hacer esto, y es primero reemplazar cada uno de  los valores de $x_i$ en las ecuaciones \eqref{ec4}, \eqref{ec5} y \eqref{ec6}, obteniendo otro nuevo conjunto de valores $x_i$, este es el planteamiento conocido como Iteración de Jacobi

Veamos un ejemplo utilizando el método de Gauss-Seidel. Obtener la solución del siguiente sistema de ecuaciones:
\begin{eqnarray}
x_1&+&10x_2&+&4x_3&-&2x_4&=&2 \nonumber\\
8x_1&&&-&10x_3&+&2x_4&=&-7 \nonumber\\
8x_1&+&3x_2&+&x_3&+&17x_4 &=&8\nonumber\\
11x_1&+&7x_2&-&3x_3&+&2x_4&=&12 \nonumber
\end{eqnarray}

Si creamos una función llamada en este caso gseidel.m en matlab, tendríamos lo siguiente:


>> A=[1 10 4 -2;8  0 -10 2;8 3 1 17;11 7 -3 2];
>> b=[2 ;-7 ;8 ;12];

>> x=gseidel(A,b,0.001,300)

x =

    2.9129
   -1.4015
    2.8661
   -0.8215

De tal forma que la solución es:

\begin{eqnarray}
x_1&=&2.9129 \nonumber\\
x_2&=&-1.4015 \nonumber\\
x_3&=&2.8661\nonumber\\
x_4&=&-0.8215 \nonumber
\end{eqnarray}




Verificamos en Matlab:



A*x

ans =

    2.0049
   -7.0011
    8.0000
   11.9898

De esta forma podemos ver que el método converge. A continuación les dejo la función implementada en Matlab

Función Gauss-Seidel

%Desarrollada por: Jorge De La Cruz
%Para un sistema de la forma 'Ax=b'
%imax: Número máxima de iteraciones en caso de
%que el método no converja
%tol: tolerancia deseada

function x=gseidel(A,b,tol,imax)
n=max(size(A));         % Obtenemos de la matriz A
                        % el número de ecuaiones
                        % y por ende de incognitas.
c=diag(A);              % Almacenamos en 'c' los elementos
                        % de la diagonal principal.
d=1;
i=1;
while (d~=0)&&(i<n)     % revisamos si en la diagonal principal 'c'
    d=c(i,1)*c(i+1,1);  % existe algún elemento igual a 0.
    i=i+1;              % En caso de que lo haya
end                     % salir del bucle y d=0.

while d==0              % Si se ha detectado algún cero en la
    B(1,:)=A(n,:);      % diagonal principal se entra en este bucle
    U(1,1)=b(n,1);      % para así cambiar las filas. Se alternarán
    for i=2:n           % las filas hasta que todos los elementos
        B(i,:)=A(i-1,:);% de la diagonal principal
        U(i,1)=b(i-1,1);% sean distintos de cero.
    end                 % Todo esto para evitar una división
    A=B;                % por cero. Aunque se recomienda
    b=U;                % siempre ordenar las ecuaciones apriori.
    c=diag(A);
    d=1;i=1;
    while (d~=0)&&(i<n)
        d=c(i,1)*c(i+1,1);
        i=i+1;
    end
end

% Se selecciona un vector X(n,1) sólo para evaluar
% el error en la primera iteración.
for i=1:n
    X(i,1)=b(i,1)/A(i,i);
end

% se sespecifica el vector x(n,1) con valores iguales a cero
% para iniciar las iteraciones.
x=zeros(n,1);

err=tol+1;
iter=1;

% Se inician las iteraciones
while (tol<err)&&(iter<imax)
    for i=1:n
        if i==1
            sum=0;
            for j=i+1:n
                sum=A(i,j)*x(j,1)+sum;
            end
            x(i,1)=(b(i,1)-sum)/A(i,i);
        end
        if i==n
            sum=0;
            for j=i-1:-1:1
                sum=A(i,j)*x(j,1)+sum;
            end
            x(i,1)=(b(i,1)-sum)/A(i,i);
        end
        if 2<=i<=n-1
            sum1=0;
            for j=i+1:n
                sum1=A(i,j)*x(j,1)+sum1;
            end
            sum2=0;
            for j=i-1:-1:1
                sum2=A(i,j)*x(j,1)+sum2;
            end
            x(i,1)=(b(i,1)-sum1-sum2)/A(i,i);
        end
    end
    err=abs((x-X)./x);
    err=max(err);
    X=x;
    iter=iter+1;
end

Método de Gauss-Seidel by Jorge De La Cruz is licensed under a Creative Commons Reconocimiento 3.0 Unported License.