[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
1. Matrices triangulares A=L.U
Referencia: Chapra 10.1 p284. Burden 6.5 p298
Una matriz A puede separarse en dos matrices triangulares que entre ellas tienen la propiedad: A = L.U:
- L de tipo triangular inferior
- U de tipo triangular superior
La matriz U se obtiene después de aplicar el proceso de eliminación hacia adelante del método de Gauss.
La matriz L contiene los factores usados en el proceso de eliminación hacia adelante del método de Gauss, escritos sobre una matriz identidad en las posiciones donde se calcularon.
[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..
2. Ejercicio
Referencia: Chapra Ejemplo 10.1 p285
Presente las matrices LU de la matriz A siguiente:
A= \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix}A = [[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]]
B = [7.85,-19.3,71.4]
El resultado del «pivoteo por filas» y «eliminación hacia adelante» se tiene:
Pivoteo parcial por filas [[ 3. -0.1 -0.2 7.85] [ 0.1 7. -0.3 -19.3 ] [ 0.3 -0.2 10. 71.4 ]]
de donde la última columna es el nuevo vector B
eliminación hacia adelante Matriz U: [[ 3. -0.1 -0.2 ] [ 0. 7.00333333 -0.29333333] [ 0. 0. 10.01204188]]
y guardando los factores del procedimiento de «eliminación hacia adelante » en una matriz L que empieza con la matriz identidad se obtiene:
matriz L: [[ 1. 0. 0. ] [ 0.03333333 1. 0. ] [ 0.1 -0.02712994 1. ]]
[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..
3. Algoritmo LU en Python
Realizado a partir del algoritmo de la sección método de Gauss y modificando las partes necesarias para el algoritmo.
Para éste algoritmo, se procede desde el bloque de pivoteo por filas«, continuando con el algoritmo de «eliminación hacia adelante» del método de Gauss. Procedimientos que dan como resultado la matriz U.
La matriz L requiere iniciar con una matriz identidad, y el procedimiento requiere que al ejecutar «eliminación hacia adelante» se almacene cada factor con el que se multiplica la fila para hacer cero. El factor se lo almacena en la matriz L, en la posición de dónde se determinó el factor.
El resultado obtenido con el algoritmo es:
Matriz aumentada [[ 3. -0.1 -0.2 7.85] [ 0.1 7. -0.3 -19.3 ] [ 0.3 -0.2 10. 71.4 ]] Pivoteo parcial: Pivoteo por filas NO requerido matriz L: [[ 1. 0. 0. ] [ 0.0333 1. 0. ] [ 0.1 -0.0271 1. ]] Matriz U: [[ 3. -0.1 -0.2 ] [ 0. 7.0033 -0.2933] [ 0. 0. 10.012 ]] >>>
Donde se puede verificar que L.U = A usando la instrucción np.dot(L.U):
>>> np.dot(L,U)
array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]])
>>>
Instrucciones en Python
# Matrices L y U # Modificando el método de Gauss import numpy as np # PROGRAMA ------------ # INGRESO A = [[ 3. , -0.1, -0.2], [ 0.1, 7. , -0.3], [ 0.3, -0.2, 10. ]] B = [7.85,-19.3,71.4] # PROCEDIMIENTO np.set_printoptions(precision=4) # 4 decimales en print casicero = 1e-15 # redondear a cero def pivoteafila(A,B,vertabla=False): ''' Pivotea parcial por filas Si hay ceros en diagonal es matriz singular, Tarea: Revisar si diagonal tiene ceros ''' A = np.array(A,dtype=float) B = np.array(B,dtype=float) # Matriz aumentada nB = len(np.shape(B)) if nB == 1: B = np.transpose([B]) AB = np.concatenate((A,B),axis=1) if vertabla==True: print('Matriz aumentada') print(AB) print('Pivoteo parcial:') # Pivoteo por filas AB tamano = np.shape(AB) n = tamano[0] m = tamano[1] # Para cada fila en AB pivoteado = 0 for i in range(0,n-1,1): # columna desde diagonal i en adelante columna = np.abs(AB[i:,i]) dondemax = np.argmax(columna) # dondemax no es en diagonal if (dondemax != 0): # intercambia filas temporal = np.copy(AB[i,:]) AB[i,:] = AB[dondemax+i,:] AB[dondemax+i,:] = temporal pivoteado = pivoteado + 1 if vertabla==True: print(' ',pivoteado, 'intercambiar filas: ',i, 'con', dondemax+i) if vertabla==True: if pivoteado==0: print(' Pivoteo por filas NO requerido') else: print(AB) return(AB) AB = pivoteafila(A,B,vertabla=True) tamano = np.shape(AB) n = tamano[0] m = tamano[1] AB0 = np.copy(AB) L = np.identity(n,dtype=float) # Inicializa L # eliminacion hacia adelante for i in range(0,n-1,1): pivote = AB[i,i] adelante = i+1 for k in range(adelante,n,1): factor = AB[k,i]/pivote AB[k,:] = AB[k,:] - AB[i,:]*factor L[k,i] = factor # llena L U = np.copy(AB[:,:n]) # SALIDA print('matriz L: ') print(L) print('Matriz U: ') print(U)
Si se requiere una respuesta unificada en una variable, se puede convertir en una arreglo de matrices [L,U].
Las matrices L y U se pueden leer como L=LU[0] y U=LU[1]
LU = np.array([L,U]) # SALIDA print('triangular inferior L') print(LU[0]) print('triangular superior U') print(LU[1])
[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
4. Algoritmo como Función de Python
El resultado a obtener es:
Matriz aumentada [[ 3. -0.1 -0.2 7.85] [ 0.1 7. -0.3 -19.3 ] [ 0.3 -0.2 10. 71.4 ]] Pivoteo parcial: Pivoteo por filas NO requerido Elimina hacia adelante: fila 0 pivote: 3.0 factor: 0.03333333333333333 para fila: 1 factor: 0.09999999999999999 para fila: 2 fila 1 pivote: 7.003333333333333 factor: -0.027129938124702525 para fila: 2 fila 2 pivote: 10.012041884816753 [[ 3. -0.1 -0.2 7.85 ] [ 0. 7.0033 -0.2933 -19.5617] [ 0. 0. 10.012 70.0843]] matriz L: [[ 1. 0. 0. ] [ 0.0333 1. 0. ] [ 0.1 -0.0271 1. ]] Matriz U: [[ 3. -0.1 -0.2 7.85 ] [ 0. 7.0033 -0.2933 -19.5617] [ 0. 0. 10.012 70.0843]] >>>
Instrucciones en Python
# Matrices L y U # Modificando el método de Gauss import numpy as np def pivoteafila(A,B,vertabla=False): ''' Pivoteo parcial por filas, entrega matriz aumentada AB Si hay ceros en diagonal es matriz singular, Tarea: Revisar si diagonal tiene ceros ''' A = np.array(A,dtype=float) B = np.array(B,dtype=float) # Matriz aumentada nB = len(np.shape(B)) if nB == 1: B = np.transpose([B]) AB = np.concatenate((A,B),axis=1) if vertabla==True: print('Matriz aumentada') print(AB) print('Pivoteo parcial:') # Pivoteo por filas AB tamano = np.shape(AB) n = tamano[0] m = tamano[1] # Para cada fila en AB pivoteado = 0 for i in range(0,n-1,1): # columna desde diagonal i en adelante columna = np.abs(AB[i:,i]) dondemax = np.argmax(columna) # dondemax no es en diagonal if (dondemax != 0): # intercambia filas temporal = np.copy(AB[i,:]) AB[i,:] = AB[dondemax+i,:] AB[dondemax+i,:] = temporal pivoteado = pivoteado + 1 if vertabla==True: print(' ',pivoteado, 'intercambiar filas: ',i, 'con', dondemax+i) if vertabla==True: if pivoteado==0: print(' Pivoteo por filas NO requerido') else: print(AB) return(AB) def gauss_eliminaAdelante(AB, vertabla=False,lu=False,casicero = 1e-15): ''' Gauss elimina hacia adelante, a partir de, matriz aumentada y pivoteada. Para respuesta en forma A=L.U usar lu=True entrega[AB,L,U] ''' tamano = np.shape(AB) n = tamano[0] m = tamano[1] L = np.identity(n,dtype=float) # Inicializa L if vertabla==True: print('Elimina hacia adelante:') for i in range(0,n,1): pivote = AB[i,i] adelante = i+1 if vertabla==True: print(' fila',i,'pivote: ', pivote) for k in range(adelante,n,1): if (np.abs(pivote)>=casicero): factor = AB[k,i]/pivote AB[k,:] = AB[k,:] - factor*AB[i,:] for j in range(0,m,1): # casicero revisa if abs(AB[k,j])<casicero: AB[k,j]=0 L[k,i] = factor # llena L if vertabla==True: print(' factor: ',factor,' para fila: ',k) else: print(' pivote:', pivote,'en fila:',i, 'genera division para cero') respuesta = AB if vertabla==True: print(AB) if lu==True: U = AB[:,:n-1] respuesta = [AB,L,U] return(respuesta) # PROGRAMA ------------------------ # INGRESO A = [[ 3. , -0.1, -0.2], [ 0.1, 7. , -0.3], [ 0.3, -0.2, 10. ]] B = [7.85,-19.3,71.4] # PROCEDIMIENTO np.set_printoptions(precision=4) # 4 decimales en print casicero = 1e-15 # redondear a cero AB = pivoteafila(A,B,vertabla=True) AB = gauss_eliminaAdelante(AB,vertabla=True, lu=True) L = AB[1] U = AB[0] # SALIDA print('matriz L: ') print(L) print('Matriz U: ') print(U)
Función en Scipy
Luego del resultado o definida la matriz a, la instrucción en la librería Scipy es:
>>> import scipy as sp
>>> [L,U] =sp.linalg.lu(A,permute_l=True)
>>> L
array([[ 1. , 0. , 0. ],
[ 0.03333333, 1. , 0. ],
[ 0.1 , -0.02712994, 1. ]])
>>> U
array([[ 3. , -0.1 , -0.2 ],
[ 0. , 7.00333333, -0.29333333],
[ 0. , 0. , 10.01204188]])
>>>
[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..
5. Solución con LY=B , UX=Y
Referencia: Chapra 10.2 p287, pdf312
Se plantea resolver el sistema de ecuaciones usando matrices triangulares
A = \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix} B = [7.85,-19.3,71.4]Para encontrar la solución al sistema de ecuaciones, se plantea resolver:
– sustitución hacia adelante de LY=B
– sustitución hacia atrás para UX=Y
Al algoritmo de la sección anterior se añade los procedimientos correspondientes con los que se obtiene la solución:
Matriz aumentada [[ 3. -0.1 -0.2 7.85] [ 0.1 7. -0.3 -19.3 ] [ 0.3 -0.2 10. 71.4 ]] Pivoteo parcial: Pivoteo por filas NO requerido matriz L: [[ 1. 0. 0. ] [ 0.03333333 1. 0. ] [ 0.1 -0.02712994 1. ]] Matriz U: [[ 3. -0.1 -0.2 ] [ 0. 7.00333333 -0.29333333] [ 0. 0. 10.01204188]] B1 : [[ 7.85] [-19.3 ] [ 71.4 ]] Y Sustitución hacia adelante [[ 7.85 ] [-19.56166667] [ 70.08429319]] X Sustitución hacia atras [ 3. -2.5 7. ] >>>
Instrucciones en Python
# Solución usando Matrices L y U # continuación de ejercicio anterior import numpy as np # PROGRAMA ------------ # INGRESO A = [[ 3. , -0.1, -0.2], [ 0.1, 7. , -0.3], [ 0.3, -0.2, 10. ]] B = [7.85,-19.3,71.4] # PROCEDIMIENTO np.set_printoptions(precision=4) # 4 decimales en print casicero = 1e-15 # redondear a cero def pivoteafila(A,B,vertabla=False): ''' Pivoteo parcial por filas, entrega matriz aumentada AB Si hay ceros en diagonal es matriz singular, Tarea: Revisar si diagonal tiene ceros ''' A = np.array(A,dtype=float) B = np.array(B,dtype=float) # Matriz aumentada nB = len(np.shape(B)) if nB == 1: B = np.transpose([B]) AB = np.concatenate((A,B),axis=1) if vertabla==True: print('Matriz aumentada') print(AB) print('Pivoteo parcial:') # Pivoteo por filas AB tamano = np.shape(AB) n = tamano[0] m = tamano[1] # Para cada fila en AB pivoteado = 0 for i in range(0,n-1,1): # columna desde diagonal i en adelante columna = np.abs(AB[i:,i]) dondemax = np.argmax(columna) # dondemax no es en diagonal if (dondemax != 0): # intercambia filas temporal = np.copy(AB[i,:]) AB[i,:] = AB[dondemax+i,:] AB[dondemax+i,:] = temporal pivoteado = pivoteado + 1 if vertabla==True: print(' ',pivoteado, 'intercambiar filas: ',i, 'con', dondemax+i) if vertabla==True: if pivoteado==0: print(' Pivoteo por filas NO requerido') else: print(AB) return(AB) AB = pivoteafila(A,B,vertabla=True) tamano = np.shape(AB) n = tamano[0] m = tamano[1] AB0 = np.copy(AB) B1 = np.copy(AB[:,m-1]) # Matrices L y U L = np.identity(n,dtype=float) # Inicializa L # eliminacion hacia adelante for i in range(0,n-1,1): pivote = AB[i,i] adelante = i+1 for k in range(adelante,n,1): factor = AB[k,i]/pivote AB[k,:] = AB[k,:] - AB[i,:]*factor L[k,i] = factor # llena L U = np.copy(AB[:,:n]) # Resolver LY = B B1 = np.transpose([B1]) AB =np.concatenate((L,B1),axis=1) # sustitución hacia adelante Y = np.zeros(n,dtype=float) Y[0] = AB[0,n] for i in range(1,n,1): suma = 0 for j in range(0,i,1): suma = suma + AB[i,j]*Y[j] b = AB[i,n] Y[i] = (b-suma)/AB[i,i] Y = np.transpose([Y]) # Resolver UX = Y AB =np.concatenate((U,Y),axis=1) # sustitución hacia atrás ultfila = n-1 ultcolumna = m-1 X = np.zeros(n,dtype=float) for i in range(ultfila,0-1,-1): suma = 0 for j in range(i+1,ultcolumna,1): suma = suma + AB[i,j]*X[j] b = AB[i,ultcolumna] X[i] = (b-suma)/AB[i,i] # SALIDA print('matriz L: ') print(L) print('Matriz U: ') print(U) print('B1 :') print(B1) print('Y Sustitución hacia adelante') print(Y) print('X Sustitución hacia atras') print(X)
[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]