Práctica 6. Sistemas de ecuaciones lineales II (Descomposición LU)
En esta práctica estudiaremos la implementación de algunos métodos de factorización LU para resolver sistemas Ax=b con A regular y factorizable en la forma LU. Concretamente se tratarán los métodos de Doolittle y de Crout.
Para ello, estudiaremos en primer lugar la implementación con MatLab de las factorizaciones de Doolittle y de Crout de una matriz A, que consisten en encontrar las matrices L y U tales que L*U = A, siendo L una matriz triangular inferior y U una matriz triangular superior, con unos en la diagonal de L en el caso de Doolittle, y con unos en la diagonal de U en el caso de Crout.
Antes de esto resultará útil saber cómo hallar con MatLab los sumatorios usando vectores.
Contents
Vectorización del sumatorio
En primer lugar, veamos como implementar los sumatorios que aparecen en la factorización de Doolittle de una matriz mediante un producto escalar de dos vectores. Los sumatorios que aparecen en la factorización de Crout son análogos.
En la factorización de Doolittle, al hallar U(k,j) debemos implementar el sumatorio de L(k,r)*U(r,j) donde r=1, 2, ..., k-1.
Es decir, debemos hallar L(k,1)*U(1,j) + L(k,2)*U(2,j) + ... + L(k, k-1)*U(k-1 ,j)
Este sumatorio puede vectorizarse como un producto escalar del vector fila L(k, 1:k-1) por el vector columna U(1:k-1 ,j).
De forma similar, al hallar L(i,k) debemos implementar el sumatorio de L(i,r)*U(r,k) donde r=1, 2, ..., k-1.
Este sumatorio puede vectorizarse como el producto escalar L(i , 1:k-1)*U(1:k-1 , k)
Factorización de Doolittle de una matriz
En la factorización de Doolitte la matriz L contiene unos en la diagonal. Comenzamos definiendo la matriz A que vamos a factorizar e inicializando variables.
A=[[6,16,21,33] ; [2,20,18,16] ; [3,18,19,21] ; [1,2,3,4]]; n=length(A); % Inicializamos las matrices U y L como matrices cuadradas de orden n % formadas por ceros. U=zeros(n); L=zeros(n);
Ahora implementamos el caso k=1, ya que en el algoritmo correspondiente aparecen sumatorios en los que el índice r recorre desde 1 hasta k-1 cuando k=2,...n. Cuando k=1, se entiende que este sumatorio no debe realizarse, y por esto en Matlab debemos hacer aparte el caso k=1. Observe cómo se asigna la primera fila de U y la primera columna de L mediante vectores.
U(1, 1:n)=A(1, 1:n); L(1,1)=1; L(2:n , 1)=( A(2:n , 1) )/U(1,1);
Hallaremos el resto de los elementos de L y de U siguiendo el siguiente orden :
Para 2<=k<=n
Para k<=j<=n U(k,j) Siguiente j
L(k,k)=1
Para k+1<=i<=n L(i,k) Siguiente i
Siguiente k
En los elementos que requieren calcular un sumatorio, estos se han calculado antes por claridad, aunque podría haberse incluido el sumatorio en la misma linea.
for k=2:n for j=k:n suma=L(k , 1:k-1)*U(1:k-1 , j); U(k,j)=A(k,j)-suma; end L(k,k)=1; for i=k+1:n suma=L(i , 1:k-1)*U(1:k-1 , k); L(i,k)=(A(i,k)-suma) / U(k,k); end end
Mostramos el resultado y comprobamos la factorización A=L*U.
disp('Resultado :') disp('L ='); disp(L); disp('U ='); disp(U); disp('Comprobación :') disp('A-L*U ='); disp(A-L*U);
Resultado :
L =
1.0000 0 0 0
0.3333 1.0000 0 0
0.5000 0.6818 1.0000 0
0.1667 -0.0455 0 1.0000
U =
6.0000 16.0000 21.0000 33.0000
0 14.6667 11.0000 5.0000
0 0 1.0000 1.0909
0 0 0 -1.2727
Comprobación :
A-L*U =
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
Ej 1. Emplee el programa proporcionado para hallar la factorización de Doolittle de la matriz A=[ [1,2,-6,2] ; [2,-1,-5,-3] ; [3,-54,77,-69] ; [4,3,-138,-104] ];
(Sol. L=[ [1,0,0,0] ; [2,1,0,0] ; [3,12,1,0] ; [4,1,-11,1] ];
U=[ [1,2,-6,2] ; [0,-5,7,-7] ; [0,0,11,9] ; [0,0,0,-6] ]; )
Factorización de Crout de una matriz
De forma similar se puede implementar la factorización de Crout, en la que aparecen unos en la diagonal de U. Después de considerar el caso k=1, para hallar la primera columna de L y la primera fila de U, podemos calcular el resto de elementos de L y de U en el siguiente orden.
Para 2<=k<=n
Para k<=i<=n L(i,k) Siguiente i
U(k,k)=1
Para k+1<=j<=n U(k,j) Siguiente j
Siguiente k
Ej 2. Implemente la factorización de Crout y compruébela factorizando la matriz A=[[2,0,-6,-4] ; [11,1,-22,-17]; [3,4,38,17]; [-1,-2,-19,-9]];
(Sol. L=[ [2,0,0,0] ; [11,1,0,0] ; [3,4,3,0] ; [-1,-2,0,-1] ]
U=[ [1,0,-3,-2] ; [0,1,11,5] ; [0,0,1,1] ; [0,0,0,1] ] )
Resolución de sistemas usando una descomposición L*U.
Si suponemos que A = L*U, el sistema A*x = b puede escribirse de la forma L*U*x = b. Para hallar x podemos introducir un vector auxiliar y = U*x, de forma que nuestro sistema se escribe como L*y = b, donde despejar y se reduce a resolver un sistema triangular. Una vez conocido el vector y, para hallar el vector x basta con resolver el sistema U*x = y, que es de nuevo un sistema triangular.
Así pues, resolver un sistema Ax=b usando la descomposición LU se reduce a hallar una factorización LU de la matriz A , y a resolver a continuación dos sistemas triangulares. Tanto la factorización LU de A como la resolución de los sistemas triangulares se realiza con pequeños programas como los que hemos descrito en esta práctica 6 y en la práctica 5.
De forma concisa repetimos los pasos que deden darse para resolver un sistema usando una descomposición L*U :
1- Descomponer A = L*U. 2- Resolver L y = b con lo que obtenemos y. 3- Resolver U x = y con lo que obtenemos x.
Ej 3. Sabiendo que A = L*U con L y U las matrices dadas a continuación y que b es el vector proporcionado, resuelva el sistema A*x = b siguiendo los pasos descritos en esta sección (Emplee, por comodidad, la división izquierda para resolver los dos sistemas triangulares). Halle A*x-b para comprobar que el vector x hallado es solución del sistema. (Sol. x=[1; 0; -2; 3]. )
L=[[1, 0, 0, 0]; [2, 1, 0, 0]; [-2, 3, 1, 0]; [-1, -1, -2, 1]] ;
U=[[-1, 4, 1, 2]; [0, 2, -2, 0]; [0, 0, -2, -1]; [0, 0, 0, 3]] ;
b=[3; 10; 7; 0] ;
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
17-Mar-2017