3.3 Gauss – Método

Referencia: Chapra 9.2 p254 pdf p278

Para la solución con Gauss procede con la matriz aumentada para aplicar :

  • eliminación hacia adelante de incógnitas:A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{a_{i,i}}
  • sustitución  hacia atras, usando la fórmula:
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, …

Ejemplo:  Se usa el ejemplo01 para matrices por lo que los datos de trabajo son:

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

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

se ejecuta el algoritmo para mostrar la matriz aumentada

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

que luego se pivotea

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

y se ejecuta el algoritmo de eliminación hacia adelante,

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

y con ello se utiliza sustitución hacia atras,

x_3 = 40.5/5 = 8.1 3.4 x_2 +6.8 x_3 = 70.38 3.4 x_2 = 70.38 -6.8 x_3 = 70.38 -6.8 (8.1) x_2 = (70.38 -6.8 (8.1))/3.4 = 4.5

y se contiua con X1

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

para verificar que el resultado es correcto, se usa A.X = B

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

El algoritmo desarrollado por partes:

# Método de Gauss,
# Recibe las matrices A,B
# presenta solucion X que cumple: A.X=B
import numpy as np

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

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

# PROCEDIMIENTO
# Matriz Aumentada
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
print('aumentada: ')
print(AB)

# pivotea por fila AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = np.abs(AB[i:,i])
    dondemax = np.argmax(columna)
    # revisa dondemax no está en la diagonal
    if (dondemax != 0):
        # intercambia fila
        temporal = np.copy(AB[i,:])
        AB[i,:] = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal
print('pivoteada: ')
print(AB)

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

# 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]
X = np.transpose([X])

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('por sustitución hacia atras')
print('el vector solución X es:')
print(X)

print('verificar que A.X = B')
print(verifica)

Tarea: convertir a funciónes cada parte.