Nos interesa examinar uno de los algoritmos más socorridos
para resolver un sistema lineal de ecuaciones y uno de los primeros en
ser analizados en su estabilidad:
La principal característica de la Eliminación Gaussian
es transformar el sistema original
| u11 u12 u13 u14 | | x1 | | y1 | | 0 u22 u23 u24 | | x2 | | y2 | | 0 0 u33 u34 | | x3 | = | y3 | | 0 0 0 u44 | | x4 | | y4 |
su solución sería:
y4 x4 = --- u44 y3 - u34*x4 x3 = ----------- u33 y2 - u24*x4 - u23 * x3 x2 = ---------------------- u22 y1 - u14*x4 - u13 * x3 - u12 * x2 x1 = --------------------------------- u11Este proceso se le conoce como "sustitución hacia atrás".
Ahora bien, ¿cómo podemos encontrar ese "sistema equivalente
Ux=y"?
Sólo contamos con la matriz A y las operaciones elementales
entre sus filas y columnas. Tomemos A una matriz de dimensión
n=4.
| a11 a12 a13 a14 | A= | a21 a22 a23 a24 | | a31 a32 a33 a34 | | a41 a42 a43 a44 |
Las operaciones disponibles son multiplicar un vector fila por un
escalar y añadirlo a otra fila.
Paso 1: Nos interesa que las posiciones 21, 31, 41 de la matriz sean ceros, esto lo podemos lograr sustrayendo un multiplo adecuado de la pimera fila de A
fila2 <- fila2 - m21 * fila1 donde m21 = a21/a11 fila3 <- fila3 - m31 * fila1 donde m31 = a31/a11 fila4 <- fila4 - m41 * fila1 donde m41 = a41/a11y entonces nuestra matriz A la hemos transformamos en
a11 a12 a13 a14 0 u22 u23 u24 0 u32 u33 u34 0 u42 u43 u44Hemos usado la notación "uij" para enfatizarr que de la segunda y hasta la cuarta fila de A han cambiado. Note, que es preferible asignar "0" sobre la columna de a11 para evitar 'ruido' que puede producir la diferencia de cantidades iguales. Los número que producirán esos "ceros", aquellos multimos de la fila primera, se conocen como multiplicadores (mij's). La parte alterada A(2:n,2:n) se le conoce como la "parte activa" del proceso.
Paso 2: Ahora, enfoquemos nuestra atención al subproblema de una matriz de 3 por 3
a11 a12 a13 a14 ------------ 0 |u22 u23 u24 0 |u32 u33 u34 0 |u42 u43 u44En forma similar, hagamos ceros los elementos por debajo de u22, para lograrlo operemos
fila3 <- fila3 - m32 * fila2 donde m32 = u32/u22 fila4 <- fila4 - m42 * fila2 donde m42 = u42/u22En forma similar al paso 1, la parte activa del proceso será u(3:n,3:n), por lo que la idea es operar con los multimplicadores m32 y m42 que hagan cero u32 y u42.
Paso 3: Luego, nuestra matriz se ha transformado en
a11 a12 a13 a14 0 u22 u23 u24 -------------- 0 0 |v33 v34 0 0 |v43 v44Se ha empleado la notación 'v' para emfatizar que el bloque que ha cambiado luego del paso 2. Ahora, sólo nos resta operar con una fila para reducir v43 a cero.
fila4 <- fila4 - l43 * fila3 donde m43 = v43/v33Lo que nos da lugar a la matriz triangular superior
a11 a12 a13 a14 0 u22 u23 u24 0 0 v33 v34 0 0 0 w44Hemos usado distintas letras (u, v, w) para ayudarnos a indicar la parte activa, aquella submatriz de A que irá cambiando conforme vayamos operando el proceso n-1 pasos.
Cómo podemos ayudarnos obtener el sistema equivalente Ux = y a partir de lo anterior?
Introduzcamos la idea de una "matriz aumentada"
Implantación del algoritmo en Matlab
Una primera idea sería concentrarnos en la parte activa
en cada paso de la eliminación gaussiana, por ejemplo, para el paso
k tendremos
Con ello, el proceso de la Eliminación Gaussiana la escribimos como:
Para
el paso k=1 y hasta el último n-1, hacer:
1. calcular los multiplicadores necesarios para eliminar la incógnita
x(k) de las
ecuaciones de la parte activa
2. actualizar la parte activa de la matriz
finalizar
Gráficamente, la evolución del algoritmo lo podemos ver
En código Matlab el algoritmo lo podemos escribir:
%
% Script 1 para el algoritmo de la
% eliminacion gaussiana
%
for k=1:n-1,
% numero de pasos
for i=k+1:n,
% recorremos las filas
m = A(i,k)/A(k,k) % calculamos
el multiplicador
% para la fila en curso
for j=k+1:n, % recorremos
las columnas
A(i,j) = A(i,j) - m*A(k,j) % actualizamos valores
end
end
end
Ahora bien, aunque hemos programado el algoritmo en Matlab,
no hemos usado la potencialidad del ambiente de trabajo: el manejo de vectores.
En el algoritmo anterior podemos observar que el proceso de actualizar
las submatrices, la parte activa en cada paso, conlleva a operar primero
todos los elementos de la fila partir de la columna k+1 y luego
todas las columnas desde k+1. En forma vectorial, el proceso de
las filas lo podemos representar en la siguiente forma:
%
% Script 2 para el algoritmo de la
% eliminacion gaussiana
%
for k=1:n-1,
% numero de pasos
for i=k+1:n,
% recorremos las filas
m = A(i,k)/A(k,k) % calculamos
el multiplicador
% para la fila en curso
A(i,k+1:n) = A(i,k+1:n) - m*A(k,k+1:n);% actualizamos
end
% toda la fila
end
Podemos calcular primero para la k columna todos los multiplicadores que serán usados y posteriormente echar mano de ellos. Si guardamos los multiplicadores en un vector v, el script anterior lo reescribimos:
%
% Script 3 para el algoritmo de la
% eliminacion gaussiana
%
for k=1:n-1,
% numero de pasos
v(k+1:n) = A(k+1:n,k)/A(k,k);
%calculamos los multiplicadores
v = v(:);
%nos interesa el vector columna
for i=k+1:n,
% recorremos las filas
A(i,k+1:n) = A(i,k+1:n) - v(i)*A(k,k+1:n);% actualizamos
end
% toda la fila
end
Aún podemos hacer algo más, por ejemplo, el ciclo for
de adentro (el loop interior) lo podemos llevar a una representación
vectorial y entonces obtener una representación más compacta
del algoritmo.
Observemos la parte activa en el paso k y veamos cómo
se actualiza y entonces representar esa actualización
%
% Script 4 para el algoritmo de la
% eliminacion gaussiana
%
for k=1:n-1,
% numero de pasos
v(k+1:n) = A(k+1:n,k)/A(k,k);
%calculamos los multiplicadores
v = v(:);
%nos interesa el vector columna
% actualizamos la parte activa
A(k+1:n,k+1:n)
= A(k+1:n,k+1:n)- v(k+1:n)*A(k,k+1:n);
end
Hasta aquí hemos descrito el algoritmo de Elimininación Gaussiana, pero sin pivoteo. En clase, hemos visto que éste algoritmo no es estable.
Ejercicio
Imprima en pantalla la matriz triangular superior obtenida al aplicar
el algoritmo de Eliminación Gaussiana a las matrices:
A = [ 2 1 1 0; 4 3 3 1; 8 7 9 5; 6 7 9 8];
A = magic(5);
A=hilb(6);
En Matlab, magic(n)
es una función que devuelve una matriz del tipo 'cuadrado mágico'
(help magic), por otra parte, hilb(n)
es una función que devuelve a la salida la matriz de hilbert de
orden n, donde h(i,j)=1/(i+j-1);.
Descomposición LU
Ahora bien, la eliminación gaussiana, (obtención de la
matriz superior U y resolver a y) en puede verse como una
transformación de la matriz A para lograr una matriz U.
En cada paso k del algoritmo, vamos transformando la matriz para hacer
ceros debajo de la columna k, para lograrlo echamos mano de los multiplicadores
por lo que es facil ver que las matrices Mk con 1's en la diagonal y, sobre
la columna k y debajo de la diagonal los multipladores con signo contrario
para el paso k esas matrices nos transforman del paso previo hacia el paso
k.
| 1 0 0 0 | M2 = | 0 1 0 0 | | 0 -m32 1 0 | | 0 -m43 0 0 |
Así, luego de n-1 pasos, cuando hemos concluido, se ha acumado
una serie de transformaciones a la matriz A que pueden ser escritas como
Mn-1 ... M2 M1 A = U
Es fácil observar que la inversa de Mk es ella misma pero con
los signos cambiados.
| 1 0 0 0 | inv(M2) = | 0 1 0 0 | | 0 m32 1 0 | | 0 m43 0 0 |
y entonces L la triangular superior de A en la descomposicón LU viene dada por
L = inv(M1) inv(M2) ... inv(Mn-1)
Con esto en mente, podemos ya construir una función en Matlab que genere la descomposición LU de A.
Si retornamos al script 4, la función podemos escribirla en la forma:%
Así bien, contando con la descomposición LU de A podemos resolver el sistema original
A x = b
en la forma
Ly = b
para luego resolver
Ux = y
al aplicarle nuestro algoritmo obtenemos:
Contando con la descomposición LU podemos resolver Ax=b
con distintas b.
Es importante hacer notar que tanto L como U son matrices
triangulares con un particularidad, la diagonal de L son unos, pero a final
de cuentas, son triangulares. En clase discutimos la sustitución
hacia adelante y hacia atrás para resolver Ax=b, para la sustitución
hacia adelante tenemos:
Describamos otra forma de representar el proceso.
Observemos que y1 aporta tanto a y2 como a y3 hasta y6 información,
luego y2 aporta a y3, a y4 y a y6 tambien información. Pero para
generar y2 es suficiente contar con y1, y desde luego, para generar y3
es suficiente contar con y2 y desde luego haber contado con y1. Esto nos
da una pauta para recorrer el procedimiento no fila por fila, sino columna
por columna aportando el valor de la variable calculada del paso anterior.
el resultado es un algoritmo muy simpático ya que en los pasos para lograr la solución se irá obteniendo una a una las incógnitas y será la última quien aporte mayor información para calcular la siguiente logrando un procedimiento 'iterativo' hacia la solución, partimos de b para lograr y, podemos incluso verlo como un proceso recursivo en el que la dimensión;n del sistema va bajando en uno.
Para el ejemplo que hemos considerado, los pasos que nos llevaría a encontrar a y los escribimos como:
de forma similar podemos interpretar a Ux=y y entonces si imprimimos lo pasos tenemos:
En Matlab la implantación del algoritmo para resolver Ly=b sería
function [y,iband] = ltrisol(L,y)
%
% Rutina que resuelve el sistema triangular
Ly = b
% haciendo uso de la orientacion columna
de la matriz
%
% Modo de uso
%
% function
y = ltrisol(L,b)
%
% Parametros de entrada
%
% L
Matriz triangular superior, no singular.
%
n x n
%
% b
Vector columna
%
n x 1
%
% Parametros de salida
%
% y
Vector solucion
%
n x 1
%
% iband
Bandera de salida
%
= 0 No hubo problemas
%
= -1 La matriz L no es cuadrada
%
= -2 No es posible resolver Ly=b puesto
%
que b no es de dimension n.
%
= -3 La matriz es singular
% Inicio
%
iband = 0;
[m ,n ] = size(L);
nb = length(y);
%
% Verificamos que L sea cuadrada
%
if (m~=n),
iband = -1;
break; break;
end
%
% Comparamos la dimension de b
%
if (nb~=n),
iband = -2
break; break;
end
%
%
%
for j=1:n-1
if (abs(L(j,j))<eps),
% se chequea la diagonal
iband=-3;
break; break; break;
end
y(j) = y(j)/L(j,j);
y(j+1:n) = y(j+1:n) -
y(j)*L(j+1:n,j);
end
y(n)=y(n)/L(n,n);
%
% Ultima instruccion
%
y de manera muy similar para Ux=y podemos describirla
function [x,iband] = utrisol(U,x)
%
% Rutina que resuelve el sistema triangular
Ux = y
% haciendo uso de la orientacion columna
de la matriz
%
% Modo de uso
%
% function
x = utrisol(U,y)
%
% Parametros de entrada
%
% U
Matriz triangular superior.
%
n x n
%
% y
Vector columna
%
n x 1
%
% Parametros de salida
%
% x
Vector solucion
%
n x 1
%
% iband
Bandera de salida
%
= 0 No hubo problemas
%
= -1 La matriz L no es cuadrada
%
= -2 No es posible resolver Ux=y puesto
%
que y no es de dimension n.
%
= -3 La matriz es singular
% Inicio
%
iband = 0;
[m ,n ] = size(U);
nb = length(x);
%
% Verificamos que L sea cuadrada
%
if (m~=n),
iband = -1;
break; break;
end
%
% Comparamos la dimension de y
%
if (nb~=n),
iband = -2
break; break;
end
%
%
%
for j=n:-1:2
if (abs(U(j,j))<eps),
% se chequea la diagonal
iband=-3;
break; break; break; % aborta el procedimiento
end
x(j) = x(j)/U(j,j);
x(1:j-1) = x(1:j-1) -
x(j)*U(1:j-1,j);
end
x(1)=x(1)/U(1,1);
%
% Ultima instruccion
%
Resuelva los siguientes sistemas usando descomposición LU con y sin pivoteo
