% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root =  ../Metodos_Numericos.tex

\chapter[Sist. Lineales]{Solución de Sistemas Lineales de Ecuaciones}

%\section{Matlab como herramienta de desarrollo}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Suma, resta, multiplicación y división de matrices}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Una matriz se define como un arreglo bidimensional de datos, que tiene n renglones y m columnas. Un elemento del arreglo puede ser identificado con $a_{ij}$

$$
A = \left(%
\begin{array}{ccccc}
   a_{1,1} & a_{1,2} & a_{1,3} & ... & a_{1,N}  \\
   a_{2,1} & a_{2,2} & a_{2,3} & ... & a_{2,N}  \\
   a_{3,1} & a_{3,2} & a_{3,3} & ... & a_{3,N}  \\
       &     &     &     &      \\
   a_{M,1} & a_{M,2} & a_{M,2} & ... & a_{M,N}  \\
\end{array}%
\right)
$$

Algunas de las operaciones básicas que pueden realizarse con matrices son suma, resta y multiplicación. La división de matrices como tal no existe y en su lugar se calcula la inversa.

\subsection{Suma de matrices}

Para que la sumar las matrices $A_{N \times M}$ y $B_{R \times S} $, se requiere que las matrices tengan el mismo número de renglones $N = R$ y de columnas $M = S$  . Si queremos encontrar la suma $C = A + B$, cada elemento de la matriz $C$ lo calculamos de la siguiente forma:

$$
C_{n,m} = A_{n,m} + B_{n,m}
$$

para todos lo $n,m$ en la matriz $C$

La implementación en Matlab es:

\begin{verbatim}
function C = suma(A, B)

  [N, M] = size(A);
  [R, S] = size(B);
  
  if N~=R || M~=S
      fprintf('Las dimensiones tienen que ser iguales\n');
      return;
  end;
  
  C = zeros(N, M);
  
  for n=1:N
      for m = 1:M
          C(n,m) = A(n,m)+B(n,m);
      end;
  end
end
\end{verbatim}

\subsection{Resta de matrices}

En este caso, se deben cumplir las mismas propiedades que la rsuma de matrices y el calculo de los elemento de la matriz $C$ se calculan como:

$$
C_{n,m} = A_{n,m} - B_{n,m}
$$

para todos lo $n,m$ en la matriz $C$

La implementación en Matlab es:

\begin{verbatim}
function C = resta(A, B)

  [N, M] = size(A);
  [R, S] = size(B);
  
  if N~=R || M~=S
      fprintf('Las dimensiones tienen que ser iguales\n');
      return;
  end;
  
  C = zeros(N, M);
  
  for n=1:N
      for m = 1:M
          C(n,m) = A(n,m)-B(n,m);
      end;
  end
end
\end{verbatim}

\subsection{Multiplicación de matrices}

Para realizar el producto $C = A*B$ tenemos que considerar que el producto existe si

1.- El número de columnas de $A_{N \times M}$ es igual al número de renglones de $B_{R \times S}$, esto es $M = R$. y

2.- Las dimensiones de la matriz resultante tendrá el mismo numero de renglones que la matriz $A$ y el mismo número de columnas que la matriz $B$.  La matriz resultante $C$ será de tamaño $N \times S$ 

3.- El cálculo de los elementos de la matriz C se lleva a cabo haciendo :

$$
C_{n,m} = \sum_{k=1}^{M} a_{n,k} * b_{k,m}
$$

para todos lo $n,m$ en la matriz $C$

La implementación en Matlab es:

\begin{verbatim}
function C = multiplica(A, B) 
   [N,M] = size(A);
   [R,S] = size(B);
   
   if M~=R 
       fprintf('El numero de renglones no es igual al numero de columnas\n');
       return;
   end;
   
   C = zeros(N, S);
   
   for n=1:N
       for m= 1:S
           suma = 0;
           for k = 1:M
               suma = suma + A(n,k)*B(k,m); 
           end;
           C(n,m) = suma;
       end;
   end;
end
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Factorización triangular LU}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Supongamos que la matriz de coeficientes $A$ de un sistema lineal $Ax=b$ admite una factorización triangular como la siguiente.

 
\begin{equation*}
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3}\\
a_{3,1} &a_{3,2} & a_{3,3}\\
\end{array}
\right]
=
\left[
\begin{array}{ccc}
1 & 0 & 0\\
l_{2,1} & 1 & 0 \\
l_{3,1} & l_{3,2} & 1 \\
\end{array}
\right]	
 *
\left[
\begin{array}{ccc}	
u_{1,1} & u_{1,2} & u_{1,3} \\
0 & u_{2,3} & u_{2,3} \\
0 & 0 & u_{3,3} \\
\end{array}
\right]
\end{equation*}

Entonces la solución la podemos plantear como

\begin{equation*}
[LU]x = b
\end{equation*}

El cual se puede solucionar resolviendo primero $Ly=b$ y luego el sistema $Ux=y$. Estos sistemas pueden resolverse de manera muy eficiente utilizando sustitución hacia adelante y sustitución hacia atrás.

\subsection{Sustitución hacia adelante}

Comencemos por dar un sistema triangular inferior

\begin{equation*}
\left[
\begin{array}{cccc}
  1         &  0         &  0        & 0  \\
  l_{2,1} &  1         &  0        & 0 \\
  l_{3,1} &  l_{3,2} & 1         & 0 \\
  l_{4,1} &  l_{4,2} & l_{4,3} & 1\\  
\end{array}
\right]
\left[
\begin{array}{c}
  y_1   \\
  y_2   \\
  y_3 \\
  y_4   
\end{array}
\right]
=
\left[
\begin{array}{c}
  b_1   \\
  b_2   \\
  b_3 \\
  b_4   
\end{array}
\right]
\end{equation*}

Si despejamos comenzando con la ecuación $I$ a $y_1$ y así sucesivamente para cada una de las ecuaciones tenemos:

\begin{eqnarray*}
y_1 &=& b_1\\
y_2 &=& b_2 - l_{2,1} y_1\\
y_3 & =& b_3 - l_{3,1} y_1 - l_{3,2} y_2\\
y_4 &=& b_4 - l_{4,1}y_1 - l_{4,2} y_2 - l_{4,3} y_3
\end{eqnarray*}

La formula para calcular la sustitución hacia adelante es:

\begin{equation*}
y_k = b_k - \sum_{n=1}^{n<k} l_{k,n} x_n
\end{equation*}

Para todos los valores $k = 1, 2, 3, \cdots, N$

La implementación en matlab es:

\begin{verbatim}
function y = sustitucion_hacia_adelante(l,b)
     N = length(l(:,1));
     y = zeros(N,1);
     for k = 1: N
          suma = 0;
          for n = 1: k-1
               suma = suma + (l(k,n)*y(n));
          end;
          y(k) = (b(k) - suma);
     end;
end
\end{verbatim}
 

\subsection{Sustitución hacia atrás}

Consideremos ahora un sistema triangular superior

\begin{equation*}
\left[
\begin{array}{cccc}
  u_{1,1} &  u_{1,2}  &  u_{1,3}  & u_{1,4}  \\
  0          &  u_{2,2}  &  u_{2,3}  & u_{2,4} \\
  0          &  0           &  u_{3,3}  & u_{3,4} \\
  0          &  0           & 0            & u_{4,4} \\  
\end{array}
\right]
\left[
\begin{array}{c}
  x_1   \\
  x_2   \\
  x_3 \\
  x_4   
\end{array}
\right]
=
\left[
\begin{array}{c}
  y_1   \\
  y_2   \\
  y_3 \\
  y_4   
\end{array}
\right]
\end{equation*}

Para este caso comenzamos despejando $x_4$ de la última ecuación y sustituimos y despejamos en las anteriores tal como:

\begin{eqnarray*}
x_4 &=& \frac{y_4}{u_{4,4}}\\
x_3 &=& \frac{y_3 - u_{3,4} x_4}{u_{3,3}}\\
x_2 &=& \frac{y_2 - u_{2,3} x_3 - u_{2,4} x_4}{u_{2,2}}\\
x_1 &=& \frac{y_1 - u_{1,2} x_2 - u_{1,3} x_3 - u_{1,4} x_4}{u_{1,1}}\\
\end{eqnarray*}

La formula para calcular la sustitución hacia atrás es:

\begin{equation*}
x_k = \frac{y_k - \sum_{n=k+1}^{N} u_{k,n} x_n}{u_{k,k}}
\end{equation*}
 

La implementación en Java es:

 
\begin{verbatim}
function x = sustitucion_hacia_atras(u, y)
    N = length(y);
    x = zeros(N,1);
    for k = N: -1: 1
        suma = 0;
        for n = k+1: N
            suma = suma +(u(k,n)*x(n));
        end;
        x(k) = (y(k)-suma)/u(k,k);
    end;
end
\end{verbatim}

\subsection{Ejemplo}

Resolver el sistema de ecuaciones

\begin{eqnarray*}
  x_1 +  2x_2 +  4x_3 +  x_4 = 21\\
2x_1 +  8x_2 +  6x_3 + 4x_4 = 52\\
3x_1 + 10x_2 +  8x_3 + 8x_4 = 79\\
4x_1 + 12x_2 + 10x_3 + 6x_4 = 82\\
\end{eqnarray*}
 
 Supongamos que tenemos la descomposición triangular siguiente
 \begin{equation*}
\left[
\begin{array}{cccc}
1 &  2 & 4 & 1 \\     
2 &  8 & 6  & 4 \\ 
3 & 10 &  8 &  8 \\
4 & 12 & 10 & 6  
\end{array}
\right]
=
\left[
\begin{array}{cccc}  
1 & 0 & 0 & 0 \\
2 & 1 & 0 & 0 \\
3 & 1 & 1 & 0 \\ 
4 & 1 & 2 & 1 \\ 
\end{array}
\right] 
\times
\left[
\begin{array}{cccc}    
1 & 2 & 4 & 1\\
0 & 4 &–2 & 2 \\
0 & 0 &–2 & 3 \\
0 & 0 & 0 &–6 \\
\end{array}
\right]
\end{equation*}

Primero resolvemos el sistema de ecuaciones

\begin{eqnarray*}
y_1                         &=& 21\\
2y_1 +  y_2              &=& 52\\
3y_1 +  y_2 +   y_3      &=& 79\\
4y_1 +  y_2 +  2y_3 +  y_4 &=& 82
\end{eqnarray*}


El cual corresponde a un sistema triangular inferior

\begin{eqnarray*}
y_1 &=&  21\\
y_2 &=&  10\\
y_3 &=&   6\\
y_4 &=& -24\
\end{eqnarray*}
 

Posteriormente solucionamos $U x = y$
\begin{eqnarray*}
x_1 + 2x_2  + 4x_3 +  x_4 &=& 21\\
      4x_2  - 2x_3 + 2x_4 &=& 10\\
           - 2x_3 + 3x_4 &=&  6\\
                 - 6x_4 &=&-24\\
\end{eqnarray*}

Utilizando sustitución hacia atrás tenemos

\begin{equation*}
X = [1,2,3,4]T
\end{equation*}

\subsection{Factorización de Crout}

Consideremos un sistema más grande, por ejemplo de 4 ecuaciones con cuatro incógnitas.

\[
\left(
\begin{array}{cccc}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4}\\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4}\\	
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4}\\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4}\\
\end{array}
\right)
=
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0\\
l_{2,1} & 1 & 0 & 0\\	
l_{3,1} & l_{3,2} & 1 & 0\\
l_{4,1} & l_{4,2} & l_{4,3} & 1\\
\end{array}
\right)
\times
\left(
\begin{array}{cccc}
u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\
0 & u_{2,2} & u_{2,3} & u_{2,4}\\	
0 & 0 & u_{3,3} & u_{3,4} \\
0 & 0 & 0 & u_{4,4}\\
\end{array}
\right)
\]
 
Si hacemos las multiplicaciones indicadas y dejamos los valores en la misma matriz tenemos:

Para el primer renglón de la matriz ($n=1$)

arriba de la diagonal

\begin{equation*}
\begin{array}{cc}
m=1  & a_{1,1} =  a_{1,1}\\
m=2 &  a_{1,2} =  a_{1,2}\\
m=3 & a_{1,3} =  a_{1,3}\\
m=4 & a_{1,4} =  a_{1,4}\\
\end{array}
\end{equation*}

Para el segundo renglón ($n=2$)

bajo la diagonal

\begin{equation*}
\begin{array}{cc}
m=1 & a_{2,1} =  a_{2,1}  / a_{1,1}\\
\end{array}
\end{equation*}

arriba de la diagonal

\begin{equation*}
\begin{array}{cc}
m=2 & a_{2,2} =  a_{2,2} - a_{2,1}a_{1,2}\\
m=3 & a_{2,3} =  a_{2,3} - a_{2,1}a_{1,3}\\
m=4 & a_{2,4} =  a_{2,4} - a_{2,1} a_{1,4}\\
\end{array}
\end{equation*}
 

Tercer renglón ($n=3$)

bajo la diagonal

\begin{equation*}
\begin{array}{cc}
m=1 & a_{3,1} = a_{3,1} /a_{1,1}\\
m=2 & a_{3,2} =  (a_{3,2} - a_{3,1}a_{1,2})/a_{2,2}\\
\end{array}
\end{equation*}
 
arriba de la diagonal

\begin{equation*}
\begin{array}{cc}
m=3 & a_{3,3} =  a_{3,3} - a_{3,1} a_{1,3} - a_{3,2}a_{2,3}\\
m=4 & a_{3,4} = a_{3,4} - a_{3,1}a_{1,4}  - a_{3,2}a_{2,4}\\
\end{array}
\end{equation*}
 
Y finalmente para el cuarto renglón ($n=4$)

bajo la diagonal

\begin{equation*}
\begin{array}{cc}
m=1 & a_{4,1} = a_{4,1} / a_{1,1}\\
m=2 & a_{4,2} = ( a_{4,2} - a_{4,1} a_{1,2} )/a_{2,2}\\
m=3 & a_{4,3} = (a_{4,3} - a_{4,1} a_{1,3} - a_{4,2} a_{2,3} ) / a_{3,3}\\
\end{array}
\end{equation*}
 
arriba de la diagonal

\begin{equation*}
\begin{array}{cc}
m=4 & a_{4,4} =  a_{4,4} - a_{4,1} a_{1,4} - a_{4,2} a_{2,4} - a_{4,3}a_{3,4}\\
\end{array}
\end{equation*}

Podemos ver, que cualquier elemento bajo la diagonal se calcula como:

\begin{equation*}
a_{n,m} = \frac{a_{n,m} - \sum_{k=1}^{m-1} a_{n,k} a_{k,m}}{a_{m,m}}
\end{equation*}
\begin{equation*}
 n=1,\cdots, N \hspace{1cm} m = 1, \cdots, n-1
\end{equation*}
 

Para los términos arriba de la diagonal

\begin{equation*}
a_{n,m} = a_{n,m} - \sum_{k=1}^{n-1} a_{k,m}
\end{equation*}

\begin{equation*}
n=1,\cdots,N \hspace{1cm} m = n,...,N.
\end{equation*}

La implementación en Matlab es:

\begin{verbatim}
function a = Factorizacion(A)
   N = length(A);
   a = A;

   for n = 1:N
      for m = 1: n-1 
        suma = 0;
        
        for k = 1:m-1
          suma = suma + a(n,k) * a(k,m);
        end;
        
        a(n,m) = (a(n,m) - suma)/ a(m,m);
      end;

      for m = n:N
        suma = 0;
        for k = 1: n-1
          suma = suma + a(n,k) * a(k,m);
        end;
        a(n,m) = a(n,m) - suma;
      end;
   end;
end
\end{verbatim}
 

\subsection{Ejemplo}

Hacer la factorización LU de la siguiente matriz, utilizando eliminación Gaussiana. Corroborar utilizando Factorización de Crout.

\[
\left(
\begin{array}{cccc}
1 & 2 & 4 & 1\\
2 & 8 & 6 & 4\\
3 & 10 & 8 & 8\\
4 & 12 & 10 & 6\\
\end{array}
\right)
\]
 
Paso 1.

\[
\left(
\begin{array}{cccc}
1 & 2 & 4 &1 \\
2 & 4 & -2 & 2\\
3 & 4 & -4 & 5\\
4 & 4 & -6 & 2\\
\end{array}
\right)
\]
 
Paso 2.
 
\[
\left(
\begin{array}{cccc}
1 & 2 & 4 & 1\\
2 & 4 & -2 & 2\\
3 & 1 & -2 & 3\\
4 & 1 & -4 & 0\\
\end{array}
\right)
\]

Paso 3.

\[
\left(
\begin{array}{cccc}
1 & 2 & 4 & 1\\
2 & 4 & -2 & 2\\
3 & 1 & -2 & 3\\
4 & 1 & -2 & 6\\
\end{array}
\right)
\]

La ejecución en Matlab queda

\begin{verbatim}
>> Factorizacion( [1  2  4  1; 2 8 6 4;3 10 8 8;4 12 10 6])

ans =

     1     2     4     1
     2     4    -2     2
     3     1    -2     3
     4     1     2    -6
\end{verbatim}

\subsection{Ejemplo}
Dado el sistema lineal de ecuaciones, calcular:
a) La factorización de Crout  y
b) Calcular la solución utilizando sustitución hacia adelante y atrás

\begin{equation*}
\left[
\begin{array}{r r r r}
    10  &  -1 &    2  &   1 \\
    -1   & 15 &   -3  &    1\\
     2   & -3  &    6  &  -3\\
     1   &  1  &   -3  &   7\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
  1 \\
  2 \\
  2 \\
  1\\
\end{array}
\right]
\end{equation*}

a) Realizamos la factorización de LU o Crout 

\begin{verbatim}
>> Af = Factorizacion([10 -1 2 1; -1 15 -3 1; 2 -3 6 -3; 1 1 -3 7])

Af =

   10.0000   -1.0000    2.0000    1.0000
   -0.1000   14.9000   -2.8000    1.1000
    0.2000   -0.1879    5.0738   -2.9933
    0.1000    0.0738   -0.5899    5.0529
\end{verbatim}

b) Par la sustitución hacia adelante hacemos

\begin{verbatim}
>> y = sustitucion_hacia_adelante(Af, [1,2,2,1]')

y =

    1.0000
    2.1000
    2.1946
    2.0397
\end{verbatim}

para luego realizar la sustitución hacia atrás

\begin{verbatim}
>> x = sustitucion_hacia_atras(Af, y)

x =

   -0.0508
    0.2372
    0.6707
    0.4037
\end{verbatim}

\subsection{Cálculo de la Matriz inversa utilizando descomposición LU}

Si consideramos un sistema lineal $Ax = b$ la solución la podemos llevar a cabo utilizando el concepto de matriz inversa tal que $x = A^{-1}b$, por lo que podemos consideran que son métodos equivalentes. Consideremos un sistema de 4 ecuaciones con 4 incógnitas como el siguiente:

\begin{equation*}
Ax = b
\end{equation*}

\begin{equation*}
\left[
\begin{array}{cccc}
a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4}\\
a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4}\\	
a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4}\\
a_{4,1} & a_{4,2} & a_{4,3} & a_{4,4}\\
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2\\	
x_3\\
x_4\\
\end{array}
\right]
=
\left[
\begin{array}{c}
0\\
0\\	
1 \\
0\\
\end{array}
\right]
\end{equation*}

Si tuviéramos la matriz  inversa,  podemos calcular los valores de $x$ haciendo:

\begin{equation*}
x = A^{-1}b
\end{equation*}

\begin{equation*}
\left[
\begin{array}{c}
x_1\\
x_2\\	
x_3\\
x_4\\
\end{array}
\right]
=
\left[
\begin{array}{cccc}
a'_{1,1} & a'_{1,2} & a'_{1,3} & a'_{1,4}\\
a'_{2,1} & a'_{2,2} & a'_{2,3} & a'_{2,4}\\	
a'_{3,1} & a'_{3,2} & a'_{3,3} & a'_{3,4}\\
a'_{4,1} & a'_{4,2} & a'_{4,3} & a'_{4,4}\\
\end{array}
\right]
\left[
\begin{array}{c}
0\\
0\\	
1 \\
0\\
\end{array}
\right]
= 
\left[
\begin{array}{c}
a'_{1,3} \\
a'_{2,3} \\	
a'_{3,3} \\
 a'_{4,3} \\
\end{array}
\right]
\end{equation*}

donde $a'_{n,m}$ son los elementos de la matriz inversa $A^{-1}$. Note que poniendo un 1 en $b_3$ obtuvimos la tercer columna de la matriz inversa, por lo tanto si solucionamos el sistema de ecuaciones utilizando descomposición LU obtendríamos la tercer columna.

Podemos generalizar el método si definimos 

\begin{equation*}
b^{(k)} =
\left\{
\begin{array}{ccc}
  b_n = 1 &  si  & n = k   \\
  b_n = 0 &  si  & n <> k 
\end{array}
\right.
\end{equation*}

Entonces si resolvemos el sistema $A x = b^{(k)}$ por cualquier método, calculamos la $k$-esima colma de la matriz inversa y la solución será $x = A^{-1}_k$. Para calcular todas las columnas de la matriz inversa solamente es necesario calcular este procedimiento para $k=1, 2, \cdots, N$. Para hacer eficiente el procedimiento solamente factorizamos las matriz una vez y utilizando sustitución hacia adelante y atrás resolvemos para cada uno de los valores de $k$.  El código para hacer esto es:

\begin{verbatim}
function Ainv = inversa(A) 
N = length(A);

Af = Factorizacion(A);

Ainv = zeros(N,N);

for k = 1: N
    b = zeros(N,1);
    b(k) = 1;
    y = sustitucion_hacia_adelante(Af, b);
    x = sustitucion_hacia_atras(Af, y);
    Ainv(:,k) = x;
end;
end
\end{verbatim}

\subsection{Ejemplo}

Calcular la inversa de la matriz
\begin{equation*}
A = 
\left[
\begin{array}{r r r r}
    10  &  -1 &    2  &   1 \\
    -1   & 15 &   -3  &    1\\
     2   & -3  &    6  &  -3\\
     1   &  1  &   -3  &   7\\
\end{array}
\right] 
\end{equation*}

Para comenzar hacemos la factorización de la matriz utilizando el metido de Crout o LU

\begin{equation*}
A = 
\left[
\begin{array}{r r r r}
    10  &  -1 &    2  &   1 \\
    -1   & 15 &   -3  &    1\\
     2   & -3  &    6  &  -3\\
     1   &  1  &   -3  &   7\\
\end{array}
\right] 
\end{equation*}

El resultado es

\begin{equation*}
L U = 
\left[
\begin{array}{rrrr}
    1.0000  &    0.0000  &  0.0000  &  0.0000 \\
   -0.1000   &   1.0000  &  0.0000  &  0.0000 \\
    0.2000   & -0.1879   &  1.0000  &  0.0000\\
    0.1000   &  0.0738   & -0.5899  &  1.0000\\
\end{array}
\right] 
\left[
\begin{array}{rrrr}
   10.0000  &  -1.0000  &  2.0000  &  1.0000 \\
    0.0000   & 14.9000  & -2.8000  &  1.1000 \\
    0.0000   &  0.0000   &  5.0738  &  -2.9933\\
    0.0000   &  0.0000   &  0.0000  &  5.0529\\

\end{array}
\right] 
\end{equation*}

Para calcular la primer columna de la matriz inversa hacemos
\begin{equation*}
\left[
\begin{array}{rrrr}
   10.0000  &  -1.0000  &  2.0000  &  1.0000 \\
    0.0000   & 14.9000  & -2.8000  &  1.1000 \\
    0.0000   &  0.0000   &  5.0738  &  -2.9933\\
    0.0000   &  0.0000   &  0.0000  &  5.0529\\
\end{array}
\right] 
\left[
\begin{array}{r}
y_1^{(1)} \\
y_2^{(1)} \\
y_3^{(1)} \\
y_4^{(1)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
1 \\
0 \\
0 \\
0 \\
\end{array}
\right] 
\end{equation*}

La solución es $y^{(1)} = [1.0000, 0.1000, -0.1812, -0.2143]^T$ 

\begin{equation*}
\left[
\begin{array}{rrrr}
    1.0000  &    0.0000  &  0.0000  &  0.0000 \\
   -0.1000   &   1.0000  &  0.0000  &  0.0000 \\
    0.2000   & -0.1879   &  1.0000  &  0.0000\\
    0.1000   &  0.0738   & -0.5899  &  1.0000\\
\end{array}
\right] 
\left[
\begin{array}{r}
x_1^{(1)} \\
x_2^{(1)} \\
x_3^{(1)} \\
x_4^{(1)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
    1.0000\\
    0.1000\\
   -0.1812\\
   -0.2143\\
\end{array}
\right] 
\end{equation*}

La primer columna de la matriz inversa es $x^{(1)} = [0.1162, -0.0016, -0.0607, -0.0424]^T$

Para la segunda columna de la matriz inversa hacemos

\begin{equation*}
\left[
\begin{array}{rrrr}
   10.0000  &  -1.0000  &  2.0000  &  1.0000 \\
    0.0000   & 14.9000  & -2.8000  &  1.1000 \\
    0.0000   &  0.0000   &  5.0738  &  -2.9933\\
    0.0000   &  0.0000   &  0.0000  &  5.0529\\
\end{array}
\right] 
\left[
\begin{array}{r}
y_1^{(2)} \\
y_2^{(2)} \\
y_3^{(2)} \\
y_4^{(2)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
0 \\
1 \\
0 \\
0 \\
\end{array}
\right] 
\end{equation*}

La solución es $y^{(2)} = [0, 1.0000, 0.1879, 0.0370]^T$ 

\begin{equation*}
\left[
\begin{array}{rrrr}
    1.0000  &    0.0000  &  0.0000  &  0.0000 \\
   -0.1000   &   1.0000  &  0.0000  &  0.0000 \\
    0.2000   & -0.1879   &  1.0000  &  0.0000\\
    0.1000   &  0.0738   & -0.5899  &  1.0000\\
\end{array}
\right] 
\left[
\begin{array}{r}
x_1^{(2)} \\
x_2^{(2)} \\
x_3^{(2)} \\
x_4^{(2)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
         0\\
    1.0000\\
    0.1879\\
    0.0370\\
\end{array}
\right] 
\end{equation*}

La segunda columna de la matriz inversa es $x^{(2)} = [   -0.0016, 0.0743, 0.0414, 0.0073]^T$

Para la tercercolumna de la matriz inversa hacemos

\begin{equation*}
\left[
\begin{array}{rrrr}
   10.0000  &  -1.0000  &  2.0000  &  1.0000 \\
    0.0000   & 14.9000  & -2.8000  &  1.1000 \\
    0.0000   &  0.0000   &  5.0738  &  -2.9933\\
    0.0000   &  0.0000   &  0.0000  &  5.0529\\
\end{array}
\right] 
\left[
\begin{array}{r}
y_1^{(3)} \\
y_2^{(3)} \\
y_3^{(3)} \\
y_4^{(3)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
0 \\
0 \\
1 \\
0 \\
\end{array}
\right] 
\end{equation*}

La solución es $y^{(3)} = [0, 0, 1.0000, 0.5899]^T$ 

\begin{equation*}
\left[
\begin{array}{rrrr}
    1.0000  &    0.0000  &  0.0000  &  0.0000 \\
   -0.1000   &   1.0000  &  0.0000  &  0.0000 \\
    0.2000   & -0.1879   &  1.0000  &  0.0000\\
    0.1000   &  0.0738   & -0.5899  &  1.0000\\
\end{array}
\right] 
\left[
\begin{array}{r}
x_1^{(3)} \\
x_2^{(3)} \\
x_3^{(3)} \\
x_4^{(3)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
         0\\
         0\\
    1.0000\\
    0.5899\\
\end{array}
\right] 
\end{equation*}

La tercer columna de la matriz inversa es $x^{(3)} = [    -0.0607, 0.0414, 0.2660, 0.1168]^T$

Para la última  columna de la matriz inversa hacemos

\begin{equation*}
\left[
\begin{array}{rrrr}
   10.0000  &  -1.0000  &  2.0000  &  1.0000 \\
    0.0000   & 14.9000  & -2.8000  &  1.1000 \\
    0.0000   &  0.0000   &  5.0738  &  -2.9933\\
    0.0000   &  0.0000   &  0.0000  &  5.0529\\
\end{array}
\right] 
\left[
\begin{array}{r}
y_1^{(4)} \\
y_2^{(4)} \\
y_3^{(4)} \\
y_4^{(4)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
0 \\
0\\
0 \\
1 \\
\end{array}
\right] 
\end{equation*}

La solución es $y^{(4)} = [0, 0, 0, 1]^T$ 

\begin{equation*}
\left[
\begin{array}{rrrr}
    1.0000  &    0.0000  &  0.0000  &  0.0000 \\
   -0.1000   &   1.0000  &  0.0000  &  0.0000 \\
    0.2000   & -0.1879   &  1.0000  &  0.0000\\
    0.1000   &  0.0738   & -0.5899  &  1.0000\\
\end{array}
\right] 
\left[
\begin{array}{r}
x_1^{(4)} \\
x_2^{(4)} \\
x_3^{(4)} \\
x_4^{(4)} \\
\end{array}
\right] 
=\left[
\begin{array}{r}
         0\\
         0\\
         0\\
         1\\
\end{array}
\right] 
\end{equation*}

La última columna de la matriz inversa es $x^{(4)} = [   -0.0424, 0.0073, 0.1168, 0.1979]^T$

Finalmente la matriz inversa es la concatenación de la soluciones

\begin{equation*}
A^{-1} = [x^{(1)},x^{(2)}, x^{(3)}, x^{(4)}]  =
\left[
\begin{array}{rrrr}
    0.1162   & -0.0016 &  -0.0607 &  -0.0424 \\
   -0.0016   & 0.0743  &  0.0414  &  0.0073\\
   -0.0607   & 0.0414  &  0.2660  &  0.1168\\
   -0.0424   & 0.0073  &  0.1168  &  0.1979\\
\end{array}
\right] 
\end{equation*}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Eliminación Gaussiana}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Consideremos que tenemos un sistema lineal $Ax=b$, donde la matriz $A$ no tiene las condiciones de ser triangular superior.

$$
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{array}%
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
$$

Comenzaremos por despejar  $x_1$ de la ecuación $I$

$$
x_1 = \frac{b_1 - a_{1,2}x_2 - a_{1,3}x_3}{a_{1,1}}
$$

y sustituimos en las ecuaciones $II$

$$
a_{2,1} \left( \frac{b_1 - a_{1,2}x_2 - a_{1,3}x_3}{a_{1,1}} \right) + a_{2,2} x_2
+ a_{2,3}x_3 = b_2
$$

agrupando términos semejantes tenemos:

$$
\left(a_{2,2}-\frac{a_{2,1} *a_{1,2}}{a_{1,1}} \right) x_2 + \left (a_{2,3}-a_{2,1}
\frac{a_{1,3}}{a_{1,1}} \right) x_3 = \left (b_2 - a_{2,1} \frac{b_1}{a_{1,1}} \right)
$$

Ahora sustituimos en la ecuación $III$

$$
\left(a_{3,1}\frac{b_1 - a_{1,2}x_2 - a_{1,3}x_3}{a_{1,1}} \right ) + a_{3,2}x_2
+ a_{3,3}x_3 = b_3
$$
$$ \left(a_{3,2}-
\frac{a_{3,1} a_{1,2}}{a_{1,1}} \right)x_2+ \left(a_{3,3}-
\frac{a_{3,1}a_{1,3}}{a_{1,1}}\right) x_3 = \left(b_3-a_{3,1} \frac{b_1}{a_{1,1}}\right)
$$

Lo cual nos da un nuevo sistema simplificado dada por

$$
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
      0 & a'_{2,2} & a'_{2,3} \\
      0 & a'_{3,2} & a'_{3,3} \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] = \left[
\begin{array}{c}
  b_1 \\
  b'_2 \\
  b'_3 \\
\end{array}
\right ]
$$

donde los valores de $a'_{n,m}$ los calculamos como:

\begin{eqnarray*}
a'_{n,m} & = & a_{n,m}-a_{n,k} \frac{a_{k,m}}{a_{k,k}}\\
b'_n & = & b_n - a_{n,k} \frac{b_k}{a_{k,k}}
\end{eqnarray*}

Si repetimos el procedimiento, ahora para $x_1$, tendremos un sistema dado como:

$$
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
      0 & a'_{2,2} & a'_{2,3} \\
      0 & 0 & a''_{3,3} \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  b_1 \\
  b'_2 \\
  b''_3 \\
\end{array}
\right]
$$

La solución del sistema Triangular superior lo podemos realizar utilizando el método de sustitución hacia atrás. La implementación en Matlab queda

\begin{verbatim}
function x = Eliminacion_Gaussiana(A0, b0)

  a = A0;
  b = b0;
  
  N = length(b);
  
  for k =1:N-1
    for n = k+1:N
        b(n) = b(n) - a(n,k)*b(k)/a(k,k);
        for m=N:-1:k
            a(n,m) = a(n,m) - a(n,k)*a(k,m)/a(k,k);
        end;
    end;
  end;
  disp(a);
  disp(b);
  
  x = sustitucion_hacia_atras(a, b);
end
\end{verbatim}

\noindent 
que corresponde a un sistema triangular superior que podemos solucionar utilizando sustitución hacia atrás.

\subsubsection{Ejemplo 1}

Dado el sistema lineal de ecuaciones, calcular el sistema Triangular Superior utilizando el algoritmo de Eliminación Gaussiana

\begin{equation*}
\left[
\begin{array}{r r r r}
    10  &  -1 &    2  &   1 \\
    -1   & 15 &   -3  &    1\\
     2   & -3  &    6  &  -3\\
     1   &  1  &   -3  &   7\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
  1 \\
  2 \\
  2 \\
  1\\
\end{array}
\right]
\end{equation*}

Primer paso $k=1$

\begin{equation*}
\left(
\begin{array}{cccc}
 10 & -1 & 2 & 1 \\
 \\
 0 & \frac{149}{10} & -\frac{14}{5} & \frac{11}{10} \\
 \\
 0 & -\frac{14}{5} & \frac{28}{5} & -\frac{16}{5} \\
 \\
 0 & \frac{11}{10} & -\frac{16}{5} & \frac{69}{10}
\end{array}
\right)\left(
\begin{array}{c}
 x_1 \\
 \\
 x_2 \\
 \\
 x_3 \\
 \\
 x_4
\end{array}
\right)\text{= } \left(
\begin{array}{c}
 1 \\
 \\
 \frac{21}{10} \\
 \\
 \frac{9}{5} \\
 \\
 \frac{9}{10}
\end{array}
\right)\end{equation*}

\begin{equation*}
\left[
\begin{array}{r r r r}
   10.0000   & -1.0000  &  2.0000  &  1.0000 \\
         0   & 14.9000 &  -2.8000  &  1.1000\\
         0   & -2.8000  &  5.6000  &  -3.2000\\
         0   & 1.1000 &  -3.2000   &  6.9000\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
    1.0000 \\
    2.1000\\
    1.8000\\
    0.9000\\
\end{array}
\right]
\end{equation*}

Segundo paso $k=2$

\begin{equation*}
\left(
\begin{array}{cccc}
 10 & -1 & 2 & 1 \\
 \\
 0 & \frac{149}{10} & -\frac{14}{5} & \frac{11}{10} \\
 \\
 0 & 0 & \frac{756}{149} & -\frac{446}{149} \\
 \\
 0 & 0 & -\frac{446}{149} & \frac{1016}{149}
\end{array}
\right)\left(
\begin{array}{c}
 x_1 \\
 \\
 x_2 \\
 \\
 x_3 \\
 \\
 x_4
\end{array}
\right)\text{= } \left(
\begin{array}{c}
 1 \\
 \\
 \frac{21}{10} \\
 \\
 \frac{327}{149} \\
 \\
 \frac{111}{149}
\end{array}
\right)
\end{equation*}


\begin{equation*}
\left[
\begin{array}{r r r r}
   10.0000  &  -1.0000  &  2.0000  &  1.0000\\
         0  & 14.9000 &  -2.8000   & 1.1000\\
         0    &     0 &   5.0738 &  -2.9933\\
         0    &     0  &  -2.9933 &   6.8188\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
    1.0000\\
    2.1000\\
    2.1946\\
    0.7450\\
\end{array}
\right]
\end{equation*}

Tercer paso $k=3$

\begin{equation*}
\left(
\begin{array}{cccc}
 10 & -1 & 2 & 1 \\
 \\
 0 & \frac{149}{10} & -\frac{14}{5} & \frac{11}{10} \\
 \\
 0 & 0 & \frac{756}{149} & -\frac{446}{149} \\
 \\
 0 & 0 & 0 & \frac{955}{189}
\end{array}
\right)\left(
\begin{array}{c}
 x_1 \\
 \\
 x_2 \\
 \\
 x_3 \\
 \\
 x_4
\end{array}
\right)\text{= } \left(
\begin{array}{c}
 1 \\
 \\
 \frac{21}{10} \\
 \\
 \frac{327}{149} \\
 \\
 \frac{257}{126}
\end{array}
\right)
\end{equation*}

\begin{equation*}
\left[
\begin{array}{r r r r}
   10.0000  & -1.0000 &   2.0000  &  1.0000\\
         0   &14.9000  & -2.8000 &  1.1000\\
         0   &      0  &  5.0738  & -2.9933\\
         0   &     0   &      0  &  5.0529\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
    1.0000\\
    2.1000\\
    2.1946\\
    2.0397\\
\end{array}
\right]
\end{equation*}

Sustitución hacia atras

Paso 1 $k=4$

\begin{eqnarray*}
x_4 & = & \frac{b_4}{a_{4,4}}\\
x_4 & = & \frac{2.0397}{5.0529} = 0.4037
\end{eqnarray*}

Paso 2 $k=3$
\begin{eqnarray*}
x_3 & = & \frac{b_3 - a_{3,4} x_4}{a_{3,3}}\\
x_3 & = & \frac{2.1946 - (-2.9933) (0.4037)}{5.0738} = 0.6707
\end{eqnarray*}

Paso 3 $k=2$

\begin{eqnarray*}
x_2 & = & \frac{b_2 - a_{2,3} x_3 - a_{2, 4} x_4}{a_{2,2}}\\
x_2 & = & \frac{2.1000 - (-2.8000 )(0.6707) - (1.1000) (0.4037)}{14.9000 } = 0.2372
\end{eqnarray*}

Paso 4 $k=1$
\begin{eqnarray*}
x_1 & = & \frac{b_1 - a_{1,2} x_2 - a_{1,3} x_3 - a_{1,4} x_4}{a_{1,1}}\\
x_1 & = & \frac{1.0000 - (-1.0000) (0.2372) - (2.0000) (0.6707) - (1.0000)(0.4037)}{10.0000} = - 0.0508
\end{eqnarray*}

La solución en Matlab es 

\begin{verbatim}
>> Eliminacion_Gaussiana([10 -1 2 1; -1 15 -3 1; 2 -3 6 -3; 1 1 -3 7], [1,2,2,1]')

   10.0000   -1.0000    2.0000    1.0000
         0   14.9000   -2.8000    1.1000
         0         0    5.0738   -2.9933
         0         0         0    5.0529

    1.0000    2.1000    2.1946    2.0397

ans =

   -0.0508
    0.2372
    0.6707
    0.4037
\end{verbatim}

\subsubsection{Ejemplo 2}

Calcular el sistema triangular superior utilizando eliminación Gaussiana.

$$
5x + 2y + 1z = 3
$$
$$
2x + 3y - 3z = -10
$$
$$
1x - 3y + 2z = 4 \\
$$

Primer paso

$$
5x +      2y +      1z = 3
$$
$$
0  + (11/5)y-(17/5)z =-(56/5)
$$
$$
0  - (17/5)y +  (9/5)z = (17/5)
$$

Segundo paso

$$
5x +      2y +      1z   = 3
$$
$$
0x + (11/5)y-( 17/5)z =-(56/5)
$$
$$
0x -      0y -(38/11)z =-(153/11)
$$

La solución en Matlab es

\begin{verbatim}
>> Eliminacion_Gaussiana([5 2 1; 2 3 -3; 1 -3 2], [3, -10, 4])
    5.0000    2.0000    1.0000
         0    2.2000   -3.4000
         0         0   -3.4545

    3.0000
  -11.2000
  -13.9091


ans =

   -0.6579
    1.1316
    4.0263
\end{verbatim}

\subsubsection{Ejemplo 3}

Resolver el sistema de ecuaciones

$$
\left[
\begin{array}{rrr}
3 &  -1 &  -1\\
-1 &  1 &  0 \\
-1 &  0 &  1\\
\end{array}
\right]
\left[
\begin{array}{r}
  x \\
  y \\
  z \\
\end{array}
\right] = \left[
\begin{array}{r}
  0 \\
  1 \\
  1 \\
\end{array}
\right ]
$$

Aplicando Eliminación Gaussiana nos queda.

$$
\left[
\begin{array}{rrr}
3 &  -1 &  -1 \\
0 &  2/3 &  -1/3\\
0 &  0 &  1/2 \\
\end{array}
\right]
\left[
\begin{array}{r}
  x \\
  y \\
  z \\
\end{array}
\right] = \left[
\begin{array}{r}
  0 \\
  1 \\
  3/2 \\
\end{array}
\right ]
$$

La solución en Matlab es

\begin{verbatim}

>> Eliminacion_Gaussiana([3 -1 -1; -1 1 0; -1 0 1], [0 1 1])
    3.0000   -1.0000   -1.0000
         0    0.6667   -0.3333
         0         0    0.5000

         0
    1.0000
    1.5000


ans =

    2.0000
    3.0000
    3.0000
\end{verbatim}

\subsubsection{Problema con la Eliminación Gaussiana}

Un ejemplo de sistema donde es necesario hacer un cambio de renglón por renglón para que tenga solución es el siguiente sistema. 

$$
\left[
\begin{array}{rrr}
1 & 2 & 6 \\
4 & 8 & -1 \\
-2 & 3 &  5\\
\end{array}
\right]
$$

Aplicando el primer paso de la eliminación gaussiana tenemos:

$$
\left[
\begin{array}{rrr}
1 & 2 & 6 \\
0 & 0 & -25 \\
0 & 7 & 17 \\
\end{array}
\right]
$$

\noindent
note, que aparece un cero en el elemento 22, lo cual nos da un sistema sin solución. Permutando los renglones 2 y 3 el sistema tiene solución.

$$
\left[
\begin{array}{rrr}
1 &  2 & 6 \\
-2 &  3 &  5 \\
4 &  8 & -1 \\
\end{array}
\right]
$$

La solución en Matlab tenemos
\begin{verbatim}
>> Eliminacion_Gaussiana([1 2 6; 4 8 -1; -2 3 5], [1 1 1])
     1     2     6
     0     0   -25
     0   NaN   Inf

     1
    -3
   Inf


ans =

   NaN
   NaN
   NaN
\end{verbatim}

Haciendo el cambio de renglones tenemos

\begin{verbatim}
>> Eliminacion_Gaussiana([1 2 6;  -2 3 5; 4 8 -1], [1 1 1])
     1     2     6
     0     7    17
     0     0   -25

     1
     3
    -3


ans =

    0.0057
    0.1371
    0.1200
\end{verbatim}


\section{Gauss-Jordan}

El método de Gauss-Jordan es una variación de la eliminación de Gauss. La principal diferencia consiste em que cuando una incógnita  se elimina en el método de Gauss-Jordan, esta es eliminada de todas las ecuaciones en lugar de hacerlo en la subsecuentes. Además todos los renglones se normalizan al dividirlos entre su elemento privote. De esta forma, el paso de eliminación genera una matriz identidad en vez de una matriz triangular. En consecuencia no es necesario usar la sustitución hacia atrás para obtener la solución. 

Comenzaremos planteando un sistema de 3 ecuaciones con 3 incógnitas dado por 

\begin{equation*}
\left[
\begin{array}{c|c|c}
a_{1,1} & a_{1,2} & a_{1,3} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} \\
\hline
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
\end{equation*}

Este sistema de ecuaciones lo podemos escribir como:

\begin{equation*}
Ax = I b
\end{equation*}

donde $I$ es la matriz identidad

\begin{equation*}
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
\end{equation*}

Lo cual da lugar a la matriz aumentada que hemos utilizado en el procedimiento de Gauss-Jordan

\begin{equation}
\left[
\begin{array}{ccc|ccc|c}
a_{1,1} & a_{1,2} & a_{1,3} & 1 & 0 & 0 & b_1 \\
a_{2,1} & a_{2,2} & a_{2,3} & 0 & 1 & 0 & b_2\\
a_{3,1} & a_{3,2} & a_{3,3} & 0 & 0 & 1 & b_3\\
\end{array}
\right]  \label{eq:Gauss_Jordan}
\end{equation}

\textbf{Primer paso}

Vamos a despejar la variable $x_1$ del sistema y la sustituimos en las otras dos ecuaciones

\begin{equation} 
\frac{b_1 - a_{1,2} x_2 - a_{1,3} x_3} {a_{1,1}} = x_1   \label{eq:eli_par_1}
\end{equation}

Sustituyendo en la ecuación 2

\begin{equation*}
a_{2,1}\left( \frac{b_1 - a_{1,2} x_2 - a_{1,3} x_3} {a_{1,1}} \right)+ a_{2,2} x_2 + a_{2,3} x_3 = b_2
\end{equation*}
\begin{equation}
0 x_1 + \left(\frac{a_{2,1}}{a_{1,1}}\right) b_1 + \left( a_{2,2} - \frac{a_{2,1}a_{1,2} }{a_{1,1}} \right) x_2 + \left( a_{2,3} \frac{a_{2,1} a_{1,3}} {a_{1,1}} \right) x_3 = b_2  \label{eq:eli_par_2}
\end{equation} 

Sustituyendo en la ecuación 3

\begin{equation*}
a_{3,1}\left( \frac{b_1 - a_{1,2} x_2 - a_{1,3} x_3} {a_{1,1}} \right)+ a_{3,2} x_2 + a_{3,3} x_3 = b_3
\end{equation*}

\begin{equation}
0 x_1 + \left(\frac{a_{3,1}}{a_{1,1}}\right) b_1 + \left( a_{3,2} - \frac{a_{3,1}a_{1,2} }{a_{1,1}} \right) x_2 + \left( a_{3,3} - \frac{a_{3,1} a_{1,3}} {a_{1,1}} \right) x_3 = b_3  \label{eq:eli_par_3}
\end{equation}

En forma matricial podemos escribir las ecuaciones \ref{eq:eli_par_1}, \ref{eq:eli_par_2} y \ref{eq:eli_par_3}

\begin{equation*}
\left[
\begin{array}{c|c|c}
1 & \frac{a_{1,2}}{a_{1,1}} & \frac{a_{1,3}}{a_{1,1}} \\
\hline
0 & \left( a_{2,2} - \frac{a_{2,1}a_{1,2} }{a_{1,1}} \right) &\left( a_{2,3} - \frac{a_{2,1} a_{1,3}} {a_{1,1}} \right) \\
\hline
0 & \left( a_{3,2} - \frac{a_{3,1}a_{1,2} }{a_{1,1}} \right) &\left( a_{3,3} - \frac{a_{3,1} a_{1,3}} {a_{1,1}} \right)  \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c|c|c}
\frac{1}{a_{1,1}} & 0 & 0 \\
\hline
-\frac{a_{2,1}}{a_{1,1}} & 1 & 0 \\
\hline
-\frac{a_{3,1}}{a_{1,1}} & 0 & 1 \\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
\end{equation*}

La matriz aumentada dada por (\ref{eq:Gauss_Jordan} ) tenemos:

\begin{equation*}
\left[
\begin{array}{c|c|c||c|c|c||c}
1 & \frac{a_{1,2}}{a_{1,1}}  & \frac{a_{1,3}}{a_{1,1}} & \frac{1}{a_{1,1}}  & 0 & 0 & \frac{1}{a_{1,1}}  b_1\\
\hline
0 & \left( a_{2,2} - \frac{a_{2,1}a_{1,2} }{a_{1,1}} \right) & \left( a_{2,3} - \frac{a_{2,1} a_{1,3}} {a_{1,1}} \right) & -\frac{a_{2,1}}{a_{1,1}}  & 1 & 0 & b_2 - \frac{a_{2,1}}{a_{1,1}} b_1\\
\hline
0 & \left( a_{3,2} - \frac{a_{3,1}a_{1,2} }{a_{1,1}} \right) &\left( a_{3,3} - \frac{a_{3,1} a_{1,3}} {a_{1,1}} \right)  & -\frac{a_{3,1}}{a_{1,1}} & 0 & 1 & b_3 - \frac{a_{3,1}}{a_{1,1}}  b_1\\
\end{array}
\right]  
\end{equation*}

\textbf{Segundo Paso}

Para un sistema equivalente 

\begin{equation}
\left[
\begin{array}{c|c|c}
1 & a'_{1,2} & a'_{1,3} \\
\hline
0 & a'_{2,2} & a'_{2,3}  \\
\hline
0 & a'_{3,2} & a'_{3,3}   \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c|c|c}
c'_{1,1} & 0 & 0 \\
\hline
c'_{2,1} & 1 & 0 \\
\hline
c'_{3,1} & 0 & 1 \\
\end{array}
\right]
\left[
\begin{array}{c}
  b'_1 \\
  b'_2 \\
  b'_3 \\
\end{array}
\right]
\end{equation}

Despejamos $x_2$ de la ecuación 2 

\begin{equation*}
x_2 = \frac{c'_{2,1} b'_1 + b'_2 - a'_{2,3} x_3}{a'_{2,2}}
\end{equation*}

Sustituyendo en la ecuación 1

\begin{equation*}
x1 + a'_{1,2} \left( \frac{c'_{2,1} b'_1 + b'_2 - a'_{2,3} x_3}{a'_{2,2}} \right)+ a'_{1,3} = c'_{1,1}b'_1
\end{equation*}


\begin{equation*}
x1 + 0 x_2 + \left (a'_{1,3}  - \frac{a'_{1,2}a'_{2,3}}{a'_{2,2}} x_3\right )= \left( c'_{1,1}b'_1 -  \frac{a'_{1,2}c'_{2,1}}{a'_{2,2}} \right) - \frac{a'_{2,1}} {a'_{2,2}} b'_2 
\end{equation*}

Sustituyendo en la ecuación 3

\begin{equation*}
0 x_1 + 0 x_2 + \left( a'_{3,3} - \frac{a'_{3,2} a'_{2,3}}{a'_{2,2}} \right) x_3 = \left ( c'_{3,1} - \frac{a'_{3,2} c'_{2,1}}{a'_{2,2}} \right) b'_1- \frac{a'_{3,2}}{a'_{2,2}} b'2 + b'3
\end{equation*}

En forma matricial

\begin{equation*}
\left[
\begin{array}{c|c|c}
1 & 0 & a'_{1,3} - \frac{a'_{1,2}a'_{2,3}}{a'_{2,2}}\\
\hline
0 & 1 & \frac{a'_{2,3}}{a'_{2,2}}  \\
\hline
0 & 0 & a'_{3,3} - \frac{a'_{3,2} a'_{2,3}}{a'_{2,2}}   \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c|c|c}
c'_{1,1} - \frac{a'_{1,2} c'_{2,1}}{a'_{2,2}} & \frac{a'_{1,2}}{a'_{2,2}} & 0 \\
\hline
\frac{c'_{2,1}}{a'_{2,2}} & \frac{1}{a'_{2,2}} & 0 \\
\hline
c'_{3,1} - \frac{a'_{3,2} c'_{2,1}}{a'_{2,2}} & -\frac{a'_{3,2}}{a'_{2,2}} & 1 \\
\end{array}
\right]
\left[
\begin{array}{c}
  b'_1 \\
  b'_2 \\
  b'_3 \\
\end{array}
\right]
\end{equation*}

\textbf{Tercer Paso}

Dado el sistema equivalente 

\begin{equation*}
\left[
\begin{array}{c|c|c}
1 & 0 & a''_{1,3}\\
\hline
0 & 1 &  a''_{2,3} \\
\hline
0 & 0 & a''_{3,3}   \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c|c|c}
c''_{1,1}  & c''_{1,2} & 0 \\
\hline
c''_{2,1} & c''_{2,2} & 0 \\
\hline
c''_{3,1}  & c''_{3,2} & 1 \\
\end{array}
\right]
\left[
\begin{array}{c}
  b''_1 \\
  b''_2 \\
  b''_3 \\
\end{array}
\right]
\end{equation*}

Despejamos $x_3$ de la ecuación 3 y sustituimos en la ecuaciones 1 y 2

\begin{equation*}
\left[
\begin{array}{c|c|c}
1 & 0 & 0\\
\hline
0 & 1 &  0 \\
\hline
0 & 0 & 1   \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c|c|c}
c''_{1,1}  & c''_{1,2} - \frac{a''_{1,3} c''_{3, 2}}{a''_{3,3}}& -\frac{a''_{1,3}}{a''_{3,3}} \\
\hline
c''_{2,1} & c''_{2,2} - \frac{a''_{2,3}}{a''_{3,3}}& -\frac{a''_{2,3}}{a''_{3,3}} \\
\hline
\frac{c''_{3,1}}{a''_{3,3}}  & \frac{c''_{3,2}}{a''_{3,3}} & \frac{1}{a''_{3,3}} \\
\end{array}
\right]
\left[
\begin{array}{c}
  b''_1 \\
  b''_2 \\
  b''_3 \\
\end{array}
\right]
\end{equation*}
Como resultado tenemos que la matrix $c''$ es la inversa de nuestro sistema y $b''$ la solución del sistema de ecuaciones

\textbf{Resumen}

Dado el sistema 

\begin{equation*}
\left[
\begin{array}{ccc|c|c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

\textbf{Paso 1}

Dividimos la primer ecuación entre $a_{1,1}$

\begin{equation*}
\begin{array}{c}
\frac{1}{a_{1,1}}\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 1 por $a_{2,1}$ y la restamos a la ecuación 2

\begin{equation*}
\begin{array}{c}
a_{2,1}\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
1 & \frac{a_{1,2}}{a_{1,1}} & \frac{a_{1,3}}{a_{1,1}}  & \cdots &  \frac{a_{1,M}}{a_{1,1}} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots  &  a_{2,M}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots  &  a_{3,M}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
1 & \frac{a_{1,2}}{a_{1,1}} & \frac{a_{1,3}}{a_{1,1}} & \cdots & \frac{a_{1,M}}{a_{1,1}} \\
\hline
0 & a_{2,2} - a_{2,1}\left( \frac{a_{1,2}}{a_{1,1}} \right )& a_{2,3} - a_{2,1}\left(\frac{a_{1,3}}{a_{1,1}} \right) & \cdots &  a_{2,M} - a_{2,1} \left( \frac{a_{1,M}}{a_{1,1}} \right )\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots  & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 1 por $a_{3,1}$ y la restamos a la ecuación 3

\begin{equation*}
\begin{array}{c}
a_{3,1}\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
1 & \frac{a_{1,2}}{a_{1,1}} & \frac{a_{1,3}}{a_{1,1}} & \cdots  & \frac{a_{1,M}}{a_{1,1}} \\
\hline
0 & a_{2,2} - a_{2,1}\left( \frac{a_{1,2}}{a_{1,1}} \right )& a_{2,3} - a_{2,1}\left(\frac{a_{1,3}}{a_{1,1}} \right) & \cdots &  a_{2,M} - a_{2,1} \left( \frac{a_{1,M}}{a_{1,1}} \right )\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}


\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
1 & \frac{a_{1,2}}{a_{1,1}} & \frac{a_{1,3}}{a_{1,1}} & \cdots & \frac{a_{1,M}}{a_{1,1}} \\
\hline
0 & a_{2,2} - a_{2,1}\left( \frac{a_{1,2}}{a_{1,1}} \right )& a_{2,3} - a_{2,1}\left(\frac{a_{1,3}}{a_{1,1}} \right) & \cdots &  a_{2,M} - a_{2,1} \left( \frac{a_{1,M}}{a_{1,1}} \right )\\
\hline
0   & a_{3,2} - a_{3,1} \left(\frac{a_{1,2}}{a_{1,1}}\right )  & a_{3,3} -a_{3,1} \left(\frac{a_{1,3}}{a_{1,1}} \right) & \cdots & a_{3,M}-  a_{3,1} \left(\frac{a_{1,M}}{a_{1,1}} \right )\\
\end{array}
\right]  
\end{equation*}

\textbf{Paso 2}

Dividimos la segunda ecuación entre $a_{2,2}$

\begin{equation*}
\begin{array}{c}
\\
\frac{1}{a_{2,2}}\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
\frac{a_{2,1}}{a_{2,2}}  & 1 & \frac{a_{2,3}}{a_{2,2}} & \cdots & \frac{a_{2,M}}{a_{2,2}}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 2 por $a_{1,2}$ y la restamos a la primer ecuación

\begin{equation*}
\begin{array}{c}
\\
a_{1,2}\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
\frac{a_{2,1}}{a_{2,2}}  & 1 & \frac{a_{2,3}}{a_{2,2}} & \cdots & \frac{a_{2,M}}{a_{2,2}}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,2} \left( \frac{a_{2,1}}{a_{2,2}} \right )& 0 & a_{1,3}- a_{1,2} \left(\frac{a_{2,3}}{a_{2,2}} \right) & \cdots & a_{1,M} - a_{1,2} \left(\frac{a_{2,M}}{a_{2,2}} \right )\\
\hline
\frac{a_{2,1}}{a_{2,2}}  & 1 & \frac{a_{2,3}}{a_{2,2}} & \cdots & \frac{a_{2,M}}{a_{2,2}}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 2 por $a_{3,2}$ y se la restamos a la ecuación 3

\begin{equation*}
\begin{array}{c}
\\
a_{3,2}\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,2} \left( \frac{a_{2,1}}{a_{2,2}} \right )& 0 & a_{1,3}- a_{1,2} \left(\frac{a_{2,3}}{a_{2,2}} \right) & \cdots & a_{1,M} - a_{1,2} \left(\frac{a_{2,M}}{a_{2,2}} \right )\\
\hline
\frac{a_{2,1}}{a_{2,2}}  & 1 & \frac{a_{2,3}}{a_{2,2}} & \cdots & \frac{a_{2,M}}{a_{2,2}}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}  \\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,2} \left( \frac{a_{2,1}}{a_{2,2}} \right )& 0 & a_{1,3}- a_{1,2} \left(\frac{a_{2,3}}{a_{2,2}} \right) & \cdots & a_{1,M} - a_{1,2} \left(\frac{a_{2,M}}{a_{2,2}} \right )\\
\hline
\frac{a_{2,1}}{a_{2,2}}  & 1 & \frac{a_{2,3}}{a_{2,2}} & \cdots & \frac{a_{2,M}}{a_{2,2}}\\
\hline
a_{3,1} - a_{3,2} \left( \frac{a_{2,1}}{a_{2,2}} \right )& 0 & a_{3,3}- a_{3,2} \left(\frac{a_{2,3}}{a_{2,2}} \right) & \cdots & a_{3,1} - a_{3,2} \left(\frac{a_{2,M}}{a_{2,2}} \right )\\
\end{array}
\right]  
\end{equation*}

\textbf{Paso 3}

Dividimos la tercera ecuación entre $a_{3,3}$

\begin{equation*}
\begin{array}{c}
\\
\\
\frac{1}{a_{3,3}}\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
a_{3,1} & a_{3,2} & a_{3,3} & \cdots & a_{3,M}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
\frac{a_{3,1}}{a_{3,3}} & \frac{a_{3,2}}{a_{3,3}} & 1 & \cdots & \frac{a_{3,M}}{a_{3,3}}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 3 por $a_{1,3}$ y la restamos a la primer ecuación

\begin{equation*}
\begin{array}{c}
\\
\\
a_{1,3}\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,M} \\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
\frac{a_{3,1}}{a_{3,3}} & \frac{a_{3,2}}{a_{3,3}} & 1 & \cdots & \frac{a_{3,M}}{a_{3,3}}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,3} \left(\frac{a_{3,1}}{a_{3,3}}\right) & a_{1,2} - a_{1,3}\left(\frac{a_{3,2}}{a_{3,3}}\right)& 0 & \cdots & a_{1,M} - a_{1,3}\left(\frac{a_{3,M}}{a_{3,3}}\right)\\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
\frac{a_{3,1}}{a_{3,3}} & \frac{a_{3,2}}{a_{3,3}} & 1 & \cdots & \frac{a_{3,M}}{a_{3,3}}\\
\end{array}
\right]  
\end{equation*}

Multiplicamos la ecuación 3 por $a_{2,3}$ y la restamos a la segunda ecuación

\begin{equation*}
\begin{array}{c}
\\
\\
a_{2,3}\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,3} \left(\frac{a_{3,1}}{a_{3,3}}\right) & a_{1,2} - a_{1,3}\left(\frac{a_{3,2}}{a_{3,3}}\right)& 0 & \cdots & a_{1,M} - a_{1,3}\left(\frac{a_{3,M}}{a_{3,3}}\right)\\
\hline
a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,M}\\
\hline
\frac{a_{3,1}}{a_{3,3}} & \frac{a_{3,2}}{a_{3,3}} & 1 & \cdots & \frac{a_{3,M}}{a_{3,3}}\\
\end{array}
\right]  
\end{equation*}

\begin{equation*}
\begin{array}{c}
\\
\\
\\
\end{array}
\left[
\begin{array}{c|c|c||c||c}
a_{1,1} - a_{1,3} \left(\frac{a_{3,1}}{a_{3,3}}\right) & a_{1,2} - a_{1,3}\left(\frac{a_{3,2}}{a_{3,3}}\right)& 0 & \cdots & a_{1,M} - a_{1,3}\left(\frac{a_{3,M}}{a_{3,3}}\right)\\
\hline
a_{2,1} - a_{2,3} \left( \frac{a_{3,1}}{a_{3,3}} \right) & a_{2,2}- a_{2,3} \left(\frac{a_{3,2}}{a_{3,3}} \right) & 0 & \cdots & a_{2,M} - a_{2,3} \left(\frac{a_{3,M}}{a_{3,3}}\right)\\
\hline
\frac{a_{3,1}}{a_{3,3}} & \frac{a_{3,2}}{a_{3,3}} & 1 & \cdots & \frac{a_{3,M}}{a_{3,3}}\\
\end{array}
\right]  
\end{equation*}

Note la similitud del método con el método de sumas y restas. No olvidar que este método es equivalente al calculo parcial de la matriz inversa. El método se ilustra mejor con un ejemplo.

\subsection{Ejemplo 1}

Utilice la técnica de Gauss-Jordan para resolver el siguiente sistema de ecuaciones.

\begin{eqnarray*}
30 x_1 -  x_2 - 2 x_3 & = & 78\\
 x_1 + 70 x_2  - 3 x_3 & = & -193\\
3 x_1 - 2 x_2 + 100 x_3 & = & 714
\end{eqnarray*}

Primero exprese los coeficientes y el lado derecho como una matriz aumentada

\[
\left(
\begin{array}{rrr|rrr|r}
 30 & -1 & -2 & 1 & 0 & 0 & 78 \\
 1 & 70 & -3 & 0 & 1 & 0 & -193 \\
 3 & -2 & 100 & 0 & 0 & 1 & 714
\end{array}
\right)
\]

\textbf{Primer iteración}

Luego se normaliza el primer renglón, al dividirlo entre el elemento pivote, 30, se obtiene

\[
\left[
\begin{array}{rrr|rrr|r}
 1 & -\frac{1}{30} & -\frac{1}{15} & \frac{1}{30} & 0 & 0 & \frac{13}{5} \\
 &  & & & & &  \\
 1 & 70 & -3 & 0 & 1 & 0 & -193 \\
 &&&&&&\\
 3 & -2 & 100 & 0 & 0 & 1 & 714
\end{array}
\right]
\]

El término $x_1$ se puede eliminar del segundo renglón restando el primer renglón multiplicado por $1$ del segundo renglón. En forma similar restando el primer renglón multiplicando por $3$ eliminará el término $x_1$ del tercer renglón

\[
\left[
\begin{array}{rrr|rrr|r}
 1 & -\frac{1}{30} & -\frac{1}{15} & \frac{1}{30} & 0 & 0 & \frac{13}{5} \\
 &&&&&&\\
 0 & \frac{2101}{30} & -\frac{44}{15} & -\frac{1}{30} & 1 & 0 & -\frac{978}{5} \\
 &&&&&&\\
 0 & -\frac{19}{10} & \frac{501}{5} & -\frac{1}{10} & 0 & 1 & \frac{3531}{5}
\end{array}
\right]
\]

\[
\left[
\begin{array}{rrr|rrr|r}
 1. & -0.0333 & -0.0667 & 0.0333 & 0. & 0. & 2.6 \\
 0. & 70.0333 & -2.9333 & -0.0333 & 1. & 0. & -195.6 \\
 0. & -1.9 & 100.2 & -0.1 & 0. & 1. & 706.2
\end{array}
\right]
\]

\textbf{Segunda Iteración}

En seguida, se normaliza el segundo renglón dividiéndolo entre $\frac{2101}{30}$

\[
\left[
\begin{array}{rrr|rrr|r}
 1 & -\frac{1}{30} & -\frac{1}{15} & \frac{1}{30} & 0 & 0 & \frac{13}{5} \\
 &&&&&&\\
 0 & 1 & -\frac{8}{191} & -\frac{1}{2101} & \frac{30}{2101} & 0 & -\frac{5868}{2101} \\
 &&&&&&\\
 0 & -\frac{19}{10} & \frac{501}{5} & -\frac{1}{10} & 0 & 1 & \frac{3531}{5}
\end{array}
\right]
\]

\[
\left[
\begin{array}{rrr|rrr|r}
 1. & -0.0333 & -0.0667 & 0.0333 & 0. & 0. & 2.6 \\
 0. & 1. & -0.0419 & -0.0005 & 0.0143 & 0. & -2.793 \\
 0. & -1.9 & 100.2 & -0.1 & 0. & 1. & 706.2
\end{array}
\right]
\]   

Al reducir los términos $x_2$ de las ecuaciones 1 y 2 tenemos
\[
\left[
\begin{array}{rrr|rrr|r}
 1 & 0 & -\frac{13}{191} & \frac{70}{2101} & \frac{1}{2101} & 0 & \frac{5267}{2101} \\
 &&&&&&\\
 0 & 1 & -\frac{8}{191} & -\frac{1}{2101} & \frac{30}{2101} & 0 & -\frac{5868}{2101} \\
 &&&&&&\\
 0 & 0 & \frac{19123}{191} & -\frac{212}{2101} & \frac{57}{2101} & 1 & \frac{1472577}{2101}
\end{array}
\right]
\]

\[
\left[
\begin{array}{rrr|rrr|r}
 1. & 0. & -0.0681 & 0.0333 & 0.0005 & 0. & 2.5069 \\
 0. & 1. & -0.0419 & -0.0005 & 0.0143 & 0. & -2.793 \\
 0. & 0. & 100.12 & -0.1009 & 0.0271 & 1. & 700.893
\end{array}
\right]
\]  


\textbf{Tercer iteración}

El tercer renglón se normaliza entonces al dividirlo entre $\frac{19123}{191}$

\[
\left[
\begin{array}{rrr|rrr|r}
 1 & 0 & -\frac{13}{191} & \frac{70}{2101} & \frac{1}{2101} & 0 & \frac{5267}{2101} \\
 &&&&&&\\
 0 & 1 & -\frac{8}{191} & -\frac{1}{2101} & \frac{30}{2101} & 0 & -\frac{5868}{2101} \\
 &&&&&&\\
 0 & 0 & 1 & -\frac{212}{210353} & \frac{57}{210353} & \frac{191}{19123} & \frac{1472577}{210353}
\end{array}
\right]
\]  

\[
\left[
\begin{array}{ccccccc}
 1. & 0. & -0.0681 & 0.0333 & 0.0005 & 0. & 2.5069 \\
 0. & 1. & -0.0419 & -0.0005 & 0.0143 & 0. & -2.793 \\
 0. & 0. & 1. & -0.001 & 0.0003 & 0.01 & 7.0005
\end{array}
\right]
\]

Por último, los términos $x_3$ se pueden reducir de la primera y segunda ecuación para obtener

 \[
\left[
\begin{array}{rrr|rrr|r}
 1 & 0 & 0 & \frac{538}{16181} & \frac{8}{16181} & \frac{1}{1471} & \frac{48274}{16181} \\
 &&&&&& \\
 0 & 1 & 0 & -\frac{109}{210353} & \frac{3006}{210353} & \frac{8}{19123} & -\frac{525828}{210353} \\
 &&&&&& \\
 0 & 0 & 1 & -\frac{212}{210353} & \frac{57}{210353} & \frac{191}{19123} & \frac{1472577}{210353}
\end{array}
\right]
\]  

\[
\left[
\begin{array}{ccccccc}
 1. & 0. & 0. & 0.0332 & 0.0005 & 0.0007 & 2.9834 \\
 0. & 1. & 0. & -0.0005 & 0.0143 & 0.0004 & -2.4997 \\
 0. & 0. & 1. & -0.001 & 0.0003 & 0.01 & 7.0005
\end{array}
\right]
\]    

\subsection{Ejemplo 2}

Dado el sistema lineal de ecuaciones calcular la solución utilizando el método de Gauss-Jordan

\begin{equation*}
\left[
\begin{array}{r r r }
     1   &   1 &    2  \\
     1   &   3 &   -6  \\
     2   &  -1  &   2  \\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right]
=
\left[
\begin{array}{r }
  -1 \\
  7 \\
  0 \\
\end{array}
\right]
\end{equation*}

La matriz aumentada para este sistema es:

\begin{equation*}
\left[
\begin{array}{r r r | r }
     1   &   1 &    2  & -1\\
     1   &   3 &   -6  &  7\\
     2   &  -1  &   2  & 0 \\
\end{array}
\right] 
\end{equation*}

\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r r r | r }
     1   &   1 &    2  & -1\\
     0   &   2 &   -8  &  8\\
     0   &  -3  &  -2  & 2 \\
\end{array}
\right] 
\end{equation*}

\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r r r | r }
     1   &   0 &    6  & -5\\
     0   &   1 &   -4  &  4\\
     0   &  0  &  -14  &  14 \\
\end{array}
\right] 
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r r r | r }
     1   &   0 &  0  & 1\\
     0   &   1 &  0  &  0\\
     0   &  0  &  1  &  -1 \\
\end{array}
\right] 
\end{equation*}


\subsection{Ejemplo 3}

Dado el sistema lineal de ecuaciones, calcular la solución y la matriz inversa utilizando el método de Gauss-Jordan

\begin{equation*}
\left[
\begin{array}{r r r r}
    10  &  -1 &    2  &   1 \\
    -1   & 15 &   -3  &    1\\
     2   & -3  &    6  &  -3\\
     1   &  1  &   -3  &   7\\
\end{array}
\right] 
\left[
\begin{array}{r }
  x_1 \\
  x_2 \\
  x_3 \\
  x_4\\
\end{array}
\right]
=
\left[
\begin{array}{r }
  1 \\
  2 \\
  2 \\
  1\\
\end{array}
\right]
\end{equation*}

El sistema inicial aumentado es
\[
\left[
\begin{array}{rrrr|rrrr|r}
    10 &   -1   &  2  &   1  &   1  &   0  &   0   &  0   &  1\\
    -1 &   15   & -3   &  1  &   0  &   1  &   0   &  0  &   2\\
     2 &   -3   &  6   & -3   &  0   &  0   &  1   &  0  &   2\\
     1  &   1  &  -3   &  7   &  0   &  0   &  0   &  1  &   1\\
 \end{array}
\right]
\]
    
\textbf{Primer iteración}

\[
\left(
\begin{array}{rrrr|rrrr|r}
 1 & -\frac{1}{10} & \frac{1}{5} & \frac{1}{10} & \frac{1}{10} & 0 & 0 & 0 & \frac{1}{10} \\
 &&&&&&&\\
 0 & \frac{149}{10} & -\frac{14}{5} & \frac{11}{10} & \frac{1}{10} & 1 & 0 & 0 & \frac{21}{10} \\
 &&&&&&&\\
 0 & -\frac{14}{5} & \frac{28}{5} & -\frac{16}{5} & -\frac{1}{5} & 0 & 1 & 0 & \frac{9}{5} \\
 &&&&&&&\\
 0 & \frac{11}{10} & -\frac{16}{5} & \frac{69}{10} & -\frac{1}{10} & 0 & 0 & 1 & \frac{9}{10}
\end{array}
\right)
\]

\[
\left[
\begin{array}{rrrr|rrrr|r}
    1.0000 &  -0.1000  &  0.2000  &  0.1000  &  0.1000  &       0   &      0  &       0  &  0.1000\\
         0     &  14.9000  & -2.8000 &   1.1000 &   0.1000 &   1.0000    &     0  &       0  &  2.1000\\
         0     &  -2.8000  &  5.6000  & -3.2000  & -0.2000  &       0   & 1.0000   &      0   & 1.8000\\
         0     &   1.1000  & -3.2000  &  6.9000  & -0.1000   &      0    &     0  &  1.0000  &  0.9000\\
\end{array}
\right]
\]

\textbf{Segunda iteración}

\[
\left[
\begin{array}{rrrr|rrrr|r}
 1 & 0 & \frac{27}{149} & \frac{16}{149} & \frac{15}{149} & \frac{1}{149} & 0 & 0 & \frac{17}{149} \\
  &&&&&&&\\
 0 & 1 & -\frac{28}{149} & \frac{11}{149} & \frac{1}{149} & \frac{10}{149} & 0 & 0 & \frac{21}{149} \\
  &&&&&&&\\
 0 & 0 & \frac{756}{149} & -\frac{446}{149} & -\frac{27}{149} & \frac{28}{149} & 1 & 0 & \frac{327}{149} \\
  &&&&&&&\\
 0 & 0 & -\frac{446}{149} & \frac{1016}{149} & -\frac{16}{149} & -\frac{11}{149} & 0 & 1 & \frac{111}{149}
\end{array}
\right]
\]

\[
\left[
\begin{array}{rrrr|rrrr|r}
    1.0000     &    0 &   0.1812 &   0.1074  &  0.1007  &  0.0067     &    0    &     0  &  0.1141\\
         0   & 1.0000 &  -0.1879 &   0.0738  &  0.0067  &  0.0671     &    0    &     0  &  0.1409\\
         0   &      0  &  5.0738  & -2.9933 &  -0.1812  &  0.1879  &  1.0000     &    0   & 2.1946\\
         0    &     0 &  -2.9933 &  6.8188  & -0.1074 &   -0.0738   &      0  &  1.0000  &  0.7450\\
\end{array}
\right]
\]

\textbf{Tercer iteración}

\[
\left(
\begin{array}{rrrr|rrrr|r}
 1 & 0 & 0 & \frac{3}{14} & \frac{3}{28} & 0 & -\frac{1}{28} & 0 & \frac{1}{28} \\
   &&&&&&&\\
 0 & 1 & 0 & -\frac{1}{27} & 0 & \frac{2}{27} & \frac{1}{27} & 0 & \frac{2}{9} \\
   &&&&&&&\\
 0 & 0 & 1 & -\frac{223}{378} & -\frac{1}{28} & \frac{1}{27} & \frac{149}{756} & 0 & \frac{109}{252} \\
   &&&&&&&\\
 0 & 0 & 0 & \frac{955}{189} & -\frac{3}{14} & \frac{1}{27} & \frac{223}{378} & 1 & \frac{257}{126}
\end{array}
\right)
\]

\[
\left[
\begin{array}{rrrr|rrrr|r}
    1.0000    &     0     &    0  &  0.2143 &   0.1071   & 0.0000 & -0.0357   &      0   & 0.0357\\
         0   & 1.0000     &    0  & -0.0370   &      0  &  0.0741   & 0.0370      &   0  &  0.2222 \\
         0   &      0   & 1.0000  & -0.5899 &  -0.0357  &  0.0370  &  0.1971   &      0  &  0.4325\\
         0    &     0   &      0  &  5.0529  & -0.2143  &  0.0370  &  0.5899   & 1.0000  &  2.0397\\
\end{array}
\right]
\]

\textbf{Cuarta iteración}

\[
\left[
\begin{array}{rrrr|rrrr|r}
 1 & 0 & 0 & 0 & \frac{111}{955} & -\frac{3}{1910} & -\frac{58}{955} & -\frac{81}{1910} & -\frac{97}{1910} \\
    &&&&&&&\\
 0 & 1 & 0 & 0 & -\frac{3}{1910} & \frac{71}{955} & \frac{79}{1910} & \frac{7}{955} & \frac{453}{1910} \\
    &&&&&&&\\
 0 & 0 & 1 & 0 & -\frac{58}{955} & \frac{79}{1910} & \frac{254}{955} & \frac{223}{1910} & \frac{1281}{1910} \\
    &&&&&&&\\
 0 & 0 & 0 & 1 & -\frac{81}{1910} & \frac{7}{955} & \frac{223}{1910} & \frac{189}{955} & \frac{771}{1910}
\end{array}
\right]
\]

\[
\left[
\begin{array}{rrrr|rrrr|r}
    1.0000    &     0    &     0    &     0  &  0.1162 &  -0.0016 &  -0.0607 &  -0.0424 &  -0.0508\\
         0  &  1.0000    &     0    &     0  & -0.0016 &   0.0743 &   0.0414 &   0.0073 &   0.2372\\
         0   &      0  &  1.0000    &     0  & -0.0607  &  0.0414  &  0.2660  &  0.1168 &   0.6707\\
         0   &      0   &      0   & 1.0000  & -0.0424  &  0.0073  &  0.1168  &  0.1979  &  0.4037\\
\end{array}
\right]
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Inversa de Shiplay}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Consideremos el siguiente sistema de ecuaciones

$$
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{array}
\right]
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
$$

de la ecuación $I$ despejamos el valor de la variable $x_1$

\begin{equation*}
x_1 = \frac{b_1 - a_{1,2}x_2 - a_{1,3} x_3}{a_{1,1}}
\end{equation*}
 
Reorganizando la ecuación tenemos
 
\begin{equation*}
\frac{b_1}{a_{1,1}} - \frac{a_{1,2}}{a_{1,1}}x_2 - \frac{a_{1,3}}{a_{1,1}}x_3 = x_1
\end{equation*}

Sustituimos el valor de $x_1$ en la ecuación $II$:

\begin{equation*}
\frac{(b_1 - a_{1,2}x_2 - a_{1,3}x_3) a_{2,1}}{a_{1,1}} + a_{2,2}x_2 + a_{2,3}x_3 = b_2
\end{equation*}

\begin{equation*}
\frac{a_{2,1}} {a_{1,1}} b_1 + \left(a_{2,2} - \frac{a_{2,1} a_{1,2}}{a_{1,1}} \right) x_2 + \left(a_{2,3} - \frac{a_{2,1} a_{1,3}}{a_{1,1}} \right)x_3 = b_2
\end{equation*}

y finalmente en $III$ tenemos :

\begin{equation*}
\frac{a_{3,1}} {a_{1,1}} b_1 + \left(a_{3,2} - \frac{a_{3,1} a_{1,2}}{a_{1,1}} \right) x_2 + \left(a_{3,3} - \frac{a_{3,1} a_{1,3}}{a_{1,1}} \right)x_3 = b_3
\end{equation*}
 

En forma matricial estas ecuaciones lucen como:

$$
\left[
\begin{array}{c|c|c}
1/a_{1,1} & a_{1,2}/a_{1,1} & a_{1,3}/a_{1,1} \\
\hline
a_{2,1}/a_{1,1} & a_{2,2} - a_{2,1} a_{1,2}/a_{1,1}  & a_{2,3} - a_{2,1} a_{1,3}/a_{1,1} \\
\hline
a_{3,1}/a_{1,1}  & a_{3,2} - a_{3,1} a_{1,2}/a_{1,1} & a_{3,3} - a_{3,1} a_{1,3}/a_{1,1} \\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  x_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
$$ 
 
Lo que nos da un sistema equivalente como el siguiente.
 
$$
\left[
\begin{array}{ccc}
a’_{1,1} & a’_{1,2} & a’_{1,3}\\
a’_{2,1} & a’_{2,2} & a’_{2,3}\\
a’_{3,1} & a’_{3,2} & a’_{3,3}\\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  x_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right]
$$

Si repetimos el procedimiento suponiendo que despejamos $x_2$ y sustituimos en las otras dos ecuaciones tendremos un sistema.

$$
\left[
\begin{array}{c|c|c}
a'_{1,1}-a'_{1,2} a'_{2,1}/a'_{2,2} & a'_{1,2}/a'_{2,2} & a'_{1,3}-a'_{1,2}a'_{2,3}/a'_{2,2}\\
\hline
-a'_{2,1}/a'_{2,2} & 1/a'_{2,2} & -a'_{2,3}/a'_{2,2} \\
\hline
a'_{3,1}-a'_{3,2}a'_{2,1}/a'_{2,2} & a'_{3,2}/a'_{2,2} & a'_{3,3}-a'_{3,2}a'_{2,3}/a'_{2,2}\\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  b_3 \\
\end{array}
\right]
$$

Lo que nos da un sistema equivalente como el siguiente.

$$
\left[
\begin{array}{ccc}
a''_{1,1} & a''_{1,2} & a''_{1,3}\\
a''_{2,1} & a''_{2,2} & a''_{2,3}\\
a''_{3,1} & a''_{3,2} & a''_{3,3}\\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  x_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  b_3 \\
\end{array}
\right]
$$ 

Finalmente si despejamos $x_3$ de la ecuación $III$  y sustituimos en las otras tenemos;

$$
\left[
\begin{array}{c|c|c}
a''_{1,1} - a''_{1,3}a''_{3,1}/a''_{3,3} & a''_{1,2} - a''_{1,3}a''_{3,2}/a''_{3,3}& a''_{1,3}/a''_{3,3} \\
\hline
a''_{2,1} - a''_{2,3}a''_{3,1}/a''_{3,3} & a''_{2,2} - a''_{2,3}a''_{3,2} /a''_{3,3} & a''_{2,3}/a''_{3,3}\\
\hline
a''_{3,1}/a''_{3,3} & a''_{3,2}/a''_{2,2}/a''_{3,3} & 1/a''_{3,3}\\
\end{array}
\right]
\left[
\begin{array}{c}
  b_1 \\
  b_2 \\
  b_3 \\
\end{array}
\right] =
\left[
\begin{array}{c}
  x_1 \\
  x_2 \\
  x_3 \\
\end{array}
\right]
$$

El procedimiento para el calculo de la inversa por el método de Shiplay es:

Dado un valor de $k= 1,2,3, ..., N$ al que llamamos pivote hacemos:
\begin{itemize}
\item $a_{n,m} = a_{n,m} – a_{n,k} a_{k,m}/a_{k,k}$ para todos los valores de $n$ y $m$ diferentes de $k$
\item $a_{n,k} = a_{n,k}/a_{k,k}$ para todos los $n$ diferentes de $k$
\item $a_{k,m} = a_{k,m}/a_{k,k}$ para todos los $m$ diferentes de $k$
\item $a_{k,k} = 1/a_{k,k}$
\end{itemize}

La implementación en Matlab queda

\begin{verbatim}
function a = Inversa_Shiplay(A)
   a=A;
   N=length(a);
   for k=1:N
       for n=1:N
           for m=1:N
               if ((m~=k)&(n~=k))
                   a(n,m)=a(n,m)-(a(n,k)*a(k,m))/a(k,k);
               end;
           end;
       end;
       for n=1:N
           if (n~=k)
               a(n,k)=a(n,k)/a(k,k);
               a(k,n)=-(a(k,n)/a(k,k));
           end;
       end;
       a(k,k)=1/a(k,k);
   end;
end
\end{verbatim}
 
\subsection{Ejemplo:}

Determinar la matriz inversa de

 \[
\left[
\begin{array}{ccc}
 3 & -1 & -1 \\
-1 &  1 &  0 \\
-1 &  0 &  1\\
\end{array}
\right]
\]
 
Primer paso:

 \[
\left[
\begin{array}{ccc}
 1/3 &  1/3 &  1/3 \\
-1/3 &  2/3 & -1/3\\
-1/3 &-1/3 &  2/3\\
 \end{array}
\right]
\]

Segundo paso:

 \[
\left[
\begin{array}{ccc}
1/2 & 1/2 & 1/2\\
1/2 & 3/2 & 1/2\\
-1/2 & -1/2 & 1/2\\
 \end{array}
\right]
\]
 
Tercer paso:

 \[
\left[
\begin{array}{ccc}
1 & 1 & 1\\
1 & 2 & 1\\
1 & 1 & 2\\
 \end{array}
\right]
\]

La solución en Matlab es:

\begin{verbatim}
>> Inversa_Shiplay([ 3 -1 -1;-1  1  0 ; -1  0  1])

ans =

    1.0000    1.0000    1.0000
    1.0000    2.0000    1.0000
    1.0000    1.0000    2.0000
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Sistemas dispersos y estrategias para conservarla}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Consideremos un sistema de ecuaciones $Ax = b$ como la siguiente

\[
A = 
\left[
\begin{array}{ccccc}
    10 &   -1  &  -1  &  -1  &   1 \\
     1  &  4   &  0   &  0   &  0\\
     1  &  0   &  4   &  0   &  0\\
     1  &   0  &   0  &   4  &   0\\
     1  &   0  &   0  &   0   &  3
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2\\
x_3\\
x_4\\
x_5
\end{array}
\right]
= 
\left[
\begin{array}{c}
1\\
2\\
3\\
4\\
5
\end{array}
\right]
\]

Para ésta matriz el total de elementos es 25 de los cuales 13 son diferentes de cero y 12 son igual cero. A este tipo de matrices o sistemas de ecuaciones se les conoce como dispersos o ralos dado que tienen una cantidad de ceros igual o mayor que los diferentes de cero.

Para este tipo de sistema cuando calculamos la solución por lo métodos dados tenemos


Descomposición LU
\begin{verbatim}
>> Factorizacion(A)

ans =

   10.0000   -1.0000   -1.0000   -1.0000    1.0000
    0.1000    4.1000    0.1000    0.1000   -0.1000
    0.1000    0.0244    4.0976    0.0976   -0.0976
    0.1000    0.0244    0.0238    4.0952   -0.0952
    0.1000    0.0244    0.0238    0.0233    2.9070
\end{verbatim}

Eliminación Gaussiana

\begin{verbatim}
>> Eliminacion_Gaussiana(A,b)
   10.0000   -1.0000   -1.0000   -1.0000    1.0000
         0    4.1000    0.1000    0.1000   -0.1000
         0         0    4.0976    0.0976   -0.0976
         0         0         0    4.0952   -0.0952
         0         0         0         0    2.9070
\end{verbatim}

Inversa de Shiplay

\begin{verbatim}
>> Inversa_Shiplay(A)

ans =

    0.0960    0.0240    0.0240    0.0240   -0.0320
   -0.0240    0.2440   -0.0060   -0.0060    0.0080
   -0.0240   -0.0060    0.2440   -0.0060    0.0080
   -0.0240   -0.0060   -0.0060    0.2440    0.0080
   -0.0320   -0.0080   -0.0080   -0.0080    0.3440
\end{verbatim}

Note que en todos los casos la dispersidad se pierde y tenemos una matriz con números diferentes de cero ahí donde existían ceros.

En este caso una opción es el manejo de Matrices Ralas del Matlab para no almacenar los valores iguales a cero y utilizar métodos iterativo como los de la siguiente sección.

La forma de representar la matriz dispersa en Matlab es

\begin{verbatim}
A = sparse(5,5);
A(1,1) = 10;
A(1,2) = -1;
A(1,3) = -1;
A(1,4) = -1;
A(1,5) = 1;
A(2,1) = 1;
A(3,1) = 1;
A(4,1) = 1;
A(5,1) = 1;
A(2,2) = 4;
A(3,3) = 4;
A(4,4) = 4;
A(5,5) = 3;
disp(A);
\end{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Método iterativo de Jacobi  y Gauss-Seidel}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Método de Jacobi}

Consideremos el siguiente sistema de ecuaciones.

\begin{equation*}
\left[
\begin{array}{ccc}
a_{1,1} & a_{1,2} & a_{1,3}\\
a_{2,1} & a_{2,2} & a_{2,3} \\	
a_{3,1} & a_{3,2} & a_{3,3}
\end{array}
\right]
\left[
\begin{array}{c}
x_1 \\
x_2 \\	
x_3
\end{array}
\right]
= 	
\left[
\begin{array}{c}
b_1 \\
b_2 \\
b_3
\end{array}
\right]
\end{equation*}

Vamos a representar cada una de las variables en términos de ellas mismas.

\begin{eqnarray*} 
x_1 & = & \frac{ b_1 - a_{1,2} x_2 - a_{1,3} x_3}{a_{1,1}} \\
x_2 & = & \frac{ b_2 - a_{2,1} x_1 - a_{2,3} x_3}{a_{2,2}} \\
x_3 & = & \frac{ b_3 - a_{3,1} x_1 - a_{3,2} x_2}{a_{3,3}} \\
\end{eqnarray*}

Lo cual nos sugiere el siguiente esquema iterativo de solución.

\begin{eqnarray*} 
x_1^{(t+1)} & = & \frac{ b_1 - a_{1,2} x_2^{(t)} - a_{1,3} x_3^{(t)}}{a_{1,1}} \\
x_2^{(t+1)} & = & \frac{ b_2 - a_{2,1} x_1^{(t)} - a_{2,3} x_3^{(t)}}{a_{2,2}} \\
x_3^{(t+1)} & = & \frac{ b_3 - a_{3,1} x_1^{(t)} - a_{3,2} x_2^{(t)}}{a_{3,3}} \\
\end{eqnarray*}
 
En general podemos escribir como

\begin{equation*} 
x_k^{(t+1)} = \frac{b_k – \sum_{l=1, l \neq k}^N a_{k,l}x_l^{(t)}}{a_{k,k}} 
\end{equation*}

La implementación en Matlab es

\begin{verbatim}
function y = Jacobi(A, x, b)

N = length(x);
y = zeros(N,1);

for iter=1:100000
    for k = 1:N
        suma =0;
        for l= 1:N
            if k ~= l
                suma = suma + A(k,l)*x(l);
            end;
        end;
        y(k) = (b(k) - suma)/A(k,k);
    end;
    if sqrt(norm(x-y)) < 1e-6 
        break;
    else
       x = y;
    end;
end;
\end{verbatim}

\subsection{Algoritmo iterativo de Gauss-Seidel}

El cambio que debemos hacer respecto al de Jacobi, es que las variables nuevas son utilizadas una vez que se realiza el cálculo de ellas así, para un sistema de tres ecuaciones tendremos:

\begin{eqnarray*} 
x_1^{(t+1)} & = & \frac{ b_1 - a_{1,2} x_2^{(t)} - a_{1,3} x_3^{(t)}}{a_{1,1}} \\
x_2^{(t+1)} & = & \frac{ b_2 - a_{2,1} x_1^{(t+1)} - a_{2,3} x_3^{(t)}}{a_{2,2}} \\
x_3^{(t+1)} & = & \frac{ b_3 - a_{3,1} x_1^{(t+1)} - a_{3,2} x_2^{(t+1)}}{a_{3,3}} \\
\end{eqnarray*}

La implementación en Matlab es

\begin{verbatim}
function y = Gauss_Seidel(A, x, b)

N = length(x);
y = zeros(N,1);

for iter=1:100000
    for k = 1:N
        suma =0;
        for l= 1:N
            if k ~= l
                suma = suma + A(k,l)*x(l);
            end;
        end;
        x(k) = (b(k) - suma)/A(k,k);
    end;
    if sqrt(norm(x-y)) < 1e-6 
        break;
    else
        y=x;
    end;
end;
\end{verbatim}

\subsubsection{Ejemplo}

Resolver el siguiente sistema de ecuaciones utilizando el método de Jacobi y comparar con el método de Gauss-Seidel.

\begin{equation*}
\left[
\begin{array}{ccc}
 4  &  -1 &  1  \\
 4  &  -8 &  1 \\
-2  &   1 &  5
\end{array}
\right]
\left[
\begin{array}{c}
x_1 \\
x_2 \\	
x_3
\end{array}
\right]
= 	
\left[
\begin{array}{c}
  7 \\
-21 \\
 15
\end{array}
\right]
\end{equation*}

Las primeras 20 iteraciones del algoritmo de Jacobi son :

\begin{center}
\begin{tabular}{|r|r|r|r|}
\hline
$k$ & $x_1$ & $x_2$ & $x_3$ \\
\hline
1 & 1.0000 & 2.0000 & 2.0000 \\
2 & 1.7500 & 3.3750 & 3.0000 \\
3 & 1.8438 & 3.8750 & 3.0250\\
4 & 1.9625 & 3.9250 & 2.9625\\
5 & 1.9906 & 3.9766 & 3.0000\\
6 & 1.9941 & 3.9953 & 3.0009\\
7 & 1.9986 & 3.9972 & 2.9986\\
8 & 1.9996 & 3.9991 & 3.0000\\
9 & 1.9998 & 3.9998 & 3.0000\\
10 & 1.9999 & 3.9999 & 2.9999 \\
11 & 2.0000 & 4.0000 & 3.0000 \\
12 & 2.0000 & 4.0000 & 3.0000 \\
\hline
\end{tabular}
\end{center}

Para correr hacer 

\begin{verbatim}
>> Jacobi([4 -1 1; 4 -8 1; -2 1 5], [1,2,2]', [7,-21, 15]')

ans =

    2.0000
    4.0000
    3.0000
\end{verbatim}

La solución utilizando Gauss-Seidel es :

\begin{center}
\begin{tabular}{|r|r|r|r|}
\hline
$k$ & $x_1$ & $x_2$ & $x_3$ \\
\hline
1 & 1.0000 & 2.0000 & 2.0000 \\
2 & 1.7500 & 3.7500 & 2.9500 \\
3 & 1.9500 & 3.9688 & 2.9863 \\
4 & 1.9956 & 3.9961 & 2.9990 \\
5 & 1.9993 & 3.9995 & 2.9998 \\
6 & 1.9999 & 3.9999 & 3.0000 \\
7 & 2.0000 & 4.0000 & 3.0000 \\
8 & 2.0000 & 4.0000 & 3.0000 \\
\hline
\end{tabular}
\end{center}
 
Note que Gauss-Seidel requiere de 7 iteraciones mientras Jacobi de 11, para convergir. Para correr hacer

\begin{verbatim}
>> Gauss_Seidel([4 -1 1; 4 -8 1; -2 1 5], [1,2,2]', [7,-21, 15]')

ans =

    2.0000
    4.0000
    3.0000
\end{verbatim}

\subsection{Ejemplo matrices dispersas} 
Resolver el sistema de ecuaciones utilizando matrices dispersas y los métodos de Jacobi y Gauss-Seidel

\[
A = 
\left[
\begin{array}{ccccc}
    10 &   -1  &  -1  &  -1  &   1 \\
     1  &  4   &  0   &  0   &  0\\
     1  &  0   &  4   &  0   &  0\\
     1  &   0  &   0  &   4  &   0\\
     1  &   0  &   0  &   0   &  3
\end{array}
\right]
\left[
\begin{array}{c}
x_1\\
x_2\\
x_3\\
x_4\\
x_5
\end{array}
\right]
= 
\left[
\begin{array}{c}
1\\
2\\
3\\
4\\
5
\end{array}
\right]
\]

\begin{verbatim}
A = sparse(5,5);
A(1,1) = 10;
A(1,2) = -1;
A(1,3) = -1;
A(1,4) = -1;
A(1,5) = 1;
A(2,1) = 1;
A(3,1) = 1;
A(4,1) = 1;
A(5,1) = 1;
A(2,2) = 4;
A(3,3) = 4;
A(4,4) = 4;
A(5,5) = 3;
disp(A);

size(A)

b = [1,2,3,4,5]';

disp(b);
Jacobi(A,[0,0,0,0,0]', b)
Gauss_Seidel(A,[0,0,0,0,0]', b)
\end{verbatim}

La solución para ambos métodos es:

\begin{verbatim}
ans =

    0.1520
    0.4620
    0.7120
    0.9620
    1.6160


ans =

    0.1520
    0.4620
    0.7120
    0.9620
    1.6160
 \end{verbatim}
