3.5 Matrices triangulares A=L.U

Referencia: Chapra 10.1 p282, pdf306

Una matriz A puede separarse en dos matrices: L y U que son de tipo triangular inferior L y triangular superior U.

Tienen la propiedad que:  A = L.U

A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

se convierte en:

triangular inferior L
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
triangular superior U
[[  3.          -0.1         -0.2       ]
 [  0.           7.00333333  -0.29333333]
 [  0.           0.          10.01204188]]

Para mantener la matriz original, se obtiene una copia de la matriz A denominada U.

Se aplica el algoritmo de eliminación hacia adelante de Gauss, determinando  el factor con el que se multiplica la fila para hacer cero. El factor se lo almacena en la matriz U, en la posición de dónde se determinó el factor.

Las matrices resultantes son L y U, que se presentan en una matriz de dos elementos LU[0] y LU[1] para entregar la respuesta unificada.

 Triangulares de una matriz LU
# Triangular inferior L = LU[0]
# Triangular superior U = LU[1]
# Referencia: Chapra 10.1.2. p.284, pdf.308
# Tarea: verificar si el pivote es cero
#        verificar diagonal dominante
import numpy as np

# Programa de prueba
# INGRESO
A = np.array([[3  , -0.1, -0.2],
              [0.1,  7  , -0.3],
              [0.3, -0.2, 10  ]])

# PROCEDIMIENTO
casicero = 1e-15 # 0
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]

U = np.copy(A)
L = np.identity(n,dtype=float)

# Gauss elimina hacia adelante
# tarea: verificar términos cero
for i in range(0,n,1):
    pivote = U[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = U[k,i]/pivote
            U[k,:] = U[k,:] - factor*U[i,:]
            L[k,i] = factor

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.

Verificando resultados: 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. ]])