Etiqueta: análisis numérico

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

1. Ejercicio

Referencia: Chapra 10.2 p292, Burden 6.3 p292, Rodríguez 4.2.5 Ejemplo 1 p118
Obtener la inversa de una matriz usando el método de Gauss-Jordan, a partir de la matriz:

\begin{pmatrix} 4 & 2 & 5 \\ 2 & 5 & 8 \\ 5 & 4 & 3 \end{pmatrix}

A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

2. Desarrollo analítico

Para el procedimiento, se crea la matriz aumentada de A con la identidad I.

AI = A|I

\begin{pmatrix} 4 & 2 & 5 & 1 & 0 & 0\\ 2 & 5 & 8 & 0 & 1 & 0 \\ 5 & 4 & 3 & 0 & 0 & 1 \end{pmatrix}

[[  4.    2.    5.   1.    0.    0. ]
 [  2.    5.    8.   0.    1.    0. ]
 [  5.    4.    3.   0.    0.    1. ]]

Con la matriz aumentada AI  se repiten los procedimientos aplicados en el método de Gauss-Jordan:

• pivoteo parcial por filas
• eliminación hacia adelante
• eliminación hacia atras

De la matriz aumentada resultante, se obtiene la inversa A^-1 en la mitad derecha de AI, lugar que originalmente correspondía a la identidad.

el resultado buscado es:

la matriz inversa es:
[[ 0.2        -0.16470588  0.10588235]
 [-0.4         0.15294118  0.25882353]
 [ 0.2         0.07058824 -0.18823529]]

\begin{pmatrix} 0.2 & -0.16470588 & 0.10588235 \\ -0.4 & 0.15294118 & 0.25882353 \\ 0.2 & 0.07058824 & -0.18823529 \end{pmatrix}

---

Verifica resultado

El resultado se verifica realizando la operación producto punto entre A y la inversa, que debe resultar la matriz identidad.

A.A^-1 = I

El resultado de la operación es una matriz identidad. Observe que los valores del orden de 10^-15 o menores se consideran como casi cero o cero.

 A.inversa = identidad
[[ 1.00000000e+00 -1.38777878e-17 -1.38777878e-16]
 [ 2.22044605e-16  1.00000000e+00 -2.22044605e-16]
 [ 5.55111512e-17 -9.71445147e-17  1.00000000e+00]]

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

3. Algoritmo en Python

El algoritmo que describe el proceso en python:

# Matriz Inversa con Gauss-Jordan
# AI es la matriz aumentada A con Identidad
# Se aplica Gauss-Jordan(AI)

import numpy as np

# INGRESO
A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]], dtype=float)

# PROCEDIMIENTO
casicero = 1e-15 # Considerar como 0

# matriz identidad
tamano = np.shape(A)
n = tamano[0]
identidad = np.identity(n)

# Matriz aumentada

AB = np.concatenate((A,identidad),axis=1)
AB0 = 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)

# 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
AB2 = np.copy(AB)

# elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    pivote = AB[i,i]
    atras = i-1 
    for k in range(atras,0-1,-1):
        factor = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor
    # diagonal a unos
    AB[i,:] = AB[i,:]/AB[i,i]

inversa = np.copy(AB[:,n:])

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminacion hacia adelante')
print(AB2)
print('eliminación hacia atrás')
print(AB)
print('Inversa de A: ')
print(inversa)

el resultado buscado es:

la matriz inversa es:
[[ 0.2        -0.16470588  0.10588235]
 [-0.4         0.15294118  0.25882353]
 [ 0.2         0.07058824 -0.18823529]]
verificando
A.inversa = identidad
[[ 1.00000000e+00 -1.38777878e-17 -1.38777878e-16]
 [ 2.22044605e-16  1.00000000e+00 -2.22044605e-16]
 [ 5.55111512e-17 -9.71445147e-17  1.00000000e+00]]

Observe que el algoritmo se pude reducir si usan los procesos de Gauss-Jordan como funciones.

Tarea: Realizar el algoritmo usando una función creada para Gauss-Jordan

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

---

4. Función en Numpy

La función en la librería Numpy es np.linalg.inv():

>>> np.linalg.inv(A)
array([[ 0.2       , -0.16470588,  0.10588235],
       [-0.4       ,  0.15294118,  0.25882353],
       [ 0.2       ,  0.07058824, -0.18823529]])
>>> 

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ]

..

---

1. Método de Gauss-Jordan

Referencia: Chapra 9.7 p277, Burden 9Ed Ex6.1.12 p370, Rodríguez 4.2 p106

El método de Gauss-Jordan presenta un procedimiento alterno al de «sustitución hacia atrás» realizado para el método de Gauss.

Se mantienen los procedimientos para:

• matriz aumentada,
• pivoteo por filas
• eliminación hacia adelante

al haber obtenido la «matriz triangular superior» , se aplica el procedimiento:

• eliminación hacia atrás

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ]
..

---

2. Ejercicio

Referencia: Rodríguez 4.0 p105,  1Eva_IT2010_T3_MN Precio artículos

\begin{cases} 4x_1+2x_2+5x_3 = 60.70 \\ 2x_1+5x_2+8x_3 = 92.90 \\ 5x_1+4x_2+3x_3 = 56.30 \end{cases}

Se continúa con el ejercicio desarrollado para el método de Gauss desde la «matriz triangular superior»y aumentada, luego de eliminación hacia adelante:

[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ]
..

---

3. Eliminación hacia atrás

El procedimiento es semejante al realizado para «eliminación hacia adelante», con la diferencia que se inicia en la última fila hacia la primera.

Las operaciones se realizan de abajo hacia arriba desde la última fila.
Para el ejercicio presentado se tiene que:

utlfila = n-1 = 3-1 = 2
ultcolumna = m-1 = 4-1 = 3

iteración 1, operación fila 3 y 2

Se aplica desde la última fila i, para las otras filas k que se encuentran hacia atrás.

i = 2
pivote = AB[2,2] = 5
k = i-1 # hacia atrás

se realizan las operaciones entre las filas i y la fila k

         [0.  3.4  6.8    70.38]
-(6.8/5)*[0.  0.   5.     40.5 ]
_______________________________
       = [0.0 3.4  8.8e-16 15.3]

para reemplazar los valores de la segunda fila en la matriz aumentada

[[5.0    4.0    3.0      56.3]
 [0.0    3.4    8.8e-16  15.3]
 [0.0    0.0    5.0      40.5]]

Observe que hay un valor muy pequeño del orden de 8.8×10^-16, que para las otras magnitudes se puede considerar como casi cero.

iteración 2, operación fila 3 y 1

Se calculan los nuevos valores de indice k

k = k-1 = 2-1 = 1 # hacia atrás

se realizan las operaciones entre las filas i y la fila k

       [5.0    4.0    3.0    56.3]
-(3/5)*[0.0    0.0    5.0    40.5]
__________________________________
     = [5.0    4.0    0.0    32.0]

que al actualizar la matriz aumentada se tiene:

[[5.0    4.0    0.0      32.0]
 [0.0    3.4    8.8e-16  15.3]
 [0.0    0.0    5.0      40.5]]}

Al haber terminado las filas hacia arriba, se puede así determinar el valor de x₃ al dividir la fila 3 para el pivote

[[5.0    4.0    0.0      32.0]
 [0.0    3.4    8.8e-16  15.3]
 [0.0    0.0    1.0       8.1]]}

iteración 3, operación fila 2 y 1

se actualizan los valores de los índices:

i = i-1 = 2-1 = 1
k = i-1 = 1-1 = 0

se pueden realizar operaciones en una sola fila hacia atrás, por lo que el resultado hasta ahora es:

[[ 5.0    0.0   -1.04e-15  14.0]
 [ 0.0    3.4    8.8e-16   15.3]
 [ 0.0    0.0    1.0        8.1]]

Se obtiene el valor de x₂, dividiendo para el valor del pivote,

[[ 5.0    0.0   -1.04e-15  14.0]
 [ 0.0    1.0    2.6e-16    4.5]
 [ 0.0    0.0    1.0        8.1]]

iteración 4, operación fila 1

No hay otras filas con las que iterar, por lo que solo se obtiene el valor de x₁ al dividir para el pivote.

[[ 1.0    0.0   -2.08e-15  2.8]
 [ 0.0    1.0    2.6e-16   4.5]
 [ 0.0    0.0    1.0       8.1]]

La solución del sistema de ecuaciones se presenta como una matriz identidad concatenada a un vector columna de constantes.

solución X: 
[2.8 4.5 8.1]

X= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}

Observación: en la matriz hay unos valores del orden de 10^-16, que corresponden a errores de operaciones en computadora (truncamiento y redondeo) que pueden ser descartados por ser casi cero. Hay que establecer entonces un parámetro para controlar los casos en que la diferencia entre los ordenes de magnitud son por ejemplo menores a 15 ordenes de magnitud 10^-15. e implementarlo en los algoritmos.

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ]
..

---

4. Algoritmo en Python

Esta sección reutiliza el algoritmo desarrollado para el Método de Gauss, por lo que los bloques de procedimiento son semejantes hasta eliminación hacia adelante. Se añade el procedimiento de eliminación hacia atrás para completar la solución al sistema de ecuaciones.

 

Se obtiene el siguiente resultado con el algoritmo:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Elimina hacia adelante:
 fila 0 pivote:  5.0
   factor:  0.4  para fila:  1
   factor:  0.8  para fila:  2
 fila 1 pivote:  3.4
   factor:  -0.3529411764705883  para fila:  2
 fila 2 pivote:  5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Elimina hacia Atras:
 fila 2 pivote:  5.0
   factor:  1.3599999999999999  para fila:  1
   factor:  0.6  para fila:  0
 fila 1 pivote:  3.4
   factor:  1.1764705882352942  para fila:  0
 fila 0 pivote:  5.0
[[ 1.00000000e+00  0.00000000e+00 -2.08983158e-16  2.80000000e+00]
 [ 0.00000000e+00  1.00000000e+00  2.61228947e-16  4.50000000e+00]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00  8.10000000e+00]]
solución X: 
[2.8 4.5 8.1]
>>> 

Instrucciones en Python

# Método de Gauss-Jordan
# 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)

def gauss_eliminaAtras(AB, vertabla=False, precision=5, casicero = 1e-15):
    ''' Gauss-Jordan elimina hacia atras
    Requiere la matriz triangular inferior
    Tarea: Verificar que sea triangular inferior
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    ultfila = n-1
    ultcolumna = m-1
    if vertabla==True:
        print('Elimina hacia Atras:')
        
    for i in range(ultfila,0-1,-1):
        pivote = AB[i,i]
        atras = i-1  # arriba de la fila i
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
            
        for k in range(atras,0-1,-1):
            if (np.abs(AB[k,i])>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
 
        AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
    X = np.copy(AB[:,ultcolumna])
    
    if vertabla==True:
        print(AB)
    return(X)

# PROGRAMA ------------------------
# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70,92.90,56.30]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_eliminaAtras(AB,vertabla=True)

# SALIDA
print('solución X: ')
print(X)

Tarea: implementar caso cuando aparecen ceros en la diagonal para dar respuesta, convertir a funciones cada parte

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ]

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 atras para UX=Y

Algoritmo en Python

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

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)
B1 = np.copy(AB[:,m-1])

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

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

---

1. Determinante de matriz

Referencia: Rodríguez 4.3.9 p129, Burden 6.4 p296, Chapra 9.1.2 p250.

El determinante de una matriz cuadrada triangular superior también puede calcularse como el producto de los coeficientes de la diagonal principal, considerando el número de cambios de fila del pivoteo k.

det(A) = (-1)^k \prod_{i=1}^n a_{i,i}

Si observamos que en las secciones anteriores se tiene desarrollado los algoritmos  para obtener la matriz triangular superior en el método de Gauss, se usan como punto de partida para obtener los resultados del cálculo del determinante.

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

---

2. Ejercicio

Referencia: Chapra Ejemplo 9.12 p277

Calcular el determinante de la matriz A.

A= \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix}

El  algoritmo para ejercicio se convierte en una extensión de los algoritmos anteriores.

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

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

---

3. Algoritmo en Python

El algoritmo parte de lo realizado en método de Gauss, indicando que la matriz a procesar es solamente A. Se mantienen los procedimientos de «pivoteo parcial por filas» y » eliminación hacia adelante»

Para contar el número de cambios de filas, se añade un contador de  pivoteado dentro del condicional para intercambio de filas.

Para el resultado del operador multiplicación, se usan todas las casillas de la diagonal al acumular las multiplicaciones.

# Determinante de una matriz A
# convirtiendo a diagonal superior 

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)

# PROCEDIMIENTO

# Matriz aumentada
AB = np.copy(A) # Para usar algoritmo previo

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
pivoteado = 0 # contador para cambio fila

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

        pivoteado = pivoteado + 1 # cuenta cambio fila
        
AB1 = np.copy(AB)

# eliminación 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

# calcula determinante
multiplica = 1
for i in range(0,n,1):
    multiplica = multiplica*AB[i,i]
determinante = ((-1)**pivoteado)*multiplica

# SALIDA
print('Pivoteo parcial por filas')
print(AB1)
print('Cambios de fila, pivoteado: ',pivoteado)
print('eliminación hacia adelante')
print(AB)
print('determinante: ')
print(determinante)

Se aplica la operación de la fórmula planteada para el método, y se presenta el resultado:

Pivoteo parcial por filas
[[ 3.  -0.1 -0.2]
 [ 0.1  7.  -0.3]
 [ 0.3 -0.2 10. ]]
Cambios de fila, pivoteado:  0
eliminación hacia adelante
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
determinante: 
210.35299999999995
>>> 

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

---

4. Algoritmo como función en Numpy

El resultado se puede verificar usando la función de Numpy:

>>> np.linalg.det(A)
210.3529999999999

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

..

---

1. Método de Gauss

Referencia: Chapra 9.2 p254, Burden 6.1 p273, Rodríguez 4.3 p119

El método de Gauss opera sobre la matriz aumentada y pivoteada por filas, añadiendo los procesos:.

• eliminación hacia adelante:

A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{a_{i,i}}

• sustitución  hacia atrás:

x_i = \frac{b_i^{(i-1)}-\sum_{j=i+1}^{n}a_{ij}^{(i-1)} x_{j}}{a_{ii}^{(i-1)}}

para i = n-1, n-2, …

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

---

2. Ejercicio

Referencia: Rodríguez 4.0 p105,  1Eva_IT2010_T3_MN Precio artículos

Se continúa a partir del resultado del tema de pivoteo parcial por filas para matrices:

\begin{cases} 4x_1+2x_2+5x_3 = 60.70 \\ 2x_1+5x_2+8x_3 = 92.90 \\ 5x_1+4x_2+3x_3 = 56.30 \end{cases}

Matriz aumentada y pivoteada por filas:

[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

---

3. Eliminación hacia adelante o eliminación Gaussiana

Consiste en simplificar la matriz a una triangular superior, con ceros debajo de la diagonal, usando operaciones entre filas, para obtener:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

Los índices de fila y columna en la matriz A[i,j] se usan de forma semejante a la nomenclatura de los textos de Álgebra Lineal. Progresivamente para cada fila, se toma como referencia o pivote el elemento de la diagonal (i=j). Luego, se realizan operaciones con las filas inferiores para convertir los elementos por debajo de la diagonal en cero. Las operaciones incluyen el vector B debido a que se trabaja sobre la matriz aumentada AB.

AB Matriz aumentada y pivoteada por filas:
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]

iteración fila 1, operación fila 1 y 2

Para la fila 1, con posición i=0, se usa el elemento a_i,i como pivote.

pivote = AB[i,i] = AB[0,0] = 5

Para las filas de que están después de la diagonal se referencian como k.Se obtiene el factor escalar de la operación entre filas de la formula

k = i+1 = 0+1 = 1

A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{pivote}

factor = AB[1,0]/pivote = 2/5

y se realiza la operación entre fila k y la fila i para actualizar la fila k,

       [ 2. 5.  8.  92.9]
-(2/5)*[ 5. 4.  3.  56.3]
__________________________
     = [ 0. 3.4 6.8 70.38]

con lo que la matriz aumentada AB se actualiza a:

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 4.    2.    5.   60.7 ]]

iteración fila 1, operación fila 1 y 3

se continúa con la siguiente fila, quedando la matriz aumentada con la columna debajo de la primera diagonal en cero:

k = i+1 = 2
factor = 4/5

        [ 4.  2.  5.   60.7] 
- (4/5)*[ 5.  4.  3.   56.3]
_____________________________
      = [ 0. -1.2 2.6  15.66]

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]

Como ya se terminaron las operaciones con la primera posición de la diagonal, el siguiente paso es usar la segunda posición, i =2.

iteración fila 2

Para la fila 2, con posición i=1, se toma el elemento de la diagonal a_i,i como pivote, la variable adelante indica la referencia de la tercera fila

pivote = A[i,i] = AB[1,1] = 3.4

Para las filas ubicadas adelante de la diagonal se referencian como k

adelante = k = i+1 = 1+1 = 2

Para aplicar la fórmula por filas, se obtiene el factor .

factor = AB[2,1]/pivote  = -1.2/3.4 = - 0,3529

            [ 0. -1.2 2.6 15.66]
-(-1.2/3.4)*[ 0.  3.4 6.8 70.38]
________________________________
         =  [ 0.  0.  5.  40.5 ]

AB =
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

Con lo que se completa el objetivo de tener ceros debajo de la diagonal.
Observe que no es necesario realizar operaciones para la última fila, por lo que k debe llegar solamente hasta la fila penúltima.

El resultado de la eliminación hacia adelante a ser usado en el próximo paso es:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

---

4. Sustitución hacia atrás

La fórmula se interpreta para facilitar el algoritmo

x_i = \frac{b_i^{(i-1)}-\sum_{j=i+1}^{n}a_{ij}^{(i-1)} x_{j}}{a_{ii}^{(i-1)}}

Para una fila i, el vector b[i] representa el valor de la constante en la fila i de la matriz aumentada, a[i] se refiere los valores de los coeficientes de la ecuación, de los que se usan los que se encuentran a la derecha de la diagonal.

Las operaciones se realizan de abajo hacia arriba desde la última fila. Para el ejercicio presentado se tiene que:

ultfila = n-1 = 3-1 = 2
ultcolumna = m-1 = 4-1 = 3

la matriz a procesar es:

Elimina hacia adelante
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]

iteración 1, fila 3, i=2

Empieza desde la última fila de la matriz,

[ 0. 0. 5. 40.5 ]

0 x_1 + 0 x_2 + 5 x_3 = 40.5

El valor de la constante es b = 40.5 y no existen elementos hacia la derecha de la diagonal. No se usa la ultima columna que es de las constantes:

5 x_3 = 40.5 x_3 = 40.5/5 = 8.1

la respuesta se interpreta en el vector X como:

X= \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 8.1 \end{pmatrix}

iteración 2, fila 2,  i = 1

De la penúltima fila se obtiene la ecuación para encontrar x₂

[ 0. 3.4 6.8 70.38]

0x_1 + 3.4 x_2 +6.8 x_3 = 70.38

se observa que b = 70.38 y  a la derecha de a diagonal se tiene un solo valor de [6.8].

3.4 x_2 = 70.38 -6.8 x_3

usa el valor de x₃ encontrado en la iteración anterior

3.4 x_2 = 70.38 -6.8 (8.1)

Muestra la ecuación para la segunda fila.

x_2 = (70.38 -6.8 (8.1))/3.4 = 4.5

la respuesta se interpreta en el vector X como:

X= \begin{pmatrix} 0 \\ 4.5 \\ 8.1 \end{pmatrix}

iteración 3 fila 1, i=0

se sigue el mismo proceso para la siguiente incógnita X₁ que se interpreta como

[ 5. 4. 3. 56.3 ]

5x_1 + 4 x_2 + 3x_3 = 56.3 5x_1 = 56.3 - ( 4 x_2 + 3x_3) x_1 = \frac{56.3 - ( 4 x_2 + 3x_3)}{5}

Se encuentra que la solución al sistema de ecuaciones es:

X= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}

por sustitución hacia atrás
el vector solución X es:
[[2.8]
 [4.5]
 [8.1]]

---

Verificar respuesta

Para verificar que el resultado es correcto, se usa el producto punto entre la matriz a y el vector resultado X. La operación A.X = B debe dar el vector B.

verificar que A.X = B
[[60.7]
 [92.9]
 [56.3]]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]
..

---

5. Algoritmo con Python

El algoritmo para el Método de Gauss, reutiliza las instrucciones para matriz aumentada y pivoteo parcial por filas.

Recordar: Asegurar que los arreglos sean de tipo Real (float), para que no se considere el vector como entero y realice operaciones entre enteros, generando errores por truncamiento.

La parte nueva a desarrollar corresponde al procedimiento de «eliminación hacia adelante» y el procedimiento de «sustitución hacia atrás».

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

import numpy as np

# INGRESO
A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])
B = np.array([[60.70],
              [92.90],
              [56.30]])
# PROCEDIMIENTO
casicero = 1e-15 # Considerar como 0

# Evitar truncamiento en operaciones
A = np.array(A,dtype=float) 

# Matriz aumentada
AB  = np.concatenate((A,B),axis=1)
AB0 = 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)

# eliminación 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

# 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]

X = np.transpose([X])

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminación hacia adelante')
print(AB)
print('solución: ')
print(X)

---

Tarea

Revisar cuando la matriz pivoteada por filas tienen un elemento cero o muy cercano a cero pues la matriz sería singular. El valor considerado como casi cero podría ser 1×10^-15

A estas alturas, por la cantidad de líneas de instrucción es recomendable reutilizar bloques de algoritmos usando funciones def-return. Por ejemplo: pivoteo por filas, eliminación hacia adelante, sustitución hacia atrás.

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

..

---

6. Algoritmo como función de Python

El resultado par el ejercicio anterior es:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Elimina hacia adelante:
 fila 0 pivote:  5.0
   factor:  0.4  para fila:  1
   factor:  0.8  para fila:  2
 fila 1 pivote:  3.4
   factor:  -0.3529411764705883  para fila:  2
 fila 2 pivote:  5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
solución: 
[2.8 4.5 8.1]
>>>  

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, casicero = 1e-15):
    ''' Gauss elimina hacia adelante
    tarea: verificar términos cero
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    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,:]
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    if vertabla==True:
        print(AB)
    return(AB)

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    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]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70,92.90,56.30]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

# SALIDA
print('solución: ')
print(X)

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ] [ Algoritmo ] [ función ]

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]

Referencia: Chapra 9.4.2 p268. Burden 6.2 p279. Rodríguez 4.0 p105

Los métodos de solución a sistemas de ecuaciones en los primeros pasos usan la matriz aumentada y pivoteada por filas. Como es un procedimiento usado en todos los métodos de la unidad, para simplificar se presenta como uno de los primeros algoritmos.

Para mostrar el desarrollo del proceso se usa como referencia un ejercicio.

..

---

1. Ejercicio

Referencia: Rodríguez 4.0 p105,  1Eva_IT2010_T3_MN Precio artículos

Ejemplo 1. Un comerciante compra tres productos A, B, C, pero en las facturas únicamente consta la cantidad comprada en Kg y el valor total de la compra.

Se necesita determinar el precio unitario de cada producto.  Dispone de solo tres facturas con los siguientes datos:

Ejemplo:

	Cantidad			Valor ($)
Factura	X1	X2	X3	Pagado
1	4	2	5	60.70
2	2	5	8	92.90
3	5	4	3	56.30

Los precios unitarios se pueden representar por las variables x₁, x₂, x₃ para escribir el sistema de ecuaciones que muestran las relaciones de cantidad, precio y valor pagado:

\begin{cases} 4x_1+2x_2+5x_3 = 60.70 \\ 2x_1+5x_2+8x_3 = 92.90 \\ 5x_1+4x_2+3x_3 = 56.30 \end{cases}

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

---

2. Matriz, Vector y la Matriz Aumentada

El sistema de ecuaciones se escribe en la forma algebraica como matrices y vectores de la forma Ax=B

\begin{pmatrix} 4 & 2 & 5 \\ 2 & 5 & 8 \\5 & 4 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 60.70 \\ 92.90 \\ 56.30 \end{pmatrix}

Para el algoritmo la matriz A y el vector B se escriben como arreglos.

A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

Observe que:

• Las matrices y vectores se ingresan como arreglos de la librería Numpy
• el vector B se escribe en forma de columna
• No se usan listas, de ser el caso se convierten hacia arreglos con np.array()

Si el vector B está como fila, se aumenta una dimensión [B] y se aplica la transpuesta

Bfila = np.array([60.70,92.90,56.30])
Bcolumna = np.transpose([Bfila])
print(Bcolumna)

>>> Bcolumna
array([[60.7],
       [92.9],
       [56.3]])

En el desarrollo de la solución, se usa la matriz aumentada, para mantener sincronía entre las operaciones entre filas de la matriz A y el vector B.

La matriz aumentada AB se forma al concatenar por columnas la matriz A con el vector B,  usando el parámetro axis=1.

AB = np.concatenate((A,B),axis=1)

el resultado AB se muestra como:

>>> AB
array([[ 4. ,  2. ,  5. , 60.7],
       [ 2. ,  5. ,  8. , 92.9],
       [ 5. ,  4. ,  3. , 56.3]])

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

---

3. Pivoteo parcial por filas

Para el pivoteo por fila de la matriz aumentada AB, tiene como primer paso revisar la primera columna desde la diagonal en adelante.

columna = [|4|,
           |2|,
           |5|]
dondemax = 2

El procedimiento de «pivoteo» se realiza si la posición dónde se encuentra el valor de  mayor magnitud no corresponde a la diagonal de la matriz (posición 0 de la columna).

En el ejercicio se encuentra que la magnitud de mayor valor está en la última fila, por lo que en AB se realiza el intercambio entre la fila 3 y la fila 1

AB = [[ 5. , 4. , 3. , 56.3],
      [ 2. , 5. , 8. , 92.9],
      [ 4. , 2. , 5. , 60.7]]

Se repite al paso anterior, pero para la segunda columna formada desde la diagonal.

columna = [|5|,
           |2|]
dondemax = 0

como la posición dondemax es la primera, índice 0, se determina que ya está en la diagonal de AB y no es necesario realizar el intercambio de filas.

Se repite el proceso para la tercera columna desde la diagonal, que resulta tener solo una casilla (columna =[5]) y no ser requiere continuar.

El resultado del pivoteo por fila se muestra como:

matriz pivoteada por fila:
AB = [[ 5. ,  4. ,  3. , 56.3],
      [ 2. ,  5. ,  8. , 92.9],
      [ 4. ,  2. ,  5. , 60.7]]

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

---

4. Algoritmo en Python

Para realizar el algoritmo, es de recordar que para realizar operaciones en una matriz sin alterar la original, se usa una copia de la matriz (np.copy). Se puede comparar y observar los cambios entre la matriz original y la copia a la que se aplicaron cambios

Si no es necesaria la comparación entre el antes y después, no se realiza la copia y se ahorra el espacio de memoria, detalle importante para matrices de «gran tamaño» y una computadora con «limitada» memoria.

# Pivoteo parcial por filas
# Solución a Sistemas de Ecuaciones

import numpy as np

# INGRESO
A = np.array([[4,2,5],
              [2,5,8],
              [5,4,3]])

B = np.array([[60.70],
              [92.90],
              [56.30]])

# PROCEDIMIENTO
# Matriz aumentada
AB  = np.concatenate((A,B),axis=1)
AB0 = 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

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB)

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]
..

---

5. Función pivoteafila(M)

Los bloques de cada procedimiento que se repiten en otros métodos se convierten a funciones def-return, empaquetando las soluciones algorítmicas a problemas resueltos.

Se usa la matriz M para generalizar y diferenciar de ‘A’ que es usada en los ejercicios en realizados en adelante.

# Pivoteo parcial por filas
# Solución a Sistemas de Ecuaciones

import numpy as np

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

# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70, 92.90, 56.30]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

# SALIDA
print('Resultado de Pivoteo parcial por filas')
print(AB)

El resultado del ejercicio es:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
Pivoteo parcial por filas
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
>>>

[ Ejercicio ] [ Matriz Aumentada ] [ Pivotea filas ] [ Algoritmo ] [ función ]