Página web del curso

Contenidos

Introducción

Métodos iterativos de resolución de sistemas lineales

Dada una aproximación inicial $\mathbf{x}^{(0)}$, un método iterativo genera una sucesión de aproximaciones

$$\mathbf{x}^{(0)}, \mathbf{x}^{(1)}, \mathbf{x}^{(2)} ,\ldots,\mathbf{x}^{(k)},\ldots \rightarrow \mathbf{x}$$

que converge a la solución.

Para generar esta sucesión, se repite el mismo esquema de operaciones hasta que se cumple un criterio de parada. Por ejemplo, hasta que hemos realizado un cierto número de iteraciones.

Los métodos iterativos clásicos (lineales) se basan en reescribir el problema

$$ A\mathbf{x}=\mathbf{b}\Longleftrightarrow\mathbf{x}=B\mathbf{x}+\mathbf{c} $$

donde $B$ es una matriz $n\times n$ y $\mathbf{c}$ es una matriz columna de $n$ elementos.

Algoritmo
  • Sea $\mathbf{x}^{(0)}$ una aproximación inicial a la solución
  • Para $k=1,2,\ldots$ $$\mathbf{x}^{(k)}=B\mathbf{x}^{(k-1)}+\mathbf{c}$$

La matriz $B$ se llama matriz de iteración y el vector $\mathbf{c}$ se llama vector de iteración.

Ventajas

  • Los métodos iterativos son, en general, más eficientes que los métodos directos para resolver sistemas de ecuaciones lineales grandes y de matriz hueca.
  • Si no se exige mucha precisión, se puede obtener una aproximación aceptable en un número pequeño de iteraciones.
  • Son menos sensibles a los errores de redondeo que los métodos directos.

Desventajas

  • En general, no es posible predecir el número de operaciones que se requieren para obtener una aproximación a la solución con una precisión determinada.
  • El tiempo de cálculo y la precisión del resultado pueden depender de la elección de ciertos parámetros (método de sobrerrelajación).
  • Generalmente no se gana tiempo por iteración si la matriz de coeficientes es simétrica.

Método de Jacobi

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

¡Atención! Aunque los nombres de las matrices son los mismos, esta descomposición no tiene nada que ver con la factorización $LU.$

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] D\mathbf{x}=-\left(L+U\right)\mathbf{x}+\mathbf{b} \quad \Rightarrow \quad \mathbf{x}=-D^{-1}\left(L+U\right)\mathbf{x}+D^{-1}\mathbf{b} \end{array} $$

Y si

$$ B_J=-D^{-1}\left(L+D\right) \quad \mathbf{c}_J=D^{-1}\mathbf{b} $$

y podemos escribir

$$ \mathbf{x}=B_J\mathbf{x}+\mathbf{c}_J $$

Y podemos hacer

$$ \mathbf{x}^{(k+1)}=B_J\mathbf{x}^{(k)}+\mathbf{c}_J $$

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

$$ \mathbf{x}^{(1)}=B_J\mathbf{x}^{(0)}+\mathbf{c}_J, \quad \mathbf{x}^{(2)}=B_J\mathbf{x}^{(1)}+\mathbf{c}_J, \quad \mathbf{x}^{(3)}=B_J\mathbf{x}^{(2)}+\mathbf{c}_J, \quad \ldots $$

Ejercicio

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_{J}$.
  • A partir de cada uno de los apartados anteriores ¿qué podemos concluir acerca de la convergencia del método de Jacobi?
  • Realizar 3 iteraciones por Jacobi comenzando con $\mathbf{x}^{(0)}=(0,0,0)^T$.

¿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|\\[0.1cm] \left|1\right|+\left|3\right|\nless\left|4\right|\\[0.1cm] \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 Jacobi

Se tiene que $B_{J}=-D^{-1}(L+U)$. O también:

  • Dividimos cada fila por el correspondiente elemento de la diagonal.
$$ \left(\begin{array}{ccc} 1 & 1/4 & 0\\ 1/4 & 1 & 3/4\\ 0 & 3/4 & 1 \end{array}\right) $$
  • Cambiamos todos los elementos de signo.
$$ \left(\begin{array}{ccc} -1 & -1/4 & 0\\ -1/4 & -1 & -3/4\\ 0 & -3/4 & -1 \end{array}\right) $$
  • Ponemos ceros en la diagonal principal.
$$ B_{J}=\left(\begin{array}{ccc} 0 & -1/4 & 0\\ -1/4 & 0 & -3/4\\ 0 & -3/4 & 0 \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}| $$

Y la norma infinito de $B_J$ viene dada por

$$ B_{J}=\left(\begin{array}{ccc} 0 & -1/4 & 0\\ -1/4 & 0 & -3/4\\ 0 & -3/4 & 0 \end{array}\right)\begin{array}{l} 0+1/4+0=1/4\\ 1/4+0+3/4=1\\ 0+3/4+0=3/4 \end{array} $$

Por lo tanto

$$ \left|\left|B_{J}\right|\right|_{\infty}=\textrm{Max} \left(1/4,1,3/4\right)=1 $$

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}| $$
$$ B_{J}^T=\left(\begin{array}{ccc} 0 & -1/4 & 0\\ -1/4 & 0 & -3/4\\ 0 & -3/4 & 0 \end{array}\right)\begin{array}{l} 0+1/4+0=1/4\\ 1/4+0+3/4=1\\ 0+3/4+0=3/4 \end{array} $$

Por lo tanto

$$ \left|\left|B_{J}\right|\right|_{1}=\textrm{Max} \left(1/4,1,3/4\right)=1 $$

Autovalores. Calculamos los autovalores de $B_J$ calculando las raíces de $|B_J-\lambda I|=0$

$$ \left|\begin{array}{ccc} 0-\lambda & -1/4 & 0\\ -1/4 & 0-\lambda & -3/4\\ 0 & -3/4 & 0-\lambda \end{array}\right|= -\lambda^3 -\left( -\frac{1}{16}\lambda-\frac{9}{16}\lambda\right) =$$
$$ = -\lambda^3 +\frac{10}{16}\lambda =\dfrac{5}{8}\lambda-\lambda^{3}=\lambda\left(\dfrac{5}{8}-\lambda^{2}\right)=0 $$

Y los autovalores son

$$ \lambda_{1}=0,\quad \lambda_{2,3}= \pm\sqrt{\frac{5}{8}} = \pm 0.79 $$

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

Es condición suficiente para la convergencia del método de Jacobi 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 Jacobi que alguna norma de la matriz de iteración $B_J$ sea menor que $1.$

  • Como $\left|\left|B_{J}\right|\right|_{\infty}\nless 1$ no podemos concluir nada.
  • Como $\left|\left|B_{J}\right|\right|_{1}\nless 1$ no podemos concluir nada.

Es condición necesaria y suficiente para la convergencia del método de Jacobi que todos los autovalores de la matriz de iteración $B_J$ en valor absoluto sean menores que $1.$

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

Realizar 3 iteraciones con el método de Jacobi

El sistema es $A\mathbf{x}=\mathbf{b}$ con

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

Si escribimos las ecuaciones

$$ \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 asumiendo que todos los valores a la derecha son los valores obtenidos en la iteración anterior

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

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

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

Segunda iteración:

$$ \begin{array}{ccl} x^{(2)} & = & (3-y^{(1)})/4=(3-0.75)/4=0.56\\ y^{(2)} & = & (3-x^{(1)}-3z^{(1)})/4=(3-0.75-3(1.25))/4=-0.38\\ z^{(2)} & = & (5-3y^{(1)})/4=(5-3(0.75))/4=0.69 \end{array} $$

Tercera iteración:

$$ \begin{array}{ccl} x^{(3)} & = & (3-y^{(2)})/4=(3-(-0.38))/4=0.84\\ y^{(3)} & = & (3-x^{(2)}-3z^{(2)})/4=(3-0.56-3(0.69))/4=0.09\\ z^{(3)} & = & (5-3y^{(2)})/4=(5-3(-0.38))/4=1.53 \end{array} $$

Para comparar, 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}^{(k)}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{(k-1)}-\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.