3.7 Método de Gauss-Seidel con Python

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


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

El método de Gauss-Seidel realiza operaciones semejantes al método de Jacobi. Gauss-Seidel itera

El método de Gauss-Sidel también usa el vector inicial X0, la diferencia consiste en que la actualización del vector X en cada iteración se realiza por cada nuevo valor del vector que se ha calculado. Por lo que las iteraciones llegan más rápido al punto cuando el método converge.

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


2. Ejercicio

Referencia: Chapra ejemplo 11.3  p311

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fue necesario desarrollar esa parte en el ejercicio.

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

Encuentre la solución usando el método de Gauss-Seidel y tolerancia de 0.00001

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


3. Desarrollo analítico

Para el desarrollo con lápiz y papel, se despeja una variable de cada ecuación, se empieza con la primera expresión para obtener x0

3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 3 x_0 = 7.85 + 0.1 x_1 + 0.2 x_2 x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3}

con la segunda ecuación se obtiene x1

0.1 x_0 +7x_1 -0.3 x_2 = -19.3 7x_1 = -19.3 - 0.1 x_0 +0.3 x_2 x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7}

usando la tercera ecuación se obtiene x2

0.3 x_0 - 0.2 x_1 + 10 x_2 = 71.4 10 x_2 = 71.4 - 0.3 x_0 + 0.2 x_1 x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}

Las ecuaciones listas para usar en el algoritmo son:

x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3} x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7} x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}

Vector inicial

Como el vector inicial no se indica en el enunciado, se considera usar el vector de ceros para iniciar  las iteraciones.

X = [0,0,0]

con tolerancia de 0.00001


Iteraciones

Con las ecuaciones obtenidas en el planteamiento, se desarrolla usando el vector inicial presentado, hasta que el |error|<tolera.

El método de Gauss -Seidel actualiza el vector de respuestas Xi+1 luego de realizar cada cálculo. es decir, aprovechas las aproximaciones de cada iteración en el momento que se realizan. Observe la diferencia con el método de Jacobi que espera a terminar con la iteración para actualizar Xi+1

itera = 0

x_0 = \frac{7.85 + 0.1 (0) + 0.2 (0)}{3} = 2.61 x_1 = \frac{ -19.3 - 0.1 (2.61) +0.3 (0)}{7} = -2.79 x_2 = \frac{71.4 - 0.3 (2.61) + 0.2 (-2.79)}{10} = 7.00 X_1 = [2.61, -2.79, 7.00] diferencias = [2.61, -2.79, 7.00] - [0,0,0] diferencias = [2.61, -2.79, 7.00] error = max |diferencias|= 7.00

itera = 1

X = [2.61, -2.79, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.79) + 0.2 (7.00)}{3} = 2.99 x_1 = \frac{ -19.3 - 0.1 (2.99) +0.3 (7.00)}{7} = -2.49 x_2 = \frac{71.4 - 0.3 (2.99) + 0.2 (-2.49)}{10} = 7.00 X_1 = [2.99, -2.49, 7.00] diferencias = [2.99, -2.49, 7.00] - [2.61, -2.79, 7.00] diferencias = [0.38, 0.3 , 0. ] error = max |diferencias| = 0.38

itera = 2

X = [2.99, -2.49, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.49) + 0.2 (7.00)}{3} = 3.00 x_1 = \frac{ -19.3 - 0.1 (3.00) +0.3 (7.00)}{7} = -2.5 x_2 = \frac{71.4 - 0.3 (3.00) + 0.2 (-2.5)}{10} = 7.00 X_3 = [3.00, -2.50, 7.00] diferencias = [3.00, -2.50, 7.00] - [2.99, -2.49, 7.00] diferencias = [ 0.01, -0.01, 0. ] error = max |diferencias|= 0.01

El error aún es mayor que tolera, por lo que es necesario continuar con las iteraciones.

Observaciones

El error disminuye en cada iteración, por lo que el método converge hacia la raíz.

El ejercicio por tener solo 3 incógnitas, la solución se puede interpretar como la intersección de planos en el espacio.

Gauss-Seidel planos

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


4. Algoritmo con Python

Con la descripción dada para el método de Gauss-Seidel, se muestra una forma básica de implementar el algoritmo.

Referencia: Chapra ejemplo 11.3  p311

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fue implementado el procedimiento. Para generalizar el algoritmo, incluya como tarea aumentar el procedimiento de pivoteo por filas.

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

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

B = np.array([7.85,-19.3,71.4])

X0  = np.array([0.,0.,0.])

tolera = 0.00001
iteramax = 100

# PROCEDIMIENTO

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)

que da como resultado:

respuesta X: 
[[ 3. ]
 [-2.5]
 [ 7. ]]
verificar A.X=B: 
[[  7.84999999]
 [-19.3       ]
 [ 71.4       ]]
>>> 

que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.


Tarea

Completar el algoritmo para usar una matriz diagonal dominante, usando al menos el pivoteo parcial por filas.

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


5. Algoritmo como función

Se agrupa las instrucciones como función. Recuerde que la matriz AB tiene que pivotearse por filas antes de usar en el algoritmo. Se obtienen los siguientes resultados:

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
Iteraciones Gauss-Seidel
itera,[X],errado
0 [0. 0. 0.] 2e-05
1 [ 2.6166667 -2.7945238  7.0056095] 7.005609523809525
2 [ 2.9905565 -2.4996247  7.0002908] 0.3738898412698415
3 [ 3.0000319 -2.499988   6.9999993] 0.00947538997430053
4 [ 3.0000004 -2.5        7.       ] 3.154544153582961e-05
5 [ 3.  -2.5  7. ] 3.544137046063156e-07
numero de condición: 3.335707415629964
respuesta con Gauss-Seidel
[ 3.  -2.5  7. ]

Instrucciones en Python como función

Se ha incorporado la función de pivoteo por filas.

# 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
    '''
    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]
    diferencia = 2*tolera*np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X],errado')
        print(itera, X, errado)
        np.set_printoptions(precision)
    while (errado>tolera and itera<iteramax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    xi = xi-A[i,j]*X[j]
            xi = xi/A[i,i]
            diferencia[i] = np.abs(xi-X[i])
            X[i] = xi
        errado = np.max(diferencia)
        itera = itera + 1
        if vertabla==True:
            print(itera, X, errado)
        
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    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 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]

X0  = [0.,0.,0.]

tolera = 0.00001
iteramax = 100
verdecimal = 7

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

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

respuesta = gauss_seidel(A,B,X0, tolera,
                         vertabla=True, precision=verdecimal)

# SALIDA
print('numero de condición:', ncond)
print('respuesta con Gauss-Seidel')
print(respuesta)

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

3.6 Método de Jacobi con Python

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


1. El método iterativo de Jacobi

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

La analogía presentadas entre la «norma como distancia 3D» y el «error de acoplamiento de aeronaves»,  

permite considerar desde un punto de partida o inicial las aproximaciones o iteraciones sucesivas hacia una solución del sistema de ecuaciones.

Las iteraciones pueden ser convergentes o no.

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.

Para describir el método iterativo de Gauss-Seidel, se usa un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante.

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

Para facilitar la escritura del algoritmo, note el uso de índices ajustado a la descripción de arreglos en Python para la primera fila i=0 y primera columna j=0.

Semejante a despejar una variable de la ecuación para representar un plano, se plantea despejar una variable de cada ecuación. Se obtiene así los valores de cada xi, por cada por cada fila i:

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

Observe que el patrón de cada fórmula para determinar x[i], tiene la forma:

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

La parte de la sumatoria se realiza para cada término de A[i,j] en la fila i, excepto para el término de la diagonal A[i,i].

Si se tiene conocimiento del problema planteado y se puede «intuir» o «suponer» una solución para el vector X. Por ejemplo, iniciando con el vector cero, es posible calcular un nuevo vector X usando las ecuaciones para cada X[i] encontradas.

X = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}

Con cada nuevo valor actualiza el el vector X, en  Xnuevo, se calcula el vector diferencias entre X y el vector Xnuevo .

El error a prestar la atención es al mayor valor de las diferencias; se toma como condición para repetir la evaluación de cada vector.

     nuevo = [ 0, 0,  0]
         X = [x0, x1, x2]
diferencia = [|0 - x0|, |0 - x1|, |0 - x2|]
errado = max(diferencia)

Se observa los resultados de errado para cada iteración, relacionados con la convergencia. Si luego de «muchas» iteraciones se encuentra que (errado>tolera),  se detiene el proceso.

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


2. Ejercicio

Referencia:  Steven C. Chapra. Numerical Methods 7th Edition,  Ejercicio 12.39 p339; 1Eva_IT2018_T3 Temperatura en nodos de placa.

La temperatura en los nodos de la malla de una placa se puede calcular con el promedio de las temperaturas de los 4 nodos vecinos de la izquierda, derecha, arriba y abajo. Placa Temperatura

Una placa cuadrada de 3 m de lado tiene la temperatura en los nodos de los bordes como se indica en la figura,

a) Plantee el sistema de ecuaciones y resuelva con el método de Jacobi  para encontrar los valores de a, b, c, d.

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


3. Desarrollo Analítico

Solución Propuesta: s1Eva_IT2018_T3 Temperatura en nodos de placa

a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:

a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

que reordenando se convierte en:

4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

simplificando:

4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

que a forma matricial se convierte en:

A = [[ 4, -1, -1, 0],
     [ 1, -4,  0, 1],
     [ 1,  0, -4, 1],
     [ 0,  1,  1,-4]]
B = [150, -80,-160,-90]

Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

El número de condición es: np.linalg.cond(A) = 3.0

que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

X_0 = [60,40,70,50]

reemplazando los valores iniciales en cada ecuación sin cambios.

iteración 0
a=\frac{50+70+100+40}{4} = 65

b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50 X1 = [65, 47.5, 67.5, 50 ] diferencia = [65-60,47.5-40, 67.5-70, 50-50] = [5, 7.5, -2.5, 0] errormax = max|diferencia| = 7.5

iteración  1
a=\frac{50+67.5+100+47.5}{4} = 66.25

b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.5+30}{4} = 51.25 X_2 = [66.25 48.75 68.75 51.25] diferencia = [66.25-65,48.75-47.5, 68.75-67.5,51.25-50] = [ 1.25, 1.25, 1.25, 1.25 ] errormax = max|diferencia| = 1.25

iteración 2
a=\frac{50+68.75+100+48.75}{4} = 66.875

b=\frac{66.25+30+50+51.25}{4} = 49.375 c=\frac{66.25+60+100+51.25}{4} = 69.375 d=\frac{48.75+60+68.75+30}{4} = 51.875 X_3 = [66.875 49.375 69.375 51.875] diferencia = [ 66.875-66.25, 49.38-48.75, 69.3875-68.75, 51.875-51.25 ] = [ 0.625, 0,63, 0.6375, 0.625 ] errormax = max|diferencia| = 0.625

siguiendo las iteraciones se debería llegar a:

>>> np.linalg.solve(A,B)
array([67.5, 50. , 70. , 52.5])

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


4. Algoritmo en Python

Tarea: Realice el algoritmo de Jacobi en Python, incluyendo el pivoteo por filas de la matriz.

Resultados con el algoritmo

Matriz aumentada
[[   4.   -1.   -1.    0.  150.]
 [   1.   -4.    0.    1.  -80.]
 [   1.    0.   -4.    1. -160.]
 [   0.    1.    1.   -4.  -90.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Jacobi
itera,[X],errado
0 [60. 40. 70. 50.] 1.0
1 [65.  47.5 67.5 50. ] 7.5
2 [66.25 48.75 68.75 51.25] 1.25
3 [66.875 49.375 69.375 51.875] 0.625
4 [67.1875 49.6875 69.6875 52.1875] 0.3125
5 [67.34375 49.84375 69.84375 52.34375] 0.15625
6 [67.42188 49.92188 69.92188 52.42188] 0.078125
7 [67.46094 49.96094 69.96094 52.46094] 0.0390625
8 [67.48047 49.98047 69.98047 52.48047] 0.01953125
9 [67.49023 49.99023 69.99023 52.49023] 0.009765625
10 [67.49512 49.99512 69.99512 52.49512] 0.0048828125
11 [67.49756 49.99756 69.99756 52.49756] 0.00244140625
12 [67.49878 49.99878 69.99878 52.49878] 0.001220703125
13 [67.49939 49.99939 69.99939 52.49939] 0.0006103515625
14 [67.49969 49.99969 69.99969 52.49969] 0.00030517578125
15 [67.49985 49.99985 69.99985 52.49985] 0.000152587890625
16 [67.49992 49.99992 69.99992 52.49992] 7.62939453125e-05
numero de condición: 3.000000000000001
respuesta X: 
[67.49992 49.99992 69.99992 52.49992]
iterado, incluyendo X0:  17

Instrucciones en Python

# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada
# Método de Jacobi

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: vertabla=True
    '''
    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]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)
    xnuevo = np.copy(X0)
    tabla = [np.copy(X0)]

    itera = 0
    if vertabla==True:
        print('Iteraciones Jacobi')
        print('itera,[X],errado')
        print(itera, xnuevo, errado)
        np.set_printoptions(precision)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        X = np.copy(xnuevo)
        
        tabla = np.concatenate((tabla,[X]),axis = 0)

        itera = itera + 1
        if vertabla==True:
            print(itera, xnuevo, errado)

    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    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

# INGRESO
A = [[ 4, -1, -1,  0],
     [ 1, -4,  0,  1],
     [ 1,  0, -4,  1],
     [ 0,  1,  1, -4]]
B = [150, -80,-160,-90]

X0 = [60,40,70,50]
tolera = 0.0001
iteramax = 100
verdecimal = 5

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

[X, puntos] = jacobi(A,B,X0,tolera,vertabla=True,precision=verdecimal)
iterado = len(puntos)

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

# SALIDA
print('numero de condición:', ncond)
print('respuesta X: ')
print(X)
print('iterado, incluyendo X0: ', iterado)

La gráfica de puntos por iteraciones en 3D mostrada al inicio se desarrolla en la propuesta de solución al ejercicio:

s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

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

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 se usa para cuantificar su nivel de mal condicionamiento.

Sea AX=B un sistema de ecuaciones lineales, entonces:

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

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

Instrucción en Python

 np.linalg.cond(A)

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] = 10c,
  • la solución [X] puede ser válida sólo para t – c dígitos
    (errores de redondeo ~ 10c–t).

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

 np.linalg.cond(A)

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 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 = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

B = np.array([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 = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

B = np.array([7.85,-19.3,71.4])

p = 2

# PROCEDIMIENTO
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.3779033047
norma euclididana: 
74.3779033047
******
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.61442432483
>>> 

compare sus resultados con las funciones numpy:

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

http://www.numpy.org/devdocs/reference/generated/numpy.linalg.norm.html

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

3.5 Normas de vector o matriz como distancias con Python

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

Normas de un vector en 3D

La norma de un vector se interpreta como una distancia entre la coordenada definida por el vector [xi, yi, zi] y el origen [0,0,0]. También se puede realizar respecto a otro punto de referencia, se conoce como Norma Ecuclidiana 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 |Xi+1 – Xi| que se interpreta mejor como una distancia.

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

    X1 =  [ 1  2  3]
    X2 =  [ 2  4 -1]
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 X1 y X2

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}

Distancia entre dos puntos en el espacio, error y Norma

Observe los siguientes videos, sobre acople de aeronaves.

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

2. Acoplamiento con estación espacial internacional ISS. RT en español . 2/Julio/2010.
El carguero ruso Progress M-06M pasó de largo la Estación Espacial Internacional fracasado en su intento de acoplarse

  • ¿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 raiz 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 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 Ecuclidiana, 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.


Algoritmo en 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:

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

X1 =  [1 2 3]
X2 =  [ 2  4 -1]
errado =  [-1 -2  4]
||errado|| =  4.58257569495584
Norma euclidiana :  4.58257569495584

las intrucciones 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 = X1 - X2
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()

3.4.1 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 = np.array([[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 atras

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
# Se aplica Gauss-Jordan(AI)

import numpy as np

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

# PROCEDIMIENTO
casicero = 1e-15 # Considerar como 0

# matriz identidad
tamano = np.shape(A)
n = tamano[0]
identidad = np.identity(n)

# Matriz aumentada

AB = np.concatenate((A,identidad),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

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

# 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
AB2 = np.copy(AB)

# elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    pivote = AB[i,i]
    atras = i-1 
    for k in range(atras,0-1,-1):
        factor = AB[k,i]/pivote
        AB[k,:] = AB[k,:] - AB[i,:]*factor
    # diagonal a unos
    AB[i,:] = AB[i,:]/AB[i,i]

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

# SALIDA
print('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminacion hacia adelante')
print(AB2)
print('eliminación hacia atrás')
print(AB)
print('Inversa de A: ')
print(inversa)

el resultado buscado es:

la matriz inversa es:
[[ 0.2        -0.16470588  0.10588235]
 [-0.4         0.15294118  0.25882353]
 [ 0.2         0.07058824 -0.18823529]]
verificando
A.inversa = identidad
[[ 1.00000000e+00 -1.38777878e-17 -1.38777878e-16]
 [ 2.22044605e-16  1.00000000e+00 -2.22044605e-16]
 [ 5.55111512e-17 -9.71445147e-17  1.00000000e+00]]

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 ]

3.4 Método de Gauss-Jordan con Python

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

..


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.

Se mantienen los procedimientos para:

  • matriz aumentada,
  • pivoteo por filas
  • eliminación hacia adelante

al haber obtenido la «matriz triangular superior» , se aplica el procedimiento:

  • eliminación hacia atrás

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


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»y aumentada, 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 ]
..


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.

Las operaciones se realizan de abajo hacia arriba desde la última fila.
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 = 2
pivote = AB[2,2] = 5
k = i-1 # hacia atrás

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

         [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 indice k

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

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

       [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 x3 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    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 x2, 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 x1 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= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}

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


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.

 

Se obtiene el siguiente resultado con el algoritmo:

Matriz aumentada
[[ 4.   2.   5.  60.7]
 [ 2.   5.   8.  92.9]
 [ 5.   4.   3.  56.3]]
Pivoteo parcial:
  1 intercambiar filas:  0 y 2
[[ 5.   4.   3.  56.3]
 [ 2.   5.   8.  92.9]
 [ 4.   2.   5.  60.7]]
Elimina hacia adelante:
 fila 0 pivote:  5.0
   factor:  0.4  para fila:  1
   factor:  0.8  para fila:  2
 fila 1 pivote:  3.4
   factor:  -0.3529411764705883  para fila:  2
 fila 2 pivote:  5.0
[[ 5.    4.    3.   56.3 ]
 [ 0.    3.4   6.8  70.38]
 [ 0.    0.    5.   40.5 ]]
Elimina hacia Atras:
 fila 2 pivote:  5.0
   factor:  1.3599999999999999  para fila:  1
   factor:  0.6  para fila:  0
 fila 1 pivote:  3.4
   factor:  1.1764705882352942  para fila:  0
 fila 0 pivote:  5.0
[[ 1.00000000e+00  0.00000000e+00 -2.08983158e-16  2.80000000e+00]
 [ 0.00000000e+00  1.00000000e+00  2.61228947e-16  4.50000000e+00]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00  8.10000000e+00]]
solución X: 
[2.8 4.5 8.1]
>>> 

Instrucciones en Python

# Método de Gauss-Jordan
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B

import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

def gauss_eliminaAdelante(AB, vertabla=False,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,:]

                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)

def gauss_eliminaAtras(AB, vertabla=False, precision=5, casicero = 1e-15):
    ''' 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,'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,:]
                
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
 
        AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
    X = np.copy(AB[:,ultcolumna])
    
    if vertabla==True:
        print(AB)
    return(X)

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

3.3.3 Método con Matrices triangulares A=L.U

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 atras para UX=Y

Algoritmo en Python

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

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 = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

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

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)

3.3.2 Matrices triangulares A=L.U con Python

[ Matriz triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]

..


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 triangular ][ Ejercicio ] [ Algoritmo ] [ Función ]
..


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 = np.array([[ 3. , -0.1, -0.2],
              [ 0.1,  7. , -0.3],
              [ 0.3, -0.2, 10. ]], dtype=float)
B = np.array([7.85,-19.3,71.4], dtype=float)

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 triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]
..


3. Algoritmo 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.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
Matriz U: 
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
>>> 

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

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 = [[ 3. , -0.1, -0.2],
     [ 0.1,  7. , -0.3],
     [ 0.3, -0.2, 10. ]]

B = [7.85,-19.3,71.4]

# PROCEDIMIENTO
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 triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]

..


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.00333333  -0.29333333 -19.56166667]
 [  0.           0.          10.01204188  70.08429319]]
matriz L: 
[[ 1.          0.          0.        ]
 [ 0.03333333  1.          0.        ]
 [ 0.1        -0.02712994  1.        ]]
Matriz U: 
[[  3.          -0.1         -0.2          7.85      ]
 [  0.           7.00333333  -0.29333333 -19.56166667]
 [  0.           0.          10.01204188  70.08429319]]
>>> 

Instrucciones en Python

# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

def gauss_eliminaAdelante(AB, vertabla=False,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,:]

                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
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 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 triangular ] [ Ejercicio ] [ Algoritmo ] [ Función ]

3.3.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: Chapra Ejemplo 9.12 p277

Calcular el determinante de la matriz A.

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

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

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

[ 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 = np.array([[3. , -0.1, -0.2],
              [0.1,  7. , -0.3],
              [0.3, -0.2, 10.  ]], dtype=float)

# PROCEDIMIENTO

# Matriz aumentada
AB = np.copy(A) # Para usar algoritmo previo

# 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)
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

        pivoteado = pivoteado + 1 # cuenta cambio fila
        
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
[[ 3.  -0.1 -0.2]
 [ 0.1  7.  -0.3]
 [ 0.3 -0.2 10. ]]
Cambios de fila, pivoteado:  0
eliminación hacia adelante
[[ 3.         -0.1        -0.2       ]
 [ 0.          7.00333333 -0.29333333]
 [ 0.          0.         10.01204188]]
determinante: 
210.35299999999995
>>> 

[ 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)
210.3529999999999

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