Establece relaciones simples entre los puntos dados que describen la función f(x). Las tabla de diferencias finitas es un elemento base para varios métodos de interpolación, por lo que se trata como un tema inicial.
Si la función f(x) de donde provienen los datos es un polinomio de grado n, entonces la n-ésima diferencia finita será una constante, y las siguientes diferencias se anularán.
Para crear la tabla de diferencias finitas, las primeras columnas requieren concatenar los valores de los índices i, xi y fi.
xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
Los índices i se crean en un vector ki, pues la variable i es usada como fila en matrices, por lo que evitamos confundirlas al usar la variable.
A la matriz con las tres columnas descritas, se le añade a la derecha una matriz de nxn para calcular las diferencias.
Se calculan las diferencias para cada columna, realizando la operación entre los valores de cada fila. Considere que no se realizan cálculos desde la diagonal hacia abajo en la tabla, los valores quedan como cero.
Al final se muestra el título y el resultado de la tabla.
# Diferencias finitas: [tabla_etiq,tabla]# Tarea: verificar tamaño de vectoresimport numpy as np
# INGRESO, Datos de prueba
xi = [0.10, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
# PROCEDIMIENTO
casicero = 1e-15 # redondear a cero# Matrices como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Tabla de Diferencias Finitas
tabla_etiq = ['i','xi','fi']
ki = np.arange(0,n,1) # filas
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)
n,m = np.shape(tabla)
# Calular tabla de Diferencias Finitas
diagonal = n-1 # derecha a izquierda
j = 3 # inicia en columna 3while (j < m): # columna
tabla_etiq.append('d'+str(j-2)+'f')
i = 0 # filawhile (i < diagonal): # antes de diagonal
tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
ifabs(tabla[i,j])<casicero: # casicero revisa
tabla[i,j]=0
i = i + 1
diagonal = diagonal - 1
j = j + 1
# SALIDAprint('Tabla Diferencia Finita: ')
print([tabla_etiq])
print(tabla)
Para tener una referencia visual sobre las primeras diferencias finitas, en una gráfica se trazan las líneas horizontales que pasan por cada punto. Para las segundas diferencias se debe incluir en la gráfica las primeras diferencias finitas vs xi repitiendo el proceso. Las líneas de distancia se añadieron con un editor de imágenes.
Las instrucciones adicionales al algoritmo anterior para añadir la gráfica son:
# Gráfica --------------import matplotlib.pyplot as plt
titulo = 'Diferencia Finita'for i inrange(0,n,1):
plt.axhline(fi[i],ls='--', color='yellow')
if i<(n-1):
plt.annotate("", xytext=(xi[i+1], fi[i]),
xy=(xi[i+1], fi[i+1]),
arrowprops=dict(arrowstyle="<->"),
color='magenta')
plt.annotate("Δ1f["+str(i)+"]",
xy=(xi[i+1], fi[i+1]),
xytext=(xi[i+1],(fi[i]+fi[i+1])/2))
plt.plot(xi,fi,'o', label = 'Puntos')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title(titulo)
plt.show()
np.diff(vector,orden)calcula la diferencia de elementos consecutivos a lo largo de un eje específico de una matriz. Incluye un parámetro de orden para la n-esima diferencia finita. Ejemplo:
Δ2xi es un valor casicero que se debe incluir en el algoritmo para redondear a cero cuando sea por ejemplo 10-15. Asunto que ya fue incluido desde Método de Gauss en la unidad de sistemas de ecuaciones
Es posible generar un sistema de ecuaciones para p(x) haciendo que pase por cada uno de los puntos o coordenadas.
Por ejemplo si se toma el primer punto con xi=0 y fi=1 se establece la ecuación:
Note que ahora las incógnitas son los coeficientes ai. Luego se continúa con los otros puntos seleccionados hasta completar las ecuaciones necesarias para el grado del polinomio seleccionado.
El sistema obtenido se resuelve con alguno de los métodos conocidos para la Solución a sistemas de ecuaciones, que requieren escribir las ecuaciones en la forma de matriz A y vector B, desarrollar la matriz aumentada, pivotear por filas, etc.
La matriz A también se conoce como Matriz Vandermonde D, que se construye observando que los coeficientes se elevan al exponente referenciado al índice columna pero de derecha a izquierda. La última columna tiene valores 1 por tener como coeficiente el valor de xi0.
Para enfocarnos en la interpolación, en la solución se propone usar un algoritmo o función en Python, obteniendo el siguiente resultado:
los coeficientes del polinomio:
[ 41.66666667 -27.5 6.83333333 1. ]
a partir del cual que se construye el polinomio con los valores obtenidos.
Para desarrollar el ejercicio, se realiza un bloque para construir la matriz A y vector B, usando los vectores que representan los puntos de muestra de la función o experimento.
xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]
Observe que en éste caso los datos se dan en forma de lista y se deben convertir hacia arreglos.
La matriz A o Matriz Vandermonde D, se construye observando que los coeficientes del vector xi se elevan al exponente referenciado al índice columna pero de derecha a izquierda.
Construir el polinomio consiste en resolver el sistema de ecuaciones indicado en la sección anterior. Para simplificar la solución del sistema, se usa Numpy, que entrega el vector solución que representan los coeficientes del polinomio de interpolación.
Para construir la expresión del polinomio, se usa la forma simbólica con Sympy, de forma semejante a la usada para construir el polinomio de Taylor.
Para mostrar el polinomio de una manera más fácil de interpretar se usa la instrucción sym.pprint(), usada al final del algoritmo.
# El polinomio de interpolación y Vandermondeimport numpy as np
import sympy as sym
# INGRESO
xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]
# PROCEDIMIENTO# Matrices como arreglo, numeros reales
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype=float)
for i inrange(0,n,1):
for j inrange(0,n,1):
potencia = (n-1)-j # Derecha a izquierda
D[i,j] = xi[i]**potencia
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,fi)
# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0*x # sym.S.Zerofor i inrange(0,n,1):
potencia = (n-1)-i # Derecha a izquierda
termino = coeficiente[i]*(x**potencia)
polinomio = polinomio + termino
# polinomio para evaluación numérica
px = sym.lambdify(x,polinomio)
# SALIDAprint('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)
Para facilitar la evaluación numérica del polinomio, se convierte el polinomio a la forma Lambda px. La gráfica se realiza con un número de muestras suficientes para suavizar la curva dentro del intervalo [a,b], por ejemplo 21, con lo que se comprueba que la curva pase por todos los puntos de la tabla xi,fi dados en el ejercicio.
Instrucciones adicionales al algoritmo para la gráfica:
# GRAFICAimport matplotlib.pyplot as plt
muestras = 21 # muestras = tramos+1
a = np.min(xi) # intervalo [a,b]
b = np.max(xi)
xk = np.linspace(a,b,muestras)
yk = px(xk)
# Usando evaluación simbólica##yk = np.zeros(muestras,dtype=float)##for k in range(0,muestras,1):## yin[k] = polinomio.subs(x,xk[k])
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xk,yk, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()
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:
# Matriz Inversa con Gauss-Jordan# AI es la matriz aumentada A con Identidadimport 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)
defpivoteafila(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 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,
'con', 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
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:',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,:]
for j inrange(0,m,1): # casicero revisaifabs(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)
defgauss_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 inrange(ultfila,0-1,-1):
pivote = AB[i,i]
atras = i-1 # arriba de la fila iif vertabla==True:
print(' fila i:',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,:]
# redondeo a cerofor j inrange(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 unosif 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 cerofor i inrange(0,n,1):
for j inrange(0,m,1):
if np.abs(A_verifica[i,j])<=casicero:
A_verifica[i,j]=0
# SALIDAprint('A_inversa: ')
print(A_inversa)
print('verifica A.A_inversa = Identidad:')
print(A_verifica)
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 = [[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 filafor i inrange(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 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:
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
# 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 cerodefpivoteafila(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,
'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 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])
# Matrices L y U# Modificando el método de Gaussimport numpy as np
defpivoteafila(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 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,
'con', 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,:]
for j inrange(0,m,1): # casicero revisaifabs(AB[k,j])<casicero:
AB[k,j]=0
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
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]
# SALIDAprint('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:
# Solución usando Matrices L y U# continuación de ejercicio anteriorimport 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 cerodefpivoteafila(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 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,
'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 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)
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 X0, 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:
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.
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:
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.
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] # numeradorfor j inrange(0,m,1): # un coeficienteif (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-Seidelimport 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 inrange(0,n,1): # una ecuacion
suma = B[i] # numeradorfor j inrange(0,m,1): # un coeficienteif (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)
# SALIDAprint('Método de Gauss-Seidel')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', itera)
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 3x3 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.162import numpy as np
defgauss_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 erradoif 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 inrange(0,n,1): # una ecuacion
suma = B[i]
for j inrange(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 convergeif (itera>iteramax):
X = np.nan
print('No converge,iteramax superado')
return(X,tabla)
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')
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)
# SALIDAprint('Metodo de Gauss-Seidel')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', n_itera)
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()
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 X0, 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.
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 xi de la fila i, que se usará en el algoritmo en Python.
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,
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:
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.
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]
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 Jacobiimport 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 inrange(0,n,1): # una ecuacion
suma = B[i] # numeradorfor j inrange(0,m,1): # un coeficienteif (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)
# SALIDAprint('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', itera)
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 3x3 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 Numericosimport numpy as np
defjacobi(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) # erradoif 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 inrange(0,n,1): # una ecuacion
suma = B[i]
for j inrange(0,m,1): # un coeficienteif (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 convergeif (itera>iteramax):
X = np.nan
print('No converge,iteramax superado')
if vertabla==True:
X = [X,tabla]
return(X)
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')
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)
# SALIDAprint('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', n_itera)
Si el sistema es de 3x3, 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)
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()
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
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] = 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
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, 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 = [[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)
# 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)
