Un ejercicio de aplicación de la factorización de Crout a la resolución de un sistema de ecuaciones lineales compatible determinado

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:

1. Resolveremos $LY=B$ para determinar el vector $Y=(y_1\,y_2\,y_3)^\top$
2. Una vez conocido $Y$, resolveremos finalmente $UX=Y$ para determinar el vector $X=(x_1\,x_2\,x_3)^\top$

Cálculo de las matrices $L$ y $U$

De acuerdo con el m. de Crout, la matriz triangular inferior es de la forma $L=\begin{pmatrix}\ell_{11}&0&0\\ \ell_{21}&\ell_{22}&0\\\ell_{31}&\ell_{32}&\ell_{33}\end{pmatrix}$ y la matriz triangular superior es $U=\begin{pmatrix}1&u_{12}&u_{13}\\ 0&1&u_{23}\\0&0&1\end{pmatrix}$. Entonces, como $$\begin{pmatrix}\ell_{11}&0&0\\ \ell_{21}&\ell_{22}&0\\\ell_{31}&\ell_{32}&\ell_{33}\end{pmatrix}\begin{pmatrix}1&u_{12}&u_{13}\\ 0&1&u_{23}\\0&0&1\end{pmatrix}=\begin{pmatrix}1&-1&0\\1&1&1\\1&0&-1\end{pmatrix}$$ por la definición de producto de matrices se tiene que
  $a_{11}=1=\ell_{11}$
  $a_{21}=1=\ell_{21}$
  $a_{31}=1=\ell_{31}$

  $a_{12}=-1=\ell_{11}\,u_{12}=1\cdot u_{12}\Rightarrow u_{12}=-1$
  $a_{22}=1=\ell_{21}\,u_{12}+1\cdot \ell_{22}=1\cdot u_{12}+\ell_{22}=-1+\ell_{22}\Rightarrow \ell_{22}=2$
  $a_{32}=0=\ell_{31}\,u_{12}=\ell_{32}\,u_{12}+1\cdot \ell_{32}=1\cdot (-1)+\ell_{32}\Rightarrow \ell_{32}=1$

  $a_{13}=0=\ell_{11}\,u_{13}=1\cdot u_{13}\Rightarrow u_{13}=0$
  $a_{23}=1=\ell_{21}\,u_{13}+\ell_{22}\,u_{23}=0+\ell_{22}\,u_{23}=0+2u_{23}\Rightarrow u_{33}=\frac{1}{2}$
  $a_{33}=-1=\ell_{31}\,u_{13}+\ell_{32}\,u_{23}+1 \cdot \ell_{33}=1 \cdot 0+1\cdot u_{23}+\ell_{33}=\frac{1}{2}+\ell_{33} \Rightarrow \ell_{33}=-\frac{3}{2}$

Por tanto $L=\begin{pmatrix}1&0&0\\ 1&2&0\\1&1&-\frac{3}{2}\end{pmatrix}$ y $U=\begin{pmatrix}1&-1&0\\ 0&1&\frac{1}{2}\\0&0&1\end{pmatrix}$

Abordamos ahora el primer paso, resolviendo $LY=B$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}y_1&&&&&=&4\\ y_1&+&y_2&&&=&3 \\ y_1&+&\frac{1}{2}\,y_2&+&y_3&=&2 \end{matrix}\right.$, encontramos fácilmente $y_1=4$, $y_2=-1$ y $y_3=-\frac{3}{2}$

Abordamos finalmente el segundo paso, resolviendo $UX=Y$

Escribiendo el sistema de ecuaciones, $\left\{ \begin{matrix}x_1&-&x_2&&&=&4\\ &&2x_2&+&x_3&=&-1 \\ &&&&-\frac{3}{2}x_3&=&-\frac{3}{2}\end{matrix}\right.$, de donde $x_1=3$, $x_2=-1$ y $x_3=1$.

Algoritmo de Crout

Arriba hemos realizado las operaciones paso a paso, si bien con un poco de paciencia podemos inducir las expresiones matemáticas que dan valor a los elementos de $L$ y $U$ en el caso general de tener que factorizar una matriz $A$ de orden $n$; esto lo podemos hacer partiendo de las regularidades que encontraremos para matrices de orden $3$. Obtendremos así lo que podemos entender como el algoritmo que cómodamente implementaremos mediante un lenguaje de programación, y, así, automatizaremos los cálculos. Se puede comprobar que:

• $\ell_{i1}=a_{i1}$ si $j=1$ para $i=1,\ldots,n$
• $u_{1j}=\dfrac{a_{1j}}{a_{11}}$ si $i=1$ para $j=2\,\ldots,n$
• $u_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{j-1}\,\ell_{ik}\,u_{kj}}{\ell_{ii}}$ si $i\lt j$ para $i=2,\ldots,n$
• $\ell_{ij}=\dfrac{a_{ij}-\displaystyle \sum_{k=1}^{i-1}\,\ell_{ik}\,u_{kj}}{u_{jj}}$ si $i\ge j$ para $j=2,\ldots,n$

$\diamond$