3.3 Método de Gauss con Python

[ 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 ai,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 ai,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 x2

[ 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 x3 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 X1 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 ]