Se quiere resolver el sistema de ecuaciones lineales $$\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ x_1&+&x_2&+&x_3&=&3 \\ x_1&&&-&x_3&=&2 \end{matrix}\right.$$, que en otro artículo ya había sido resuelto por el método de Doolittle. Ahora vamos a resolverlo por el método de Crout.
Escribamos el sistema en forma matricial $$\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}4\\3\\2\end{pmatrix}$$ Expresado así, $AX=B$, donde $X=(x_1\,x_2\,x_3)^\top$, $b=(4\,3\,2)^\top$ y $A=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$, vamos a factorizar la matriz $A$ (que es regular) de la forma $A=LU$ por el método de Crout, donde $L$ es una matriz triangular inferior, y $U$ es una matriz triangular superior con unos en la diagonal principal —recordemos que los unos en la digonal principal estaban en la matriz $L$ en el método de Doolittle—.
Al obtener $A=LU$, el sistema de ecuacione puede escribirse de la forma $LUX=B$, esto es, $L(UX)=b$. Denotando $UX=Y$, la resolución constará de los siguientes pasos:
Las matrices de permutación permiten formalizar los procedimientos de pivotaje en la resolución de sistemas de ecuaciones lineales. Tal como ya se ha comentado en otros artículos, el pivotaje puede ser parcial (p. maximal por columnas) —en el que solamente se permutan filas de una matriz o de un sistema de ecuaciones lineales— o bien p. total (o p. completo) —en el que se permutan filas (entre el conjunto de las mismas), y también columnas (entre el conjunto de columnas)—. En este breve artículo hablaré de las operaciones de permutación maximal por columnas, que ilustraré con un sencillo ejemplo, aprovechándolo para mostrar cómo podemos utilizar esta operación para formalizar la resolución de un sistema de ecuaciones lineales en el que intercambiaremos (permutaremos) dos de sus filas.
Denominamos matriz de permutación a una matriz cuadrada $P$ formada a partir de la alteración del orden de las filas de la matriz identidad $I$ tal que al multiplicar por la izquierda a una matriz cuadrada del mismo orden $A$ permuta sus filas convenientemente. Así, por ejemplo, dada la matriz $A=\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ para intercambiar la primera y la tecera filas —siendo $I=\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}$ la matriz identidad del orden de la matriz dada— realizaremos el producto de $P=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}$ —notemos que intercambiamos las filas primera y tercera de la matriz identidad; las mismas que queremos inteercambiar en la matriz $A$— por la matriz $A$ por la izquierda: $PA=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ obteniendo como resultado la matriz pedida (con las filas primera y tercera permutadas): $\begin{pmatrix}4&-1&1\\1&2&3\\0&1&2\end{pmatrix}$.
Consideremos ahora el sistema de ecuacions lineales $$\left\{\begin{matrix} &&x_2&+&2x_3&=&1 \\ x_1&+&2x_2&+&3x_3&=&0 \\ 4x_1&-&x_2&+&x_3&=&0 \\ \end{matrix}\right.$$ que podemos expresar matricialmente de la forma $AX=B \quad (1)$, donde $A=\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}$ es la matriz de los coeficientes del sistema; $X=(x_1\,x_2\,x_3)^\top$ es la matriz columna de las incógnitas, y $B=(1\,0\,0)^\top$ es la matriz columna de los términos independientes.
Para reducir el sistema por el método de Gauss, nos interesa permutar la primera ecuación con alguna de las otras dos, habida cuenta de que el coeficiente de la primera incógnita de la primera ecuación es $0$. Para ello, debemos elegir el elemento pivote de la primera columna, que, como norma general, escogemos el elemento máximo de la misma; por lo tanto, en el caso que nos ocupa, intercambiaremos (permutaremos) la primera y la tercera filas. Expresado en forma matricial esto se va a traducir en multiplicar por la matriz de permutación correspondiente por la izquierda a ambos miembros de (1), así que escribiremos el sistema equivalente $PAX=PB$, esto es $$\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}0&1&2\\1&2&3\\4&-1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0&0&1\\0&1&0\\1&0&0\end{pmatrix}\begin{pmatrix}1\\0\\0\end{pmatrix}$$ y realizando las operaciones de multiplicación de matrices llegamos a $$\begin{pmatrix}4&-1&1\\1&2&3\\0&1&2\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}0\\0\\1\end{pmatrix}$$ Es decir, las ecuaciones del sistema equivalente quedan ahora en el siguiente orden $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ x_1&+&2x_2&+&3x_3&=&0 \\ &&x_2&+&2x_3&=&1 \\ \end{matrix}\right.$$ Podemos ahora realizar ya los pasos de reducción por Gauss: de $-\frac{1}{4}\,e_1+e_2 \rightarrow e_2$ obtenemos el sistema equivalente $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ &&9x_2&+&11x_3&=&0 \\ &&x_2&+&2x_3&=&1 \\ \end{matrix}\right.$$ y, finalmente, mediante la operación elemental (entre la segunda fila y la tercera del paso anterior) $-\frac{1}{9}\,e_1+e_3 \rightarrow e_3$ llegamos a $$\left\{\begin{matrix} 4x_1&-&x_2&+&x_3&=&0 \\ &&9x_2&+&11x_3&=&0 \\ &&&&\frac{7}{9}x_3&=&1 \\ \end{matrix}\right.$$ por lo que, despejando la tercera incóngita de la tercera ecuación, y sustituyendo regresivamente a las ecuaciones segunda y primera, obtenemos la solución: $x_1=-\dfrac{5}{7}$, $x_2=-\dfrac{11}{7}$ y $x_3=\dfrac{9}{7}$.
En ulteriores artículos mostraré cómo intervienen las matrices de permutación en la resolución de sistemas de ecuaciones lineales en el caso de optar por los métodos de factorización $LU$ de la matriz de los coeficientes del sistema.
Conviene tener en cuenta que toda matriz de permutación $P$ es simétrica: $P=P^\top$; y ortogonal: $P^{-1}=P^\top$, como puede demostrarse sin dificultad. Estas propiedades de las matrices de permutación facilitan la realización de determinados cálculos, tal y como veremos más adelante en este contexto de los métodos de reducción con pivotaje.$\diamond$
Se quiere resolver el sistema de ecuaciones lineales $$\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ x_1&+&x_2&+&x_3&=&3 \\ x_1&&&-&x_3&=&2 \end{matrix}\right.$$
Escribamos el sistema en forma matricial $$\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}4\\3\\2\end{pmatrix}$$ Expresado así, $Ax=b$, donde $x=(x_1,x_2,x_3)^\top$, $b=(4,3,2)^\top$ y $A=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$, vamos a factorizar la matriz $A$ (que es regular) de la forma $A=LU$ por el método de Doolittle, donde $L$ es una matriz triangular inferior con unos en la diagonal principal, y $U$ es una matriz triangular superior.
Al obtener $A=LU$, el sistema de ecuacione puede escribirse de la forma $LUX=B$, esto es, $L(UX)=b$. Denotando $UX=Y$, la resolución constará de los siguientes pasos:
Toda matriz de orden $n$, regular, $A$, cuyos determinantes principales sean distintos de cero, puede expresarse de la forma $A=LU$ mediante operaciones de reducción gaussiana sin necesidad de pivotaje, donde $L$ es una matriz triangular inferior de orden $n$ con unos en la diagonal principal y $U$ es una matriz triangular superior de orden $n$, siendo ambas matrices regulares (sus determinantes son no nulos).
En estas condiciones, las matrices $L$ y $U$ son únicas. En efecto, supongamos que no lo sean, esto es, si podemos escribir $A=L_1\,U_1=L_2\,U_2$; entonces, como todas estas matrices son inversibles, podemos escribir $L_{2}^{-1}\,L_1\,U_1=L_{2}^{-1}\,L_2\,U_2$, luego $L_{2}^{-1}\,L_1\,U_1=I\,U_2$ y por tanto $L_{2}^{-1}\,L_1\,U_1=U_2$; multiplicando ahora por $U_{2}^{-1}$ por la derecha en ambos miembros, tenemos $L_{2}^{-1}\,L_1\,U_1\,U_{2}^{-1}=U_2\,U_{1}^{-1}$, esto es $L_{2}^{-1}\,L_1\,U_1\,U_{1}^{-1}=U_2\,U_{1}^{-1}$ y por tanto, $L_{2}^{-1}\,L_1\,I=U_2\,U_{1}^{-1}$, que es lo mismo que $L_{2}^{-1}\,L_1=U_2\,U_{1}^{-1}$. La matriz inversa de una matriz triangular inferior es otra matriz triangular inferior, y lo mismo sucede con la inversa de una matriz triangular superior, así que, como el producto de dos matrices triangulares inferiores es otra matriz triangular inferior (primer miembro) y el producto de dos matrices triangulares superiores es una matriz triangular superior, la única manera de que esto sea posible es que ambos miembros de la igualdad sean la matriz identidad, por consiguiente: $L_{2}^{-1}\,L_1=I \Leftrightarrow L_2=L_1$ y $U_2\,U_{1}^{-1}\Leftrightarrow U_2=U_1$. $\square$
Observaciones:
Más adelante, extenderemos lo que vamos a deducir ahora al caso en que sí se requiere pivotamiento.
Empezaremos pues con el proceso de reducción. Siendo esta matriz de orden $3$, tendremos que llevar a cabo dos etapas.
Obtenemos ceros por debajo de la columna del elemento $a_{11}$ haciendo las siguientes transformaciones $\begin{matrix}f_1 \rightarrow f_1 \\ -\frac{1}{3}f_1+f_2\rightarrow f_2\\-\frac{2}{3}f_1+f_3\rightarrow f_3\end{matrix}$, obteniendo la matriz $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&\frac{5}{3}&-\frac{2}{3}\end{pmatrix}$. Finalmente, anularemos el elemento de la segunda fila y la segunda columna mendiante $\begin{matrix}f_1 \rightarrow f_1 \\f_2\rightarrow f_2\\-\frac{5}{4}f_2+f_3\rightarrow f_3\end{matrix}$ lo cual da lugar a la matriz reducida por Gauss $\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$, que, obviamente es una matriz triangular superior, y a la que denominaremos $U$.
La factorización de una matriz regular como producto de una matriz triangular inferior y una matriz triangular superior facilita la resolución de un sistema de ecuaciones lineales, así como diversos cálculos matriciales. Veamos tres utilidades importantes:
Notemos ahora que, de los multiplicadores de la transformación, podemos escribir esta otra matriz triangular inferior y a la que denominaremos $L$ y escribiremos $L=\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}$, donde hemos cambiado el signo de los números distintos de cero y de uno que hemos utilizado para efectuar las transformaciones. Pues bien, puede comprobarse que haciendo el producto $LU$ se obtiene la matriz $A$; en efecto $$\begin{pmatrix}1&0&0\\\frac{1}{3}&1&0\\\frac{2}{3}&\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}=\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}$$ Esta es pues una manera sencilla de efectuar la factorización de una matriz regular dada $A$ como producto de una matriz triangular inferior por una matriz triangular superior, $A=LU$, siempre que no tengamos que utilizar técnicas de pivotamiento parcial —por encontrarnos con ceros en los elementos de la diagonal principal en las etapas del proceso de reducción, en cuyo caso deberemos extender esta técnica, tal y como ya trataré más adelante artícul—. Veremos también las utilidades de dicha factorización en el próximo artículo.
Observación: Conociendo $L$, para calcular $U$, fijémonos en que de la factorización que buscamos, $A=LU$, tenemos que, multiplicando ambos miembros por $L^{-1}$ por la izquierda, se tiene que $U=L^{-1}\,A=(E(1)^{-1}\,E(2)^{-1})^{-1}\,A=(E(2)^{-1})^{-1}\,(E(1)^{-1})^{-1}\,A=E(2)\,E(1)\,A$, con lo cual la matriz triangular superior $U$ es, como debe ser —valga como comprobación—, $$U=\begin{pmatrix}1&0&0\\0&1&0\\0&-\frac{5}{4}&1\end{pmatrix}\begin{pmatrix}1&0&0\\-\frac{1}{3}&1&0\\-\frac{2}{3}&0&1\end{pmatrix}\begin{pmatrix}3&-1&1\\1&1&1\\2&1&0\end{pmatrix}=\begin{pmatrix}3&-1&1\\0&\frac{4}{3}&\frac{2}{3}\\0&0&-\frac{3}{2}\end{pmatrix}$$
Así, en el caso de factorizar una matriz de orden $4$, la matriz triangular inferior haremos el siguiente cálculo, $L=E(1)\,E(2)\,E(3)$; y para calcular la matriz triangular superior $U$ tendríamos tres pasos en lugar de dos: $\left\{\begin{matrix} A(2)=E(1)\,A(1)\\ A(3)=E(2)\,A(2) \\ U\equiv A(4)=E(3)\,A(3) \end{matrix}\right\}$ donde $A(1) \equiv A$, con las respectivas matrices de multiplicadores (ahora hay tres): $E(1)=\begin{pmatrix}1&0&0&0 \\ -\frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ -\frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ -\frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -\frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & -\frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -\frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$, así como sus respectivas inversas: $E(1)^{-1}=\begin{pmatrix}1&0&0&0 \\ \frac{a_{21}^{(1)}}{a_{11}^{(1)}} & 1 & 0 & 0 \\ \frac{a_{31}^{(1)}}{a_{11}^{(1)}} & 0 & 1 & 0 \\ \frac{a_{41}^{(1)}}{a_{11}^{(1)}} & 0 & 0 & 1 \end{pmatrix}$, $E(2)^{-1}=\begin{pmatrix}1& 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & \frac{a_{32}^{(2)}}{a_{22}^{(2)}} & 1 & 0 \\ 0 & \frac{a_{42}^{(2)}}{a_{22}^{(2)}} & 0 & 1 \end{pmatrix}$ y $E(3)^{-1}=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & \frac{a_{43}^{(3)}}{a_{33}^{(3)}} & 1\end{pmatrix}$
Referencias:
[1] A. Aubanell, A. Benseny y A. Deslshams, Útiles de cálculo numérico (Labor, Barcelona, 1993).
[2] C. Moreno, Introducción al cálculo numérico (UNED, Madrid, 2011).
