3.1 Solución de sistema de 3×3

Solución de un sistema de 3×3 como un punto de intersección de planos

Referencia del ejercicio: http://blog.espol.edu.ec/matg1013/1eva_iit2011_t2-sistema-de-ecuaciones/

Considere el siguiente sistema de ecuaciones:

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}

Se pueden representar como planos en el espacio despejando la variable z para cada ecuación, de la forma:

\begin{cases} z=(1+ 2x - 5y)/(9)\\z=(6 -7x-y)/(1)\\z=(-26+3x-7y)/(-1)\end{cases}

Se visualizan los resultados de cada ecuación como los planos mostrados en la gráfica:

La intersección de los planos genera un punto cuyas coordenadas corresponden a la solución del sistema.

Para observar mejor del resultado, ejecute las intrucciones en python propuestas y rote el gráfico resultante con el cursor.


Intrucciones en python

Para generalizar el ingreso de las ecuaciones se usa la forma matricial Ax=B, usando solo los coeficientes.

# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

Para visualizar los planos se requiere cada ecuación como funciones(x,y), por simplicidad se usan formulas en formato lambda. Con la forma matricial se generan las ecuaciones de cada plano usando los coeficientes al despejar la variable z.

# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

Para evaluar las funciones se requieren muestras en cada punto del plano X,Y, semejante a las gráficas en 2D que solo requerian muestra en el eje X.

– Primero se generan las muestras para cada eje en los vectores xi, yi.

a = -5
b = 5
muestras = 21

xi = np.linspace(a,b, muestras)
yi = np.linspace(a,b, muestras)

= Las combinaciones entre las muestras de cada eje se obtienen generando las matrices Xi, Yi que representan la malla de muestras en el plano para evaluar cada punto.

Xi, Yi = np.meshgrid(xi,yi)

– Se evaluan los puntos Xi,Yi en cada ecuación, generando las matrices Zi

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

Se grafica cada ecuación como una «malla de alambre» o wireframe, usando librerias 3D, usando Xi,Yi, Z0, luego con Z1 y Z2

El Valor de la solución se grafica como un punto usando «dispersión» o scatter.

El algoritmo completo se muestra a continuación:

# 1ra Evaluación II Término 2011
# Tema 2. Sistema de ecuaciones 3x3
# Desarrollo conceptual como interseccion de Planos

import numpy as np

# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

a = -5
b = 5
muestras = 21

# PROCEDIMIENTO --------
# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

xi = np.linspace(a,b, muestras)
yi = np.linspace(a,b, muestras)
Xi, Yi = np.meshgrid(xi,yi)

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

# solución al sistema
punto = np.linalg.solve(A,B)

# SALIDA
print('respuesta de A.X=B : ')
print(punto)

# GRAFICA de planos
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

figura = plt.figure()
grafica = figura.add_subplot(111, projection='3d')

grafica.plot_wireframe(Xi,Yi,Z0, color ='blue',   label='Ecuación 1')
grafica.plot_wireframe(Xi,Yi,Z1, color ='green',  label='Ecuación 2')
grafica.plot_wireframe(Xi,Yi,Z2, color ='orange', label='Ecuación 3')

grafica.scatter(punto[0],punto[1],punto[2],
                color = 'red', marker='o', label ='punto')

grafica.set_title('Sistema de ecuaciones 3x3')
grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()
plt.show()

3.4 Gauss-Jordan Método

Referencia: Chapra 9.7 p277 pdf301

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

  • eliminación hacia adelante de incógnitas
  • eliminación hacia atras de incógnitas

sigue el ejemplo01 para matrices, usando la forma de función para el algoritmo.

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

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

Se crea la matriz aumentada AB.
Para el procedo de eliminación hacia adelante, en cada fila el elemento diagonal será el pivote, con el que se crea el coeficiente para eliminación de términos hacia adelante.

Para la eliminación hacia atras, se normaliza la fila haciendo el elemento diagonal uno, para calcular el coeficiente.

matriz aumentada:
[[  4.    2.    5.   60.7]
 [  2.    5.    8.   92.9]
 [  5.    4.    3.   56.3]]
 *** Gauss elimina hacia adelante ***
coeficiente:  2.0
[[   4.     2.     5.    60.7]
 [   0.     8.    11.   125.1]
 [   5.     4.     3.    56.3]]
coeficiente:  0.8
[[   4.      2.      5.     60.7 ]
 [   0.      8.     11.    125.1 ]
 [   0.      1.2    -2.6   -15.66]]
coeficiente:  6.66666666667
[[   4.            2.            5.           60.7       ]
 [   0.            8.           11.          125.1       ]
 [   0.            0.          -28.33333333 -229.5       ]]
 *** Gauss-Jordan elimina hacia atras *** 
coeficiente:  0.0909090909091
[[  4.           2.           5.          60.7       ]
 [  0.           0.72727273   0.           3.27272727]
 [ -0.          -0.           1.           8.1       ]]
coeficiente:  0.2
[[ 0.8         0.4         0.          4.04      ]
 [ 0.          0.72727273  0.          3.27272727]
 [-0.         -0.          1.          8.1       ]]
coeficiente:  2.5
[[ 2.   0.   0.   5.6]
 [ 0.   1.   0.   4.5]
 [-0.  -0.   1.   8.1]]
aumentada: 
[[ 1.   0.   0.   2.8]
 [ 0.   1.   0.   4.5]
 [-0.  -0.   1.   8.1]]
X: 
[[ 2.8]
 [ 4.5]
 [ 8.1]]
>>> 

El algoritmo desarrollado por partes:

# Gauss Jordan
# Supone que los elementos de diagonal no son cero
# Tarea: verificar diagonal no tiene ceros
# Muestra proceso
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)
print('matriz aumentada:')
print(AB)

print(' *** Gauss elimina hacia adelante ***')
casicero = 0 # 1e-15
# Gauss elimina hacia adelante
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(AB[k,i])>=casicero):
            coeficiente = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*coeficiente - AB[i,:]
        else:
            coeficiente= 'division para cero'
        print('coeficiente: ',coeficiente)
        print(AB)

print(' *** Gauss-Jordan elimina hacia atras *** ')
# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            coeficiente = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*coeficiente - AB[i,:]
        else:
            coeficiente= 'division para cero'
        print('coeficiente: ', coeficiente)
        print(AB)
X = AB[:,ultcolumna]
X = np.transpose([X])

# SALIDA
print('aumentada: ')
print(AB)
print('X: ')
print(X)

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

3.3.1 Gauss – Determinante

Determinante calculado de la forma triangular de la matriz, multiplicando todos los coeficientes de la diagonal.

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

calculando el determinante se llega al resultado de:

determinante de A: 
210.353

Primero se pivotea la fila, se obtiene la forma LU de la matriz, usando solo la triangular superior, se multiplican todos los coeficientes de la diagonal, con lo que se obtiene el determinante.

# determinante de A
# Referencia: Rodriguez 4.3.9 p.129
import numpy as np

def determinante(A):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    C = np.copy(A)
    C, cuenta = pivoteafila(C, contar=1)
    C_LU = triang_LU(C)
    U = C_LU[1]
    detA = 1
    for i in range(0,n,1):
        detA = detA*U[i,i]
    detA = detA*((-1)**cuenta) 
    return(detA)

# se modifica el pivoteo para que cuente los cambios de fila.
def pivoteafila(A, contar = 0):
    # Si contar=1, cuenta los cambios de fila
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    C = np.copy(A)
    cambio =0
    # mayor en diagonal
    for i in range(0,n-1,1):
        imax = np.argmax(np.abs(C[i:,i]))
        # intercambia si imax no es la diagonal
        if ((imax+i)!= i):
            # copia a temporal para intercambiar fila
            temporal = np.copy(C[i,:])
            C[i,:] = C[imax+i,:]
            C[imax+i,:] = temporal
            cambio = cambio+1
    if (contar==1):
        C=[C,cambio]
    return(C)

# Triangulares de una matriz
def triang_LU(A, casicero = 1e-15):
    # Triangulares de una matriz A=L.U
    # Triangular inferior L = LU[0]
    # Triangular superior U = LU[1]
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]

    U = np.copy(A)
    L = np.identity(n,dtype=float)

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

    LU = np.array([L,U])
    return(LU)

# Programa de prueba #######
# INGRESO
A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])


# PROCEDIMIENTO
detA = determinante(A)

# SALIDA
print('determinante de A: ')
print(detA)

verifique usando la función de numpy:

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

3.5 Matrices triangulares A=L.U

Referencia: Chapra 10.1 p282, pdf306

Una matriz A puede separarse en dos matrices: L y U que son de tipo triangular inferior L y triangular superior U.

Tienen la propiedad que:  A = L.U

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

se convierte en:

triangular inferior L
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
triangular superior U
[[  3.          -0.1         -0.2       ]
 [  0.           7.00333333  -0.29333333]
 [  0.           0.          10.01204188]]

Para mantener la matriz original, se obtiene una copia de la matriz A denominada U.

Se aplica el algoritmo de eliminación hacia adelante de Gauss, determinando  el factor con el que se multiplica la fila para hacer cero. El factor se lo almacena en la matriz U, en la posición de dónde se determinó el factor.

Las matrices resultantes son L y U, que se presentan en una matriz de dos elementos LU[0] y LU[1] para entregar la respuesta unificada.

 Triangulares de una matriz LU
# Triangular inferior L = LU[0]
# Triangular superior U = LU[1]
# Referencia: Chapra 10.1.2. p.284, pdf.308
# Tarea: verificar si el pivote es cero
#        verificar diagonal dominante
import numpy as np

# Programa de prueba
# INGRESO
A = np.array([[3  , -0.1, -0.2],
              [0.1,  7  , -0.3],
              [0.3, -0.2, 10  ]])

# PROCEDIMIENTO
casicero = 1e-15 # 0
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]

U = np.copy(A)
L = np.identity(n,dtype=float)

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

LU = np.array([L,U])

# SALIDA
print('triangular inferior L')
print(LU[0])
print('triangular superior U')
print(LU[1])

Tarea: Verificar los resultados, y considerar las divisiones para cero.

Verificando resultados: A=L.U

>>> np.dot(LU[0],LU[1])
array([[  3. ,  -0.1,  -0.2],
       [  0.1,   7. ,  -0.3],
       [  0.3,  -0.2,  10. ]])

3.2 Pivoteo parcial por filas

Referencia: Chapra 9.4.2 p268, pdf 292.

Para minimizar errores en la solución de matrices, se prefiere tener el valor mayor de cada fila en la diagonal, denominado «pivote».

Para los algoritmos, el pivoteo se realiza sobre la matriz aumentada, para así ordenar también los valores del vector B.

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


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

Matriz Aumentada, se realiza al concatenar la matriz A con el vector B por medio de columnas. Observe que el vector B se utilza en forma de columna.

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

el resultado de AB corresponde a:

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

El pivoteo parcial corresponde al procedimiento realizado por filas en la matriz AB, por ejemplo:

Para el pivoteo por fila, para la primera fila se inicia revisando la primera columna, desde la diagonal en adelante y se observa dónde se encuentra la mayor magnitud escalar (usando valor absoluto).

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

El mayor debe encontrarse en la primera posición de la columna, que corresponde a la diagonal de la matriz.  de no se así, debe intercambiar las filas de la matriz.

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

Para la segunda fila , se toma la segunda columna desde la diagonal en adelante y se repite el proceso.

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

que repetido para todas las filas de la matriz, tiene como resultado:

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

En el algoritmo propuesto, se usa una copia de la matriz, de tal manera que se conserve la matriz original para comparar con el resultado.

Si no se requiere la matriz original, no se realiza la copia y se ahorra el espacio de memoria; detalle importante para matrices de «gran» tamaño.

# Pivotea por filas, Pivoteo parcial
# Recibe la matriz A, copia en C para no modificar A
# Si hay ceros en diagonal, la matriz es singular,
# Tarea: Revisar si diagonal tiene ceros
import numpy as np

def pivoteafila(M):
    tamano = np.shape(M)
    n = tamano[0]
    m = tamano[1]
    C = np.copy(M)
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(C[i:,i])
        dondemax = np.argmax(columna)
        # revisa dondemax no está en la diagonal
        if (dondemax != 0):
            # intercambia fila
            temporal = np.copy(C[i,:])
            C[i,:] = C[dondemax+i,:]
            C[dondemax+i,:] = temporal
    return(C)

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

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

# PROCEDIMIENTO
AB = np.concatenate((A,B),axis=1)
respuesta = pivoteafila(AB)
# SALIDA
print('matriz pivoteada por fila:')
print(respuesta)

3.4.2 Gauss-Jordan Inversa Ejemplo01

Referencia: Chapra 10.2 p292, pdf 316; Rodriguez Cap.4.2.5 Ejemplo 1 pdf.118

Obtener la inversa de una matriz usando el método de Gauss-Jordan, a partir de la matriz:

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

Para el proceso se utiliza la matriz aumentada de A con la identidad.

AI = A|I

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

Con la matriz aumentada AI  se repite el proceso de Gauss-Jordan.

La inversa se encuentra en la mitad derecha de AI, lugar que originalmente correspondía a la identidad .

El algoritmo que describe el proceso en python:

# Método de Gauss para inversa de A
# 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]])

# PROCEDIMIENTO
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]

# Añade la matriz identidad
identidad = np.identity(n)
AI = np.concatenate((A,identidad), axis=1)
m = 2*m

# Gauss elimina hacia adelante
# tarea: verificar términos cero
casicero = 1e-15
for i in range(0,n,1):
    pivote = AI[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = AI[k,i]/pivote
            AI[k,:] = AI[k,:] - factor*AI[i,:]
        else:
            factor = 'division para cero'
# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AI[i,:] = AI[i,:]/AI[i,i]
    pivote = AI[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(pivote)>=casicero):
            factor = AI[k,i]/pivote
            AI[k,:] = AI[k,:] - factor*AI[i,:]
        else:
            factor= 'division para cero'

inversa = AI[:,n:]

# Para verificar el resultado
verifica = np.dot(A,inversa)

# SALIDA
print('la matriz inversa es:')
print(inversa)

print('verificando A.inversa = identidad')
print(verifica)

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   2.77555756e-17   0.00000000e+00]
 [  0.00000000e+00   1.00000000e+00   0.00000000e+00]
 [ -5.55111512e-17   8.32667268e-17   1.00000000e+00]]

Observe que el algoritmo se pude reducir si usa el proceso de Gauss-Jordan como una función.

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

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 y la solución obtenida: es

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

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

por sustitución hacia atras
el vector solución X es:
[[2.8]
 [4.5]
 [8.1]]
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]])

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

# PROCEDIMIENTO
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('aumentada: ')
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.

3.4.1 Gauss-Jordan Ejemplo01

Ref: Rodriguez Cap.4 Ejemplo 1 pdf.105

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 A B C 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 x1, x2, x3 para generar el sistema de ecuaciones lineales con tres variables:

4 x1 + 2 x2 + 5 x3 = 60.70
2 x1 + 5 x2 + 8 x3 = 92.90
5 x1 + 4 x2 + 3 x3 = 56.30

para resolver:
\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}

que para la computadora será:

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

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

algoritmo propuesto, tiene como resultado:

matriz aumentada: 
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Elimina hacia adelante
[[  4.        2.        5.       60.7    ]
 [  0.        4.        5.5      62.55   ]
 [  0.        0.       -5.3125  -43.03125]]
Elimina hacia atras
[[ 1.   0.   0.   2.8]
 [ 0.   1.   0.   4.5]
 [-0.  -0.   1.   8.1]]
el vector solución X es:
[[2.8]
 [4.5]
 [8.1]]
verificar que A.X = B
[[60.7]
 [92.9]
 [56.3]]
# Método de Gauss-Jordan ,
# 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]])

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

# PROCEDIMIENTO
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('aumentada: ')
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)

# Gauss elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            factor = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*factor - AB[i,:]
        else:
            factor= 'division para cero'
print('Elimina hacia atras')
print(AB)

X = AB[:,ultcolumna]

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

# SALIDA
print('el vector solución X es:')
print(np.transpose([X]))

print('verificar que A.X = B')
print(np.transpose([verifica]))

Finalmente, si verificamos los resultados multiplicando la matriz de cantidad por el vector resultante deberá obtener el vector de valor pagado:

>>> np.dot(cantidad,precio)
array([[ 60.7],
       [ 92.9],
       [ 56.3]])

una función de numpy que resuelve el sistema de ecuaciones:

>>> np.linalg.solve(cantidad, pagado)
array([[ 2.8],
       [ 4.5],
       [ 8.1]])