4  Sistemas de ecuaciones lineales

Definición 4.1 (Sistema de ecuaciones lineales) Un sistema de \(m\) ecuaciones lineales con \(n\) incógnitas es un conjunto de ecuaciones sobre un cuerpo \(\mathbb K\) como sigue: \[\left\{\begin{array}{ccccccccc} a_{11}x_1& + &a_{12}x_2& +& \cdots &+& a_{1n}x_n & = & b_1\\ a_{21}x_1& + &a_{22}x_2& +& \cdots &+& a_{2n}x_n & = & b_2\\ \vdots& &\vdots& & \ddots& & \vdots & & \vdots\\ a_{m1}x_1& + &a_{m2}x_2& +& \cdots &+& a_{mn}x_n & = & b_m \end{array}\right.\] donde \(a_{i,j}, b_i, x_j\in\mathbb K\) para todo \(i\in\{1,\ldots,m\}\), \(j\in\{1,\ldots,n\}\).

A los \(a_{i,j}\) se les denomina coeficientes; y a los términos \(b_i\) que aparecen en el segundo miembro, términos independientes.

Comúnmente, nos referiremos a los \(x_i\) como incógnitas.

Lo más interesante de un sistema de ecuaciones es conocer si existen valores que, asignados convenientemente a las incógnitas \(x_j\) en las ecuaciones anteriores, hacen que todas las igualdades sean ciertas. A esos valores los denominaremos solución.

Definición 4.2 (Solución) Llamamos solución de un sistema de ecuaciones lineales a todo elemento \(\symbfit{x} = (x_1,x_2,\ldots,x_n)\in\mathbb{K}^n\) que haga ciertas todas las igualdades presentes en el sistema. Asimismo, llamamos resolver el sistema a encontrar todas las soluciones.

En cuanto a la existencia de soluciones, hay tres tipos de sistemas de ecuaciones lineales:

• Si un sistema no tiene ninguna solución, decimos que es incompatible.

• En caso de que tenga, al menos, una solución, diremos que el sistema es compatible. Además, si tiene:

  • Una única solución, decimos que es compatible determinado.
  • Más de una solución, decimos que es un sistema compatible indeterminado.

Observación 4.1. Un sistema homogéneo (\(b_i = 0\)) siempre es compatible.

4.1 Sistemas equivalentes

Comenzamos por una definición:

Definición 4.3 (Sistemas equivalentes) Dos sistemas de ecuaciones con el mismo número de incógnitas se dicen equivalentes si tienen el mismo conjunto de soluciones.

Dado un sistema de ecuaciones lineales podemos pasar a otro sistema equivalente realizando alguna de las siguientes manipulaciones:

Tipo 1. Intercambiar el orden en el que figuran las ecuaciones en el sistema.

Tipo 2. Multiplicar una de las ecuaciones por cualquier escalar no nulo.

Tipo 3. Sumar a una ecuación un múltiplo de otra.

A cada una de estas estas operaciones, que trivialmente mantienen el conjunto de soluciones, se le llama operación o transformación elemental.

Aplicar sucesivamente estas transformaciones para pasar a un sistema equivalente más fácil de resolver.

Observación 4.2. Éste es el método de Gauss: transformar un sistema de ecuaciones en uno triangular fácil de resolver (sustitución regresiva).

4.2 Forma matricial

Dado el sistema lineal de \(m\) ecuaciones y \(n\) incógnitas \[ \left\{\begin{array}{ccccccccc} a_{11}x_1& + &a_{12}x_2& +& \cdots &+& a_{1n}x_n & = & b_1\\ a_{21}x_1& + &a_{22}x_2& +& \cdots &+& a_{2n}x_n & = & b_2\\ \vdots& &\vdots& & \ddots & & \vdots & & \vdots\\ a_{m1}x_1& + &a_{m2}x_2& +& \cdots &+& a_{mn}x_n & = & b_m \end{array}\right. \] se definen las siguientes matrices \[ A = \begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots& \ddots &\vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn} \end{pmatrix}, \qquad \symbfit{x} = \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{pmatrix}, \qquad \symbfit{b} = \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{pmatrix} \] Entonces, el sistema se puede escribir como \(A\,\symbfit{x}={b}\).

La matriz \(A\in\mathcal{M}_{m\times n}(\mathbb K)\) se denomina matriz de coeficientes del sistema lineal, \(\symbfit{b}\in\mathcal{M}_{m, 1}(\mathbb K)\) es el vector de términos independientes.

\(A^* = (A \mid \symbfit{b} )\) es la matriz aumentada del sistema lineal: \[A^* = ( A \mid \symbfit{b} ) = \left(\begin{array}{cccc|c} a_{11}&a_{12}&\cdots&a_{1n}& b_1\\ a_{21}&a_{22}&\cdots&a_{2n}& b_2\\ \vdots&\vdots& \ddots &\vdots & \vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn}& b_m \end{array}\right) \]

Las operaciones elementales que vimos anteriormente sobre un sistema se traducen de manera obvia a operaciones elementales por filas en la matriz aumentada.

Ejercicio 4.1 Consideremos el sistema dado por: \[\left\{\begin{array}{rrrcr} x & & -z & = & 3\\ & y & + 2z & = & 1\\ 2x & + y & & = & 7\\ \end{array}\right.\] Apliquemos el método de Gauss para obtener una versión equivalente más simple.

¿Qué podemos concluir?

Solución 4.1.

1. La tercera ecuación ha desaparecido, es combinación lineal de las demás: \(E_3 = 2E_1 + E_2\).
2. La última fila es completa de 0s. ¿Qué habría pasado si no fuera así?.

4.3 Estudio de sistemas de ecuaciones lineales

Comenzamos por expresar formalmente las operaciones permitidas entre las ecuaciones de un sistema: aquellas que no cambian el conjunto de soluciones y que, a la postre, nos servirán para conocer el carácter del sistema.

Definición 4.4 (Combinación lineal) Decimos que una fila \(F\) de una matriz es combinación lineal de otras filas \(F_1, \ldots, F_k\) si existen escalares \(\lambda_1, \ldots, \lambda_k\in\mathbb{K}\) tales que \[F = \lambda_1 F_1 + \cdots + \lambda_k F_k\]

Un conjunto de filas es linealmente dependiente si alguna de ellas se puede poner como combinación lineal del resto. Se llamará linealmente independiente si no es así.

Definición 4.5 (Rango de una matriz) Llamamos rango de la matriz \(A\) al tamaño del mayor conjunto de filas de ella que es linealmente independiente. Lo denotaremos por \(\mathrm{Rg}(A)\).

Teorema 4.1 (Rouché-Frobenius) Sea \(A\) la matriz de coeficientes de un sistema de ecuaciones lineales, \(A^*\) su matriz ampliada y sea \(n\) el número de incógnitas del sistema. Entonces:

• El sistema es incompatible si, y sólo si, \(\mathrm{Rg}(A) < \mathrm{Rg}(A^*)\).
• El sistema es compatible determinado si, y sólo si, \(\mathrm{Rg}(A) = \mathrm{Rg}(A^*) = n\).
• El sistema es compatible indeterminado si, y sólo si, \(\mathrm{Rg}(A) = \mathrm{Rg}(A^*) < n\).

Observación 4.3. En el caso de un sistema compatible indeterminado, el número de incógnitas libres (a las que asignamos parámetros) es \(n − \mathrm{Rg}(A)\).

Ejercicio 4.2 Estudiar el sistema \[\left\{\begin{array}{rrrcr} x & & -z & = & 3\\ & y & + 2z & = & 1\\ 2x & + y & & = & 7\\ \end{array}\right.\] y encontrar su conjunto de soluciones.

4.4 Matrices y sistemas escalonados

La idea de esta sección es mostrar cómo reducir la complejidad del estudio de un sistema de ecuaciones lineales mediante su transformación en un sistema equivalente más sencillo.

Definición 4.6 (Forma escalonada) Una matriz \(m\times n\) está en forma escalonada por filas si verifica:

1. Si hay filas que solo contienen ceros, están en la parte inferior de la matriz.

2. En cada fila, al leer de izquierda a derecha, la primera entrada distinta de cero, llamada entrada principal o pivote de su fila, es un 1.

3. Si las filas \(i\) e \(i+1\) no constan solo de ceros, entonces la entrada principal de la fila \(i+1\) está a la derecha de la entrada principal de la fila \(i\).

Una matriz escalonada por filas se dice que está en forma escalonada reducida por filas si además verifica:

4. Si una columna contiene la entrada principal de alguna fila, entonces el resto de los elementos de esa columna son iguales a cero.

Ejemplo 4.1 Estas matrices están en forma escalonada por filas: \[ A = \begin{pmatrix} \color{ocre}{1}&2&3&5&4\\ 0&0&\color{ocre}{1}&7&2\\ 0&0&0&\color{ocre}{1}&0 \end{pmatrix}, \ B = \begin{pmatrix} \color{ocre}{1}&0&0&4\\ 0&\color{ocre}{1}&0&5\\ 0&0&\color{ocre}{1}&2 \end{pmatrix}, \ C= \begin{pmatrix} \color{ocre}{1}&2&0&3&0\\ 0&0&\color{ocre}{1}&4&0\\ 0&0&0&0&\color{ocre}{1}\\ 0&0&0&0&0\\ 0&0&0&0&0 \end{pmatrix}\]

Teorema 4.2 Toda matriz se puede transformar, usando las operaciones elementales, en una única matriz escalonada reducida por filas.

Definición 4.7 (Matrices equivalentes por filas) Dos matrices \(A\) y \(B\) se llaman equivalentes por filas si podemos pasar de una a otra mediante las transformaciones elementales por filas ya vistas. Se denota por \(A \approx_f B\).

Este teorema se puede reescribir de varias formas:

1. “Toda matriz es equivalente por filas a una única matriz escalonada reducida por filas”.
2. “Todo sistema de ecuaciones lineales se puede transformar en uno equivalente escalonado reducido usando únicamente transformaciones elementales”.

El Teorema 4.2 da lugar al conocido como método de Gauss-Jordan (Algoritmo 1).

Carl Friedrich Gauss (Gauß en alemán) (Braunschweig, 30 de abril de 1777 - Gotinga, 23 de febrero de 1855) fue un prolífico matemático, al que llegaron a llamar el príncipe de las Matemáticas.

\begin{algorithm} \caption{Método de Gauss-Jordan} \begin{algorithmic} \REQUIRE $A$ y $\symbfit{b}$, matriz de coeficientes y segundo término de un sistema de ecuaciones. \ENSURE $\symbfit{x}$, la solución del sistema de ecuaciones. \STATE Formar la matriz aumentada $(A \mid \symbfit{b} )$ \STATE Mediante operaciones elementales de filas, transformar la matriz aumentada $(A \mid \symbfit{b})$ a su forma escalonada reducida por filas $(U \mid \symbfit{c})$. \STATE En cada fila distinta de cero de la matriz $(U \mid \symbfit{c})$, despejar la incógnita correspondiente a la entrada principal de la fila, teniendo en cuenta las posibles variables libres. \end{algorithmic} \end{algorithm}

Resolver el sistema \[\left\{\begin{array}{rrrrcr} &2y&+3z&-4t&=&1\\ &&2z&+3t&=&4\\ 2x&+2y&-5z&+2t&=&4\\ 2x&&-6z&+9t&=&7 \end{array}\right.\] mediante Gauss-Jordan.

Sea \(I_n\) la matriz identidad en \(\mathcal{M}_n(\mathbb{K})\), \(1 \leq i, j \leq n\) y \(c\in\mathbb{K}\). Definimos las siguientes matrices elementales fila:

• \(\mathcal{E}_{ij}\), como la matriz obtenida aplicando \(F_i \leftrightarrow F_j\) a la matriz \(I_n\).
• \(\mathcal{E}_i(c)\) con \(c\ne 0\), como la matriz obtenida aplicando \(F_i \leftarrow cF_i\) a la matriz \(I_n\).
• \(\mathcal{E}_{ij}(c)\), como la matriz obtenida aplicando \(F_i \leftarrow F_i + cF_j\) a la matriz \(I_n\).

Multiplicar una matriz \(A\) por la izquierda por una matriz elemental fila equivale a realizar en \(A\) la transformación elemental correspondiente.

\[\mathcal{E}_{12} = \begin{pmatrix} 0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1\end{pmatrix},\quad \mathcal{E}_{2}(3) = \begin{pmatrix} 1 & 0 & 0\\ 0 & 3 & 0\\ 0 & 0 & 1\end{pmatrix},\quad \mathcal{E}_{12}(-5) = \begin{pmatrix} 1 & -5 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{pmatrix}\]

Nos podemos preguntar si estas matrices son regulares y si sus inversas tienen alguna propiedad.

Las matrices elementales fila son regulares y verifican: \[\mathcal{E}_{ij}^{-1} = \mathcal{E}_{ji},\qquad\mathcal{E}_{i}(c)^{-1} = \mathcal{E}_{i}(c^{-1}),\qquad\mathcal{E}_{ij}(c)^{-1} = \mathcal{E}_{ij}(-c)\]

4.5 Equivalencia de matrices

Todos los conceptos y resultados que hemos explicado por filas se pueden definir también por columnas:

• Transformaciones elementales por columnas.
• Matrices escalonadas (reducidas) por columnas.
• Matrices equivalentes por columnas (\(A\approx_c B\)).
• Matrices elementales columna.
• \(\ldots\).

Diremos que dos matrices \(A\) y \(B\) son equivalentes si se puede transformar una en la otra mediante transformaciones elementales por filas y columnas.

Hemos de notar que realizar operaciones por filas equivale a multiplicar por una matriz por la izquierda, mientras que operar por columnas equivale a multiplicar por una matriz a la derecha.

Sean \(A,B\in\mathcal{M}_{m,n}(\mathbb{K})\).

• \(A\) y \(B\) son equivalentes por filas si, y sólo si, existe \(P\in\mathcal{M}_m(\mathbb{K})\) regular tal que \(B = PA\).
• \(A\) y \(B\) son equivalentes por columnas si, y sólo si, existe \(Q\in\mathcal{M}_n(\mathbb{K})\) regular tal que \(B = AQ\).
• \(A\) y \(B\) son equivalentes si, y sólo si, existen \(P\in\mathcal{M}_m(\mathbb{K}), Q\in\mathcal{M}_n(\mathbb{K})\) regulares tales que \(B = PAQ\).

La idea de la demostración (para la equivalencia por filas) es como sigue:

• Para que \(A\approx_f B\), se necesita que haya una serie de transformaciones elementales fila que aplicar a \(A\) para obtener \(B\): \[A\rightarrow A_1\rightarrow A_2\rightarrow\cdots\rightarrow A_p = B\]
• Cada transformación elemental fila lleva asociada una matriz elemental fila, de forma que multiplicar por la izquierda por ella es hacer la transformación. Si llamamos \(T_k\) a la matriz elemental fila aplicada en el paso \(k\), nos queda: \[A\,\rightarrow\,T_1A\,\rightarrow\,T_2T_1A\,\rightarrow\,\cdots \,\rightarrow\,T_pT_{p-1}\cdots T_2T_1A = B\]
• Basta llamar \(P = T_pT_{p-1}\cdots T_2T_1\).

¿Cómo podemos encontrar las \(P\) y \(Q\) del teorema que hacen que \(B = PAQ\)?

Realizar transformaciones elementales fila y columna para obtener \(B\), almacenando los cambios que se han ido haciendo: \[ \begin{array}{c|c} A & I_m \\ \hline I_n & \end{array} \longrightarrow \begin{array}{c|c} B & P \\ \hline Q & \end{array} \]

Encontrar una matriz triangular superior \(U\) equivalente a \[A = \begin{pmatrix} 2 & -1 & 3 \\ -4 & 2 & -6 \\ -2 & 2 & 4 \end{pmatrix},\] así como las matrices \(P, Q\) tales que \(U = PAQ\).

4.6 Aplicación al cálculo del rango y de la inversa

Veamos cómo aplicar lo estudiado a dos situaciones que, generalmente, se han venido estudiando mediante determinantes: el cálculo del rango y de la inversa de una matriz dada.

Una matriz \(A\in\mathcal{M}_{m,n}(\mathbb{K})\) tiene rango \(r\) si, y sólo si, es equivalente a la matriz \[\left(\begin{array}{c|c} I_r & 0 \\ \hline 0 & 0 \end{array}\right)\]

Por tanto, nos bastará con aplicar transformaciones elementales por filas y por columnas a la matriz \(A\) para hallar una matriz como la del teorema, y así determinar su rango.

Calcular el rango de las matrices \[ \begin{pmatrix} 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & 2 & -3 \\ 0 & 1 & 2 & -1 & 0 \\ \end{pmatrix}\quad\text{y}\quad \begin{pmatrix} 1 & 2 & 1 & 3 & 1\\ 0 & 1 & 2 & 3 & 0\\ 1 & 3 & 3 & 6 & 2\\ 2 & 6 & 6 & 12 & 3 \end{pmatrix} \]

Un caso particular de lo anterior sirve para calcular la matriz inversa de una matriz regular \(A\).

\(A\in\mathcal{M}_n(\mathbb{K})\) es regular si, y sólo si, \(A\) es equivalente por filas (o por columnas) a \(I_n\).

Como \(A\approx_f I_n\), existe \(P\) regular tal que \(I_n = PA\). Para calcular \(P\), podremos usar el método de Gauss-Jordan para transformar \[(A | I_n) \approx_f (I_n | P)\]

Pero \(I_n = PA\) significa que \(P = A^{-1}\).

Calcular la inversa de \[\begin{pmatrix} 1 & 2 & 3 & 4 \\ -2 & 3 & 4 & 5 \\ 3 & -4 & 5 & 6 \\ -4 & 5 & 6 & 7 \\ \end{pmatrix}\]

4.7 Métodos de factorización

Recordemos que:

1. Los sistemas de ecuaciones fáciles de resolver son los diagonales, y los triangulares superiores e inferiores.
2. Las matrices escalonadas son generalizaciones del concepto de matrices triangulares.
3. Toda matriz es equivalente por filas a una matriz escalonada por filas (única si es escalonada reducida).

¿Podemos usar estos hechos para transformar un sistema de ecuaciones \(A\symbfit{x} = \symbfit{b}\) en otro sistema equivalente que tenga estructura escalonada?

Esto lo hacía el método de Gauss. Ahora veremos otro tipo de métodos directos de resolución de sistemas de ecuaciones lineales.

4.7.1 Factorización LU

Los métodos de factorización buscan descomponer la matriz de coeficientes \(A\) como producto de dos o más matrices cuya estructura sea simple.

En particular, el método LU tiene como objetivo expresar \(A = LU\), donde \(L\) es una matriz triangular inferior y \(U\) es una matriz escalonada por filas.

La idea es reescribir el sistema \(A\symbfit{x} = \symbfit{b}\) como \(LU\symbfit{x} = \symbfit{b}\). Entonces:

1. Construimos dos sistemas de ecuaciones: \[\left\{\begin{array}{l} U\symbfit{x} = \symbfit{y}\\ L\symbfit{y} = \symbfit{b} \end{array}\right.\]

Es fácil ver que la solución de \(A\symbfit{x}=\symbfit{b}\) y de \(U\symbfit{x} = \symbfit{y}\) es la misma. 2. Resolvemos \(L\symbfit{y} = \symbfit{b}\) mediante sustitución progresiva (aquí la \(\symbfit{y}\) es un vector incógnita cuyo valor encontramos). 3. Con el valor calculado de \(\symbfit{y}\), resolvemos \(U\symbfit{x}=\symbfit{y}\) mediante sustitución regresiva: ahora el vector \(\symbfit{y}\) es el término independiente en este nuevo sistema. Por tanto, hallamos \(\symbfit{x}\), que es la solución de \(A\symbfit{x}=\symbfit{b}\).

Como hemos visto antes, hallar matrices escalonadas está muy relacionado con el método de Gauss. Así, el procedimiento para encontrar \(L\) y \(U\) a partir de \(A\) se basa en dicho método, usando únicamente transformaciones de tipo 3 (en particular, ).

\[\begin{array}{ccc} \text{Transformación} & L & U \\ \hline & I_m & A \\ T_1 & T_1^{-1} & T_1A \\ T_2 & T_2^{-1}T_1^{-1} & T_2T_1A \\ \vdots & \vdots & \vdots \\ & L & U\\\hline \end{array}\]

Para nosotros, la principal ventaja de los métodos de factorización es que la factorización se puede hacer con la matriz de coeficientes, antes incluso de conocer el vector de términos independientes.

Además, una vez hecha la factorización, se puede usar para resolver tantos sistemas (con la misma matriz de coeficientes) como se desee, sin cálculos adicionales.

Resolvamos \[\left\{\begin{array}{rrrl} x_1 & +2x_2 & +3x_3 & = 2\\ 3x_1 & +2x_2 &+x_3 & = 2\\ x_1 &+x_2 &+x_3 &=1 \end{array}\right.\] mediante factorización LU.

Paso 1: Calcular la factorización.

Paso 2: Resolver \(L\symbfit{y} = \symbfit{b}\).

Paso 3: Resolver \(U\symbfit{x} = \symbfit{y}\).

4.7.2 Condiciones para la factorización

¿Toda matriz admite factorización LU?

Podemos comprobar que no podemos aplicar el método anterior a la matriz \[\begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix}\] aunque esta matriz es cuadrada y regular.

Hay que buscar algún criterio más para garantizar la existencia de la factorización.

Sea \(A = (a_{i,j})\in\mathcal{M}_n(\mathbb{K})\). Denominamos menor principal de orden \(k\), con \(1\leq k\leq n\), al determinante siguiente: \[A_k = \left|\begin{array}{ccc} a_{1,1} & \ldots & a_{1,k}\\ \vdots & \ddots & \vdots\\ a_{k,1} & \ldots & a_{k,k} \end{array}\right|\]

Una matriz regular \(A \in\mathcal{M}_n(\mathbb{K})\) admite factorización LU si, y sólo si \(A_k\ne 0\) para todo \(1 \leq k \leq n\).

Para algunas aplicaciones de la factorización LU interesa que los elementos de la diagonal de \(L\) sean todos unos. Denotaremos esta factorización por L\(_1\)U.

En el ejemplo anterior, habíamos hecho realmente la factorización L\(_1\)U: \[\begin{pmatrix}1 & 2 & 3\\ 3 & 2 & 1 \\ 1 & 1 & 1\end{pmatrix} = \begin{pmatrix}\textcolor{ocre}{1} & 0 & 0 \\ 3 & \textcolor{ocre}{1} & 0 \\ 1 & \frac{1}{4} & \textcolor{ocre}{1}\end{pmatrix}\cdot \begin{pmatrix}1 & 2 & 3\\ 0 & -4 & -8 \\ 0 & 0 & 0\end{pmatrix}\]

Esta factorización, si existe, es única.

¿Podremos extender todo esto a matrices donde no se dé la condición sobre los menores principales? ¿Puede bastar que la matriz de coeficientes sea regular?

4.7.3 Factorización LU con permutación

Recordemos que podemos conseguir sistemas equivalentes aplicando transformaciones de tipo 1 (intercambiar ecuaciones de orden).

Llamamos factorización LU con permutación, también conocida como factorización LU con pivotes parciales, de una matriz \(A\) a una igualdad de tipo \(PA = LU\) donde \(P\) es una matriz de permutación (producto de \(\mathcal{E}_{ij}\)), \(L\) es una matriz triangular inferior y \(U\) es una matriz triangular superior. Si los elementos de la diagonal de \(L\) son todos unos, decimos que es una factorización L\(_1\)U con permutación.

Toda matriz cuadrada \(A\) admite una factorización LU con permutación.

¿Cómo usamos la factorización LU con permutación?

Supongamos que tenemos \(P, L, U\) tales que \(PA = LU\). Entonces:

1. Partiendo de \(A\symbfit{x} = \symbfit{b}\), multiplicamos por \(P\) por la izquierda ambos términos: \(PA\symbfit{x} = P\symbfit{b}\).
2. Como \(PA = LU\), realmente tenemos \(LU\symbfit{x} = P\symbfit{b}\).
3. Al igual que la factorización LU:
   1. Resolvemos \(L\symbfit{y} = P\symbfit{b}\): hallamos \(y\) mediante sustitución progresiva.
   2. Finalmente, resolvemos \(U\symbfit{x} = \symbfit{y}\) mediante sustitución regresiva.

Esta estrategia se puede usar para matrices no cuadradas.

Sabemos cómo usar la factorización, pero ¿cómo la calculamos?

Se hace igual que la normal, llevando cuenta de las permutaciones.

\[\begin{array}{cccc} \text{Transformación} & P & L & U \\ \hline &I_m & I_m & A \\ T_1 & P_1 & T_1^{-1} & T_1A \\ T_2 & P_2P_1& T_2^{-1}T_1^{-1} & T_2T_1A \\ \vdots & \vdots & \vdots \\ & P & L & U\\\hline \end{array}\]

Resolvamos \[\left\{\begin{array}{rrrrl} x_1 &+ x_2 &+ 2 x_3 &+ x_4 &= 3\\ 2x_1 &+2x_2 &+x_3 &+x_4 &= 9\\ 3x_1 &+2x_2 &+x_3 &+x_4 &= 6 \end{array}\right.\] mediante factorización LU.

Paso 1: Calcular la factorización.

Recordar que hace falta sumar la identidad para conseguir \[L = \begin{pmatrix}\textcolor{ocre}{1} & 0 & 0 \\ 3 & \textcolor{ocre}{1} & 0 \\ 2 & 0 & \textcolor{ocre}{1}\end{pmatrix}\]

Con la \(P\) hallada, podemos obtener el sistema equivalente que podrá ser resuelto usando la factorización LU: \[A\symbfit{x} = \symbfit{b} \Leftrightarrow PA\symbfit{x}=P\symbfit{b} \Leftrightarrow LU\symbfit{x} = P\symbfit{b} = \begin{pmatrix}3\\6\\9\end{pmatrix}\]

Paso 2: Resolver \(L\symbfit{y} = P\symbfit{b}\).

\[L\symbfit{y} = \symbfit{b} \Longleftrightarrow \begin{pmatrix}1 & 0 & 0 \\ 3 & 1 & 0 \\ 2 & 0 & 1\end{pmatrix} \begin{pmatrix}y_1\\y_2\\y_3\end{pmatrix} = \begin{pmatrix}3\\6\\9\end{pmatrix} \Longleftrightarrow \left\{\begin{array}{rrrl} y_1 & & & =3 \\ 3y_1 &+y_2 & & = 6 \\ 2y_1 & &+y_3 &= 9 \end{array}\right.\]

Paso 3: Resolver \(U\symbfit{x} = \symbfit{y}\). \[\begin{align*} U\symbfit{x} = \symbfit{y} & \Leftrightarrow \begin{pmatrix}1&1&2&1\\ 0&-1&-5&-2 \\ 0&0&-3&-1\end{pmatrix} \begin{pmatrix}x_1\\x_2\\x_3\\x_4\end{pmatrix} = \begin{pmatrix}3\\-3\\3\end{pmatrix}\Longleftrightarrow \\ & \Leftrightarrow\left\{\begin{array}{rrrrr} x_1 & +x_2 &+2x_3 &+x_4 & = 3 \\ & −x_2 &−5x_3 & -2x_4 &= −3 \\ & & -3x_3 & - x_4&= 3 \end{array}\right. \end{align*}\]

Usar la factorización LU de \(A\) para resolver el sistema \(A^2\symbfit{x} = \symbfit{b}\) (sin calcular \(A^2\)), donde \[A = \begin{pmatrix}−1 & 0 & 1 \\3 & 1 & 0\\-2 & 0 & 4\end{pmatrix}\qquad b = \begin{pmatrix}20\\40\\80\end{pmatrix}\]

4.8 Métodos numéricos iterativos

4.8.1 ¿Qué problemas hay con los métodos directos?

Los métodos directos resultan, en general, inservibles para \(n>50\) incógnitas. El problema es que los errores que se van cometiendo se propagan, y al aumentar el número de operaciones llegan a invalidar el resultado obtenido.

Otro problema es que los métodos directos necesitan almacenar la matriz \(A\) en memoria (pensar la memoria necesaria para almacenar una matriz de n=2500 a 8 bytes por coeficiente en doble precisión).

Los grandes sistemas de ecuaciones que surgen en la práctica, tienen la matriz \(A\) dispersa (muchos coeficientes son cero) y, aunque existen métodos directos especiales, usualmente se resuelven por métodos iterativos.

Los métodos iterativos tienen la ventaja de no propagar el error.

4.8.2 ¿En qué se basan los métodos iterativos?

En vez de determinar la solución de un sistema de forma directa, construiremos una sucesión de aproximaciones a la solución.

Es decir, en lugar de calcular directamente \(\symbfit{x}\) (la solución del sistema \(A\symbfit{x}=\symbfit{b}\)), construimos en un primer paso un vector \(\symbfit{x^{(1)}}\); a partir de él, determinamos un \(\symbfit{x^{(2)}}\); etc., con la esperanza de que, en el “infinito”, los términos de esta sucesión se aproximen muy bien a \(\symbfit{x}\).

Ya hemos dicho que esto tiene ventajas (menor uso de memoria, mejor representación de las matrices, menor propagación de errores).

Pero…

1. ¿Cómo se usan las sucesiones de vectores? ¿Qué relación tienen con las sucesiones numéricas?
2. Una vez que sabemos manejar sucesiones de vectores, ¿cómo construimos la sucesión que necesitamos para resolver un sistema de ecuaciones?

4.8.2.1 Breves nociones de normas y sucesiones

Sea \(A\in\mathcal{M}_{m,n}(\mathbb{K})\) una matriz y \(\symbfit{v}\in\mathcal{M}_{m,1}(\mathbb{K})\cong\mathbb{K}^m\) un vector (matriz columna). Definimos:

• La norma del máximo de \(\symbfit{v}\) como \(\|\symbfit{v}\|_{\infty} = \max\{|v_1|,|v_2|,\ldots,|v_m|\}\).
• La norma del máximo de la matriz \(A\) como \(\|A\|_{\infty} = \max_i\sum_j|a_{i,j}|\).

A esta norma también se le denomina norma infinito.

\[\begin{align*} \left\|\begin{pmatrix}-3\\2\\-6\end{pmatrix}\right\|_{\infty} & = \max\{|-3|, |2|, |-6|\} = 6\\ \left\|\begin{pmatrix}1 & -4 & 3\\-1 & 2 & -3\end{pmatrix}\right\|_{\infty} & = \max\{|1| + |-4| + |3|, |-1| + |2| + |-3|\} = \\ & = \max\{8, 6\} = 8 \end{align*}\]

Con la ayuda de una norma, podemos medir distancias entre objetos.

Podemos definir la convergencia de una sucesión vectorial en términos de normas:

Consideremos una sucesión de vectores \(\{\symbfit{x^{(1)}}, \symbfit{x^{(2)}}, \ldots,\symbfit{x^{(k)}},\ldots\}\) de vectores columna en \(\mathbb{K}^n\). Diremos que \(\symbfit{x}\in\mathbb{K}^n\) es el límite de la sucesión \(\{\symbfit{x^{(k)}}\}\)} si \[\lim_{k\to\infty}\|\symbfit{x^{(k)}} - \symbfit{x}\|_{\infty} = 0\]

En este caso, diremos que \(\{\symbfit{x^{(k)}}\}\) es **convergente y denotaremos \(\lim \symbfit{x^{(k)}} = \symbfit{x}\) o \(\symbfit{x^{(k)}}\longrightarrow \symbfit{x}\).

Toda esta discusión se puede escribir con cualquier norma matricial (hay muchas). Por ejemplo, la \(\|\cdot\|_1\), que también es sencilla de calcular: \(\|A\|_1 = \|A^T\|_{\infty}\).

4.8.2.2 Formulación de punto fijo

Sea \(F\colon \mathbb{K}^n\to\mathbb{K}^n\) una función (campo vectorial). Un vector \(\symbfit{x}\in\mathbb{K}^n\) se denomina punto fijo de \(F\)} si \(F(\symbfit{x}) = \symbfit{x}\).

El punto \(x = 0\) es un punto fijo de la función seno, puesto que \(\sin 0 = 0\).

Vamos a reescribir el problema \(A\symbfit{x}=\symbfit{b}\) como un problema de punto fijo con la misma solución. Es decir, transformaremos \(A\symbfit{x}=\symbfit{b}\) en otro problema del tipo \(\symbfit{x} = F(\symbfit{x})\), de forma que \(\symbfit{x}\) es solución de \(A\symbfit{x}=\symbfit{b}\) si, y sólo si, \(\symbfit{x} = F(\symbfit{x})\).

¿Hay alguna forma de determinar los puntos fijos de \(F\)? ¿Hay que imponer alguna condición sobre \(F\)?

Sea \(F\colon \mathbb{K}^n\to\mathbb{K}^n\) una función tal que existe \(L<1\) para el que \(\|F(\symbfit{a}) - F(\symbfit{b})\| \le L \|\symbfit{a}-\symbfit{b}\|\) para todo \(\symbfit{a},\symbfit{b}\in\mathbb{K}^n\). Entonces:

1. El campo \(F\) tiene un único punto fijo \(\symbfit{x}\).
2. Si construimos una sucesión \(\{\symbfit{x^{(k)}}\}\) de la siguiente forma: \[\symbfit{x^{(0)}}\in\mathbb{K}^n\text{ cualquiera}, \symbfit{x^{(1)}} = F(\symbfit{x^{(0)}}), \symbfit{x^{(2)}} = F(\symbfit{x^{(1)}}), \symbfit{x^{(3)}} = F(\symbfit{x^{(2)}}), \ldots\] entonces \(\{\symbfit{x^{(k)}}\}\longrightarrow x\).

A una función \(F\) con la propiedad del enunciado, la llamaremos contractiva.

Calcular el punto fijo de la función seno en \(\left(-\pi, \pi\right)\) con la calculadora.

Definir un método iterativo para resolver sistemas de ecuaciones lineales es, esencialmente, definir la función \(F\) de forma adecuada para que la solución de \(A\symbfit{x}=\symbfit{b}\) sea punto fijo de \(F\), es decir, para que \(F(\symbfit{x}) = \symbfit{x}\).

Nos centraremos en funciones \(F\) definidas de forma sencilla: \[F(\symbfit{x}) = T\symbfit{x} + \symbfit{c}\] donde \(T\) es una matriz del mismo tamaño que \(A\) y \(\symbfit{c}\) un vector columna del tamaño de \(\symbfit{b}\).

Entonces, el método iterativo vendrá dado de esta forma: \[\symbfit{x^{(0)}}\in\mathbb{K}^n\text{ cualquiera}, \symbfit{x^{(1)}} = T\symbfit{x^{(0)}}+\symbfit{c}, \symbfit{x^{(2)}} = T\symbfit{x^{(1)}}+\symbfit{c}, \symbfit{x^{(3)}} = T\symbfit{x^{(2)}}+\symbfit{c}, \ldots\]

¿Cuándo podremos asegurar que la solución de \(A\symbfit{x}=\symbfit{b}\) y el punto fijo de \(F(\symbfit{x}) = T\symbfit{x} + \symbfit{c}\) son el mismo?

Cuando esto sucede, el método iterativo se denomina consistente.

Sea \(A\in\mathcal{M}_n(\mathbb{K})\) una matriz regular. El método iterativo dado por \[\symbfit{x^{(0)}}\in\mathbb{K}^n\text{ cualquiera}, \symbfit{x^{(1)}} = T\symbfit{x^{(0)}}+\symbfit{c}, \symbfit{x^{(2)}} = T\symbfit{x^{(1)}}+\symbfit{c}, \symbfit{x^{(3)}} = T\symbfit{x^{(2)}}+\symbfit{c}, \ldots\] es consistente si, y sólo si, \(I-T\) es regular y \(\symbfit{c}=(I-T)A^{-1}\symbfit{b}.\)

Definiendo \(T\) y \(\symbfit{c}\) de distintas formas (consistentemente) aparecen distintos métodos. En concreto, estudiaremos el método de Jacobi y el método de Gauss-Seidel.

4.8.3 Método de Jacobi

Observemos que la matriz \(A\) se puede descomponer de la forma \(A = D+L+R\): \[\underbrace{ \begin{pmatrix} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \ldots & a_{nn} \end{pmatrix}}_{A} =\underbrace{ \left( \begin{array}{c@{}c@{}c@{}c} a_{11} & 0 & \ldots & 0 \\ 0 & a_{22} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & a_{nn} \end{array} \right)}_{D}+ \underbrace{ \left( \begin{array}{c@{}c@{}c@{}c} 0 & 0 & \ldots & 0 \\ a_{21} & 0 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \ldots & 0 \end{array} \right)}_{L}+ \underbrace{ \left( \begin{array}{c@{}c@{}c@{}c} 0 & a_{12} & \ldots & a_{1n} \\ 0 & 0 & \ldots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & 0 \end{array} \right).}_{R}\]

Podemos entonces reescribir el sistema de la siguiente forma \[A\symbfit{x}=\symbfit{b} \Leftrightarrow (D+L+R)\symbfit{x}=\symbfit{b} \Leftrightarrow D\symbfit{x}=-(L+R)\symbfit{x}+\symbfit{b}\Leftrightarrow \symbfit{x}=\underbrace{-D^{-1}(L+R)}_{T}\symbfit{x}+\underbrace{D^{-1}\symbfit{b}}_{\symbfit{c}}\]

Por tanto: \(\symbfit{c}=D^{-1}\symbfit{b}\) y la matriz del método de Jacobi es: \(T=-D^{-1}(L+R)\).

En la práctica, dado el sistema: \[\left. \begin{array}{cccccc} a_{1,1}x_{1} & +a_{1,2}x_{2} & +a_{1,3}x_{3} & \ldots & +a_{1,n}x_{n} & =b_{1} \\ a_{2,1}x_{1} & +a_{2,2}x_{2} & +a_{2,3}x_{3} & \ldots & +a_{2,n}x_{n} & =b_{2} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ a_{n,1}x_{1} & +a_{n,2}x_{2} & +a_{n,3}x_{3} & \ldots & +a_{n,n}x_{n} & =b_{n} \\ \end{array} \right\}\] despejamos la incógnita \(x_{i}\) de la ecuación i-ésima: \[\hspace*{-1cm}\left. \begin{array}{ll} x_{1} = &\frac{1}{a_{1,1}}\left( b_{1} -a_{1,2}x_{2} -a_{1,3}x_{3} - \ldots -a_{1,n}x_{n}\right) \\ x_{2} = &\frac{1}{a_{2,2}}\left( b_{2} -a_{2,1}x_{2} -a_{2,3}x_{3} - \ldots -a_{2,n}x_{n}\right) \\ \multicolumn{1}{c}{\vdots} & \multicolumn{1}{c}{\vdots} \\ x_{i} = &\frac{1}{a_{i,i}}\left( b_{i} -a_{i,1}x_{1} -\ldots -a_{i,i-1}x_{i-1}-a_{i,i+1}x_{i+1}-\ldots-a_{i,n}x_{n}\right)\\ \multicolumn{1}{c}{\vdots}&\multicolumn{1}{c}{\vdots}\\ x_{n} = &\frac{1}{a_{n,n}}\left( b_{n} -a_{n,1}x_{1} -a_{n,2}x_{2} - \ldots -a_{n,n-1}x_{n-1}\right) \\ \end{array} \right\}\]

El método iterativo es:

\[\left. \begin{array}{ll} x_{1}^{(k+1)} = &\frac{1}{a_{1,1}}\left( b_{1} -a_{1,2}x_{2}^{(k)} -a_{1,3}x_{3}^{(k)} - \ldots -a_{1,n}x_{n}^{(k)}\right) \\ x_{2}^{(k+1)} = &\frac{1}{a_{2,2}}\left( b_{2} -a_{2,1}x_{1}^{(k)} -a_{2,3}x_{3}^{(k)} - \ldots -a_{2,n}x_{n}^{(k)}\right) \\ \multicolumn{1}{c}{\vdots} & \multicolumn{1}{c}{\vdots} \\ x_{i}^{(k+1)} = &\frac{1}{a_{i,i}}\left( b_{i} -a_{i,1}x_{1}^{(k)} -\ldots -a_{i,i-1}x_{i-1}^{(k)}-a_{i,i+1}x_{i+1}^{(k)}-\ldots-a_{i,n}x_{n}^{(k)}\right)\\ \multicolumn{1}{c}{\vdots}&\multicolumn{1}{c}{\vdots}\\ x_{n}^{(k+1)} = &\frac{1}{a_{n,n}}\left( b_{n} -a_{n,1}x_{1}^{(k)} -a_{n,2}x_{2}^{(k)} - \ldots -a_{n,n-1}x_{n-1}^{(k)}\right) \\ \end{array} \right\}\]

Resolver el siguiente sistema utilizando el método de Jacobi, dando tres iteraciones y partiendo del vector nulo: \[ \left\{ \begin{array}{cccccc} 10x & -y & +2z & & = & 6 \\ -x & +11y & -z & +3t & = & 25 \\ 2x & -y &+10z & -t & = & -11 \\ & 3y & -z & +8t & = & 15 \end{array} \right. \]

Partiremos del vector inicial \[ \begin{pmatrix} x^{(0)} \\ y^{(0)} \\ z^{(0)} \\ t^{(0)} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \] y realizaremos tres iteraciones.

En el caso del método de Jacobi iteramos usando las igualdades \[ \left\{ \begin{array}{ccccccc} x^{(k+1)} & = & & \frac{1}{10}y^{(k)} & -\frac{2}{10}z^{(k)}& &+\frac{6}{10} \\ \noalign{\vspace{1mm}} y^{(k+1)} & = & \frac{1}{11} x^{(k)} & & +\frac{1}{11}z^{(k)} & -\frac{3}{11}t^{(k)} & + \frac{25}{11} \\ \noalign{\vspace{1mm}} z^{(k+1)} & = & -\frac{2}{10} x^{(k)} & + \frac{1}{10}y^{(k)} & & +\frac{1}{10} t^{(k)} & -\frac{11}{10} \\ \noalign{\vspace{1mm}} t^{(k+1)} & = & & -\frac{3}{8} y^{(k)} & + \frac{1}{8}z^{(k)} & & +\frac{15}{8} \end{array} \right. \] para obtener: \[ \begin{array}{l} \left( \begin{array}{c} x^{(0)} \\ y^{(0)} \\ z^{(0)} \\ t^{(0)} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \end{array} \right) \quad \Rightarrow \quad \left( \begin{array}{c} x^{(1)} \\ y^{(1)} \\ z^{(1)} \\ t^{(1)} \end{array} \right) = \left( \begin{array}{c} 0.6000 \\ 2.2727 \\ -1.100 \\ 1.8750 \end{array} \right) \\ \quad \Rightarrow \quad \left( \begin{array}{c} x^{(2)} \\ y^{(2)} \\ z^{(2)} \\ t^{(2)} \end{array} \right) = \left( \begin{array}{c} 1.0473 \\ 1.7159 \\ -0.8052 \\ 0.8852 \end{array} \right) \quad \Rightarrow \quad \left( \begin{array}{c} x^{(3)} \\ y^{(3)} \\ z^{(3)} \\ t^{(3)} \end{array} \right) = \left( \begin{array}{c} 0.9326 \\ 2.0522 \\ -1.0493 \\ 1.1309 \end{array} \right). \end{array} \]

4.8.4 Método de Gauss-Seidel

En el ejemplo de Jacobi anterior, podemos ver cómo se calcula \(y^{(1)}\): \[y^{(1)} = \frac{1}{11} x^{(0)} +\frac{1}{11}z^{(0)} -\frac{3}{11}t^{(0)} + \frac{25}{11}\]

Si nos fijamos en que justo antes habíamos calculado \(x^{(1)}\), podemos pensar que este valor está más próximo al valor real de la componente \(x\) de la solución real. Así que para calcular \(y^{(1)}\) será mejor emplear \(x^{(1)}\) en lugar de \(x^{(0)}\).

Por tanto, la expresión quedaría \[y^{(1)} = \frac{1}{11} \fcolorbox{ocre}{ocre!15}{\color{ocre}\ensuremath{x^{(1)}}} +\frac{1}{11}z^{(0)} -\frac{3}{11}t^{(0)} + \frac{25}{11}\]

Siguiendo este razonamiento, podríamos calcular \(z^{(1)}\) usando \(x^{(1)}\) e \(y^{(1)}\), y así con todas las iteraciones, usando siempre el valor más actualizado posible para cada variable.

Ésta es la expresión del método: \[\left. \begin{array}{cc} x_{1}^{(k+1)} = &\frac{1}{a_{1,1}}\left( b_{1} -a_{1,2}x_{2}^{(k)} -a_{1,3}x_{3}^{(k)} - \ldots -a_{1,n}x_{n}^{(k)}\right) \\ x_{2}^{(k+1)} = &\frac{1}{a_{2,2}}\left( b_{2} -a_{2,1}\textcolor{ocre}{x_{1}^{(k+1)}} -a_{2,3}x_{3}^{(k)} - \ldots -a_{2,n}x_{n}^{(k)}\right) \\ \multicolumn{1}{c}{\vdots} & \multicolumn{1}{c}{\vdots} \\ x_{i}^{(k+1)} = &\frac{1}{a_{i,i}}\left( b_{i} -a_{i,1}\textcolor{ocre}{x_{1}^{(k+1)}} -\ldots -a_{i,i-1}\textcolor{ocre}{x_{i-1}^{(k+1)}}-a_{i,i+1}x_{i+1}^{(k)}-\ldots-a_{i,n}x_{n}^{(k)}\right)\\ \multicolumn{1}{c}{\vdots}&\multicolumn{1}{c}{\vdots}\\ x_{n}^{(k+1)} = &\frac{1}{a_{n,n}}\left( b_{n} -a_{n,1}\textcolor{ocre}{x_{1}^{(k+1)}} -a_{n,2}\textcolor{ocre}{x_{2}^{(k+1)}} - \ldots -a_{n,n-1}\textcolor{ocre}{x_{n-1}^{(k+1)}}\right) \\ \end{array} \right\}\]

• Al calcular cada incógnita usamos la estimación más reciente del resto de variables, lo cual proporciona mayor velocidad al método.
• La expresión se puede deducir fácilmente como se hizo con el método de Jacobi, despejando \(x_i\) de la ecuación \(i\)-ésima, pero haciendo las actualizaciones mencionadas.
• En la expresión matricial del método de Gauss-Seidel, \(\symbfit{c}=(D+L)^{-1}\symbfit{b}\) y la matriz del método es \(T=-(D+L)^{-1}R\).

Resolver el siguiente sistema utilizando el método de Gauss-Seidel, dando tres iteraciones y partiendo del vector nulo: \[ \left\{ \begin{array}{cccccc} 10x & -y & +2z & & = & 6 \\ -x & +11y & -z & +3t & = & 25 \\ 2x & -y &+10z & -t & = & -11 \\ & 3y & -z & +8t & = & 15 \end{array} \right. \]

\[ \left\{ \begin{array}{ccccccc} x^{(k+1)} & = & & \frac{1}{10}y^{(k)} & -\frac{2}{10}z^{(k)}& &+\frac{6}{10} \\ \noalign{\vspace{1mm}} y^{(k+1)} & = & \frac{1}{11} x^{(k+1)} & & +\frac{1}{11}z^{(k)} & -\frac{3}{11}t^{(k)} & + \frac{25}{11} \\ \noalign{\vspace{1mm}} z^{(k+1)} & = & -\frac{2}{10} x^{(k+1)} & + \frac{1}{10}y^{(k+1)} & & +\frac{1}{10} t^{(k)} & -\frac{11}{10} \\ \noalign{\vspace{1mm}} t^{(k+1)} & = & & -\frac{3}{8} y^{(k+1)} & + \frac{1}{8}z^{(k+1)} & & +\frac{15}{8} \end{array} \right. \]

Dar tres iteraciones de Gauss-Seidel.

\[ \begin{array}{l} \left( \begin{array}{c} x^{(0)} \\ y^{(0)} \\ z^{(0)} \\ t^{(0)} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \end{array} \right) \quad \Rightarrow \quad \left( \begin{array}{c} x^{(1)} \\ y^{(1)} \\ z^{(1)} \\ t^{(1)} \end{array} \right) = \left( \begin{array}{c} 0.6000 \\ 2.3272 \\ -0.9873\\ 0.8789 \end{array} \right) \\ \quad \Rightarrow \quad \left( \begin{array}{c} x^{(2)} \\ y^{(2)} \\ z^{(2)} \\ t^{(2)} \end{array} \right) = \left( \begin{array}{c} 1.0300 \\ 2.0370 \\ -1.0140 \\ 0.9844 \end{array} \right) \quad \Rightarrow \quad \left( \begin{array}{c} x^{(3)} \\ y^{(3)} \\ z^{(3)} \\ t^{(3)} \end{array} \right) = \left( \begin{array}{c} 1.0065 \\ 2.0036 \\ -1.0025 \\ 0.9983 \end{array} \right). \end{array} \]

4.8.5 Análisis de la convergencia y del error

¿Cuándo podremos asegurar que un método iterativo es convergente? ¿Podemos estimar el error que cometemos con las aproximaciones?

Recordemos que los métodos iterativos se construyen a partir de una matriz \(T\) y un vector \(c\) que, bajo ciertas condiciones, nos aseguran la consistencia: \[\symbfit{x^{(0)}}\text{ cualquiera},\quad \symbfit{x^{(1)}} = T\symbfit{x^{0)}}+\symbfit{c}, \ldots,\symbfit{x^{(k+1)}} = T\symbfit{x^{(k)}}+\symbfit{c},\ldots\]

Vamos a analizar el error que se comete en la aproximación \(k\)-ésima: la diferencia o distancia entre \(\symbfit{x^{(k)}}\) y la solución real \(\symbfit{x}\), es decir, \[e_k = \|\symbfit{x^{(k)}} - \symbfit{x}\|_{\infty}\]

Lo ideal es que este error \(e_k\) tienda a 0 (es la definición de método convergente).

Sea \(\{\symbfit{x^{(k)}}\}\) la sucesión de vectores construida usando un método iterativo consistente como el expresado en la diapositiva anterior, y sea \(\symbfit{x}\in\mathbb{K}^n\) el punto fijo (solución de \(A\symbfit{x}=\symbfit{b}\)). Entonces \[\symbfit{x^{(k)}} - \symbfit{x} = T^k \left(\symbfit{x^{(0)}} - \symbfit{x}\right)\]

Como consecuencia, \(\symbfit{x^{(k)}} \longrightarrow \symbfit{x}\) si, y sólo si, o bien \(\symbfit{x^{(0)}} = \symbfit{x}\) o \(T^k\) tiende a la matriz nula cuando \(k\longrightarrow \infty\).

Esto se debe a que: \[e_k = \|\symbfit{x^{(k)}} - \symbfit{x}\|_{\infty} = \|T^k \left(\symbfit{x^{(0)}} - \symbfit{x}\right)\|_{\infty}\le \left(\|T\|_{\infty}\right)^k\cdot\|\symbfit{x^{(0)}} - \symbfit{x}\|_{\infty}\] donde hemos usado una propiedad de las normas matriciales: \(\|XY\|_{\infty} \le \|X\|_{\infty}\cdot\|Y\|_{\infty}\) para todo par de matrices \(X,Y\) con dimensiones compatibles con su producto.

Sea \(\{\symbfit{x^{(k)}}\}\) la sucesión de vectores construida usando un método iterativo consistente como el expresado an la diapositiva anterior, y sea \(\symbfit{x}\in\mathbb{K}^n\) el punto fijo (solución de \(A\symbfit{x}=\symbfit{b}\)).

Si \(\|T\|_{\infty} < 1\) o \(\|T\|_1 < 1\), entonces \(\symbfit{x^{(k)}} \longrightarrow \symbfit{x}\).

Pero la matriz \(T\) no siempre es sencilla de calcular (especialmente en el método de Gauss-Seidel). ¿Hay algún criterio sobre la matriz de coeficientes \(A\) para asegurar la convergencia de estos métodos?

Si nos dijamos cómo se construye la matriz \(T\) en el método de Jacobi, podemos comprobar que \(\|T\|_{\infty} < 1\) si, y sólo si, para cada fila \(i=1,\ldots, n\), se tiene

Es decir, si, y sólo si, \(\displaystyle |a_{i,i}| > \sum_{1\le j\le n, j\ne i}|a_{i,j}|\) para todo \(i\).

Una matriz \(A=(a_{i,j})\in\mathcal{M}_n(\mathbb{K})\) se denomina de diagonal estrictamente dominante si verifica \[ |a_{i,i}| > \sum_{1\le j\le n, j\ne i}|a_{i,j}|\] para todo \(i\in\{1,\ldots,n\}\).

Sea \(A\symbfit{x}=\symbfit{b}\) un sistema de ecuaciones lineales con \(A\) una matriz regular. Si \(A\) es de diagonal estrictamente dominante, entonces los métodos de Jacobi y Gauss-Seidel proporcionan sucesiones convergentes a la solución del sistema, independientemente de la aproximación inicial \(\symbfit{x^{(0)}}\) usada.

¿Cómo sabemos cuántas iteraciones tenemos que dar por alguno de estos métodos? Lo usual es dejar de iterar cuando el error cometido en la iteración actual es menor que una cota de error tolerable, o cuando se alcanza un número de iteraciones prefijado.

Hasta ahora, sabíamos que \[e_k \le \left(\|T\|_{\infty}\right)^k\cdot\|\symbfit{x^{(0)}} - \symbfit{x}\|_{\infty}\] pero esto no nos ayuda a determinar el error en la iteración \(k\)-ésima pues desconocemos el valor real de \(\symbfit{x}\).

Por tanto, habremos de buscar alguna forma de expresar el error en términos únicamente de cantidades que podamos calcular.

Sea \(\{\symbfit{x^{(k)}}\}\) la sucesión de vectores construida usando un método iterativo consistente como el expresado anteriormente, y sea \(\symbfit{x}\in\mathbb{K}^n\) el punto fijo (solución de \(A\symbfit{x}=\symbfit{b}\)). Entonces, para cada \(k\in\mathbb{N}\setminus\{0\}\), se tiene:

1. \(\displaystyle e_k \leq \frac{\|T\|_{\infty}}{1 - \|T\|_{\infty}}\left\|\symbfit{x^{(k)}} - \symbfit{x^{(k-1)}}\right\|_{\infty}\)
2. \(\displaystyle e_k \leq \frac{\left(\|T\|_{\infty}\right)^k}{1 - \|T\|_{\infty}}\left\|\symbfit{x^{(1)}} - \symbfit{x^{(0)}}\right\|_{\infty}\)

1. Podemos usar mientras iteramos, llevando cuenta de la aproximación actual y de la anterior, para obtener una cota del error cometido en la iteración actual.
2. El punto se puede usar para, una vez realizada la primera iteración (es decir, una vez conocemos \(\symbfit{x^{(1)}}\)), estimar el número de iteraciones necesarias para asegurar un error menor que el umbral de tolerancia prefijado.

Consideremos, en \(\mathbb{R}\), el sistema de ecuaciones \[\left\{\begin{array}{cc} 2x_1+x_2 &= 4 \\ x_1 + 3x_2 & = 7\end{array}\right.\]

1. Partiendo del vector \((4, 1)^T\), ¿cuántas iteraciones del método de Jacobi estimamos a priori que se necesitan para obtener un error menor que \(10^{-5}\)?
2. Después de 4 iteraciones por Jacobi, ¿qué cota de error se obtiene para esa iteración?

La matriz es \(T = \begin{pmatrix}0 & \sfrac{-1}{2}\\ \sfrac{-1}{3} & 0\end{pmatrix}\), de donde \(\|T\|_{\infty} = \sfrac{1}{2} < 1\), luego el método es convergente. Se podría haber deducido también observando que la matriz de coeficientes es de diagonal estrictamente dominante.

1. Para calcular el número de iteraciones a priori, haremos uso de la expresión \[\displaystyle e_k \leq \frac{\left(\|T\|_{\infty}\right)^k}{1 - \|T\|_{\infty}}\left\|\symbfit{x^{(1)}} - \symbfit{x^{(0)}}\right\|_{\infty}\]

Ya conocemos \(\|T\|_{\infty} = \sfrac{1}{2}\) y \(\symbfit{x^{(0)}} = \begin{pmatrix}4 \\ -1\end{pmatrix}\). Nos falta conocer \(\symbfit{x^{(1)}}\), que se have mediante Jacobi.

Resulta \(\symbfit{x^{(1)}} = \begin{pmatrix}2.5 \\ 1\end{pmatrix}\), luego, para la iteración \(k\)-ésima, es \[e_k \leq \frac{\left(\sfrac{1}{2}\right)^k}{1 - \sfrac{1}{2}}\left\|\begin{pmatrix}2.5 \\ 1\end{pmatrix} - \begin{pmatrix}4 \\ -1\end{pmatrix}\right\|_{\infty} = \frac{\sfrac{1}{2^k}}{\sfrac{1}{2}}\left\|\begin{pmatrix}-1.5 \\ 2\end{pmatrix} \right\|_{\infty} = \frac{\sfrac{1}{2^k}}{\sfrac{1}{2}}\cdot 2 = \frac{1}{2^{k-2}}\]

Como queremos que el error sea menor que \(10^{-5}\), basta con imponer que \(\frac{1}{2^{k-2}} < 10^{-5}\), pues así el error en la iteración \(k\) será \(e_k \le \frac{1}{2^{k-2}} < 10^{-5}\).

Despejamos entonces \(k\) de la expresión anterior: \[k > 2+5\log_2 10 \approx 18.61\] por tanto, a partir de la iteración 19, podemos asegurar que la diferencia entre la aproximación obtenida y la solución real es menor que \(10^{-5}\).

2. Para calcular una cota del error tras la cuarta iteración, emplearemos la expresión:

luego, debemos calcular hasta \(\symbfit{x^{(4)}}\). \end{enumerate}

Se tiene que: \[\symbfit{x^{(2)}} = \begin{pmatrix}1.5\\1.5\end{pmatrix}, \symbfit{x^{(3)}} = \begin{pmatrix}1.25\\1.8333\end{pmatrix}, \symbfit{x^{(4)}} = \begin{pmatrix}1.0833\\1.9167\end{pmatrix}\] Luego \[e_k \leq \left\|\begin{pmatrix}1.0833\\1.9167\end{pmatrix} - \begin{pmatrix}1.25\\1.8333\end{pmatrix}\right\|_{\infty} = \left\|\begin{pmatrix}-0.1667\\0.0834\end{pmatrix}\right\|_{\infty} = 0.1667\]

La solución real de este sistema es \(\symbfit{x} = \begin{pmatrix}1\\2\end{pmatrix}\). Si calculamos el error real tras 4 iteraciones, tenemos:

Como era de esperar, por lo visto en teoría, efectivamente el error real \(e_k = 0.0833\) es menor que la cota calculada, \(0.1667\).

4.9 Ejercicios

4.10 FALTA:

• Descomposición de Cholesky (aquí o mejor en el de formas bilineales definidas positivas)
• Método numérico de relajación
• Vector residuo y número de condición.