Métodos de Jacobi y Gauss-Seidel

Describimos los métodos de Jacobi y Gauss-Seidel para el cálculo aproximado de la solución de un sistema lineal compatible y determinado.

Enunciado
Sea $A\in\mathbb{C}^{n\times n}$, $\lambda_1,\ldots,\lambda_n$ sus valores propios (repetidos o no) y sea $$\rho (A)=\max \left\{\left|\lambda_1\right|,\ldots,\left|\lambda_n\right|\right\}.$$ su radio espectral. Según sabemos, $$\displaystyle \lim_{k\to +\infty}A^k=0\Leftrightarrow \rho (A)<1\qquad (1)$$ para cualquier norma en $\mathbb{C}^{n\times n}.$
1. Sea ahora $A$ invertible y consideremos el sistema lineal $Ax=b$ con $b\in\mathbb{C}^{n\times 1}$ (sistema que por supuesto es compatible y determinado). Se considera la sucesión definida por $$x^{(0)} \in\mathbb{C}^{n\times 1},\quad x^{(m)}=(I-A)x^{(m-1)}+b.\qquad (2)$$ Demostrar que la sucesión $x^{(m)}$ converge a la única solución $x$ del sistema $Ax=b$ si y sólo si, $\rho (I-A)<1.$
2. Método de Jacobi. Supongamos que la matriz $A$ se puede descomponer en la forma $A=D+T$ con $D$ diagonal y todos los términos de su diagonal principal no nulos. Demostrar que la sucesión iterativa $$x^{(m)}=-D^{-1}Tx^{(m-1)}+D^{-1}b\qquad (3)$$ tiende a la única solución $x$ del sistema $Ax=b$ para toda elección de $x^{(0)}$ si y sólo si, $\rho (-D^{-1}T)<1.$
3. Se considera el sistema lineal $\begin{pmatrix}{2}&{-1}\\{-1}&{2}\end{pmatrix}\begin{pmatrix}{x_1}\\{x_2}\end{pmatrix}=\begin{pmatrix}{8}\\{-7}\end{pmatrix}.$
   Hallar su solución exacta y una solución aproximada mediante el método de Jacobi para la elección $x^{(0)}=(1,0)^T$ y hasta la cuarta iteración.
4. Método de Gauss-Seidel. Supongamos que la matriz $A$ se puede descomponer en la forma $A=S+T$ con $S$ triangular inferior invertible y $T$ triangular superior con ceros en la diagonal principal. Demostrar que la sucesión iterativa $$x^{(m)}=-S^{-1}Tx^{(m-1)}+S^{-1}b\qquad (4)$$ tiende a la única solución $x$ del sistema $Ax=b$ para toda elección de $x^{(0)}$ si y sólo si, $\rho (-S^{-1}T)<1.$
   (e) Para el sistema del apartado (c) hallar una solución aproximada mediante el método de Gauss-Seidel para la elección $x^{(0)}=(1,0)^T$ y hasta la cuarta iteración.

Solución
1. Llamemos $\epsilon^{(m)}=x^{(m)}-x $. Demostrar que $x^{(m)}\to x$ equivale a demostrar que $\epsilon^{(m)}\to 0$. Llamemos $B=I-A.$ Tenemos $$B\epsilon^{(m-1)}=(I-A)\left(x^{(m-1)}-x \right)=(I-A)x^{(m-1)}-x+\underbrace{Ax}_{b}=x^{(m)}-x=\epsilon^{(m)}.$$ De la relación $\epsilon^{(m)}=B\epsilon^{(m-1)}$ se deduce que $\epsilon^{(m)}=B^m\epsilon^{(0)}.$ Entonces, si $\rho (B)<1$ $$\lim_{m\to +\infty}\epsilon^{(m)}=\lim_{m\to +\infty}\left(B^m\epsilon^{(0)}\right)=\left(\lim_{m\to +\infty}B^m\right)\epsilon^{(0)}\underbrace{=}_{\text{por }(1)}0\epsilon^{(0)}=0.$$ Recíprocamente, si $\epsilon^{(m)}\to 0$ para cualquier elección de $x^{(0)}$, eligiendo $\epsilon^{(0)}\ne0$ $$0=\lim_{m\to +\infty}\epsilon^{(m)}=\lim_{m\to +\infty}\left(B^m\epsilon^{(0)}\right)=\left(\lim_{m\to +\infty}B^m\right)\underbrace{\epsilon^{(0)}}_{\ne 0}$$ $$\Rightarrow \lim_{m\to +\infty}B^m=0\underbrace{\Rightarrow}_{\text{por }(1)} \rho (B)<1.$$
2. Tenemos las equivalencias $$Ax=b\Leftrightarrow (D+T)x=b\Leftrightarrow Dx+Tx=b$$ $$\Leftrightarrow x=-D^{-1}Tx+D^{-1}b\Leftrightarrow \left(I+D^{-1}T\right)x=D^{-1}b.$$ Nuestra matriz $A$ es ahora $I+D^{-1}T$ y nuestro $b$ es $D^{-1}b$, en consecuencia la sucesión iterativa $(2)$ se convierte en $x^{(m)}=-D^{-1}Tx^{(m-1)}+D^{-1}b$ que converge a la única solución $x$ del sistema $Ax=b$ para toda elección de $x^{(0)}$ si y sólo si, $\rho (-D^{-1}T)<1.$
3. Resolviendo el sistema por el método de Gauss obtenemos su solución exacta: $x=(3,-2)^T$. Apliquemos ahora el método de Jacobi. Tenemos $$A=\begin{pmatrix}{2}&{-1}\\{-1}&{2}\end{pmatrix}=\begin{pmatrix}{2}&{0}\\{0}&{2}\end{pmatrix}+\begin{pmatrix}{0}&{-1}\\{-1}&{0}\end{pmatrix}=D+T,$$ $$-D^{-1}T=-\dfrac{1}{2}\begin{pmatrix}{1}&{0}\\{0}&{1}\end{pmatrix}\begin{pmatrix}{0}&{-1}\\{-1}&{0}\end{pmatrix}=\frac{1}{2}\begin{pmatrix}{0}&{1}\\{1}&{0}\end{pmatrix},$$ $$D^{-1}b=\dfrac{1}{2}\begin{pmatrix}{1}&{0}\\{0}&{1}\end{pmatrix}\begin{pmatrix}{8}\\{-7}\end{pmatrix}=\dfrac{1}{2}\begin{pmatrix}{8}\\{-7}\end{pmatrix}.$$ Los valores propios de $-D^{-1}T$ son $\lambda =\pm 1/2$ luego $\rho(-D^{-1}T)=1/2<1$, lo cual garantiza la convergencia del método de Jacobi. Llamando $x^{(m)}=(x_1^{(m)},x_2^{(m)})^T$ la sucesión iterativa $(3)$ es $$x^{(m)}=\begin{pmatrix}{x_1^{(m)}}\\{x_2^{(m)}}\end{pmatrix}=\frac{1}{2}\begin{pmatrix}{0}&{1}\\{1}&{0}\end{pmatrix}\begin{pmatrix}{x_1^{(m-1)}}\\{x_2^{(m-1)}}\end{pmatrix}+\dfrac{1}{2}\begin{pmatrix}{8}\\{-7}\end{pmatrix}=\dfrac{1}{2}\begin{pmatrix}{8+x_2^{(m-1)}}\\{-7+x_1^{(m-1)}}\end{pmatrix}.$$ Para la elección $x^{(0)}=(1,0)^T$ vamos obteniendo $$x^{(1)}=\dfrac{1}{2}\begin{pmatrix}{8+0}\\{-7+1}\end{pmatrix}=\begin{pmatrix}{4}\\{-3}\end{pmatrix},$$ $$x^{(2)}=\dfrac{1}{2}\begin{pmatrix}{8-3}\\{-7+4}\end{pmatrix}=\begin{pmatrix}{5/2}\\{-3/2}\end{pmatrix}=\begin{pmatrix}{2.5}\\{-1.5}\end{pmatrix},$$ $$x^{(3)}=\dfrac{1}{2}\begin{pmatrix}{8-3/2}\\{-7+5/2}\end{pmatrix}=\begin{pmatrix}{13/4}\\{-9/4}\end{pmatrix}=\begin{pmatrix}{3.25}\\{-2.25}\end{pmatrix},$$ $$x^{(4)}=\dfrac{1}{2}\begin{pmatrix}{8-9/4}\\{-7+13/4}\end{pmatrix}=\begin{pmatrix}{23/8}\\{-15/8}\end{pmatrix}=\begin{pmatrix}{2.875}\\{-1.875}\end{pmatrix}.$$
4. La demostración es totalmente análoga a la del apartado 2.
5. Tenemos $$A=\begin{pmatrix}{2}&{-1}\\{-1}&{2}\end{pmatrix}=\begin{pmatrix}{2}&{0}\\{-1}&{2}\end{pmatrix}+\begin{pmatrix}{0}&{-1}\\{0}&{0}\end{pmatrix}=S+T,$$ $$-S^{-1}T=-\dfrac{1}{4}\begin{pmatrix}{2}&{0}\\{1}&{2}\end{pmatrix}\begin{pmatrix}{0}&{-1}\\{0}&{0}\end{pmatrix}=\frac{1}{4}\begin{pmatrix}{0}&{2}\\{0}&{1}\end{pmatrix},$$ $$S^{-1}b=\dfrac{1}{4}\begin{pmatrix}{2}&{0}\\{1}&{2}\end{pmatrix}\begin{pmatrix}{8}\\{-7}\end{pmatrix}=\dfrac{1}{4}\begin{pmatrix}{16}\\{-6}\end{pmatrix}.$$ Los valores propios de $-S^{-1}T$ son $\lambda =0$ y $\lambda =1/4$ luego $\rho(-S^{-1}T)=1/4<1$, lo cual garantiza la convergencia del método de Gauss-Seidel. Llamando $x^{(m)}=(x_1^{(m)},x_2^{(m)})^T$ la sucesión iterativa $(4)$ es $$x^{(m)}=\begin{pmatrix}{x_1^{(m)}}\\{x_2^{(m)}}\end{pmatrix}=\frac{1}{4}\begin{pmatrix}{0}&{2}\\{0}&{1}\end{pmatrix}\begin{pmatrix}{x_1^{(m-1)}}\\{x_2^{(m-1)}}\end{pmatrix}+\dfrac{1}{4}\begin{pmatrix}{16}\\{-6}\end{pmatrix}=\dfrac{1}{4}\begin{pmatrix}{16+2x_2^{(m-1)}}\\{-6+x_2^{(m-1)}}\end{pmatrix}.$$ Para la elección $x^{(0)}=(1,0)^T$ vamos obteniendo $$x^{(1)}=\dfrac{1}{4}\begin{pmatrix}{16+2\cdot 0}\\{-6+0}\end{pmatrix}=\begin{pmatrix}{4}\\{-3/2}\end{pmatrix}=\begin{pmatrix}{4}\\{-1.5}\end{pmatrix},$$ $$x^{(2)}=\dfrac{1}{4}\begin{pmatrix}{16+2\cdot (-3/2)}\\{-6-3/2}\end{pmatrix}=\begin{pmatrix}{13/4}\\{-15/8}\end{pmatrix}=\begin{pmatrix}{3.25}\\{-1.875}\end{pmatrix},$$ $$x^{(3)}=\dfrac{1}{4}\begin{pmatrix}{16+2\cdot (-15/8)}\\{-6-15/8}\end{pmatrix}=\begin{pmatrix}{49/16}\\{-63/32}\end{pmatrix}=\begin{pmatrix}{3.0625}\\{-1.96875}\end{pmatrix},$$ $$x^{(4)}=\dfrac{1}{4}\begin{pmatrix}{16+2\cdot (-63/32)}\\{-6-63/32}\end{pmatrix}=\begin{pmatrix}{193/64}\\{-255/128}\end{pmatrix}=\begin{pmatrix}{3.015625}\\{-1.9921875}\end{pmatrix}.$$