3.4 Método de Gauss-Jordan con Python

Referencia: Chapra 9.7 p277 pdf301, Burden 9Ed Ex6.1.12 p370, Rodriguez 4.2 p106

El método de Gauss-Jordan es semejante al método de Gauss en los procedimientos para obtener:

  • la matriz aumentada,
  • pivoteada por filas
  • eliminación hacia adelante

El cambio se presenta a partir de la matriz triangular superior, donde se aplica el procedimiento para obtener la solución:

  • eliminación hacia atrás

Método de Gauss – Jordan

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

A partir de haber terminado el procedimiento de «eliminación hacia adelante» y haber obtenido la «matriz triangular superior» aumentada, se aplica el procedimiento de.  «eliminación hacia atrás».

Se continúa con el ejercicio desde la «matriz triangular superior» aumentada:

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

Eliminación hacia atras

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

de 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 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 x3 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 x2, 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]]

iteracion 4, operacion fila 1

No hay otras filas con las que iterar, por lo que solo se obtiene el valor de x1 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 columa de constantes.

solución X es:
[[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.


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 atras para completar la solución al sistema de ecuaciones.

El algoritmo desarrollado por partes:

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

# 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]
X = np.copy(AB[:,ultcolumna])
X = np.transpose([X])


# 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('solución de X: ')
print(X)

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