Unidad 3 Sistemas de ecuaciones – Métodos numéricos

2017-06-132025-08-16

3.9.2 Método de Gauss-Jordan para matriz Inversa con Python

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

1. Ejercicio

Referencia: Chapra 10.2 p292, Burden 6.3 p292, Rodríguez 4.2.5 Ejemplo 1 p118
Obtener la inversa de una matriz usando el método de Gauss-Jordan, a partir de la matriz:

\begin{pmatrix} 4 & 2 & 5 \\ 2 & 5 & 8 \\ 5 & 4 & 3 \end{pmatrix}

A = [[4,2,5],
     [2,5,8],
     [5,4,3]]

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

2. Desarrollo analítico

Para el procedimiento, se crea la matriz aumentada de A con la identidad I.

AI = A|I

\begin{pmatrix} 4 & 2 & 5 & 1 & 0 & 0\\ 2 & 5 & 8 & 0 & 1 & 0 \\ 5 & 4 & 3 & 0 & 0 & 1 \end{pmatrix}

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

Con la matriz aumentada AI se repiten los procedimientos aplicados en el método de Gauss-Jordan:

pivoteo parcial por filas
eliminación hacia adelante
eliminación hacia atrás

De la matriz aumentada resultante, se obtiene la inversa A^-1 en la mitad derecha de AI, lugar que originalmente correspondía a la identidad.

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

\begin{pmatrix} 0.2 & -0.16470588 & 0.10588235 \\ -0.4 & 0.15294118 & 0.25882353 \\ 0.2 & 0.07058824 & -0.18823529 \end{pmatrix}

Verifica resultado

El resultado se verifica realizando la operación producto punto entre A y la inversa, que debe resultar la matriz identidad.

A.A^-1 = I

El resultado de la operación es una matriz identidad. Observe que los valores del orden de 10^-15 o menores se consideran como casi cero o cero.

 A.inversa = identidad
[[ 1.00000000e+00 -1.38777878e-17 -1.38777878e-16]
 [ 2.22044605e-16  1.00000000e+00 -2.22044605e-16]
 [ 5.55111512e-17 -9.71445147e-17  1.00000000e+00]]

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

3. Algoritmo en Python

El algoritmo que describe el proceso en python:

# Matriz Inversa con Gauss-Jordan
# AI es la matriz aumentada A con Identidad
import numpy as np

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

# PROCEDIMIENTO
# B = matriz_identidad
n,m = np.shape(A)
identidad = np.identity(n)

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(identidad,dtype=float)

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

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 gauss_eliminaAtras(AB, vertabla=False,
                       inversa=False,
                       casicero = 1e-14):
    ''' Gauss-Jordan elimina hacia atras
    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 Atras:')
        
    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)

AB = pivoteafila(A,identidad,vertabla=True)
AB = gauss_eliminaAdelante(AB,vertabla=True)
A_inversa = gauss_eliminaAtras(AB,inversa=True,
                               vertabla=True)
A_verifica = np.dot(A,A_inversa)
# redondeo a cero
for i in range(0,n,1):
    for j in range(0,m,1): 
        if np.abs(A_verifica[i,j])<=casicero:
            A_verifica[i,j]=0

# SALIDA
print('A_inversa: ')
print(A_inversa)
print('verifica A.A_inversa = Identidad:')
print(A_verifica)

el resultado buscado es:

Matriz aumentada
[[4. 2. 5. 1. 0. 0.]
 [2. 5. 8. 0. 1. 0.]
 [5. 4. 3. 0. 0. 1.]]
Pivoteo parcial:
  1 intercambiar filas:  0 con 2
[[5. 4. 3. 0. 0. 1.]
 [2. 5. 8. 0. 1. 0.]
 [4. 2. 5. 1. 0. 0.]]
Elimina hacia adelante:
 fila i: 0  pivote: 5.0
  fila k: 1  factor: 0.4
  fila k: 2  factor: 0.8
[[ 5.   4.   3.   0.   0.   1. ]
 [ 0.   3.4  6.8  0.   1.  -0.4]
 [ 0.  -1.2  2.6  1.   0.  -0.8]]
 fila i: 1  pivote: 3.4
  fila k: 2  factor: -0.3529411764705883
[[ 5.      4.      3.      0.      0.      1.    ]
 [ 0.      3.4     6.8     0.      1.     -0.4   ]
 [ 0.      0.      5.      1.      0.3529 -0.9412]]
 fila i: 2  pivote: 5.0
[[ 5.      4.      3.      0.      0.      1.    ]
 [ 0.      3.4     6.8     0.      1.     -0.4   ]
 [ 0.      0.      5.      1.      0.3529 -0.9412]]
Elimina hacia Atras:
 fila i: 2  pivote: 5.0
  fila k: 1  factor: 1.3599999999999999
  fila k: 0  factor: 0.6
[[ 5.      4.      0.     -0.6    -0.2118  1.5647]
 [ 0.      3.4     0.     -1.36    0.52    0.88  ]
 [ 0.      0.      1.      0.2     0.0706 -0.1882]]
 fila i: 1  pivote: 3.4
  fila k: 0  factor: 1.1764705882352942
[[ 5.      0.      0.      1.     -0.8235  0.5294]
 [ 0.      1.      0.     -0.4     0.1529  0.2588]
 [ 0.      0.      1.      0.2     0.0706 -0.1882]]
 fila i: 0  pivote: 5.0
[[ 1.      0.      0.      0.2    -0.1647  0.1059]
 [ 0.      1.      0.     -0.4     0.1529  0.2588]
 [ 0.      0.      1.      0.2     0.0706 -0.1882]]
A_inversa: 
[[ 0.2    -0.1647  0.1059]
 [-0.4     0.1529  0.2588]
 [ 0.2     0.0706 -0.1882]]
verifica A.A_inversa = Identidad:
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Observe que el algoritmo se pude reducir si usan los procesos de Gauss-Jordan como funciones.

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

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..

4. Función en Numpy

La función en la librería Numpy es np.linalg.inv():

>>> np.linalg.inv(A)
array([[ 0.2       , -0.16470588,  0.10588235],
       [-0.4       ,  0.15294118,  0.25882353],
       [ 0.2       ,  0.07058824, -0.18823529]])
>>>

matriz inversa: [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]

2017-06-122025-08-16

3.9.1 Método de Gauss – determinante de matriz con Python

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

1. Determinante de matriz

Referencia: Rodríguez 4.3.9 p129, Burden 6.4 p296, Chapra 9.1.2 p250.

El determinante de una matriz cuadrada triangular superior también puede calcularse como el producto de los coeficientes de la diagonal principal, considerando el número de cambios de fila del pivoteo k.

det(A) = (-1)^k \prod_{i=1}^n a_{i,i}

Si observamos que en las secciones anteriores se tiene desarrollado los algoritmos para obtener la matriz triangular superior en el método de Gauss, se usan como punto de partida para obtener los resultados del cálculo del determinante.

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

2. Ejercicio

Referencia: 1Eva_IIT2010_T2 Sistema ecuaciones, X0 = unos, Chapra Ejemplo 9.12 p277

Calcular el determinante de la matriz A del sistema A.X=B dado por

A= \begin{pmatrix} 0.4 & 1.1 & 3.1 \\ 4 & 0.15 & 0.25 \\2 & 5.6 & 3.1 \end{pmatrix}

B= [7.5, 4.45, 0.1]

El algoritmo para ejercicio se convierte en una extensión de los algoritmos anteriores.

A = [[0.4, 1.1 ,  3.1],
     [4.0, 0.15, 0.25],
     [2.0, 5.6 , 3.1]]
B = [7.5, 4.45, 0.1]

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

3. Algoritmo en Python

El algoritmo parte de lo realizado en método de Gauss, indicando que la matriz a procesar es solamente A. Se mantienen los procedimientos de «pivoteo parcial por filas» y » eliminación hacia adelante»

Para contar el número de cambios de filas, se añade un contador de pivoteado dentro del condicional para intercambio de filas.

Para el resultado del operador multiplicación, se usan todas las casillas de la diagonal al acumular las multiplicaciones.

# Determinante de una matriz A
# convirtiendo a diagonal superior 

import numpy as np

# INGRESO
A = [[0.4, 1.1 ,  3.1],
     [4.0, 0.15, 0.25],
     [2.0, 5.6 , 3.1]]
B = [7.5, 4.45, 0.1]

# 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
pivoteado = 0 # contador para cambio fila

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
        
        pivoteado = pivoteado + 1 # cuenta cambio fila
        
# Actualiza A y B pivoteado
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

# calcula determinante
multiplica = 1
for i in range(0,n,1):
    multiplica = multiplica*AB[i,i]
determinante = ((-1)**pivoteado)*multiplica

# SALIDA
print('Pivoteo parcial por filas')
print(AB1)
print('Cambios de fila, pivoteado: ',pivoteado)
print('eliminación hacia adelante')
print(AB)
print('determinante: ')
print(determinante)

Se aplica la operación de la fórmula planteada para el método, y se presenta el resultado:

Pivoteo parcial por filas
[[4.   0.15 0.25 4.45]
 [2.   5.6  3.1  0.1 ]
 [0.4  1.1  3.1  7.5 ]]
Cambios de fila, pivoteado:  2
eliminación hacia adelante
[[ 4.          0.15        0.25        4.45      ]
 [ 0.          5.525       2.975      -2.125     ]
 [ 0.          0.          2.49076923  7.47230769]]
determinante: 
55.04599999999999

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..

4. Algoritmo como función en Numpy

El resultado se puede verificar usando la función de Numpy:

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

[ Determinante ] [ Ejercicio ] [ Algoritmo ] [ Función ]

2017-06-112025-08-16

3.8 Matrices triangulares A=L.U con Python

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]

1. Matrices triangulares A=L.U

Referencia: Chapra 10.1 p284. Burden 6.5 p298

Una matriz A puede separarse en dos matrices triangulares que entre ellas tienen la propiedad: A = L.U:

L de tipo triangular inferior
U de tipo triangular superior

La matriz U se obtiene después de aplicar el proceso de eliminación hacia adelante del método de Gauss.

La matriz L contiene los factores usados en el proceso de eliminación hacia adelante del método de Gauss, escritos sobre una matriz identidad en las posiciones donde se calcularon.

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..

2. Ejercicio

Referencia: Chapra Ejemplo 10.1 p285

Presente las matrices LU de la matriz A siguiente:

A= \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix}

A = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

El resultado del «pivoteo por filas» y «eliminación hacia adelante» se tiene:

Pivoteo parcial por filas
[[  3.    -0.1   -0.2    7.85]
 [  0.1    7.    -0.3  -19.3 ]
 [  0.3   -0.2   10.    71.4 ]]

de donde la última columna es el nuevo vector B

eliminación hacia adelante
Matriz U: 
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]

y guardando los factores del procedimiento de «eliminación hacia adelante » en una matriz L que empieza con la matriz identidad se obtiene:

matriz L: 
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..

3. Algoritmo LU en Python

Realizado a partir del algoritmo de la sección método de Gauss y modificando las partes necesarias para el algoritmo.

Para éste algoritmo, se procede desde el bloque de pivoteo por filas«, continuando con el algoritmo de «eliminación hacia adelante» del método de Gauss. Procedimientos que dan como resultado la matriz U.

La matriz L requiere iniciar con una matriz identidad, y el procedimiento requiere que al ejecutar «eliminación hacia adelante» se almacene cada factor con el que se multiplica la fila para hacer cero. El factor se lo almacena en la matriz L, en la posición de dónde se determinó el factor.

El resultado obtenido con el algoritmo es:

Matriz aumentada
[[  3.    -0.1   -0.2    7.85]
 [  0.1    7.    -0.3  -19.3 ]
 [  0.3   -0.2   10.    71.4 ]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
matriz L: 
[[ 1.      0.      0.    ]
 [ 0.0333  1.      0.    ]
 [ 0.1    -0.0271  1.    ]]
Matriz U: 
[[ 3.     -0.1    -0.2   ]
 [ 0.      7.0033 -0.2933]
 [ 0.      0.     10.012 ]]
>>>

Donde se puede verificar que L.U = A usando la instrucción np.dot(L.U):

>>> np.dot(L,U)
array([[ 3. , -0.1, -0.2],
       [ 0.1,  7. , -0.3],
       [ 0.3, -0.2, 10. ]])
>>>

Instrucciones en Python

# Matrices L y U
# Modificando el método de Gauss
import numpy as np

# PROGRAMA ------------
# INGRESO
A = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

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

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

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

tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

AB0 = np.copy(AB)

L = np.identity(n,dtype=float) # Inicializa L

# eliminacion 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
        
        L[k,i] = factor # llena L

U = np.copy(AB[:,:n])

# SALIDA
print('matriz L: ')
print(L)
print('Matriz U: ')
print(U)

Si se requiere una respuesta unificada en una variable, se puede convertir en una arreglo de matrices [L,U].
Las matrices L y U se pueden leer como L=LU[0] y U=LU[1]

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

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

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]

4. Algoritmo como Función de Python

El resultado a obtener es:

Matriz aumentada
[[  3.    -0.1   -0.2    7.85]
 [  0.1    7.    -0.3  -19.3 ]
 [  0.3   -0.2   10.    71.4 ]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Elimina hacia adelante:
 fila 0 pivote:  3.0
   factor:  0.03333333333333333  para fila:  1
   factor:  0.09999999999999999  para fila:  2
 fila 1 pivote:  7.003333333333333
   factor:  -0.027129938124702525  para fila:  2
 fila 2 pivote:  10.012041884816753
[[  3.      -0.1     -0.2      7.85  ]
 [  0.       7.0033  -0.2933 -19.5617]
 [  0.       0.      10.012   70.0843]]
matriz L: 
[[ 1.      0.      0.    ]
 [ 0.0333  1.      0.    ]
 [ 0.1    -0.0271  1.    ]]
Matriz U: 
[[  3.      -0.1     -0.2      7.85  ]
 [  0.       7.0033  -0.2933 -19.5617]
 [  0.       0.      10.012   70.0843]]
>>>

Instrucciones en Python

# Matrices L y U
# Modificando el método de Gauss
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)
    return(AB)

def gauss_eliminaAdelante(AB, vertabla=False,lu=False,casicero = 1e-15):
    ''' Gauss elimina hacia adelante, a partir de,
    matriz aumentada y pivoteada.
    Para respuesta en forma A=L.U usar lu=True entrega[AB,L,U]
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    L = np.identity(n,dtype=float) # Inicializa L
    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,:]
                for j in range(0,m,1): # casicero revisa
                    if abs(AB[k,j])<casicero:
                        AB[k,j]=0

                L[k,i] = factor # llena L
                
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    respuesta = AB
    if vertabla==True:
        print(AB)
    if lu==True:
        U = AB[:,:n-1]
        respuesta = [AB,L,U]
    return(respuesta)

# PROGRAMA ------------------------
# INGRESO
A = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

# 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, lu=True)
L = AB[1]
U = AB[0]

# SALIDA
print('matriz L: ')
print(L)
print('Matriz U: ')
print(U)

Función en Scipy

Luego del resultado o definida la matriz a, la instrucción en la librería Scipy es:

>>> import scipy as sp
>>> [L,U] =sp.linalg.lu(A,permute_l=True)
>>> L
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.03333333,  1.        ,  0.        ],
       [ 0.1       , -0.02712994,  1.        ]])
>>> U
array([[ 3.        , -0.1       , -0.2       ],
       [ 0.        ,  7.00333333, -0.29333333],
       [ 0.        ,  0.        , 10.01204188]])
>>>

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]
..

5. Solución con LY=B , UX=Y

Referencia: Chapra 10.2 p287, pdf312

Se plantea resolver el sistema de ecuaciones usando matrices triangulares

A = \begin{pmatrix} 3 & -0.1 & -0.2 \\ 0.1 & 7 & -0.3 \\0.3 & -0.2 & 10 \end{pmatrix}

B = [7.85,-19.3,71.4]

Para encontrar la solución al sistema de ecuaciones, se plantea resolver:
– sustitución hacia adelante de LY=B
– sustitución hacia atrás para UX=Y

Al algoritmo de la sección anterior se añade los procedimientos correspondientes con los que se obtiene la solución:

Matriz aumentada
[[  3.    -0.1   -0.2    7.85]
 [  0.1    7.    -0.3  -19.3 ]
 [  0.3   -0.2   10.    71.4 ]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
matriz L: 
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
Matriz U: 
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
B1 :
[[  7.85]
 [-19.3 ]
 [ 71.4 ]]
Y Sustitución hacia adelante
[[  7.85      ]
 [-19.56166667]
 [ 70.08429319]]
X Sustitución hacia atras
[ 3.  -2.5  7. ]
>>>

Instrucciones en Python

# Solución usando Matrices L y U
# continuación de ejercicio anterior
import numpy as np

# PROGRAMA ------------
# INGRESO
A = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

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

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

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

tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

AB0 = np.copy(AB)
B1 = np.copy(AB[:,m-1])

# Matrices L y U
L = np.identity(n,dtype=float) # Inicializa L
# eliminacion 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
        
        L[k,i] = factor # llena L

U = np.copy(AB[:,:n])

# Resolver LY = B
B1  = np.transpose([B1])
AB =np.concatenate((L,B1),axis=1)

# sustitución hacia adelante
Y = np.zeros(n,dtype=float)
Y[0] = AB[0,n]
for i in range(1,n,1):
    suma = 0
    for j in range(0,i,1):
        suma = suma + AB[i,j]*Y[j]
    b = AB[i,n]
    Y[i] = (b-suma)/AB[i,i]

Y = np.transpose([Y])

# Resolver UX = Y
AB =np.concatenate((U,Y),axis=1)

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

# SALIDA
print('matriz L: ')
print(L)
print('Matriz U: ')
print(U)
print('B1 :')
print(B1)
print('Y Sustitución hacia adelante')
print(Y)
print('X Sustitución hacia atras')
print(X)

[ Matriz A=L.U ] [ Ejercicio ] [ Algoritmo LU ] [ Función LU ] ||
[ Solución LY=B , UX=Y ]

2017-06-102025-08-16

3.7 Método de Gauss-Seidel con Python

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

Referencia: Chapra 11.2 p310, Burden 7.3 p337, Rodríguez 5.2 p162

El método de Gauss-Seidel, reescribe las expresiones del sistema de ecuaciones, despejando la incógnita de la diagonal en cada ecuación. Usa el vector inicial X₀, actualizando el vector X por cada nuevo valor calculado del vector. La diferencia principal con el método de Jacobi es la actualización de X en cada cálculo, por lo que las iteraciones llegan más rápido al punto cuando el método converge.

Las iteraciones pueden ser o no, convergentes.

La expresión para encontrar nuevos valores usando la matriz de coeficientes A, el vector de constantes B y como incógnitas X es la misma que Jacobi:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n a_{i,j}x_j\bigg) \frac{1}{a_{ii}}

La convergencia del método iterativo se puede revisar usando el número de condición de la matriz de coeficientes A. Un resultado cercano a 1, implica que el sistema será convergente.

cond(A) = ||A|| ||A^-1||

np.linalg.cond(A)

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

2. 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}

Luego del pivoteo parcial por filas, el sistema de ecuaciones a usar para el método de Gauss-Seidel es:

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

Desarrolle el ejercicio usando el método de Gauss-Seidel, tomando como vector inicial X₀=[0,0,0] y tolerancia de 0.0001

Compare los resultados con el método de Jacobi.

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

3. Desarrollo Analítico

Para el método de Gauss-Seidel, se usan las ecuaciones reordenas con el pivoteo parcial por filas, para que en lo posible sea diagonal dominante que mejora la convergencia. El resultado del algoritmo a partir de las ecuaciones presentadas al inicio es:

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

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

se reescriben despejando la incógnita de la diagonal en cada ecuación.

x=\frac{6-y-z}{7}

y= \frac{-26+3x+z}{7}

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

El vector inicial X₀=[0,0,0] representa un punto de partida, valores estimados o semilla para buscar la solución a las variables X₀=[x₀,y₀,z₀].

Como las ecuaciones son el comportamiento del ‘sistema‘, se usan con X₀ para obtener un nuevo y mejor valor estimado para la solución.

itera = 0

X_0=[x_0,y_0,z_0]

X_1 = X_0=[0,0,0]

x_1=\frac{6-(0)-(0)}{7} =\frac{6}{7}=0.8571

Se actualiza el vector X₁,

X_1 = [0.8571,0,0]

y_1= \frac{-26+3(0.8571)+(0)}{7} =-3.3469

Se actualiza el vector X₁,

X_1 = [0.8571,-3.3469,0]

z_1=\frac{1+2(0.8571)-5(-3.3469)}{9} =\frac{1}{9}=2.161

X_1=[0.8571, -3.3469, 2.161]

Para revisar si se está más cerca o no de la solución se calculan los desplazamientos o diferencias entre los valores anteriores y los nuevos:

\text{diferencias}=X_1-X_0

=[0.8571-0, -3.3469-0, 2.161-0]

 [0.8571, -3.3469, 2.161]
-[0.    ,  0.    , 0.    ]
__________________________
 [0.8571, -3.3469, 2.161]

\text{diferencias}=[0.8571, -3.3469, 2.161]

Para estimar el valor del error, se presenta un valor simplificado o escalar de las diferencias (norma del error). Por simplicidad se usa el valor de magnitud máxima. Se usa como variable ‘errado‘, para evitar el conflicto de nombre de variable error usada como palabra reservada en Python:

\text{errado}= max \Big|\text{diferencias}\Big|

\text{errado}= max \Big|[0.8571, -3.3469, 2.161]\Big| = 3.3469

Al comparar el valor errado con la tolerancia, se observa que el error aún es mayor que lo tolerado. Por lo que se sigue con la siguiente iteración.

itera = 1

X_1=[x_1,y_1,z_1]

X_2 = X_1=[0.8571, -3.3469, 2.161]

x_2=\frac{6-(-3.3469)-(2.161)}{7} = 1.0266

Se actualiza el vector X₂,

X_2 = [1.0266, -3.3469, 2.161]

y_2= \frac{-26+3(1.0266)+(2.161)}{7} =-2.9656

Se actualiza el vector X₂,

X_2 = [1.0266, -2.9656, 2.161]

z_2=\frac{1+2(1.0266)-5(-2.9656)}{9} =1.9868

X_2 = [1.0266, -2.9656, 1.9868]

\text{diferencias}=|X_2-X_1|

=|[1.0266-0.8571, -2.9656-(-3.3469), 1.9868-2.161]|

 [1.0266, -2.9656, 1.9868]
-[0.8571, -3.3469, 2.161 ] 
__________________________
 [0.1694,  0.3813, 0.1742]

\text{diferencias}= |[0.1694, 0.3813, 0.1742]|

\text{errado}=max |X_2-X_1|= 0.3813

itera = 2

X_3 = X_2=[1.0266, -2.9656, 1.9868]

x_3=\frac{6-(-2.9656)-(1.9868)}{7} = 0.997

X_3 =[0.997, -2.9656, 1.9868]

y_3= \frac{-26+3(0.997)+(1.9868)}{7} = -3.0032

X_3 =[0.997, -3.0032, 1.9868]

z_3=\frac{1+2(0.997)-5(-3.0032)}{9} = 2.0011

X_3 =[0.997, -3.0032, 2.0011]

\text{diferencias}=|X_3-X_2|

=|[0.997-1.0266, -3.0032-(-2.9656), 2.0011-1.9868]|

 [ 0.997 , -3.0032, 2.0011]
-[ 1.0266, -2.9656, 1.9868] 
__________________________
 [-0.0296,  0.0376, 0.0143]

\text{diferencias}=|[-0.0296, 0.0376, 0.0143]|

\text{errado}=max |X_2-X_1|= 0.0376

En cada iteración el valor de error disminuye, por lo que se estima que el método converge.

Al usar el algoritmo, de los resultados se observa que se detiene luego de 6 iteraciones, los errores son menores a la tolerancia:

Iteraciones Gauss-Seidel
itera,[X]
   errado,[diferencia]
0 [0. 0. 0.] 0.0002
1 [ 0.8571 -3.3469  2.161 ]
  3.3469387755102042 [0.8571 3.3469 2.161 ]
2 [ 1.0266 -2.9656  1.9868]
  0.381322597066037 [0.1694 0.3813 0.1742]
3 [ 0.997  -3.0032  2.0011]
  0.03756644188210334 [0.0296 0.0376 0.0143]
4 [ 1.0003 -2.9997  1.9999]
  0.00346691102194141 [0.0033 0.0035 0.0012]
5 [ 1. -3.  2.]
  0.0003256615309844557 [3.2566e-04 3.0918e-04 9.9398e-05]
6 [ 1. -3.  2.]
  2.996898191776065e-05 [2.9969e-05 2.7044e-05 8.3644e-06]
Metodo de Gauss-Seidel
numero de condición: 1.9508402675105447
X:  [ 1. -3.  2.]
errado: 2.996898191776065e-05
iteraciones: 6

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

4. Algoritmo en Python

El algoritmo para Gauss-Seidel es semejante al realizado para Jacobi, se debe incorporar las instrucciones para actualizar los valores de X[i] en cada uso de ecuación o fila.

Para el algoritmo se usa el sistema de ecuaciones en su forma matricial A.X=B.
Se asume que se ha aplicado el pivoteo parcial por filas y matriz aumentada.

i=0 # una ecuacion
suma = B[i] # numerador
for j in range(0,m,1): # un coeficiente
    if (i!=j): # excepto diagonal de A
        suma = suma-A[i,j]*X[j]
                
nuevo = suma/A[i,i]
diferencia[i] = abs(nuevo-X[i])
X[i] = nuevo
errado = np.max(diferencia)

La actualización para el algoritmo es:

# Método de Gauss-Seidel
import numpy as np

# INGRESO
# NOTA: Para AB, se ha usado pivoteo parcial por filas
A = [[ 7,  1,  1],
     [-3,  7, -1],
     [-2,  5,  9]]
B = [6, -26, 1]

X0 = [0,0,0]
tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A) # tamaño A
n = tamano[0]
m = tamano[1]

# valores iniciales
diferencia = np.ones(n, dtype=float)
errado = 2*tolera # np.max(diferencia)

itera = 0 
X = np.copy(X0)
while errado>tolera and itera<=iteramax:
    
    for i in range(0,n,1): # una ecuacion
        suma = B[i] # numerador
        for j in range(0,m,1): # un coeficiente
            if (i!=j): # excepto diagonal de A
                suma = suma-A[i,j]*X[j]
                
        nuevo = suma/A[i,i]
        diferencia[i] = abs(nuevo-X[i])
        X[i] = nuevo
    errado = np.max(diferencia)
    
    print(itera, X)
    print('  ',errado,diferencia)
    itera = itera + 1

if (itera>iteramax): # No converge
    X = np.nan
    print('No converge,iteramax superado')

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('Método de Gauss-Seidel')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', itera)

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

5. Algoritmo en Python como función

Las instrucciones del algoritmo, empaquetadas como función, añaden algunas facilidades para el desarrollo del ejercicio. Antes de usar el método de Gauss-Seidel, se aplica el pivoteo parcial por filas para mejorar de ser posible la convergencia.

En la función gauss_seidel() se añaden las instrucciones para mostrar los resultados parciales en cada iteración. Se añade una tabla para revisión de valores al final del algoritmo, o para realizar la gráfica de errores por iteración.

Si el sistema es de 3×3 se puede añadir una gráfica de los resultados por iteración, mostrando la trayectoria descrita por los resultados parciales en 3D.

# Algoritmo Gauss-Seidel
# solución de matrices
# métodos iterativos
# Referencia: Chapra 11.2, p.310,
#      Rodriguez 5.2 p.162
import numpy as np

def gauss_seidel(A,B,X0,tolera, iteramax=100, vertabla=False, precision=4):
    ''' Método de Gauss Seidel, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    # Matrices como arreglo, numeros reales
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float)
    X0 = np.array(X0, dtype=float)
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]

    # valores iniciales
    diferencia = 2*tolera*np.ones(n, dtype=float)
    errado = 2*tolera # np.max(diferencia)
   
    tabla = [np.copy(X0)] # tabla de iteraciones
    tabla = np.concatenate((tabla,[[np.nan]]),
                            axis=1) # añade errado

    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X]')
        print('   errado,[diferencia]')
        print(0,X0)
        print(' ',np.nan)
        np.set_printoptions(precision)

    itera = 0 # Gauss-Sediel
    X = np.copy(X0)
    while (errado>tolera and itera<iteramax):
        
        for i in range(0,n,1): # una ecuacion
            suma = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    suma = suma-A[i,j]*X[j]
            nuevo = suma/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)

        Xfila= np.concatenate((X,[errado]),axis=0)
        tabla = np.concatenate((tabla,[Xfila]),axis = 0)
        itera = itera + 1
        if vertabla==True:
            print(itera, X)
            print(' ', errado,diferencia)
        
    # No converge
    if (itera>iteramax):
        X = np.nan
        print('No converge,iteramax superado')
    return(X,tabla)

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)

# PROGRAMA Búsqueda de solucion  --------
# INGRESO
A = [[ 7,  1,  1],
     [-3,  7, -1],
     [-2,  5,  9]]
B = [6, -26, 1]

X0 = [0,0,0]
tolera = 0.0001
iteramax = 100
verdecimal = 4

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
n,m = np.shape(AB)

A = AB[:,:n] # separa en A y B
B = AB[:,n]

[X, tabla] = gauss_seidel(A,B,X0, tolera,
                           vertabla=True,
                           precision=verdecimal)
n_itera = len(tabla)-1
errado = tabla[-1,-1]

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('Metodo de Gauss-Seidel')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', n_itera)

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]
..

4.1 Gráficas de iteraciones

Se muestra la gráfica de errores vs iteraciones.

# GRAFICA de iteraciones
errados = tabla[:,n]

import matplotlib.pyplot as plt
# grafica de error por iteracion
fig2D = plt.figure()
graf = fig2D.add_subplot(111)
graf.plot(errados,color='purple')
graf.set_xlabel('itera')
graf.set_ylabel('|error|')
graf.set_title('Método de Gauss-Seidel: error[itera]')
graf.grid()

plt.show()

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]

4.2 Gráfica de iteraciones para sistema de 3×3

Si el sistema es de 3×3, se muestra una gráfica de la trayectoria de resultados parciales, semejante a la presentada para la descripción del concepto. Permite observar si la espiral es convergente o no.

# grafica de iteraciones x,y,z
# solo para 3x3
xp = tabla[:,0]
yp = tabla[:,1]
zp = tabla[:,2]

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

graf3D.plot(xp,yp,zp,color='purple')
graf3D.plot(xp[0],yp[0],zp[0],'o',label='X[0]')
graf3D.plot(xp[n_itera],yp[n_itera],zp[n_itera],
            'o',label='X['+str(n_itera)+']')

graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_title('Método de Gauss-Seidel: iteraciones X')

# graf3D.view_init(45,45)
plt.legend()
plt.show()

[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ] [Gráfica]

2017-06-092025-08-16

3.6 Método de Jacobi con Python

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]

1. El método iterativo de Jacobi

Referencia: Burden 7.3 p334, Rodríguez 5.1 p154, Chapra 11.2.1p312

Los métodos iterativos para sistemas de ecuaciones, son semejantes al método de punto fijo para búsqueda de raíces, requieren un punto inicial para la búsqueda «de la raíz o solución» que satisface el sistema. Se plantea considerar aproximaciones o iteraciones de forma sucesiva desde un punto de partida o inicial X₀, hacia una solución del sistema de ecuaciones.

Las iteraciones pueden ser o no, convergentes.

Para describir el método iterativo de Jacobi, se usa como ejemplo: un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante, planteado en su forma matricial.

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

La expresión para encontrar nuevos valores usando la matriz de coeficientes A, el vector de constantes B y como incógnitas X es:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n a_{i,j}x_j\bigg) \frac{1}{a_{ii}}

El método de Jacobi, reescribe las expresiones despejando la incógnita de la diagonal en cada ecuación. El patrón que se observa en las ecuaciones despejadas, se describe en la fórmula para x_i de la fila i, que se usará en el algoritmo en Python.

x_0 = \frac{b_{0} -a_{0,1}x_1 -a_{0,2} x_2 }{a_{0,0}}

x_1 = \frac{b_{1} -a_{1,0} x_0 -a_{1,2} x_2}{a_{1,1}}

x_2 = \frac{b_{2} -a_{2,0} x_0 - a_{2,1} x_1}{a_{2,2}}

La convergencia del método iterativo se puede revisar usando el número de condición de la matriz de coeficientes A. Un resultado cercano a 1, implica que el sistema será convergente.

cond(A) = ||A|| ||A^-1||

np.linalg.cond(A)

Usando de forma experimental los ejercicios de las evaluaciones de la sección matriciales iterativos, el número de condición no debería superar el valor de 7 para que el sistema sea convergente. Puede revisar que si es mayor, el método tiende a no ser convergente. Aunque los textos del libro no lo indican directamente,

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]

2. Ejercicio

Referencia: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

Considere el siguiente sistema de ecuaciones A.X=B dado por:

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

Luego del pivoteo parcial por filas, el sistema de ecuaciones a usar para el método de Jacobi es:

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

Desarrolle el ejercicio usando el método de Jacobi, tomando como vector inicial X₀=[0,0,0] y tolerancia de 0.0001.

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]

3. Desarrollo Analítico

Para el método de Jacobi, se usan las ecuaciones reordenas con el pivoteo parcial por filas, para que en lo posible sea diagonal dominante que mejora la convergencia. El resultado del algoritmo a partir de las ecuaciones presentadas al inicio es:

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

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

Se despejan la incógnita de la diagonal en cada ecuación.

x=\frac{6-y-z}{7}

y= \frac{-26+3x+z}{7}

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

El vector inicial X₀=[0,0,0] representa un punto de partida, valores estimados o semilla para buscar la solución a las variables x₀,y₀,z₀.

Como las ecuaciones son el comportamiento del ‘sistema‘, se usan con X₀ para obtener un nuevo y mejor valor estimado para la solución.

itera = 0

X_0=[x_0,y_0,z_0]

X_0=[0,0,0]

x_1=\frac{6-(0)-(0)}{7} =\frac{6}{7}=0.8571

y_1= \frac{-26+3(0)+(0)}{7} =\frac{-26}{7}=-3.7142

z_1=\frac{1+2(0)-5(0)}{9} =\frac{1}{9}=0.1111

X_1=[0.8571, -3.7142, 0.1111]

Para revisar si se está más cerca de la solución, se calculan las diferencias, tramos o desplazamientos entre los valores nuevos y anteriores:

\text{diferencias}=X_1-X_0

=[0.8571-0, -3.7142-0, 0.1111-0]

 [0.8571, -3.7142, 0.1111]
-[0.    ,  0.    , 0.    ]
__________________________
 [0.8571, -3.7142, 0.1111]

\text{diferencias}=[0.8571, -3.7142, 0.1111]

Para estimar el valor del error, se presenta un valor simplificado o escalar de las diferencias (norma del error). Por simplicidad se usa el valor de magnitud máxima de los errores. Se usa como variable ‘errado‘, para evitar el conflicto de nombre de variable error usada como palabra reservada en Python:

\text{errado}= max \Big|\text{diferencias}\Big|

\text{errado}= max \Big|[0.8571, -3.7142, 0.1111]\Big| = 3.7142

Al comparar el valor errado con la tolerancia, se observa que el error aún es mayor que lo tolerado. Por lo que se continúa con otra iteración.

itera = 1

X_1=[x_1,y_1,z_1]

X_1=[0.8571, -3.7142, 0.1111]

x_2=\frac{6-(-3.7142)-(0.1111)}{7} = 1.3718

y_2= \frac{-26+3(0.8571)+(0.1111)}{7} =-3.3310

z_2=\frac{1+2(0.8571)-5(-3.7142)}{9} =2.3650

X_2=[1.3718, -3.3310, 2.3650]

\text{diferencias}=|X_2-X_1|

=|[0.8571-1.3718, -3.3310-(-3.7142), 2.3650-0.1111]|

 [1.3718, -3.3310, 2.3650]
-[0.8571, -3.7142, 0.1111] 
__________________________
[0.51474, -0.38322, 2.25397]

\text{diferencias}= |[0.51474, -0.38322, 2.25397]|

\text{errado}=max |X_2-X_1|= 2.2539

itera = 2

X_2=[1.3718, -3.3310, 2.3650]

x_3=\frac{6-(-3.3310)-(2.3650)}{7} = 0.9951

y_3= \frac{-26+3(1.3718)+(2.3650)}{7} = -2.7884

z_3=\frac{1+2(1.3718)-5(-3.3310)}{9} = 2.2665

X_3=[0.99514, -2.7884, 2.2665]

\text{diferencias}=|X_3-X_2|

=|[0.99514-1.3718, -2.7884-(-3.3310), 2.2665-2.3650]|

 [ 0.9951, -2.7884, 2.2665]
-[ 1.3718, -3.3310, 2.3650] 
__________________________
 [-0.3767,  0.5426, 0.0985]

\text{diferencias}=|[-0.3767, 0.5426, 0.0985]|

\text{errado}=max |X_2-X_1|= 0.5425

En cada iteración el valor de error disminuye, por lo que se estima que el método converge.

Al usar el algoritmo, de los resultados se observa que se detiene luego de 14 iteraciones, los errores son menores a la tolerancia:

Iteraciones Jacobi
itera,[X]
     ,errado,|diferencia|
0 [0. 0. 0.]
  nan
1 [ 0.85714 -3.71429  0.11111]
  3.7142857142857144 [0.85714 3.71429 0.11111]
2 [ 1.37188 -3.33107  2.36508]
  2.2539682539682544 [0.51474 0.38322 2.25397]
3 [ 0.99514 -2.78847  2.26657]
  0.5425979915775843 [0.37674 0.5426  0.09851]
4 [ 0.9317 -2.964   1.8814]
  0.3851635892452223 [0.06344 0.17553 0.38516]
...
  0.0003598675094789 [2.30116e-04 3.59868e-04 6.47473e-05]
13 [ 1.00004 -3.00002  2.00007]
  0.0002510632800638 [4.21600e-05 1.07871e-04 2.51063e-04]
14 [ 0.99999 -2.99997  2.00002]
  5.393476706316e-05 [5.12763e-05 5.39348e-05 5.05593e-05]
Metodo de Jacobi
numero de condición: 1.9508402675105447
X:  [ 0.99999 -2.99997  2.00002]
errado: 5.3934767063168465e-05
iteraciones: 14

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..

4. Algoritmo en Python

Para el algoritmo se usa el sistema de ecuaciones en su forma matricial A.X=B. Se asume que se ha aplicado el pivoteo parcial por filas y matriz aumentada (luego se añade como una función). Se asegura que A, X y B sean arreglos de números reales.

Las instrucciones inician seleccionando la primera ecuación, primera fila, i=0. El numerador de la expresión se acumula en suma, que inicia con B[i]. Luego se restan todas las combinaciones de A[i,j]*X[j], menos en la diagonal. El nuevo valor de X[i], o Xnuevo, es suma dividido para el coeficiente de la diagonal A[i,i]

X_i = \bigg(B_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}

i=0 # una ecuacion
suma = B[i] # numerador
for j in range(0,m,1): # un coeficiente
    if (i!=j): # excepto diagonal de A
        suma = suma-A[i,j]*X[j]
Xnuevo[i] = suma/A[i,i]

Se repite el proceso para las siguientes ecuaciones, cambiando el valor de ‘i’

En cada iteración se determinan las diferencias |X[i+1]-X[i]| como los errores de iteración para cada variable.

diferencia = abs(Xnuevo-X)
errado = np.max(diferencia)

Como las diferencias es un vector y la tolerancia solo un número (escalar), se simplifica usando la norma np.max(diferencia).
De ser necesario, se repite la iteración, actualizando el vector X.

# Método de Jacobi
import numpy as np

# INGRESO
# NOTA: Para AB, se ha usado pivoteo parcial por filas
A = [[ 7,  1,  1],
     [-3,  7, -1],
     [-2,  5,  9]]
B = [6, -26, 1]

X0 = [0,0,0]
tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A) # tamaño A
n = tamano[0]
m = tamano[1]

# valores iniciales
diferencia = np.ones(n, dtype=float)
errado = 2*tolera # np.max(diferencia)

itera = 0 
X = np.copy(X0)
Xnuevo = np.copy(X0)
while errado>tolera and itera<=iteramax:
    
    for i in range(0,n,1): # una ecuacion
        suma = B[i] # numerador
        for j in range(0,m,1): # un coeficiente
            if (i!=j): # excepto diagonal de A
                suma = suma-A[i,j]*X[j]
        Xnuevo[i] = suma/A[i,i]
    
    diferencia = abs(Xnuevo-X)
    errado = np.max(diferencia)
    
    print(itera, X)
    print('  ',errado,diferencia)
    
    X = np.copy(Xnuevo) # actualiza X
    itera = itera + 1

if (itera>iteramax): # No converge
    X = np.nan
    print('No converge,iteramax superado')

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', itera)

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..

5. Algoritmo en Python como función

Las instrucciones del algoritmo, empaquetadas como función, añaden algunas facilidades para el desarrollo del ejercicio. Antes de usar el método de Jacobi, se aplica el pivoteo parcial por filas para mejorar de ser posible la convergencia.

En la función jacobi() se añaden las instrucciones para mostrar los resultados parciales en cada iteración. Se añade una tabla para revisión de valores al final del algoritmo, o para realizar la gráfica de errores por iteración.

Si el sistema es de 3×3 se puede añadir una gráfica de los resultados por iteración, mostrando la trayectoria descrita por los resultados parciales en 3D.

# 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
# Metodos Numericos
import numpy as np

def jacobi(A,B,X0, tolera, iteramax=100,
           vertabla=False, precision=4):
    ''' Método de Jacobi, tolerancia, vector inicial X0
        para mostrar iteraciones y tabla: vertabla=True
    '''
    # Matrices como arreglo, numeros reales
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    X0 = np.array(X0,dtype=float)
    tamano = np.shape(A) # tamaño A
    n = tamano[0]
    m = tamano[1]
    
    # valores iniciales
    diferencia = 2*tolera*np.ones(n, dtype=float)
    errado = 2*tolera     # np.max(diferencia)
    tabla = [np.copy(X0)] # tabla de iteraciones
    tabla = np.concatenate((tabla,[[np.nan]]),
                            axis=1) # errado

    if vertabla==True:
        print('Iteraciones Jacobi')
        print('itera,[X]')
        print('     ,errado,|diferencia|')
        print(0,X0)
        print(' ',np.nan)
        np.set_printoptions(precision)

    itera = 0 # Jacobi
    X = np.copy(X0)
    Xnuevo = np.copy(X0)
    while errado>tolera and itera<=iteramax:
        
        for i in range(0,n,1): # una ecuacion
            suma = B[i]
            for j in range(0,m,1): # un coeficiente
                if (i!=j): # excepto diagonal de A
                    suma = suma-A[i,j]*X[j]
            Xnuevo[i] = suma/A[i,i]
        diferencia = abs(Xnuevo-X)
        errado = np.max(diferencia)
        X = np.copy(Xnuevo)

        Xfila = np.concatenate((Xnuevo,[errado]),axis=0)
        tabla = np.concatenate((tabla,[Xfila]),axis=0)

        itera = itera + 1

        if vertabla==True:
            print(itera, Xnuevo)
            print(' ',errado,diferencia)

    # No converge
    if (itera>iteramax):
        X = np.nan
        print('No converge,iteramax superado')

    if vertabla==True:
        X = [X,tabla]
    return(X)

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)

# PROGRAMA Búsqueda de solucion  --------
# INGRESO
A = [[ 7,  1,  1],
     [-3,  7, -1],
     [-2,  5,  9]]
B = [6, -26, 1]

X0 = [0,0,0]
tolera = 0.0001
iteramax = 100
verdecimal = 5

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
n,m = np.shape(AB)
A = AB[:,:n] # separa en A y B
B = AB[:,n]

[X, tabla] = jacobi(A,B,X0,tolera,iteramax,
                    vertabla=True,
                    precision=verdecimal)
n_itera = len(tabla)-1
errado = tabla[-1,-1]

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', n_itera)

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..

5.1 Gráficas de iteraciones

Se muestra la gráfica de errores vs iteraciones.

# GRAFICA de iteraciones
errados = tabla[:,n]

import matplotlib.pyplot as plt
# grafica de error por iteracion
fig2D = plt.figure()
graf = fig2D.add_subplot(111)
graf.plot(errados)
graf.set_xlabel('itera')
graf.set_ylabel('|error|')
graf.set_title('Método de Jacobi: error[itera]')
graf.grid()

plt.show()

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]

5.2 Gráfica de iteraciones para sistema de 3×3

# Grafica de iteraciones x,y,z
# solo para 3x3
xp = tabla[:,0]
yp = tabla[:,1]
zp = tabla[:,2]

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

graf3D.plot(xp,yp,zp)
graf3D.plot(xp[0],yp[0],zp[0],'o',label='X[0]')
graf3D.plot(xp[n_itera],yp[n_itera],zp[n_itera],
            'o',label='X['+str(n_itera)+']')

graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_title('Método de Jacobi: iteraciones X')

# graf3D.view_init(45,45)
plt.legend()

plt.show()

[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]

2017-06-082025-08-16

3.5.2 Número de Condición

Referencia: Chapra 10.3.2 p300/pdf324, Burden 9Ed 7.5 p470, Rodríguez 4.4.3 p133

El número de condición de una matriz A, es una forma de medir la sensibilidad del resultado del sistema de ecuaciones ante pequeños cambios en los valores de la entrada X. El número de condición de una matriz se usa para cuantificar su nivel de mal condicionamiento.

Sea A.X=B un sistema de ecuaciones lineales, entonces:

cond(A) = ||A|| ||A^-1||

es el número de condición de la matriz A.

Los pequeños errores, por ejemplo por truncamiento, pueden amplificarse y afectar la precisión de la solución. Si el número de condición es cercano a 1 o ‘bajo’, implica que la matriz está bien condicionada, errores pequeños en los valores tienen poco impacto en la solución.

Instrucción en Python

 np.linalg.cond(A)

Ejemplo: Si la matriz A es la identidad, el número de condición es 1.

A = [[1,0,0],
     [0,1,0],
     [0,0,1]]
>>> np.linalg.cond(A)
1.0

Otro ejemplo: 1Eva_IT2010_T3_MN Precio artículos

A = [[2,2,4,1],
     [2,2,5,2],
     [4,1,1,2],
     [2,5,2,1]]
np.linalg.cond(A)
18.46408777611575

realizando el pivoteo parcial por filas a la matriz A, no afecta el número de condición.

A = [[4,1,1,2],
     [2,5,2,1],
     [2,2,5,2],
     [2,2,4,1]]
np.linalg.cond(A)
18.464087776115733

Se observará que si el número de condición está alejado de 1, es ‘alto’, al aplicar los métodos iterativos de Jacobi o Gauss-Seidel, resultan NO convergentes. Lo que puede ser un factor para seleccionar el método a usar para encontrar la solución al sistema de ecuaciones.

Tarea

Usando como base los procedimientos desarrollados en Python, elabore un algoritmo para encontrar el número de condición de una matriz.

«el error relativo de la norma de la solución calculada puede ser tan grande como el error relativo de la norma de los coeficientes de [A], multiplicada por el número de condición.»

Por ejemplo,

si los coeficientes de [A] se encuentran a t dígitos de precisión
(esto es, los errores de redondeo son del orden de 10^–t) y
Cond [A] = 10^c,
la solución [X] puede ser válida sólo para t – c dígitos
(errores de redondeo ~ 10^c–t).

verifique el resultado obtenido con el algoritmo, comparando con usar la instrucción

 np.linalg.cond(A)

2017-06-072025-08-16

3.5.1 Normas de Vector o Matriz con Python

Referencia: Chapra 10.3.1 p298 pdf322, Burden 7.1 p320, Rodríguez 4.4.1 p132, MATG1049 Algebra Lineal – Norma, distancias y ángulos

Es una manera de expresar la magnitud de sus componentes:

Sea X un vector de n componentes:

||X|| = \sum_{i=1}^{n}|X_i|

||X|| = max|X_i| , i=1,2, ...,n

||X|| = \left( \sum_{i=1}^{n}X_i^2 \right) ^{1/2}

Sea una matriz A de nxn componentes:

||A|| = max(\sum_{j=1}^{n}|a_{i,j}|, i = 1,2,...,n)

||A|| = max(\sum_{i=1}^{n}|a_{i,j}|, j = 1,2,...,n)

||A|| = \left( \sum_{i=1}^{n}\sum_{j=1}^{n} a_{i,j}^2 \right) ^{1/2}

Ejercicios

Ejercicio 1. Usando los conceptos de normas mostradas, para el siguiente vector:

 x= [5, -3, 2]

a) calcule las normas mostradas (en papel),
b) Realice los respectivos algoritmos en python,
c) Determine los tiempos de ejecución de cada algoritmo. ¿Cúál es el más rápido?

Ejercicio 2. Usando los conceptos de normas mostradas, para la siguiente matriz:

A = [[5, -3, 2],
     [4,  8,-4],
     [2,  6, 1]]

Repita los literales del ejercicio anterior.

Nota: para convertir una lista X en arreglo use: np.array(X)

Normas con Python y Numpy

Algunas normas vectoriales y matriciales con Python. Cálculo del número de condición.

Se presenta un ejemplo usando la matriz A y el vector B en un programa de prueba.

A = [[3,-0.1,-0.2],
     [0.1,7,-0.3],
     [0.3,-0.2,10]]
B = [7.85,-19.3,71.4]

Instrucciones en Python:

# Normas vectoriales y matriciales
# Referencia: Chapra 10.3, p299, pdf323
import numpy as np

def norma_p(X,p):
    Xp = (np.abs(X))**p
    suma = np.sum(Xp)
    norma = suma**(1/p)
    return(norma)

def norma_euclidiana(X):
    norma = norma_p(X,2)
    return(norma)

def norma_filasum(X):
    sfila = np.sum(np.abs(X),axis=1)
    norma = np.max(sfila)
    return(norma)

def norma_frobenius(X):
    tamano = np.shape(X)
    n = tamano[0]
    m = tamano[1]
    norma = 0
    for i in range(0,n,1):
        for j in range(0,m,1):
            norma =  norma + np.abs(X[i,j])**2
    norma = np.sqrt(norma)
    return(norma)

def num_condicion(X):
    M = np.copy(X)
    Mi = np.linalg.inv(M)
    nM = norma_filasum(M)
    nMi= norma_filasum(Mi)
    ncondicion = nM*nMi
    return(ncondicion)

# Programa de prueba #######
# INGRESO
A = [[3,-0.1,-0.2],
     [0.1,7,-0.3],
     [0.3,-0.2,10]]
B = [7.85,-19.3,71.4]
p = 2

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

normap    = norma_p(B, p)
normaeucl = norma_euclidiana(B)
normafilasuma   = norma_filasum(A)
numerocondicion = num_condicion(A)

# SALIDA
print('vector:',B)
print('norma p: ',2)
print(normap)

print('norma euclididana: ')
print(normaeucl)

print('******')
print('matriz: ')
print(A)
print('norma suma fila: ',normafilasuma)

print('n mero de condici n:')
print(numerocondicion)

cuyos resultados del ejercicio serán:

vector: [  7.85 -19.3   71.4 ]
norma p:  2
74.3779033046778
norma euclididana: 
74.3779033046778
******
matriz: 
[[ 3.  -0.1 -0.2]
 [ 0.1  7.  -0.3]
 [ 0.3 -0.2 10. ]]
norma suma fila:  10.5
número de condición:
3.6144243248254124
>>>

compare sus resultados con las funciones numpy:

np.linalg.norm(A)
np.linalg.cond(A)

Referencias: [1] http://www.numpy.org/devdocs/reference/generated/numpy.linalg.norm.html

[2] http://www.numpy.org/devdocs/reference/generated/numpy.linalg.cond.html

2017-06-062025-09-12

3.5 Normas de vector o matriz como distancias con Python

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

2017-06-052025-08-16

3.4 Método de Gauss-Jordan con Python

[ 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.8×10^-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)
    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 ]

2017-06-042025-08-16

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:

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

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