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 ]

3.2 Pivoteo parcial por filas con Python

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

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

Matriz PivoteadaLos 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ículosvende frutas

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

3.1 Sistema de ecuaciones 3×3, planos 3D con Python

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]

Referencia: Chapra 9.1 p247

Como punto de partida para la unidad de sistema de ecuaciones se usa un ejemplo para ilustrar con gráficos la solución para : 3 ecuaciones y 3 incógnitas.

La solución es un punto de intersección de los planos en el espacio dados por cada ecuación en el sistema.

La intersección entre dos planos es una recta, que al añadir un plano adicional la intersección es un punto. Las coordenadas del punto es la solución buscada.

Para generalizar y usar en otros casos, se realiza el desarrollo analítico del tema, junto al algoritmo y las gráficas en 3D en Python y Matplotlib.

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]
..


1. Ejercicio

Referencia: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

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}

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]
..


2. Desarrollo analítico. Ecuaciones como Planos en el espacio

Si observamos la primera ecuación en una gráfica, encontramos que para graficar el plano se requiere definir los intervalos para los ejes X, Y, por ejemplo:

z=\frac{1+ 2x - 5y}{9} -5 \leq x \leq 5 -7 \leq y \leq 7

Al incorporar el plano de la segunda ecuación encontramos que la solución al sistema es la recta intersección de ambos planos.

Usando la tercera ecuación, la intersección de los planos genera un punto cuyas coordenadas corresponden a la solución del sistema.

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]
..


3. Algoritmo en Python

Para observar mejor del resultado en tres dimensiones, se plantean algunas instrucciones básicas para realizar la gráfica y obtener una mejor vista rotando el gráfico con el cursor.

Las ecuaciones del problema se pueden escribir en forma matricial Ax=B, usando solamente los coeficientes.

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}
# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

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

Con la forma matricial se generan las ecuaciones de cada plano usando los coeficientes al despejar la variable z.

Para visualizar los planos generados por cada ecuación se escribe cada ecuación como funciones(x,y) de dos variables. Para simplificar la escritura se usa el formato lambda.

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

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]
..


4. Gráficos. Intervalos de observación

Para la gráfica se requiere usar intervalos y muestras para las variables independientes x,y de forma semejante a las gráficas en 2D.

Las muestras en cada eje generan puntos en el plano X,Y y luego se usan todas las combinaciones que se generan.

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

ax = -5    # Intervalo [a,b] para eje x
bx = 5
ay = ax-2  # Intervalo [a,b] para eje y
by = bx+2
muestras = 11

xi = np.linspace(ax,bx, muestras)
yi = np.linspace(ay,by, muestras)

con lo que se obtiene como resultado:

>>> xi
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.])
>>> yi
array([-7. , -5.6, -4.2, -2.8, -1.4,  0. ,  1.4,  2.8,  4.2,  5.6,  7. ])
>>> 

– 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)

el resultado son las matrices con los puntos de evaluación en la malla. Por ejemplo el primer punto [x0,y0] = [-5,-7]:

>>> Xi
array([[-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.],
       [-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.]])
>>> Yi
array([[-7. , -7. , -7. , -7. , -7. , -7. , -7. , -7. , -7. , -7. , -7. ],
       [-5.6, -5.6, -5.6, -5.6, -5.6, -5.6, -5.6, -5.6, -5.6, -5.6, -5.6],
       [-4.2, -4.2, -4.2, -4.2, -4.2, -4.2, -4.2, -4.2, -4.2, -4.2, -4.2],
       [-2.8, -2.8, -2.8, -2.8, -2.8, -2.8, -2.8, -2.8, -2.8, -2.8, -2.8],
       [-1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4, -1.4],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
       [ 1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4],
       [ 2.8,  2.8,  2.8,  2.8,  2.8,  2.8,  2.8,  2.8,  2.8,  2.8,  2.8],
       [ 4.2,  4.2,  4.2,  4.2,  4.2,  4.2,  4.2,  4.2,  4.2,  4.2,  4.2],
       [ 5.6,  5.6,  5.6,  5.6,  5.6,  5.6,  5.6,  5.6,  5.6,  5.6,  5.6],
       [ 7. ,  7. ,  7. ,  7. ,  7. ,  7. ,  7. ,  7. ,  7. ,  7. ,  7. ]])
>>> 

– Se evalúan los puntos Xi,Yi en cada ecuación, que dan como resultado las matrices Zi

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

Nota: Si es necesario, revise la sección de Gráficas 3D wireframe en la sección de Actividades y Recursos/Resumen Python.

El gráfico para cada ecuación ser realiza con la instrucción wireframe (malla de alambre) de las librerías 3D, se usa como entrada para el primer plano Xi,Yi,Z0. Luego se repite el procedimiento para Z1 y Z2.

Para simplicidad del algoritmo, pues la revisión del concepto es sobre intersección de planos, la solución se obtiene con las funciones de Numpy.

punto = np.linalg.solve(A,B)

El punto se incorpora a la gráfica como un punto usando scatter (dispersión).

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]
..


5. Instrucciones en Python

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

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

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

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

ax = -5     # Intervalo X
bx = 5
ay = ax-2   # Intervalo Y
by = bx+2
muestras = 11

# 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(ax,bx, muestras)
yi = np.linspace(ay,by, 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)

# Interseccion entre ecuación 1 y 2
# PlanoXZ, extremo inferior de y
Aa  = np.copy(A[0:2,[0,2]])
Ba  = np.copy(B[0:2])
Ba  = Ba-ay*A[0:2,1]
pta = np.linalg.solve(Aa,Ba)
pa  = np.array([ay])
pxza = np.array([pta[0],ay,pta[1]])

# PlanoXZ, extremo superior de y
Ba  = Ba-by*A[0:2,1]
ptb = np.linalg.solve(Aa,Ba)
pb  = np.array([by])
pxzb = np.array([ptb[0],by,ptb[1]])

# GRAFICA de planos
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111, projection='3d')

graf3D.plot_wireframe(Xi,Yi,Z0,
                       color ='blue',
                       label='Ecuación 1')
graf3D.plot_wireframe(Xi,Yi,Z1,
                       color ='green',
                       label='Ecuación 2')
graf3D.plot_wireframe(Xi,Yi,Z2,
                       color ='orange',
                       label='Ecuación 3')
# recta intersección planos 1 y 2
graf3D.plot([pxza[0],pxzb[0]],
             [pxza[1],pxzb[1]],
             [pxza[2],pxzb[2]],
             label='Sol 1y2',
             color = 'violet',
             linewidth = 4)
# Punto solución del sistema 3x3
graf3D.scatter(punto[0],punto[1],punto[2],
                color = 'red',
                marker='o',
                label ='punto',
                linewidth = 6)

graf3D.set_title('Sistema de ecuaciones 3x3')
graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_xlim(ax, bx)
graf3D.set_ylim(ay, by)
graf3D.legend()
graf3D.view_init(45, 45)
# rotacion de ejes
for angulo in range(45, 360+45, 5 ):
    graf3D.view_init(45, angulo)
    plt.draw()
    plt.pause(.001)
plt.show()

Ejemplo: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfico ] [ Python ]