Referencia: Chapra 10.1 p284, pdf306. Burden 9Ed p400
Una matriz A puede separarse en dos matrices triangulares:
- L de tipo triangular inferior
- U de tipo triangular superior
que entre ellas tienen la propiedad que: A = L.U
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.
Ejercicio
Ejemplo Chapra 10.1 p285, pdf309
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 = np.array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]], dtype=float)
B = np.array([7.85,-19.3,71.4], dtype=float)
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. ]]
Algoritmo 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.
# Matrices L y U # Modificando el método de Gauss import numpy as np # INGRESO A = np.array([[ 3. , -0.1, -0.2], [ 0.1, 7. , -0.3], [ 0.3, -0.2, 10. ]], dtype=float) B = np.array([7.85,-19.3,71.4], dtype=float) # PROCEDIMIENTO B = np.transpose([B]) AB = np.concatenate((A,B),axis=1) AB = np.copy(AB) # Pivoteo parcial por filas tamano = np.shape(AB) n = tamano[0] m = tamano[1] # Para cada fila en AB for i in range(0,n-1,1): # columna desde diagonal i en adelante columna = abs(AB[i:,i]) dondemax = np.argmax(columna) # dondemax no está en diagonal if (dondemax !=0): # intercambia filas temporal = np.copy(AB[i,:]) AB[i,:] = AB[dondemax+i,:] AB[dondemax+i,:] = temporal AB1 = np.copy(AB) A1 = np.copy(AB[:,:m-1]) B1 = np.copy(AB[:,m-1]) # eliminacion hacia adelante # se inicializa L L = np.identity(n,dtype=float) 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 U = np.copy(AB[:,:m-1]) # SALIDA print('Pivoteo parcial por filas') print(AB1) print('eliminación hacia adelante') print('Matriz U: ') print(U) print('matriz L: ') print(L)
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])
Tarea
Verificar los resultados, y considerar las divisiones para cero y «casicero».
Verificar resultados con la operación A=L.U
>>> np.dot(LU[0],LU[1])
array([[ 3. , -0.1, -0.2],
[ 0.1, 7. , -0.3],
[ 0.3, -0.2, 10. ]])
Para encontrar la solución al sistema de ecuaciones, se plantea resolver:
– sustitución hacia adelante de LY=B1
– sustitución hacia atras para UX=Y
Referencia: Chapra 10.2 p287, pdf312.