Página web del curso

Sea el sistema $A\mathbf{x}=\mathbf{b}$ donde

$$ A=\left(\begin{array}{ccc} 4 & 1 & 0\\ 1 & 4 & 3\\ 0 & 3 & 4 \end{array}\right)\quad\mathbf{b}=\left(\begin{array}{c} 3\\ 3\\ 5 \end{array}\right) $$
  • ¿Es $A$ diagonal dominante por filas?
  • Calcular la norma infinito, la norma uno y los autovalores de $B_{G-S}$.
  • A partir de cada uno de los apartados anteriores ¿qué podemos concluir acerca de la convergencia del método de Gauss-Seidel?
  • Realizar 3 iteraciones por Gauss-Seidel comenzando con $\mathbf{x}^{(0)}=(0,0,0)^T$.

Introducción

Método de Gauss-Seidel

Se basa en la descomposición

$$ A=L+D+U $$

donde si $$ A= \left(\begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1n}\\ a_{21} & a_{22} & \ldots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n1} & \ldots & a_{nn} \end{array}\right) $$

entonces

$$ L=\left(\begin{array}{cccc} 0 & 0 & \cdots & 0\\ a_{21} & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & \cdots & 0 \end{array}\right)\;D=\left(\begin{array}{cccc} a_{11} & 0 & \cdots & 0\\ 0 & a_{22} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & a_{nn} \end{array}\right)\;U=\left(\begin{array}{cccc} 0 & a_{12} & \cdots & a_{1n}\\ 0 & 0 & \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 0 \end{array}\right) $$

Tenemos que

$$ \begin{array}{l} A\mathbf{x}=\mathbf{b} \quad \Rightarrow \quad \left(L+D+U\right)\mathbf{x}=\mathbf{b} \quad \Rightarrow \quad \\[0.5cm] (L+D)\,\mathbf{x}=-U\mathbf{x}+\mathbf{b} \quad \Rightarrow \quad \mathbf{x}=-\left(L+D\right)^{-1}U\mathbf{x}+\left(L+D\right)^{-1}\mathbf{b} \end{array} $$

con

$$ B_{G-S}=-\left(L+D\right)^{-1}U \quad \mathbf{c}_{G-S}=\left(L+D\right)^{-1}\mathbf{b} $$

y podemos escribir

$$ \mathbf{x}=B_{G-S}\mathbf{x}+\mathbf{c}_{G-S} $$

Y el método vendrá dado por

$$ \mathbf{x}^{(k+1)}=B_{G-S}\mathbf{x}^{(k)}+\mathbf{c}_{G-S} $$

Y si comenzamos con un $\mathbf{x}^{(0)}$

$$ \mathbf{x}^{(1)}=B_{G-S}\mathbf{x}^{(0)}+\mathbf{c}_{G-S}, \quad \mathbf{x}^{(2)}=B_{G-S}\mathbf{x}^{(1)}+\mathbf{c}_{G-S}, \quad \mathbf{x}^{(3)}=B_{G-S}\mathbf{x}^{(2)}+\mathbf{c}_{G-S}, \quad \ldots $$

Ejercicio

¿Es A diagonal dominante por filas?

Se dice que una matriz $A$ de $n$ filas y $n$ columnas es diagonal dominante por filas si

$$ \sum_{j=1\\i\ne j}^n|a_{ij}|\lt |a_{ii}|\quad i=1,\ldots,n $$

Para nuestra matriz $A$

$$A=\left(\begin{array}{crr} \fbox{$4$} & 1 & 0\\ 1 & \fbox{$4$} & 3\\ 0 & 3 & \fbox{$4$} \end{array}\right)\begin{array}{r} \left|1\right|+\left|0\right|\lt\left|4\right|\\ \left|1\right|+\left|3\right|\nless\left|4\right|\\ \left|0\right|+\left|3\right|\lt\left|4\right| \\ \end{array}$$

Y no es diagonal dominante por filas puesto que en la segunda fila no se verifica la desigualdad ya que $\left|4\right|=\left|1\right|+\left|3\right|$

Calcular la norma infinito, la norma uno y los autovalores de la matriz de iteración de Gauss-Seidel

Se tiene que $B_{G-S}=-(L+D)^{-1}U$. Si

$$ A=\left(\begin{array}{ccc} 4 & 1 & 0\\ 1 & 4 & 3\\ 0 & 3 & 4 \end{array}\right) $$

se tiene que

$$ L=\left(\begin{array}{ccc} 0 & 0 & 0\\ 1 & 0 & 0\\ 0 & 3 & 0 \end{array}\right) \quad D=\left(\begin{array}{ccc} 4 & 0 & 0\\ 0 & 4 & 0\\ 0 & 0 & 4 \end{array}\right) \quad U=\left(\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 3\\ 0 & 0 & 0 \end{array}\right) $$$$ L+D=\left(\begin{array}{ccc} 4 & 0 & 0\\ 1 & 4 & 0\\ 0 & 3 & 4 \end{array}\right) $$

Calculemos su inversa por Gauss-Jordan:

Escribimos la matriz $\left[(L+D)|I\right]$. Como el pivote es $\mathbf{{\color{red}4}}$ dividimos toda la fila por 4 para convertirlo en 1.

$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} {\mathbf{{\color{red}4}}} & 0 &0 & 1 & 0 & 0\\ \mathbf{{\color{blue}1}} & 4 & 0 & 0 & 1 & 0\\ \mathbf{{\color{ForestGreen}0}} & 3 & 4 & 0 & 0 & 1 \end{array}\right) \begin{array}{rrl} f_{1} & = & f_{1}/(\mathbf{\color{red}{4}})\\ & &\\ & & \end{array} $$
$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} {\mathbf{{\color{red}1}}} & 0 & 0 & 1/4 & 0 & 0\\ \mathbf{{\color{blue}1}} & 4 & 0 & 0 & 1 & 0\\ \mathbf{{\color{ForestGreen}0}} & 3 & 4 & 0 & 0 & 1 \end{array}\right) \begin{array}{rrl} f_{1}\\ f_{2} & = & f_{2}-\mathbf{{\color{blue}1}}f_{1}\\ f_3 & = & f_{3}-\mathbf{{\color{ForestGreen}0}}f_{1} \end{array} $$

Hacemos cero por debajo del pivote.

$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} 1 & \mathbf{{\color{blue}0}} & 0 & 1/4 & 0 & 0\\ 0 & {\mathbf{{\color{red}4}}} & 0 & -1/4 & 1 & 0\\ 0 & \mathbf{{\color{ForestGreen}3}} & 4 & 0 & 0 & 1 \end{array}\right) \begin{array}{rrl} & &\\ f_{2} & = & f_{2}/(\mathbf{\color{red}{4}})\\ & & \end{array} $$

Como el pivote es $\mathbf{{\color{red}4}}$ dividimos toda la fila por 4 para convertirlo en 1.

$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} 1 & \mathbf{{\color{blue}0}} & 0 & 1/4 & 0 & 0\\ 0 & {\mathbf{{\color{red}1}}} & 0 & -1/16 & 1/4 & 0\\ 0 & \mathbf{{\color{ForestGreen}3}} & 4 & 0 & 0 & 1 \end{array}\right) \begin{array}{rrl} f_{1} & = & f_{1}-\mathbf{{\color{blue}0}}\,f_{2}\\ f_{2}\\ f_3 & = & f_{3}-\mathbf{{\color{ForestGreen}3}}\,f_{2} \end{array} $$

Hacemos ceros por encima y por debajo del pivote (aunque por encima ya hay un cero y no hay que hacer nada)

$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} 1 & 0 & \mathbf{{\color{blue}0}} & 1/4 & 0 & 0\\ 0 & 1 & \mathbf{{\color{ForestGreen}0}} & -1/16 & 1/4 & 0\\ 0 & 0 & {\mathbf{{\color{red}4}}} & 3/16 & -3/4 & 1 \end{array}\right) \begin{array}{rrl} & &\\ & &\\ f_{3} & = & f_{3}/(\mathbf{\color{red}{4}}) \end{array} $$

Dividimos por el pivote la fila del pivote, y como por encima de él sólo hay ceros, ya hemos acabado.

$$ \begin{array}{c} f_{1}\\ f_{2}\\ f_{3} \end{array}\left(\begin{array}{rrr|rrr} 1 & 0 & 0 & 1/4 & 0 & 0\\ 0 & 1 & 0 & -1/16 & 1/4 & 0\\ 0 & 0 & 1 & 3/64 & -3/16 & 1/4 \end{array}\right) $$

Como a la izquierda ya tenemos la matriz $I$ la matriz de la derecha será $(L+D)^{-1}$.

$$ (L+D)^{-1} = \left(\begin{array}{rrr} 1/4 & 0 & 0\\ -1/16 & 1/4 & 0\\ 3/64 & -3/16 & 1/4 \end{array}\right) $$

La multiplicamos por $U$

$$ B_{G-S} = -(L+D)^{-1}\,U= -\left(\begin{array}{rrr} 1/4 & 0 & 0\\ -1/16 & 1/4 & 0\\ 3/64 & -3/16 & 1/4 \end{array}\right) \left(\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 3\\ 0 & 0 & 0 \end{array}\right)= \left(\begin{array}{rrr} 0 & -1/4 & 0\\ 0 & 1/16 & -3/4\\ 0 & -3/64 & 9/16 \end{array}\right) $$

Norma infinito. Si $A$ es una matriz $m\times n$ su norma infinito viene dada por:

$$ \|A\|_{\infty} = \max_{1\leq i \leq m}\sum_{j=1}^n|a_{ij}| $$

Pasamos a coma flotante para comparar mejor cantidades y estudiamos la norma infinito

$$ B_{G-S}=\left(\begin{array}{rrr} 0 & -0.25 & 0\\ 0 & 0.06 & -0.75\\ 0 & -0.05 & 0.56 \end{array}\right)\begin{array}{l} 0+0.25+0=0.25\\ 0+0.06+0.75=0.81\\ 0+0.05+0.56=0.61 \end{array} $$

Por lo tanto

$$ \left|\left|B_{G-S}\right|\right|_{\infty}=\textrm{Max} \left(0.25,0.81,0.61\right)=0.81 $$

Norma uno. Si $A$ es una matriz $m\times n$ su norma uno viene dada por:

$$ \|A\|_1 = \max_{1\leq j \leq n}\sum_{i=1}^m|a_{ij}| $$

Calculamos la norma uno

$$ B_{G-S}^T=\left(\begin{array}{ccc} 0 & 0 & 0\\ -0.25 & 0.06 & -0.05\\ 0 & -0.75 & 0.56 \end{array}\right)\begin{array}{l} 0+0+0=0\\ 0.25+0.06+0.05=0.36\\ 0+0.75+0.56=1.31 \end{array} $$

Por lo tanto

$$ \left|\left|B_{G-S}\right|\right|_{1}=\textrm{Max} \left(0,0.36,1.31\right)=1.31 $$

Autovalores. Calculamos los autovalores de $B_{G-S}$

$$ \left(\begin{array}{rrr} 0 & -1/4 & 0\\ 0 & 1/16 & -3/4\\ 0 & -3/64 & 9/16 \end{array}\right) $$

Calculamos las raíces de $|B_{G-S}-\lambda I|=0$

$$ \left|\begin{array}{ccc} 0-\lambda & -1/4 & 0\\ 0 & (1/16)-\lambda & -3/4\\ 0 & -3/64 & (9/16)-\lambda \end{array}\right|= $$

. $$ =-\lambda\left|\begin{array}{cc} (1/16)-\lambda & -3/4\\ -3/64 & (9/16)-\lambda \end{array}\right|= $$ . $$ =-\lambda \left[ \left(\frac{1}{16}-\lambda\right)\left(\frac{9}{16}-\lambda\right)-\left(-\frac{3}{4}\right)\left(-\frac{3}{64}\right)\right]= $$ . $$ =-\lambda\left[\left(\frac{1}{16}-\lambda\right)\left(\frac{9}{16}-\lambda\right)-\frac{9}{4\times 64} \right]= $$ . $$ =-\lambda\left[\frac{1}{16}\times\frac{9}{16}-\frac{1}{16}\lambda-\frac{9}{16}\lambda+\lambda^2-\frac{9}{2^2\times 2^6} \right]= $$ .

$$ =-\lambda\left[\frac{9}{2^8}-\frac{10}{16}\lambda+\lambda^2-\frac{9}{2^2\times 2^6} \right]=-\lambda^2\left(-\frac{10}{16}+\lambda\right)=-\lambda^2\left(\lambda-\frac{5}{8}\right) $$

Y los autovalores son

$$ \lambda_{1,2}=0,\quad\lambda_{3}=\frac{5}{8}=0.625 $$

A partir de cada uno de los apartados anteriores ¿qué podemos concluir acerca de la convergencia del método de Gauss-Seidel?

Es condición suficiente para la convergencia del método de Gauss-Seidel que la matriz de coeficientes $A$ sea diagonal dominante por filas:

  • Como $A$ no es diagonal dominante, no podemos concluir nada.

Es condición suficiente para la convergencia del método de Gauss-Seidel que alguna norma de la matriz de iteración $B_{G-S}$ sea menor que $1.$

  • Como $\left|\left|B_{G-S}\right|\right|_{\infty}\lt 1$ podemos concluir que el método de Gauss-Seidel converge para este sistema. (Y no haría falta seguir investigando y obtener la otra norma ni los autovalores, aunque lo hemos hecho para comprobar que los resultados son coherentes).
  • De $\left|\left|B_{G-S}\right|\right|_{1}\gt 1$ no podríamos concluir nada.

Es condición necesaria y suficiente para la convergencia del método de Gauss-Seidel que todos los autovalores de la matriz de iteración $B_{G-S}$ en valor absoluto sean menores que $1.$

  • Como todos los autovalores de $B_{G-S}$ son menores que uno en valor absoluto podemos concluir que el Método de Gauss-Seidel será convergente para cualquier valor inicial.

Realizar 3 iteraciones con Gauss-Seidel

El sistema es

$$ \begin{array}{rrrrrr} 4x & +&y & & & = & 3\\ x & +&4y & +&3z & = & 3\\ & &3y & +&4z & = & 5 \end{array} $$

En la primera ecuación despejamos la primera incógnita, $x$, en la segunda, la segunda incógnita, $y$, y finalmente, $z.$

$$ \begin{array}{ccl} x & = & \dfrac{3-y}{4}\\ y & = & \dfrac{3-x-3z}{4}\\ z & = & \dfrac{5-3y}{4} \end{array} $$

Realizamos las iteraciones tomando los valores de la iteración anterior, pero los valores nuevos que vayamos calculando, los usamos ya

$$ \begin{array}{ccl} x^{\color{red}{(1)}} & = & \dfrac{3-y^{(0)}}{4}\\ y^{\color{red}{(1)}} & = & \dfrac{3-x^{\color{red}{(1)}}-3z^{(0)}}{4}\\ z^{\color{red}{(1)}} & = & \dfrac{5-3y^{\color{red}{(1)}}}{4} \end{array} $$

Realizamos una iteración, tomando como valor inicial, el vector nulo

$$ \begin{array}{ccl} x^{\color{red}{(1)}} & = & (3-y^{(0)})/4=(3-0)/4=3/4= 0.75\\ y^{\color{red}{(1)}} & = & (3-x^{\color{red}{(1)}}-3z^{(0)})/4=(3-0.75-0)/4=0.56\\ z^{\color{red}{(1)}} & = & (5-3y^{\color{red}{(1)}})/4= (5-3(0.56))/4= 0.83 \end{array} $$

Realizamos la iteración 2

$$ \begin{array}{ccl} x^{\color{red}{(2)}} & = & (3-y^{(1)})/4=(3-0.56)/4=0.61\\ y^{\color{red}{(2)}} & = & (3-x^{\color{red}{(2)}}-3z^{(1)})/4=(3-0.61-3(0.83))/4=-0.02\\ z^{\color{red}{(2)}} & = & (5-3y^{\color{red}{(2)}})/4= (5-3(-0.02))/4= 1.27 \end{array} $$

Realizamos la iteración 3

$$ \begin{array}{ccl} x^{\color{red}{(3)}} & = & (3-y^{(2)})/4=(3-(-0.02))/4=0.76\\ y^{\color{red}{(3)}} & = & (3-x^{\color{red}{(3)}}-3z^{(2)})/4=(3-0.76-3(1.27))/4=-0.39\\ z^{\color{red}{(3)}} & = & (5-3y^{\color{red}{(3)}})/4= (5-3(-0.39))/4= 1.54 \end{array} $$

La solución exacta es

$$x = 1\quad y = -1 \quad z = 2$$

Algoritmo

El algoritmo que hemos aplicado se expresa formalmente

  • Elegir una aproximación inicial $\mathbf{x}^{(0)}$
  • Para $k=1,2,\ldots,$MaxIter
    • Para $i=1,2,\ldots,n$, calcular $$x_{i}^{\color{red}{(k)}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
    • Si se cumple el criterio de parada, tomar $\mathbf{x}^{(k)}$ como aproximación a la solución.