3.3.2 Matrices triangulares A=L.U con Python

[ Matriz triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]

..


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 triangular ][ Ejercicio ] [ Algoritmo ] [ Función ]
..


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 = 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.        ]]

[ Matriz triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..


3. 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.

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.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
Matriz U: 
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
>>> 

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

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,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

# 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
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 triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]

..


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.00333333  -0.29333333 -19.56166667]
 [  0.           0.          10.01204188  70.08429319]]
matriz L: 
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
Matriz U: 
[[  3.          -0.1         -0.2          7.85      ]
 [  0.           7.00333333  -0.29333333 -19.56166667]
 [  0.           0.          10.01204188  70.08429319]]
>>> 

Instrucciones en Python

# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

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,'y', 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,:]

                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
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 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 triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]