El polinomio de interpolación de Lagrange reformula el polinomio de interpolación de Newton evitando el cálculo de la tabla de diferencias divididas. El método de Lagrange tolera las diferencias entre distancias x de los puntos de muestra.
El polinomio de Lagrange se construye a partir de las fórmulas:
Las expresiones del polinomio contiene los binomios en su forma básica, para resolver y simplificar las ecuaciones se usa polinomio.expand().
Para realizar la gráfica del polinomio es conveniente convertirlo a la forma lambda con Numpy, de esta forma se evalúa en una sola linea todos los puntos para el intervalo [a,b] en x.
# Interpolacion de Lagrange# divisoresL solo para mostrar valoresimport numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO , Datos de prueba
xi = [0, 0.2, 0.3, 0.4]
fi = [1, 1.6, 1.7, 2.0]
# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i inrange(0,n,1):
# Termino de Lagrange
numerador = 1
denominador = 1
for j inrange(0,n,1):
if (j!=i):
numerador = numerador*(x-xi[j])
denominador = denominador*(xi[i]-xi[j])
terminoLi = numerador/denominador
polinomio = polinomio + terminoLi*fi[i]
divisorL[i] = denominador
# simplifica el polinomio
polisimple = polinomio.expand()
# para evaluación numérica
px = sym.lambdify(x,polisimple)
# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)
# SALIDAprint(' valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)
# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()
El método se usa en el caso que los puntos en el «eje x» se encuentran espaciados de forma arbitraria y provienen de una función desconocida pero supuestamente diferenciable.
Para lo cual se debe interpretar la tabla de diferencias divididas, también como una aproximación a una derivada:
i
xi
f[xi]
Primero
Segundo
Tercero
0
x0
f[x0]
f[x1,x0]
f[x2,x1,x0]
f[x3,x2,x1,x0]
1
x1
f[x1]
f[x2,x1]
f[x3,x2,x1]
2
x2
f[x2]
f[x3,x2]
3
x2
f[x3]
En la fórmula del polinomio, las diferencias divididas sirven para evaluar los coeficientes de cada término adicional para aumentar el grado. Semejante al proceso realizado para «diferencias finitas divididas»:
Se dispone de los datos (x, f(x)), en donde x es un valor de inversión y f(x) es un valor de ganancia, ambos en miles de dólares:
inversión
ganancia
3.2
5.12
3.8
6.42
4.2
7.25
4.5
6.85
Considere que los valores invertidos en materia prima para producción, dependen de la demanda de del producto en el mercado, motivo por el que los valores de inversión no guardan distancias equidistantes entre si.
Se toman los datos de la tabla como arreglos para el algoritmo
xi = [3.2 , 3.8 , 4.2 , 4.5 ]
fi = [5.12, 6.42, 7.25, 6.85]
Con los datos se llena la tabla de diferencias divididas, donde por simplicidad, se escriben las operaciones en cada casilla. La última columna o cuarta diferencia dividida es cero por no disponer de datos para hacer el cálculo.
i
xi
f[xi]
Primero
Segundo
Tercero
0
3.2
5.12
\frac{6.42-5.12}{3.8-3.2} =2.1667
\frac{2.075-2.1667}{4.2-3.2} =-0.0917
\frac{-4.869-(-0.0917)}{4.5-3.2} =-3.6749
1
3.8
6.42
\frac{7.25-6.42}{4.2-3.8} =2.075
\frac{-1.3333-2.075}{4.5-3.8} =-4.869
2
4.2
7.25
\frac{6.85-7.25}{4.5-4.2} =-1.3333
3
4.5
6.85
Las diferencias divididas de la primera fila son los valores usados para la expresión del polinomio de interpolación:
El algoritmo para interpolación de Diferencias Divididas o Newton, considera reutilizar el procedimiento de cálculo de diferencias finitas incorporando la parte del denominador.
La creación de la expresión del polinomio también es semejante a la usada para diferencias finitas avanzadas.
Se incorpora la parte gráfica para observar los resultados en el intervalo xi, con el número de muestras = 101 para tener una buena resolución de la línea del polinomio.
# Polinomio interpolación# Diferencias Divididas de Newton# Tarea: Verificar tamaño de vectores,# verificar puntos equidistantes en ximport numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO , Datos de prueba
xi = [3.2, 3.8, 4.2, 4.5]
fi = [5.12, 6.42, 7.25, 6.85]
# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
# Tabla de Diferencias Divididas
titulo = ['i ','xi ','fi ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)
# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
# Añade título para cada columna
titulo.append('F['+str(j-2)+']')
# cada fila de columna
i = 0
paso = j-2 # inicia en 1while (i < diagonal):
denominador = (xi[i+paso]-xi[i])
numerador = tabla[i+1,j-1]-tabla[i,j-1]
tabla[i,j] = numerador/denominador
i = i+1
diagonal = diagonal - 1
j = j+1
# POLINOMIO con diferencias Divididas# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)
# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j inrange(1,n,1):
factor = dDividida[j-1]
termino = 1
for k inrange(0,j,1):
termino = termino*(x-xi[k])
polinomio = polinomio + termino*factor
# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()
# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)
# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)
# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)
# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):## plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()
El resultado del algoritmo se muestra a continuación:
1. Interpolación por Diferencias finitas avanzadas
Referencia: Rodríguez 6.6.4 p221, Burden 9Ed p129
Se usa en interpolación cuando los puntos en el «eje x» se encuentran igualmente espaciados, la diferencia entre puntos consecutivos xi es una constante denominada h.
h = xi+1 – xi
Una relación entre derivadas y diferencias finitas se establece mediante:
f^{(n)}(z) = \frac{\Delta ^{n} f_{0}}{h^{n}} para algún k en el intervalo [x0,xn]
\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]
El polinomio de interpolación se puede construir por medio de diferencias finitas avanzadas con las siguiente fórmula:
Se toman los datos del ejercicio de diferencias finitas , observando que se requiere que el tamaño de paso h sea constante entre los puntos consecutivos xi.
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
El polinomio se construye usando el ejercicio de diferencias finitas.
Para construir la expresión del polinomio añadiendo los términos de la fórmula, se define la variable simbólica x con Sympy.
Para simplificar el polinomio resultante con las expresiones de multiplicación, se utiliza la instrucción sym.expand().
En caso de requerir evaluar la fórmula con un vector de datos se la convierte a la forma lambda.
Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores xi, definiendo el número de muestras = 101 que son suficientes para que la gráfica se observe sin distorsión.
# Polinomio interpolación# Diferencias finitas avanzadas# Tarea: Verificar tamaño de vectores,# verificar puntos equidistantes en ximport numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO , Datos de prueba
xi = [0.1, 0.2, 0.3, 0.4]
fi = [1.45, 1.6, 1.7, 2.0]
# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
fi = np.array(xi,dtype=float)
# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)
# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
# Añade título para cada columna
titulo.append('df'+str(j-2))
# cada fila de columna
i = 0
while (i < diagonal):
tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
i = i+1
diagonal = diagonal - 1
j = j+1
# POLINOMIO con diferencias Finitas avanzadas# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)
# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j inrange(1,n,1):
denominador = np.math.factorial(j)*(h**j)
factor = dfinita[j-1]/denominador
termino = 1
for k inrange(0,j,1):
termino = termino*(x-xi[k])
polinomio = polinomio + termino*factor
# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()
# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)
# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)
# SALIDAprint('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)
# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):## plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()
Tarea: se recomienda realizar las gráficas comparativas entre métodos, debe mostrar la diferencia con los métodos que requieren el tamaño de paso equidistante h, y los que no lo requieren. Permite detectar errores de selección de método para interpolación.
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.
Para un ejemplo se toman los siguientes puntos:
xi
0.10
0.2
0.3
0.4
fi
1.45
1.6
1.7
2.0
La tabla de diferencias finitas se construye tomando los datos y sus índices como parte de las primeras columnas.
i
xi
fi
Δ1fi
Δ2fi
Δ3fi
Δ4fi
0
0.1
1.45
0.15
-0.05
0.25
0.
1
0.2
1.6
0.1
0.2
0.
0.
2
0.3
1.7
0.3
0.
0.
0.
3
0.4
2.0
0.
0.
0.
0.
Cada casilla de diferencia finita Δjfi se obtiene restando los dos valores consecutivos de la columna anterior. Por ejemplo, para la primera columna:
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 índice i, xi y fi.
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([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.
# Tabla de Diferencias finitas# resultado en: [título,tabla]# Tarea: verificar tamaño de vectoresimport numpy as np
# INGRESO, Datos de prueba
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])
# PROCEDIMIENTO# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)
# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)
# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
# Añade título para cada columna
titulo.append('df'+str(j-2))
# cada fila de columna
i = 0
while (i < diagonal):
tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
i = i + 1
diagonal = diagonal - 1
j = j + 1
# SALIDAprint('Tabla Diferencia Finita: ')
print([titulo])
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 graficar 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áficaimport matplotlib.pyplot as plt
for i inrange(0,n,1):
plt.axhline(fi[i],ls='--', color='yellow')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencia Finita')
plt.show()
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ónimport 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# Convierte a arreglos numpy
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
B = fi
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
# Aplicar métodos Unidad03. Tarea# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)
# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i inrange(0,n,1):
potencia = (n-1)-i # Derecha a izquierda
termino = coeficiente[i]*(x**potencia)
polinomio = polinomio + termino
# 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
# Polinomio a forma Lambda x:# para evaluación con vectores de datos xin# muestras = tramos+1
muestras = 21
px = sym.lambdify(x,polinomio)
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)
# Usando evaluación simbólica##yin = np.zeros(muestras,dtype=float)##for j in range(0,muestras,1):## yin[j] = polinomio.subs(x,xin[j])
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()
El método de Gauss-Seidel realiza operaciones semejantes al método de Jacobi.
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.
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
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 iterativosimport 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
whilenot(errado<=tolera or itera>iteramax):
# por filafor i inrange(0,n,1):
# por columna
suma = 0
for j inrange(0,m,1):
# excepto diagonal de Aif (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 convergeif (itera>iteramax):
X=0
# revisa respuesta
verifica = np.dot(A,X)
# SALIDAprint('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)
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:
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.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
'''
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]')
print(' diferencia,errado')
print(itera, X, errado)
np.set_printoptions(precision)
while (errado>tolera and itera<iteramax):
for i inrange(0,n,1):
xi = B[i]
for j inrange(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)
print(' ', diferencia, errado)
# No convergeif (itera>iteramax):
X=itera
print('iteramax superado, No converge')
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')
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)
# SALIDAprint('numero de condición:', ncond)
print('respuesta con Gauss-Seidel')
print(respuesta)
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.
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:
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.
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.
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.
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.
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.
# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada# Método de Jacobiimport 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: 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)
whilenot(errado<=tolera or itera>iteramax):
for i inrange(0,n,1):
nuevo = B[i]
for j inrange(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 convergeif (itera>iteramax):
X=itera
print('iteramax superado, No converge')
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')
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)
# SALIDAprint('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:
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
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:
