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, pdf323import numpy as np
defnorma_p(X,p):
Xp = (np.abs(X))**p
suma = np.sum(Xp)
norma = suma**(1/p)
return(norma)
defnorma_euclidiana(X):
norma = norma_p(X,2)
return(norma)
defnorma_filasum(X):
sfila = np.sum(np.abs(X),axis=1)
norma = np.max(sfila)
return(norma)
defnorma_frobenius(X):
tamano = np.shape(X)
n = tamano[0]
m = tamano[1]
norma = 0
for i inrange(0,n,1):
for j inrange(0,m,1):
norma = norma + np.abs(X[i,j])**2
norma = np.sqrt(norma)
return(norma)
defnum_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)
# SALIDAprint('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)
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.
¿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:
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:
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:
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.
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:
# Método de Gauss-Jordan# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (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)
defgauss_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 Lif vertabla==True:
print('Elimina hacia adelante:')
for i inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
if vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(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 Lif 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)
defgauss_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 inrange(ultfila,0-1,-1):
pivote = AB[i,i]
atras = i-1 # arriba de la fila iif vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(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)
# SALIDAprint('solución X: ')
print(X)
Tarea: implementar caso cuando aparecen ceros en la diagonal para dar respuesta, convertir a funciones cada parte
# Solución usando Matrices L y U# continuación de ejercicio anteriorimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (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 adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i+1
for k inrange(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 inrange(1,n,1):
suma = 0
for j inrange(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 inrange(ultfila,0-1,-1):
suma = 0
for j inrange(i+1,ultcolumna,1):
suma = suma + AB[i,j]*X[j]
b = AB[i,ultcolumna]
X[i] = (b-suma)/AB[i,i]
# SALIDAprint('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)
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.
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.
# Matrices L y U# Modificando el método de Gaussimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (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 adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i+1
for k inrange(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])
# SALIDAprint('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])
# SALIDAprint('triangular inferior L')
print(LU[0])
print('triangular superior U')
print(LU[1])
# Método de Gauss# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (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)
defgauss_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 Lif vertabla==True:
print('Elimina hacia adelante:')
for i inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
if vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(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 Lif 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:
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.
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 filafor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (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 adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i + 1
for k inrange(adelante,n,1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
# calcula determinante
multiplica = 1
for i inrange(0,n,1):
multiplica = multiplica*AB[i,i]
determinante = ((-1)**pivoteado)*multiplica
# SALIDAprint('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:
Los índices de fila y columna en la matriz A[i,j] se usan de forma semejante a la nomenclatura de los textos de Álgebra Lineal. Progresivamente para cada fila, se toma como referencia o pivote el elemento de la diagonal (i=j). Luego, se realizan operaciones con las filas inferiores para convertir los elementos por debajo de la diagonal en cero. Las operaciones incluyen el vector B debido a que se trabaja sobre la matriz aumentada AB.
AB Matriz aumentada y pivoteada por filas:
[[ 5. 4. 3. 56.3]
[ 2. 5. 8. 92.9]
[ 4. 2. 5. 60.7]]
iteración fila 1, operación fila 1 y 2
Para la fila 1, con posición i=0, se usa el elemento ai,i como pivote.
pivote = AB[i,i] = AB[0,0] = 5
Para las filas de que están después de la diagonal se referencian como k.Se obtiene el factor escalar de la operación entre filas de la formula
k = i+1 = 0+1 = 1
A_{k} = A_{k} -A_{i}\frac{a_{k,i}}{pivote}
factor = AB[1,0]/pivote = 2/5
y se realiza la operación entre fila k y la fila i para actualizar la fila k,
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:
Para una fila i, el vector b[i] representa el valor de la constante en la fila i de la matriz aumentada, a[i] se refiere los valores de los coeficientes de la ecuación, de los que se usan los que se encuentran a la derecha de la diagonal.
Las operaciones se realizan de abajo hacia arriba desde la última fila. Para el ejercicio presentado se tiene que:
Se encuentra que la solución al sistema de ecuaciones es:
X= \begin{pmatrix} 2.8\\ 4.5 \\ 8.1 \end{pmatrix}
por sustitución hacia atrás
el vector solución X es:
[[2.8]
[4.5]
[8.1]]
Verificar respuesta
Para verificar que el resultado es correcto, se usa el producto punto entre la matriz a y el vector resultado X. La operación A.X = B debe dar el vector B.
El algoritmo para el Método de Gauss, reutiliza las instrucciones para matriz aumentada y pivoteo parcial por filas.
Recordar: Asegurar que los arreglos sean de tipo Real (float), para que no se considere el vector como entero y realice operaciones entre enteros, generando errores por truncamiento.
La parte nueva a desarrollar corresponde al procedimiento de «eliminación hacia adelante» y el procedimiento de «sustitución hacia atrás».
# Método de Gauss# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
# INGRESO
A = np.array([[4,2,5],
[2,5,8],
[5,4,3]])
B = np.array([[60.70],
[92.90],
[56.30]])
# PROCEDIMIENTO
casicero = 1e-15 # Considerar como 0# Evitar truncamiento en operaciones
A = np.array(A,dtype=float)
# Matriz aumentada
AB = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)
# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
AB1 = np.copy(AB)
# eliminación hacia adelantefor i inrange(0,n-1,1):
pivote = AB[i,i]
adelante = i + 1
for k inrange(adelante,n,1):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - AB[i,:]*factor
# sustitución hacia atrás
ultfila = n-1
ultcolumna = m-1
X = np.zeros(n,dtype=float)
for i inrange(ultfila,0-1,-1):
suma = 0
for j inrange(i+1,ultcolumna,1):
suma = suma + AB[i,j]*X[j]
b = AB[i,ultcolumna]
X[i] = (b-suma)/AB[i,i]
X = np.transpose([X])
# SALIDAprint('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB1)
print('eliminación hacia adelante')
print(AB)
print('solución: ')
print(X)
Tarea
Revisar cuando la matriz pivoteada por filas tienen un elemento cero o muy cercano a cero pues la matriz sería singular. El valor considerado como casi cero podría ser 1×10-15
A estas alturas, por la cantidad de líneas de instrucción es recomendable reutilizar bloques de algoritmos usando funciones def-return. Por ejemplo: pivoteo por filas, eliminación hacia adelante, sustitución hacia atrás.
# Método de Gauss# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (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)
defgauss_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 inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
if vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(adelante,n,1):
if (np.abs(pivote)>=casicero):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - factor*AB[i,:]
if vertabla==True:
print(' factor: ',factor,' para fila: ',k)
else:
print(' pivote:', pivote,'en fila:',i,
'genera division para cero')
if vertabla==True:
print(AB)
return(AB)
defgauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
''' Gauss sustituye hacia atras
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Sustitución hacia atras
X = np.zeros(n,dtype=float)
ultfila = n-1
ultcolumna = m-1
for i inrange(ultfila,0-1,-1):
suma = 0
for j inrange(i+1,ultcolumna,1):
suma = suma + AB[i,j]*X[j]
X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
return(X)
# INGRESO
A = [[4,2,5],
[2,5,8],
[5,4,3]]
B = [60.70,92.90,56.30]
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
AB = gauss_eliminaAdelante(AB,vertabla=True)
X = gauss_sustituyeAtras(AB,vertabla=True)
# SALIDAprint('solución: ')
print(X)
Los métodos de solución a sistemas de ecuaciones en los primeros pasos usan la matriz aumentada y pivoteada por filas. Como es un procedimiento usado en todos los métodos de la unidad, para simplificar se presenta como uno de los primeros algoritmos.
Para mostrar el desarrollo del proceso se usa como referencia un ejercicio.
Ejemplo 1. Un comerciante compra tres productos A, B, C, pero en las facturas únicamente consta la cantidad comprada en Kg y el valor total de la compra.
Se necesita determinar el precio unitario de cada producto. Dispone de solo tres facturas con los siguientes datos:
Ejemplo:
Cantidad
Valor ($)
Factura
X1
X2
X3
Pagado
1
4
2
5
60.70
2
2
5
8
92.90
3
5
4
3
56.30
Los precios unitarios se pueden representar por las variables x1, x2, x3 para escribir el sistema de ecuaciones que muestran las relaciones de cantidad, precio y valor pagado:
Para el pivoteo por fila de la matriz aumentada AB, tiene como primer paso revisar la primera columna desde la diagonal en adelante.
columna = [|4|,
|2|,
|5|]
dondemax = 2
El procedimiento de «pivoteo» se realiza si la posición dónde se encuentra el valor de mayor magnitud no corresponde a la diagonal de la matriz (posición 0 de la columna).
En el ejercicio se encuentra que la magnitud de mayor valor está en la última fila, por lo que en AB se realiza el intercambio entre la fila 3 y la fila 1
Para realizar el algoritmo, es de recordar que para realizar operaciones en una matriz sin alterar la original, se usa una copia de la matriz (np.copy). Se puede comparar y observar los cambios entre la matriz original y la copia a la que se aplicaron cambios
Si no es necesaria la comparación entre el antes y después, no se realiza la copia y se ahorra el espacio de memoria, detalle importante para matrices de «gran tamaño» y una computadora con «limitada» memoria.
# Pivoteo parcial por filas# Solución a Sistemas de Ecuacionesimport numpy as np
# INGRESO
A = np.array([[4,2,5],
[2,5,8],
[5,4,3]])
B = np.array([[60.70],
[92.90],
[56.30]])
# PROCEDIMIENTO# Matriz aumentada
AB = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)
# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en ABfor i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no está en diagonalif (dondemax !=0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
# SALIDAprint('Matriz aumentada:')
print(AB0)
print('Pivoteo parcial por filas')
print(AB)
Los bloques de cada procedimiento que se repiten en otros métodos se convierten a funciones def-return, empaquetando las soluciones algorítmicas a problemas resueltos.
Se usa la matriz M para generalizar y diferenciar de ‘A’ que es usada en los ejercicios en realizados en adelante.
# Pivoteo parcial por filas# Solución a Sistemas de Ecuacionesimport numpy as np
defpivoteafila(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 inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print(AB)
return(AB)
# INGRESO
A = [[4,2,5],
[2,5,8],
[5,4,3]]
B = [60.70, 92.90, 56.30]
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
# SALIDAprint('Resultado de Pivoteo parcial por filas')
print(AB)
