[ norma 3D ] [ acople aeronave ] [ algoritmo ]

..

---

Normas de un vector en 3D

Referencia: Chapra 10.3.1 p298, Burden 7.1 p320, Rodríguez 4.4.1 p132

La norma de un vector se interpreta como una distancia entre la coordenada definida por el vector [x_i, y_i, z_i] y el origen [0,0,0]. También se puede realizar respecto a otro punto de referencia, se conoce como Norma Euclidiana o Norma p=2.

Previamente, en búsqueda de raíces, el error de aproximación se considera como la diferencia entre dos iteraciones sucesivas:.

error = |x_{i+1} - x_{i}|

Si el concepto se extiende a vectores en tres dimensiones, se observa como el error entre vectores |X_i+1 - X_i| que se interpreta mejor como una distancia.

Por ejemplo, si se usan los puntos y se visualizan en un gráfico:

    X2:  [ 2  4 -1]
   -X1: -[ 1  2  3]
______________________
errado =|[ 1  2 -4]|

La diferencia entre los vectores |[ 1, 2, -4]|  es más simple de observar como un número escalar. Una forma de convertir el punto a un escalar es usar la distancia al origen.

errado = \sqrt{1^2 + 2^2 + (-4)^2} = 4.58

El resultado también se interpreta como la distancia entre los puntos X₁ y X₂.

De ampliar el concepto a n dimensiones se conoce como norma de un vector o matriz.

||x|| = \sqrt{x^2+y^2+z^2} ||x|| = \Big[ \sum_{i=0}^{n} x_i ^2 \Big] ^{1/2}

[ norma 3D ] [ acople aeronave ] [ algoritmo ]
..

---

2. Distancia entre dos puntos en el espacio, error y Norma

Observe los siguientes videos, sobre acople de aeronaves.

Acoplamiento de aviones para recarga de combustible . www.AiirSource.com. 30/diciembre/2015.
KC-135 Stratotanker in Action - Aircraft Air Refueling

El carguero ruso Progress M-06M pasó de largo la Estación Espacial Internacional fracasado en su intento de acoplarse

Soyuz MS-22 docking. 21 sept 2022.

• ¿Considera importante que exista acoplamiento para iniciar la carga el combustible? o ¿para abrir la escotilla del transbordador?

Si el error de acoplamiento entre artefactos se calcula como la Norma entre los puntos de "contacto",

• ¿es necesario calcular la raíz cuadrada de los cuadrados de las diferencias? o,
• ¿Solo toma en cuenta el mínimo o el máximo de las diferencias entre las coordenadas?,
• ¿cuál de las formas tiene menos operaciones, es más rápida de realizar?

considere sus respuestas para el siguiente concepto.

[ norma 3D ] [ acople aeronave ] [ algoritmo ]

---

3. Norma infinito

Se determina como el valor mas grande entre los elementos del vector o matriz.

||x|| = max\Big[ |x_i| \Big]

Es más sencilla de calcular que la Norma Euclidiana, Norma P=2, mostrada al principio.

Se interpreta como si alguna de las diferencia entre las coordenadas de los puntos de acople entre aeronaves es mayor que lo tolerado, no se debe permitir abrir la válvula de combustible o la escotilla del transbordador. No es prioritario calcular la suma de los cuadrados de las diferencias para saber que no existe acople aún.

Existen otras formas de calcular la Norma, que se presentan en la siguiente página web.

[ norma 3D ] [ acople aeronave ] [ algoritmo ]

..

---

4. Algoritmo con Python

Principalmente se usa para generar las gráficas 3D, que ayudan a la interpretación del concepto con los puntos de coordenadas:

X0 = np.array([0.0, 0, 0])
X1 = np.array([1.0, 2, 3])
X2 = np.array([2.0, 4,-1])

Las instrucciones gráficas principales son:

• Para puntos, gráficos de dispersión o scatter()
• Para las lineas, gráficos de línea o plot()

El resultado de la parte numérica se obtiene con pocas instrucciones en el bloque de procedimiento

     X2:  [ 2 4 -1]
    -X1: -[ 1 2  3]

____________________
errado =  [ 1 2 -4]
||errado|| =  4.58257569495584
Norma euclidiana :  4.58257569495584

las instrucciones en Python son:

# Norma como error
# o distancia entre dos puntos
# caso 3D
import numpy as np

# INGRESO
X0 = np.array([0.0, 0, 0])
X1 = np.array([1.0, 2, 3])
X2 = np.array([2.0, 4,-1])

# PROCEDIMIENTO
errado = X2 - X1
distancia = np.sqrt(np.sum(errado**2))
# funciones numpy
Nerrado = np.linalg.norm(errado)

# SALIDA
print('X1 = ', X1)
print('X2 = ', X2)
print('errado = ', errado)
print('||errado|| = ', distancia)
print('Norma euclidiana : ',Nerrado)


# Grafica
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
figura  = plt.figure()
grafica = figura.add_subplot(111,projection = '3d')

# puntos en el espacio
[x, y , z] = X0
grafica.scatter(x,y,z, c = 'blue',
                marker='o', label = 'X0')

[x, y , z] = X1
grafica.scatter(x,y,z, c = 'red',
                marker='o', label = 'X1')

[x, y , z] = X2
grafica.scatter(x,y,z, c = 'green',
                marker='o', label = 'X2')

# líneas entre puntos
linea = np.concatenate(([X0],[X1]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||X1||')

linea = np.concatenate(([X0],[X2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||X2||')

linea = np.concatenate(([X1],[X2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
grafica.plot(x,y,z, label = '||error||')

grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()

grafica.view_init(35, 25)
plt.show()

---

[ norma 3D ] [ acople aeronave ] [ algoritmo ]

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]
..

---

1. Método de Gauss-Jordan

Referencia: Chapra 9.7 p277, Burden 9Ed Ex6.1.12 p370, Rodríguez 4.2 p106

El método de Gauss-Jordan presenta un procedimiento alterno al de "sustitución hacia atrás" realizado para el método de Gauss. Luego de obtener la "matriz triangular superior", aplica el procedimiento de "Eliminación hacia atrás".

\left( \begin{array}{rrr|r} a_{0,0} & a_{0,1} & a_{0,2} & b_{0} \\ 0 & a_{1,1} & a_{1,2} & b_{1} \\ 0 & 0 & a_{2,2} & b_{2} \end{array} \right)

• Eliminación hacia atrás

Desde la última fila,
para i = n-1, n-2, ...
sobre las filas k = i-1 superiores

a_{k} = a_{k} -a_{i}\frac{a_{ki}}{a_{ii}}

Se mantienen al inicio los procedimientos para:
matriz aumentada,
pivoteo parcial por filas,
eliminación hacia adelante desarrollados anteriormente.

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]
..

---

2. Ejercicio

Referencia: Rodríguez 4.0 p105,  1Eva_IT2010_T3_MN Precio artículos

\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}

Se continúa con el ejercicio desarrollado para el método de Gauss desde la "matriz triangular superior",
luego de eliminación hacia adelante:

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

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]
..

---

3. Eliminación hacia atrás

El procedimiento es semejante al realizado para "eliminación hacia adelante", con la diferencia que se inicia en la última fila hacia la primera. 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

Se aplica desde la última fila i, para las otras filas k que se encuentran hacia atrás.

i = ultfila = 2
pivote = AB[2,2] = 5
k = i-1 # hacia atrás

se realizan las operaciones entre las filas i y la fila j

         [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 8.8x10^-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 índice k

k = k-1 = 2-1 = 1 # hacia atrás

se realizan las operaciones entre las filas i y la fila j

       [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 x₃ 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    4.0    0.0     32.0]
-(4/3.4)*[0.0    3.4    8.8e-16 15.3]
__________________________________
       = [5.0    0.0   -1.04e-15 14.0]

[[ 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 x₂,
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]]

iteración 4, operación fila 1

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

solución X: 
[2.8, 4.5, 8.1]

X= [2.8, 4.5 , 8.1 ]

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.

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]
..

---

4. 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 atrás para completar la solución al sistema de ecuaciones.

Instrucciones en Python

# Método de Gauss-Jordan
# Sistemas de Ecuaciones A.X=B
import numpy as np

# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]
B = [60.70, 92.90, 56.30]

# PROCEDIMIENTO
np.set_printoptions(precision=4) # 4 decimales en print
casicero = 1e-15  # redondear a cero
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)

# Matriz aumentada AB
B_columna = np.transpose([B])
AB  = np.concatenate((A,B_columna),axis=1)

# 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)
    
    if (dondemax !=0): # NO en diagonal
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal
print('Pivoteo parcial por filas AB: ----')
print(AB)

# Eliminación hacia adelante
print('Eliminación hacia adelante: ------')
for i in range(0,n-1,1): # una fila
    pivote   = AB[i,i]
    adelante = i + 1
    for k in range(adelante,n,1): # diagonal adelante
        factor  = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor

        print('i:',i,'k:',k, 'factor:',factor)
    print(AB)

# Eliminación hacia atras
print('Eliminación hacia atras: ------')
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1): # una fila
    pivote   = AB[i,i]
    print('i:',i,' pivote:',pivote)
    
    atras = i - 1
    for k in range(atras,0-1,-1): # diagonal adelante
        factor  = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor
        for j in range(0,m,1): # casicero revisa
            if abs(AB[k,j])<casicero:
                AB[k,j]=0
        
        print(' k:',k, 'factor:',factor)
    print(AB)
    
    AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
    
print('AB:')
print(AB)
X = np.copy(AB[:,ultcolumna])

# SALIDA
print('Método de Gauss-Jordan')
print('solución X: ')
print(X)

El resultado del algoritmo es:

Pivoteo parcial por filas AB: ----
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Eliminación hacia adelante: ------
i: 0 k: 1 factor: 0.4
i: 0 k: 2 factor: 0.8
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]
i: 1 k: 2 factor: -0.3529411764705883
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Eliminación hacia atras: ------
i: 2  pivote: 5.0
 k: 1 factor: 1.3599999999999999
 k: 0 factor: 0.6
[[ 5.   4.   0.  32. ]
 [ 0.   3.4  0.  15.3]
 [ 0.   0.   5.  40.5]]
i: 1  pivote: 3.4
 k: 0 factor: 1.1764705882352942
[[ 5.   0.   0.  14. ]
 [ 0.   3.4  0.  15.3]
 [ 0.   0.   1.   8.1]]
i: 0  pivote: 5.0
[[ 5.   0.   0.  14. ]
 [ 0.   1.   0.   4.5]
 [ 0.   0.   1.   8.1]]
AB:
[[1.  0.  0.  2.8]
 [0.  1.  0.  4.5]
 [0.  0.  1.  8.1]]
Método de Gauss-Jordan
solución X: 
[2.8 4.5 8.1]

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]
..

---

5. Algoritmo como función en Python

El algoritmo considera que valores con magnitud menores a "casicero" se redondean a 0. Se aplica recorriendo la fila de referencia k, por cada elemento j.

# redondeo a cero
for j in range(0,m,1): 
    if np.abs(AB[k,j])<=casicero:
        AB[k,j]=0

Adicionalmente se añade el parámetro "inversa", que se usa para encontrar A^-1 como se indica en Matriz Inversa.

Con el algoritmo se obtiene el siguiente resultado:

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 i: 0  pivote: 5.0
  fila k: 1  factor: 0.4
  fila k: 2  factor: 0.8
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]
 fila i: 1  pivote: 3.4
  fila k: 2  factor: -0.3529411764705883
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
 fila i: 2  pivote: 5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Elimina hacia Atrás:
 fila i: 2  pivote: 5.0
  fila k: 1  factor: 1.3599999999999999
  fila k: 0  factor: 0.6
[[ 5.   4.   0.  32. ]
 [ 0.   3.4  0.  15.3]
 [ 0.   0.   1.   8.1]]
 fila i: 1  pivote: 3.4
  fila k: 0  factor: 1.1764705882352942
[[ 5.   0.   0.  14. ]
 [ 0.   1.   0.   4.5]
 [ 0.   0.   1.   8.1]]
 fila i: 0  pivote: 5.0
[[1.  0.  0.  2.8]
 [0.  1.  0.  4.5]
 [0.  0.  1.  8.1]]
solución X: 
[2.8 4.5 8.1]

Instrucciones en Python

# Método de Gauss-Jordan
# Sistemas de Ecuaciones A.X=B
import numpy as np

def gauss_eliminaAtras(AB, vertabla=False,
                       inversa=False,
                       casicero = 1e-15):
    ''' Gauss-Jordan elimina hacia atrás
    Requiere la matriz triangular inferior
    Tarea: Verificar que sea triangular inferior
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    ultfila = n-1
    ultcolumna = m-1
    if vertabla==True:
        print('Elimina hacia Atrás:')
        
    for i in range(ultfila,0-1,-1):
        pivote = AB[i,i]
        atras = i-1  # arriba de la fila i
        if vertabla==True:
            print(' fila i:',i,' pivote:', pivote)
            
        for k in range(atras,0-1,-1):
            if np.abs(AB[k,i])>=casicero:
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                
                # redondeo a cero
                for j in range(0,m,1): 
                    if np.abs(AB[k,j])<=casicero:
                        AB[k,j]=0
                if vertabla==True:
                    print('  fila k:',k,
                          ' factor:',factor)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
        AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
        
        if vertabla==True:
            print(AB)
    
    respuesta = np.copy(AB[:,ultcolumna])
    if inversa==True: # matriz inversa
        respuesta = np.copy(AB[:,n:])
    return(respuesta)

def gauss_eliminaAdelante(AB,vertabla=False,
                          lu=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:',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,:]
                for j in range(0,m,1): # casicero revisa
                    if abs(AB[k,j])<casicero:
                        AB[k,j]=0
                if vertabla==True:
                    print('  fila k:',k,
                          ' factor:',factor)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
        if vertabla==True:
            print(AB)
    respuesta = np.copy(AB)
    if lu==True: # matriz triangular A=L.U
        U = AB[:,:n-1]
        respuesta = [AB,L,U]
    return(respuesta)

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

# PROGRAMA ------------------------
# INGRESO
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

B = [60.70,92.90,56.30]

# PROCEDIMIENTO
np.set_printoptions(precision=4) # 4 decimales en print
casicero = 1e-15  # redondear a cero

AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_eliminaAtras(AB,vertabla=True)

# SALIDA
print('solución X: ')
print(X)

---

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

[ Gauss-Jordan ] [ Ejercicio ] [ Eliminación atrás ] [ Algoritmo ] [ función ]

[ 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:

\begin{cases} a_{0,0} x_0+a_{0,1}x_1+a_{0,2} x_2 = b_{0} \\ a_{1,0} x_0+a_{1,1}x_1+a_{1,2} x_2 = b_{1} \\ a_{2,0} x_0+a_{2,1}x_1+a_{2,2} x_2 = b_{2} \end{cases}

1.1 Eliminación hacia adelante

Desde la primera fila i ,
sobre las filas inferiores

k = i+1, i+2, ...

a_{k} = a_{k} -a_{i}\frac{a_{ki}}{a_{ii}}

1.2 Sustitución hacia atrás

Desde la última fila, para

i = n-1, n-2, ...

x_i = \Bigg( b_i-\sum_{j=i+1}^{n} a_{ij} x_{j} \Bigg) \frac{1}{a_{ii}}

y las columnas j luego de la diagonal.

[ 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

Para el sistema de ecuaciones mostrado, desarrolle usando el método de Gauss.

\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}

Se aplica el pivoteo parcial por filas, que tiene como resultado:

\left( \begin{array}{rrr|r} 5 & 4 & 3 & 56.3 \\ 2 & 5 & 8 & 92.9 \\ 4 & 2 & 5 & 60.7 \end{array} \right)

[[ 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:

\left( \begin{array}{rrr|r} 5 & 4 & 3 & 56.3 \\ 0 & 3.4 & 6.8 & 70.38 \\ 0 & 0 & 5 & 40.5 \end{array} \right)

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 ya 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 a_i,i como pivote.

pivote = AB[i,i] = AB[0,0] = 5

Para las filas 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 fórmula.

k = i+1 = 0+1 = 1
factor = AB[1,0]/pivote = 2/5

a_{k} = a_{k} -a_{i}\frac{a_{k,i}}{pivote}

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 = k+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 =1.

iteración fila 2
Para la fila 2, con posición i=1, se toma el elemento de la diagonal a_i,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. Para aplicar la fórmula por filas, se obtiene el factor .

k = i+1 = 1+1 = 2
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:

\left( \begin{array}{rrr|r} 5 & 4 & 3 & 56.3 \\ 0 & 3.4 & 6.8 & 70.38 \\ 0 & 0 & 5 & 40.5 \end{array} \right)

[ 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.
Desde la última fila, para i = n-1, n-2, ...
Y las columnas j luego de la diagonal

x_i = \Bigg( b_i-\sum_{j=i+1}^{n} a_{ij} x_{j} \Bigg) \frac{1}{a_{ii}}

Para una fila i, el vector B[i] representa el valor de la constante en la fila i de la matriz aumentada, A[i,j] se refiere los valores de los coeficientes de la ecuación, de 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 la resultante de eliminación hacia adelante:

\left( \begin{array}{rrr|r} 5 & 4 & 3 & 56.3 \\ 0 & 3.4 & 6.8 & 70.38 \\ 0 & 0 & 5 & 40.5 \end{array} \right)

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[i] = 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 = \frac{40.5}{5} = 8.1

la respuesta se interpreta en el vector X como:

X= [ x_1 , x_2 , x_3] = [ 0 , 0 , 8.1 ]

iteración 2, fila 2,  i = 1
De la penúltima fila se obtiene la ecuación para encontrar x₂

[ 0. 3.4 6.8 70.38]

0x_1 + 3.4 x_2 +6.8 x_3 = 70.38

se observa que B[i] = 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 3.4 x_2 = 70.38 -6.8 (8.1)

Usa el valor de x₃ encontrado en la iteración anterior. Muestra la ecuación para la segunda fila.

x_2 = \frac{70.38 -6.8 (8.1)}{3.4} = 4.5

la respuesta se interpreta en el vector X como:

X= [ 0 , 4.5 , 8.1 ]

iteración 3 fila 1, i=0
Se sigue el mismo proceso para la siguiente incógnita x₁ 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 (4.5) + 3(8.1))}{5} = 2.8

Se encuentra que la solución al sistema de ecuaciones es:

X= [ 2.8, 4.5 , 8.1 ]

Método de Gauss
solución X: 
[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.

>>> A
array([[4., 2., 5.],
       [2., 5., 8.],
       [5., 4., 3.]])
>>> X
array([2.8, 4.5, 8.1])
>>> np.dot(A,X)
array([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 número Real (float), para que no se considere el vector como entero y realice operaciones entre enteros, generando errores por truncamiento.

El resultado del algoritmo muestra los siguientes resultados:

Pivoteo parcial por filas AB: ----
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Eliminación hacia adelante: ------
i: 0 k: 1 factor: 0.4
i: 0 k: 2 factor: 0.8
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]
i: 1 k: 2 factor: -0.3529411764705883
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Método de Gauss
solución X: 
[2.8 4.5 8.1]

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
# Sistemas de Ecuaciones A.X=B
import numpy as np

# INGRESO
# NOTA: Para AB, se ha usado pivoteo parcial por filas
A = [[4,2,5],
     [2,5,8],
     [5,4,3]]
B = [60.70, 92.90, 56.30]

# PROCEDIMIENTO
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)

# Matriz aumentada AB
B_columna = np.transpose([B])
AB  = np.concatenate((A,B_columna),axis=1)

# 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)
    
    if (dondemax !=0): # NO en diagonal
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal
print('Pivoteo parcial por filas AB: ----')
print(AB)

# Eliminación hacia adelante
print('Eliminación hacia adelante: ------')
for i in range(0,n-1,1): # una fila
    pivote   = AB[i,i]
    adelante = i + 1
    for k in range(adelante,n,1): # diagonal adelante
        factor  = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor

        print('i:',i,'k:',k, 'factor:',factor)
    print(AB)

# 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 = AB[i,ultcolumna]
    for j in range(i+1,ultcolumna,1):
        suma = suma - AB[i,j]*X[j]
    X[i] = suma/AB[i,i]

# SALIDA
print('Método de Gauss')
print('solución X: ')
print(X)

---

Tarea

Revisar cuando la matriz pivoteada por filas tienen en la diagonal un elemento cero o muy cercano a 0, la matriz sería singular. El valor considerado como casi cero podría ser 1x10^-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

Se integra la función pivoteafila(), se crean las funciones para los procedimientos gauss_eliminaAdelante() y gauss_sustituyeAtras(), usando la variable vertabla para observar los resultados parciales.

# Método de Gauss
# Sistemas de Ecuaciones A.X=B
import numpy as np

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:',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,:]
                for j in range(0,m,1): # casicero revisa
                    if abs(AB[k,j])<casicero:
                        AB[k,j]=0
                if vertabla==True:
                    print('  fila k:',k,
                          ' factor:',factor)
            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
    '''
    n,m = np.shape(AB)
    # Sustitución hacia atras
    if vertabla==True:
        print('Sustitute hacia atras:')
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        suma = AB[i,ultcolumna]   # B[i]
        for j in range(i+1,ultcolumna,1):
            suma = suma - AB[i,j]*X[j]
        X[i] = suma/AB[i,i]    # suma/pivote
        if vertabla==True:
            print(' fila i:',i,
                  '   X['+str(i)+']:',X[i])
    return(X)

def pivoteafila(A,B,vertabla=False):
    '''
    Pivoteo 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,
                      'con', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print('AB:')
            print(AB)
    return(AB)

# PROGRAMA Búsqueda de solucion  --------
# 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('Método de Gauss')
print('solución X: ')
print(X)

El resultado para 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 con 2
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Elimina hacia adelante:
 fila i: 0  pivote: 5.0
  fila k: 1  factor: 0.4
  fila k: 2  factor: 0.8
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.   -1.2   2.6  15.66]]
 fila i: 1  pivote: 3.4
  fila k: 2  factor: -0.3529411764705883
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
 fila i: 2  pivote: 5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Sustitute hacia atras:
 fila i: 2    X[2]: 8.100000000000003
 fila i: 1    X[1]: 4.499999999999997
 fila i: 0    X[0]: 2.8
Método de Gauss
solución X: 
[2.8 4.5 8.1]

[ Gauss ] [ Ejercicio ] [ Eliminación adelante ] [ Sustitución atrás ]
[ Algoritmo ] [ función ]

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]

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

Los 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.

El objetivo del pivoteo es ordenar la matriz para que en lo posible sea diagonal dominante.

|a_{i,i}| \geq \sum_{j=0 ; i\neq j}^n |a_{i,j}|

..

---

1. Ejercicio

Referencia: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

Considere el siguiente sistema de ecuaciones AX=B dado por:

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

Realice el pivoteo parcial por filas, de tal manera que la diagonal de A sea estrictamente dominante.

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]
..

---

2. Matriz Aumentada

El sistema de ecuaciones se escribe en la forma algebraica como matrices y vectores de la forma Ax=B

\begin{pmatrix} -2 & 5 & 9 \\ 7 & 1 & 1 \\-3 & 7 & -1 \end{pmatrix} \begin{pmatrix} x_0 \\ x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 6 \\ -26 \end{pmatrix}

Para el algoritmo, las ecuaciones se escriben en forma de matriz A y vector B.

A = [[-2, 5, 9],
     [ 7, 1, 1],
     [-3, 7,-1]]

B = [1,6,-26]

Las matrices y vectores se procesan como arreglos np.array() de números reales. No se usan listas.

A = np.array(A,dtype=float)
B = np.array(B,dtype=float)

El vector B para la matriz aumentada se usa en forma de columna. Se aumenta una dimensión [B] y se aplica la transpuesta.

B_columna = np.transpose([B])

>>> B_columna
array([[  1],
       [  6],
       [-26]])

La matriz aumentada AB se forma al concatenar por columnas la matriz A con el vector B_columna,  usando el parámetro axis=1.

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

el resultado AB se muestra como:

>>> AB
array([[ -2. ,  5. ,   9. ,   1.],
       [  7. ,  1. ,   1. ,   6.],
       [ -3. ,  7. ,  -1. , -26.]])

\left( \begin{array}{rrr|r} -2 & 5 & 9 & 1 \\ 7 & 1 & 1 & 6 \\ -3 & 7 & -1 & -26 \end{array} \right)

En varios algoritmos se usa la matriz aumentada AB, para mantener el orden entre filas de A y B al aplicar cambios de fila.

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]
..

---

3. Pivoteo parcial por filas

Para el pivoteo por filas de la matriz aumentada AB, tiene como primer paso revisar la primera columna, desde la posición diagonal hacia todas las filas por debajo, usando solo la magnitud de los coeficientes |A|.

La representación gráfica de la matriz A, permite visualizar los intercambios de fila.

>>> AB
array([[ -2. ,  5. ,   9. ,   1.],
       [  7. ,  1. ,   1. ,   6.],
       [ -3. ,  7. ,  -1. , -26.]])

Las instrucciones para seleccionar la columna, usar el valor absoluto y dónde se encuentra el máximo son:

columna = abs(AB[i:,i])
dondemax = np.argmax(columna)

Los valores de la columna se reducen a:

columna = [|-2.|,   fila=0
           | 7.|,   fila=1
           |-3.|]   fila=2
dondemax = 1

Si la magnitud de mayor valor no está en la diagonal (dondemax≠0), se realiza el intercambio entre la fila dondemax y la fila de la diagonal i en AB.

AB = [[  7. ,  1. ,   1. ,   6.],
      [ -2. ,  5. ,   9. ,   1.],
      [ -3. ,  7. ,  -1. , -26.]]

En éste paso se asegura que en la casilla AB[0,0] se encuentre el mayor valor de la columna 1.

El intercambio de filas requiere tres pasos: realizar una copia de la fila que se esta analizando en un vector temporal. Luego se copian los valores en la fila AB[dondemax+i] en la fila AB[i]. Finalmente actualizando la fila AB[dondemax+i] con los valores de temporal.

    if (dondemax !=0): # NO en diagonal
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

Se repite el proceso anterior para la siguiente columna (j=1), pero siempre formada desde la diagonal (i=j).

La columna en éste paso se reduce el tamaño en una posición. Se consideran solo las magnitudes para la selección de dondemax.

AB = [[  7. ,  1. ,   1. ,   6.],
      [ -2. ,  5. ,   9. ,   1.],
      [ -3. ,  7. ,  -1. , -26.]]

Los valores de la columna se reducen a:

columna = [|5|,  fila=0
           |7|]  fila=1
dondemax = 1

En el ejercicio se encuentra que: la magnitud de mayor valor está en la fila=1, no en la diagonal. Se realiza el intercambio entre la fila 'i' y la fila (dondemax+i) de AB.

AB = ([[  7. ,  1. ,   1. ,   6.],
       [ -3. ,  7. ,  -1. , -26.],
       [ -2. ,  5. ,   9. ,   1.]]

Se repite el proceso para la tercera columna desde la diagonal, que resulta tener solo una casilla (columna =[9]) y no ser requiere continuar.

En la gráfica se destaca el resultado del pivoteo por filas, acercando en lo posible a una matriz diagonal dominante.

Pivoteo parcial por filas AB:
[[  7. , 1. ,  1. ,   6.],
 [ -3. , 7. , -1. , -26.],
 [ -2. , 5. ,  9. ,   1.]]
A:
[[ 7. , 1. , 1.],
[ -3. , 7. , -1.],
[ -2. , 5. , 9.]]
B:
[6. , -26. , 1.]

Tarea: Revisar cuando en una columna se tiene dos valores mayores de igual valor. Ampliar el algoritmo. Considere como ejemplo del caso planteado: 1Eva_IT2011_T2_MN Alimentos para animales

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo 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).  Con la copia de la matriz, 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
import numpy as np

# INGRESO
A = [[-2, 5, 9],
     [ 7, 1, 1],
     [-3, 7,-1]]
B = [1,6,-26]

# PROCEDIMIENTO
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)

# Matriz aumentada AB
B_columna = np.transpose([B])
AB  = np.concatenate((A,B_columna),axis=1)
AB0 = np.copy(AB) # copia de 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)
    
    if (dondemax !=0): # NO en diagonal
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal
# Actualiza A y B pivoteado
A = AB[:,:n]
B = AB[:,n]

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas AB:')
print(AB)
print('A:')
print(A)
print('B:')
print(B)

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo 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.

# Pivoteo parcial por filas, funcion
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivoteo 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,
                      'con', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print('AB:')
            print(AB)
    return(AB)

# INGRESO
A = [[-2, 5, 9],
     [ 7, 1, 1],
     [-3, 7,-1]]
B = [1,6,-26]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
n = len(AB)

# Actualiza A y B pivoteado
A = AB[:,:n]
B = AB[:,n]

# SALIDA
print('Pivoteo parcial por filas AB:')
print(AB)
print('A:')
print(A)
print('B:')
print(B)

El resultado del ejercicio es:

Matriz aumentada
[[ -2.   5.   9.   1.]
 [  7.   1.   1.   6.]
 [ -3.   7.  -1. -26.]]
Pivoteo parcial:
  1 intercambiar filas:  0 con 1
  2 intercambiar filas:  1 con 2
[[  7.   1.   1.   6.]
 [ -3.   7.  -1. -26.]
 [ -2.   5.   9.   1.]]
Pivoteo parcial por filas AB:
[[  7.   1.   1.   6.]
 [ -3.   7.  -1. -26.]
 [ -2.   5.   9.   1.]]
A:
[[ 7.  1.  1.]
 [-3.  7. -1.]
 [-2.  5.  9.]]
B:
[  6. -26.   1.]

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]

---

6. Gráfica de Matriz A en 3D

Instrucciones para la gráfica en 3D de una matriz A. Las alturas de las barras son relativas a los valores de la casilla A[i,j]

# Matriz A como Barras3D
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
##A = [[-2, 5, 9],
##     [ 7, 1, 1],
##     [-3, 7,-1]]
##titulo = 'matriz A'

##A = [[ 7, 1, 1],
##     [-2, 5, 9],
##     [-3, 7,-1]]
##titulo = 'A pivotea fila0 y fila1'

##A = [[ 7, 1, 1],
##     [-3, 7,-1],
##     [-2, 5, 9]]
##titulo = 'A pivotea fila1 y fila2'

A = [[ 7, 1, 1],
     [-3, 7,-1],
     [-2, 5, 9]]
titulo='A: Pivoteo por Filas'

# PROCEDIMIENTO
A = np.array(A, dtype=float)
n,m = np.shape(A)

# para grafica 3D
xi = np.arange(0,n,1)
yj = np.arange(0,m,1)
Xi,Yj = np.meshgrid(xi,yj)

A_min = np.min([np.min(A),-3])
A_max = np.max([np.max(A),9])

# SALIDA
# Grafica barras 3D
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111,projection='3d')

AT = np.transpose(A) # filas en x, columnas en y
h = 0.5 # Barras ancho y profundidad
color_A = AT.ravel()/A_max
color_ij = plt.cm.rainbow(color_A)

graf3D.bar3d(Xi.ravel(),Yj.ravel(),
             np.zeros(n*m),h,h,
             AT.ravel(),color=color_ij)

# escala eje z
graf3D.set_zlim(A_min,A_max)
# etiqueta
graf3D.set_xlabel('filas i')
graf3D.set_ylabel('columnas j')
graf3D.set_zlabel('A[i,j]')
graf3D.set_title(titulo)
# etiquetas ejes
graf3D.xaxis.set_ticks(xi + h/2)
graf3D.xaxis.set_ticklabels(xi)
graf3D.yaxis.set_ticks(yj + h/2)
graf3D.yaxis.set_ticklabels(yj)
# Barra de color
color_escala = plt.Normalize(A_min,A_max)
color_barra = plt.cm.ScalarMappable(norm=color_escala,
                cmap=plt.colormaps["rainbow"])
fig3D.colorbar(color_barra,ax=graf3D,label="A[i,j]")

graf3D.view_init(35,35)
plt.show()

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]

---

7. Ejercicio 2

Referencia: Rodríguez 4.0 p105,  1Eva_IT2010_T3_MN Precio artículos

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	X0	X1	X2	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 x₀, x₁, x₂ para escribir el sistema de ecuaciones que muestran las relaciones de cantidad, precio y valor pagado:

\begin{cases} 4x_0+2x_1+5x_2 = 60.70 \\ 2x_0+5x_1+8x_2 = 92.90 \\ 5x_0+4x_1+3x_2 = 56.30 \end{cases}

Desarrollo

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_0 \\ x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 60.70 \\ 92.90 \\ 56.30 \end{pmatrix}

La matriz aumentada usando el algoritmo tiene como resultado:

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

En varios algoritmos se usa la matriz aumentada AB, para mantener el orden entre filas al aplicar cambios de fila.

Para el pivoteo por fila de la matriz aumentada AB, tiene como primer paso revisar la primera columna desde la diagonal en adelante.

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

Los valores de la columna se reducen a:

columna = [|4|,   fila=0
           |2|,   fila=1
           |5|]   fila=2
dondemax = 2

En el ejercicio se encuentra que: la magnitud de mayor valor está en la última fila, no en la diagonal. Se realiza el intercambio entre la fila 2 y la fila 0 de AB.

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

Se repite el proceso anterior para la siguiente columna, pero siempre formada desde la diagonal.

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

Los valores de la columna se reducen a:

columna = [|5|,  fila=0
           |2|]  fila=1
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 a continuación:

Pivoteo parcial por filas AB:
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
A:
[[5. 4. 3.]
 [2. 5. 8.]
 [4. 2. 5.]]
B:
[56.3 92.9 60.7]

[ Ejercicio ] [ Matriz Aumentada ] [ Pivoteo filas ] [ algoritmo ] [ función ]

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]

Referencia: Chapra 9.1 p247
Como punto de partida para la unidad de sistema de ecuaciones, se usa un ejemplo para ilustrar la solución del sistema en gráficos 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 representar lo indicado, se realiza el desarrollo analítico del ejercicio junto al algoritmo y las gráficas en 3D con Python y Matplotlib.

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]
..

---

1. Ejercicio

Referencia: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

Considere el siguiente sistema de ecuaciones AX=B dado por:

\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 ] [ Python ]
..

---

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

Para diferenciar las ecuaciones, se añade el índice por filas i a cada una:

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

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

-5 ≤ x ≤ 5
-7 ≤ y ≤ 7

Al combinar los planos Z₀ y Z₁, 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.

Recta intersección de planos Z0 y Z1

Usando solo las dos primeras ecuaciones y considerando para el eje 'y' una constante, por ejemplo el lado izquierdo del intervalo y_a.

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

Se reordenan las ecuaciones y reemplazan valores:

\begin{cases} -2x+9z=1-5y_a\\7x+z=6-y_a\end{cases} \begin{cases} -2x+9z=1-5(-7)\\7x+z=6-(-7)\end{cases}

Al resolver las ecuaciones se encuentra un punto de la recta de intersección, del lado izquierdo del intervalo para el eje y. Se aplica el mismo proceso para el lado derecho del intervalo del eje 'y' y se tienen las coordenadas:

array([[ 1.24615385, -7.        ,  4.27692308],
       [ 0.38461538,  7.        , -3.69230769]])

Con las que se puede trazar la línea de intersección entre los planos Z0 y Z1

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]
..

---

3. Algoritmo con Python

Para observar el sistema de ecuaciones en un gráfico de tres dimensiones, se utilizan gráficas 3D  y observar los resultados con varios puntos de vista al rotar el gráfico con el cursor.

Para el algoritmo, las ecuaciones del sistema se escriben en forma matricial Ax=B.

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

A es una matriz de coeficientes cuadrada, el número de filas es igual al de columnas. B es un vector de constantes, con el mismo número de elementos que filas de A.

# INGRESO Ax=B
A = [[-2, 5, 9],
     [ 7, 1, 1],
     [-3, 7,-1]]

B = [1,6,-26]

Las ecuaciones de cada plano escriben usando los coeficientes de A y B, según lo realizado en el desarrollo analítico al despejar la variable z. Cada ecuación es una expresión de dos variables (x,y), para simplificar se usa el formato lambda x,y.

# Ecuaciones de planos
f0 = lambda x,y: (B[0]-A[0,0]*x-A[0,1]*y)/A[0,2]
f1 = lambda x,y: (B[1]-A[1,0]*x-A[1,1]*y)/A[1,2]
f2 = lambda x,y: (B[2]-A[2,0]*x-A[2,1]*y)/A[2,2]

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]

..

3.1 Intervalos, muestras para cada eje en los vectores xi, yj

Para la gráfica 3D se requiere usar intervalos para las variables independientes x,y de forma semejante a las gráficas en 2D. En cada intervalo se realizaran muestras en varios puntos del intervalo. La combinación de muestras en cada eje, generan puntos en el plano X,Y.

Se generan las muestras para cada eje en los vectores xi, yj

xa = -5    # Intervalo [xa,xb] para eje x
xb = 5
ya = 7     # Intervalo [ya,yb] para eje y
yb = -7
muestras = 11

xi = np.linspace(xa,xb, muestras)
yj = np.linspace(ya,yb, muestras)

el resultado de las muestras en cada eje para el ejercicio es:

>>> xi
array([-5., -4., -3., -2., -1.,  0.,  1.,  2.,  3.,  4.,  5.])
>>> yj
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 realizan en las matrices X_i,Y_j, (en mayúsculas).

Xi, Yj = np.meshgrid(xi,yj)

Cada matriz representa una malla de puntos en el plano, usando luego cada punto para evaluar el valor de Z. El resultado son las matrices con los puntos de evaluación en la malla. Por ejemplo el primer punto [x₀,y₀] = [-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.]])
>>> Yj
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. ]])
>>> 

La matriz de resultados Z_i, se genera al evaluar una ecuación en cada punto (i, j). Como hay una ecuación por cada fila, se usa el índice para diferenciar cada resultado Z_i.

Z0 = f0(Xi,Yj)
Z1 = f1(Xi,Yj)
Z2 = f2(Xi,Yj)

El gráfico del plano Z₀, se realiza con las matrices X_i,Y_j, con la instrucción plot_wireframe(Xi,Yj,Z0) de las librerías 3D. Luego se repite el procedimiento para Z₁ y Z₂.

graf3D.plot_wireframe(Xi,Yj,Z0,
                      color ='blue',
                      label='Z0')

Si es necesario, revise la sección de Gráficas 3D para wireframe.

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

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

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]
..

---

4. Instrucciones con Python

# Sistema de ecuaciones 3x3
# Interseccion de Planos
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO Ax=B
A = [[-2, 5, 9],
     [ 7, 1, 1],
     [-3, 7,-1]]
B = [1,6,-26]

xa = -5     # Intervalo Xi ,[ax,bx]
xb = 5
ya = -7     # Intervalo Yj ,[ay,by]
yb = 7
muestras = 11 # en cada intervalo

# PROCEDIMIENTO
# matriz y vector como Numpy
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)

# solucion al sistema
punto = np.linalg.solve(A,B)

# Interseccion entre ecuacion Z0 y Z1
A01  = np.copy(A[0:2,[0,2]])
# usando ay, extremo izquierdo eje y
B01  = B[0:2] - ya*A[0:2,1]
x01 = np.linalg.solve(A01,B01)
recta01 = [x01[0],ya,x01[1]]  # tabla coordenadas
# usando by, extremo derecho eje y
B01 = B[0:2] - yb*A[0:2,1]
x01 = np.linalg.solve(A01,B01)
recta01 = np.concatenate(([recta01 ],
                  [[x01[0],yb,x01[1]]]),
                  axis=0)

# ecuaciones de planos
f0 = lambda x,y: (B[0]-A[0,0]*x-A[0,1]*y)/A[0,2]
f1 = lambda x,y: (B[1]-A[1,0]*x-A[1,1]*y)/A[1,2]
f2 = lambda x,y: (B[2]-A[2,0]*x-A[2,1]*y)/A[2,2]

# muestras
xi = np.linspace(xa,xb, muestras)
yj = np.linspace(ya,yb, muestras)
Xi, Yj = np.meshgrid(xi,yj) # en plano XY

# evalua planos Zi
Z0 = f0(Xi,Yj)
Z1 = f1(Xi,Yj)
Z2 = f2(Xi,Yj)

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

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

# Planos
graf3D.plot_wireframe(Xi,Yj,Z0,
                       color ='blue',
                       label='Z0')
graf3D.plot_wireframe(Xi,Yj,Z1,
                       color ='orange',
                       label='Z1')
graf3D.plot_wireframe(Xi,Yj,Z2,
                       color ='green',
                       label='Z2')

# recta entre planos Z0 y Z1
graf3D.plot(recta01[:,0],recta01[:,1], recta01[:,2],
            color='purple',
            label='recta01',linewidth=4)

# Punto solucion Z0,Z1,Z2
graf3D.plot(punto[0],punto[1],punto[2],
            'o',color='red',
            label='Punto',linewidth=6)

# etiquetas y entorno gráfico
graf3D.set_title('Sistema de ecuaciones 3x3')
graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_xlim(xa, xb)
graf3D.set_ylim(ya, yb)
graf3D.legend()
graf3D.view_init(45, 45) # elevacion, azimut

### rotacion de ejes
##for angulo in range(45, 360+45, 5 ):
##    graf3D.view_init(45, angulo)
##    plt.draw()
##    plt.pause(.1)

plt.show()

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]

---

Referencia: Dear linear algebra students, This is what matrices (and matrix manipulation) really look like. Zach Star. 5 mar 2020.

Ejemplo: [ ejercicio ] [ analítico ] [ algoritmo ] [ Python ]