Consideramos sistemas de ecuaciones lineales $n \times n$ de la forma $Ax=b$ con incógnita $x \in \mathbb{R}^n$.
Se recomienda el uso de los métodos iterativos para resolver sistemas grandes (mayores a $10^{8}$ entradas, por ejemplo para matrices mayores a un número de renglones y columnas que $A \in \mathbb{R}^{10^4 \times 10^4}$) y también para matrices sparse o ralas (con muchos ceros) pues tales métodos son eficientes en el uso de memoria y cómputo.
Todos los métodos iterativos inician con una aproximación inicial $x^{(0)}$ a la solución $x$ para generar una secuencia de vectores $x^{(k)}$, $k=0,1,\dots$. Bajo condiciones descritas en el método después de un número de iteraciones $x^{(k)}$ converge a $x$. La elección de la aproximación incial y la convergencia (o no convergencia) es dependiente del método. Hay métodos que convergen bajo condiciones simples y con cualquier elección de aproximación inicial, otros métodos son más restrictivos en estas elecciones.
En esta nota revisarmos de una forma resumida los métodos iterativos de Jacobi y de Gauss-Seidel. Para profundizar en éstos métodos se recomienda consultar las referencias al final de la nota.
En una primera descripción del método se asume que $a_{ii} \neq 0$, $\forall i=1,2,\dots,n$.
El método iterativo de Jacobi se obtiene al resolver la $i$-ésima equación en $Ax=b$ para $x_i$ resultando en la ecuación:
$$x_i = \displaystyle \sum_{\overset{j=1}{j \neq i}} ^n \left ( - \frac{a_{ij}x_j}{a_{ii}} \right ) + \frac{b_i}{a_{ii}}, i=1,2,\dots,n$$Una vez realizado lo anterior, para cada $k \geq 1$ las componentes $x_i^{(k)}$ de $x^{(k)}$ se generan a partir de las componentes de $x^{(k-1)}$ con la expresión:
$$ x_i^{(k)} = \frac{1}{a_{ii}} \left [ \displaystyle \sum_{\overset{j=1}{j \neq i}} ^n \left ( -a_{ij}x_{j}^{(k-1)} \right ) + b_i \right ], i=1,2,\dots,n $$Ejemplo:
Considérese el sistema lineal $Ax=b$ dado por:
$$ \begin{array}{ccc} 10x_1-x_2+2x_3 &=& 6 \\ -x_1+11x_2-x_3+3x_4 &=& 25 \\ 2x_1-x_2+10x_3-x_4 &=& -11 \\ 3x_2-x_3+8x_4 &=& 15 \end{array} $$el cual tiene solución $x=(1,2,-1,1)^t$. Usar el método iterativo de Jacobi para encontrar aproximaciones $x^{(k)}$ a $x$ iniciando con $x^{(0)} = (0,0,0,0)^t$ hasta tener $\frac{||x^{(k)}-x^{(k-1)}||_\infty}{||x^{(k)}||_\infty} < 10^{-3}$.
Solución:
Resolvemos la $i$-ésima ecuación para cada $x_i$, $i=1,2,3,4$ para obtener:
$$ \begin{array}{ccc} x_1 &=& \frac{1}{10}x_2 - \frac{1}{5}x_3 + \frac{3}{5} \\ x_2 &=& \frac{1}{11}x_1 + \frac{1}{11}x_3-\frac{3}{11}x_4+\frac{25}{11}\\ x_3 &=& -\frac{1}{5}x_1 + \frac{1}{10}x_2 + \frac{1}{10}x_4 - \frac{11}{10}\\ x_4 &=& -\frac{3}{8}x_2 + \frac{1}{8}x_3 + \frac{15}{8} \end{array} $$Utilizando la aproximación inicial $x^{(0)} = (0,0,0,0)^t$ se obtiene $x^{(1)}$ dado por:
$$ \begin{array}{ccc} x_1^{(1)} &=& \frac{1}{10}x_2^{(0)} - \frac{1}{5}x_3^{(0)} + \frac{3}{5} &=& 0.6 \\ x_2^{(1)} &=& \frac{1}{11}x_1^{(0)} + \frac{1}{11}x_3^{(0)} - \frac{3}{11}x_4^{(0)} + \frac{25}{11} &\approx& 2.2727 \\ x_3^{(1)} &=& -\frac{1}{5}x_1^{(0)}+\frac{1}{10}x_2^{(0)}+\frac{1}{10}x_4^{(0)}-\frac{11}{10}&=&-1.1 \\ x_4^{(1)} &=& -\frac{3}{8}x_2^{(0)}+\frac{1}{8}x_3^{(0)} + \frac{15}{8} &=& 1.875 \end{array} $$Iteraciones $x^{(k)} = (x_1^{(k)},x_2^{(k)}, x_3^{(k)}, x_4^{(k)})$ son generadas de una forma similar y los resultados de las mismas son presentadas a continuación:
| k | 0 | 1 | 2 | ... | 10 |
| $x_1$(k) | 0 | 0.6 | 1.0473 | ... | 1.0001 |
| $x_2$(k) | 0 | 2.2727 | 1.7159 | ... | 1.998 |
| $x_3$(k) | 0 | -1.1 | -0.8052 | ... | -0.9998 |
| $x_4$(k) | 0 | 1.875 | 0.8852 | ... | 0.9998 |
Se cumple que para la iteración $10$: $\frac{||x^{(10)}-x^{(9)}||_\infty}{||x^{(10)}||_\infty}= \frac{.8 \times 10^{-3}}{1.9998} < 10^{-3}$ y $||x^{(10)}-x||_\infty = .2 \times 10^{-3}$.
Ejercicio: verificar la iteración $2$ de la tabla anterior.
Comentario:
Algoritmo
Entrada:
Salida:
Paso 1 $k=1$
Paso 2 Mientras $k \leq maxiter$ realizar pasos 3 a 6:
Paso 3 Para $i=1,\dots,n$. Establecer:
$$x_i = \frac{1}{a_{ii}} \left [ - \displaystyle \sum_{\overset{j=1}{j \neq i}}^n(a_{ij}x_j^{(0)})+b_i \right ]$$para las entradas del vector $x \in \mathbb{R}^n$.
Paso 4 Revisar si se cumple la condición $||x-x^{(k-1)}|| < (1+||x||)*TOL$. Si se cumple lo anterior entonces devolver $x$ y un mensaje del tipo "procedimiento exitoso. El error entre iteraciones es de: ... ". Detener algoritmo
Paso 5 Establecer $k=k+1$.
Paso 6 Para $i=1,\dots, n$ establecer $x_i^{(0)} = x_i$.
Paso 7 Devolver mensaje del tipo "máximo número de iteraciones alcanzadas".
Nota: los pasos en el algoritmo anterior representan una guía para la implementación. Al describirse los pasos de un algoritmo no implica que se tengan que implementar uno a continuación del otro como se describe. Si una implementación respeta la lógica y al mismo método, entonces pueden seguirse los pasos de una forma distinta.
Este método se basa en las iteraciones de Jacobi pero se aprovechan los cálculos realizados de una mejor forma. Para el ejemplo anterior con Jacobi en la primera iteración se escribió:
$$x_1^{(1)} = \frac{1}{10}x_2^{(0)} - \frac{1}{5}x_3^{(0)} + \frac{3}{5} = 0.6$$Y para la actualización $x_2^{(1)}$ se escribió:
En el método de Gauss-Seidel se hace el ajuste de utilizar la $x_1^{(0)}$ calculada para actualizar $x_2^{(1)}$:
$$x_2^{(1)} = \frac{1}{11}x_1^{(1)} + \frac{1}{11}x_3^{(0)} - \frac{3}{11}x_4^{(0)} + \frac{25}{11}$$De modo que $x_2^{(1)}$ queda:
$$x_2^{(1)} = \frac{1}{11}x_1^{(1)} + \frac{1}{11}x_3^{(0)} - \frac{3}{11}x_4^{(0)} + \frac{25}{11} \approx 2.3272$$Para $x_3^{(1)}$ se tendrá:
$$x_3^{(1)} = -\frac{1}{5}x_1^{(1)}+\frac{1}{10}x_2^{(1)}+\frac{1}{10}x_4^{(0)}-\frac{11}{10}$$En general para la actualización $x^{(k)}$ se tendrá:
$$x_i^{(k)} = \frac{1}{a_{ii}} \left [-\displaystyle \sum_{j=1}^{i-1}(a_{ij}x_j^{(k)})-\sum_{j=i+1}^n(a_{ij}x_j^{(k-1)}) + b_i \right ], i=1,2,\dots,n$$Comentario:
Ejemplo:
Considérese el mismo sistema lineal $Ax=b$ del ejemplo del método de Jacobi dado por:
$$ \begin{array}{ccc} 10x_1-x_2+2x_3 &=& 6 \\ -x_1+11x_2-x_3+3x_4 &=& 25 \\ 2x_1-x_2+10x_3-x_4 &=& -11 \\ 3x_2-x_3+8x_4 &=& 15 \end{array} $$el cual tiene solución $x=(1,2,-1,1)^t$. Usar el método iterativo de Gauss-Seidel para encontrar aproximaciones $x^{(k)}$ a $x$ iniciando con $x^{(0)} = (0,0,0,0)^t$ hasta tener $\frac{||x^{(k)}-x^{(k-1)}||_\infty}{||x^{(k)}||_\infty} < 10^{-3}$.
La actualización en el método de Gauss-Seidel se describe para cada $k=1,2,\dots$ como:
$$ \begin{array}{ccc} x_1^{(k)} &=& \frac{1}{10}x_2^{(k-1)} - \frac{1}{5}x_3^{(k-1)} + \frac{3}{5} \\ x_2^{(k)} &=& \frac{1}{11}x_1^{(k)} + \frac{1}{11}x_3^{(k-1)} - \frac{3}{11}x_4^{(k-1)} + \frac{25}{11} \\ x_3^{(k)} &=& -\frac{1}{5}x_1^{(k)}+\frac{1}{10}x_2^{(k)}+\frac{1}{10}x_4^{(k-1)}-\frac{11}{10}\\ x_4^{(k)} &=& -\frac{3}{8}x_2^{(k)}+\frac{1}{8}x_3^{(k)} + \frac{15}{8} \end{array} $$Con las actualizaciones anteriores se obtiene una tabla de la forma:
| k | 0 | 1 | 2 | ... | 5 |
| $x_1$(k) | 0 | 0.6 | 1.030 | ... | 1.0001 |
| $x_2$(k) | 0 | 2.3272 | 2.037 | ... | 2.0 |
| $x_3$(k) | 0 | -0.9873 | -1.014 | ... | -1.0 |
| $x_4$(k) | 0 | 0.8789 | 0.9844 | ... | 1.0 |
Se cumple que para la iteración $5$: $\frac{||x^{(5)}-x^{(4)}||_\infty}{||x^{(5)}||_\infty}= \frac{.8 \times 10^{-3}}{2.0} = 4 \times 10^{-4} < 10^{-3}$ y $||x^{(5)}-x||_\infty \approx 4.9 \times 10^{-5}$.
Ejercicio: verificar la iteración $2$ de la tabla anterior.
Comentario: Para este ejemplo el método de Gauss-Seidel realiza menos iteraciones que el de Jacobi. Esto normalmente se cumple pero existen sistemas en los que el método de Jacobi converge y el de Gauss-Seidel no o viceversa (Gauss-Seidel converge y Jacobi no).
Algoritmo
Entrada:
Salida:
Paso 1 $k=1$
Paso 2 Mientras $k \leq maxiter$ realizar pasos 3 a 6:
Paso 3 Para $i=1,\dots,n$. Establecer:
$$x_i = \frac{1}{a_{ii}} \left [-\displaystyle \sum_{j=1}^{i-1}(a_{ij}x_j)-\sum_{j=i+1}^n(a_{ij}x_j^{(0)}) + b_i \right ], i=1,2,\dots,n$$para las entradas del vector $x \in \mathbb{R}^n$.
Paso 4 Revisar si se cumple la condición $||x-x^{(k-1)}|| < (1+||x||)*TOL$. Si se cumple lo anterior entonces devolver $x$ y un mensaje del tipo "procedimiento exitoso. El error entre iteraciones es de: ... ". Detener algoritmo
Paso 5 Establecer $k=k+1$.
Paso 6 Para $i=1,\dots, n$ establecer $x_i^{(0)} = x_i$.
Paso 7 Devolver mensaje del tipo "máximo número de iteraciones alcanzadas".
Nota: los pasos en el algoritmo anterior representan una guía para la implementación. Al describirse los pasos de un algoritmo no implica que se tengan que implementar uno a continuación del otro como se describe. Si una implementación respeta la lógica y al mismo método, entonces pueden seguirse los pasos de una forma distinta.
Comentario: Tanto para el método de Jacobi como para el de Gauss-Seidel se requiere encontrar alguna $a_{ii} \neq 0$ para la actualización de la variable $x_i$. Si se cumple tal requisito y además condiciones que aseguran convergencia de los métodos, entonces el sistema $Ax=b$ tiene solución.
Las siguientes condiciones para asegurar convergencia de los métodos iterativos requieren que sus actualizaciones se escriban en una forma:
$$x^{(k)} = Tx^{(k-1)} + c, \text{ para } k=1,2,\dots$$con $x^{(0)}$ aproximación inicial, $T \in \mathbb{R}^{n \times n}, c \in \mathbb{R}^n$. Los métodos de Jacobi y Gauss-Seidel sí podemos escribirlos en la forma anterior. Por ejemplo $T_j = D^{-1}(L+U)$, $T_g = (D-L)^{-1}U)$ para Jacobi y Gauss-Seidel respectivamente. La matriz $D$ es la diagonal de $A$, $L$ es la porción triangular inferior sin la diagonal de $A$ y $U$ es la porción triangular superior sin la diagonal de $A$. Ver la referencia $1$ para más detalles.
Teorema
Para cualquier $x^{(0)} \in \mathbb{R}^n$ la secuencia $x^{(k)}$ definida por $x^{(k)} = Tx^{(k-1)} + c, \text{ para } k=1,2,\dots$ con $c \in \mathbb{R}^n$, converge a la solución única de $Ax=b$ si y sólo si $\rho(T) < 1$.
Obs: $\rho(T)$ es el radio espectral de $T$ definido como $\rho(T) = \max\{|\lambda_1|, |\lambda_2|, \dots, |\lambda_n|\}$ para $\lambda_i$ $i$-ésimo eigenvalor de $T$.
Teorema
Si $A$ es una matriz dominante esctricta en la diagonal por renglones* entonces para cualquier elección $x^{(0)}$ los métodos de Jacobi y de Gauss-Seidel con sus actualizaciones respectivas convergen a la solución única de $Ax=b$.
*Una matriz es dominante estricta en la diagonal por renglones si $|a_{ii}| > \displaystyle \sum_{\overset{j=1}{j\neq i}}^n |a_{ij}|$ $\forall i=1,\dots,n$.
Comentario: se tiene la propiedad: si $A$ es diagonal dominante estricta entonces $A$ es no singular.
Teorema
Si $||T|| < 1$ para cualquier norma matricial y $c$ es un vector, entonces la secuencia $x^{(k)}$ definida por $x^{(k)} = Tx^{(k-1)} + c, \text{ para } k=1,2,\dots$ converge para cualquier $x^{(0)}$ al vector $x \in \mathbb{R}^{n}$ solución de $Ax=b$. Además $||x^{(k)}-x|| \approx \rho(T)^{k} ||x^{(0)}-x||$
Comentario: por el teorema anterior la velocidad de convergencia depende del valor $\rho(T)$. No hay resultado general que defina cuál de los dos métodos, Jacobi o Gauss-Seidel, converja más rápido para un sistema de ecuaciones lineales arbitrario aunque hay casos especiales. Uno de los casos especiales es si $a_{ij} \leq 0$ para cada $i \neq j$ y $a_{ii}>0$ para $i=1,2,\dots,n$ entonces si uno de los dos métodos converge entonces el otro también converge y el método de Gauss-Seidel lo hace de forma más rápida.
Referencias: