Los errores disminuyen en cada iteración, por lo que el método converge,
si se analiza en número de condición de la matriz A es 2, lo que es muy cercano a 1, por lo tanto el método converge.
Revisión de resultados
Resultados usando algoritmos desarrollados en clase:
matriz aumentada:
[[-4. 2. 0. 3. ]
[ 0. 3. 1. 10. ]
[ 2. 0. -3. -6. ]]
Elimina hacia adelante
[[-4. 2. 0. 3. ]
[ 0. 3. 1. 10. ]
[ 0. 1. -3. -4.5]]
Elimina hacia adelante
[[-4. 2. 0. 3. ]
[ 0. 3. 1. 10. ]
[ 0. 0. -3.333333 -7.833333]]
Elimina hacia adelante
[[-4. 2. 0. 3. ]
[ 0. 3. 1. 10. ]
[ 0. 0. -3.333333 -7.833333]]
Elimina hacia atras
[[ 1. -0. -0. 0.525]
[ 0. 1. 0. 2.55 ]
[-0. -0. 1. 2.35 ]]
el vector solución X es:
[[0.525]
[2.55 ]
[2.35 ]]
verificar que A.X = B
[[ 3.]
[10.]
[-6.]]
Número de condición de A: 2.005894
los resultados para [a,b,c]:
[0.525 2.55 2.35 ]
producto punto entre vectores:
v12: 0.0
v13: 1.3322676295501878e-15
v23: 2.0
Algoritmos en Python:
# 1Eva_IT2019_T3 Vector perpendicular a planoimport numpy as np
# INGRESO
A = np.array([[-4.,2,0],
[ 0., 3,1],
[ 2.,0,-3]])
B = np.array([3.,10,-6])
# PROCEDIMIENTO
B = np.transpose([B])
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
print('matriz aumentada: ')
print(AB)
# Gauss elimina hacia adelante# tarea: verificar términos cerofor i inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
for k inrange(adelante,n,1):
if (np.abs(pivote)>=casicero):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - factor*AB[i,:]
print('Elimina hacia adelante')
print(AB)
# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i inrange(ultfila,0-1,-1):
# Normaliza a 1 elemento diagonal
AB[i,:] = AB[i,:]/AB[i,i]
pivote = AB[i,i] # uno# arriba de la fila i
atras = i-1
for k inrange(atras,0-1,-1):
if (np.abs(AB[k,i])>=casicero):
factor = pivote/AB[k,i]
AB[k,:] = AB[k,:]*factor - AB[i,:]
else:
factor= 'division para cero'print('Elimina hacia atras')
print(AB)
X = AB[:,ultcolumna]
# Verifica resultado
verifica = np.dot(A,X)
# SALIDAprint('el vector solución X es:')
print(np.transpose([X]))
print('verificar que A.X = B')
print(np.transpose([verifica]))
numcond = np.linalg.cond(A)
# para comprobar respuestas
v1 = [2,-3,X[0]]
v2 = [X[1],1,-4]
v3 = [3,X[2],2]
v12 = np.dot(v1,v2)
v13 = np.dot(v1,v3)
v23 = np.dot(v2,v3)
# SALIDAprint('\n Número de condición de A: ', numcond)
print('\n los resultados para [a,b,c]:')
print(X)
print('\n productos punto entre vectores:')
print('v12: ',v12)
print('v13: ',v13)
print('v23: ',v23)
Tarea, usar el algoritmo de Jacobi usado en el taller correspondiente.
Por no disponer de valor inicial para TA, considere que el cable colgado no debería tener tensión TA=0 N, pues en la forma x=12/TA se crea una indeterminación. Si no dispone de algún criterio para seleccionar el valor de TA puede iniciar un valor positivo, por ejemplo 120 con lo que el valor de x0=12/120=0.1
Se requiere un polinomio de grado 3 siendo el eje x correspondiente a temperatura. Son necesarios 4 puntos de referencia alrededor de 15 grados, dos a la izquierda y dos a la derecha.
Se observa que los datos en el eje x son equidistantes, h=8, y ordenados en forma ascendente, se cumple con los requisitos para usar diferencias finitas avanzadas. que tiene la forma de:
Resolviendo y simplificando el polinomio, se puede observar que al aumentar el grado, la constante del término disminuye.
p_3(x)=12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3
para el cálculo del error se puede usar un término adicional del polinomio, añadiendo un punto más a la tabla de diferencia finitas. Se evalúa éste término y se estima el error que dado que el término de grado 3 es del orden de 10-5, el error será menor. (Tarea)
El valor de oxígeno usado como referencia es 9, cuyos valores de temperatura se encuentran entre 16 y 24 que se toman como rango inicial de búsqueda [a,b]. Por lo que el polinomio se iguala a 9 y se crea la forma estandarizada del problema:
Para mostrar el procedimiento se realizan solo tres iteraciones,
1ra Iteración
a=16 , b = 24, c = (16+24)/2 = 20
f(a) = 0.9, f(b) = -0.6, f(c) = 0.011
error = |24-16| = 8
como f(c) es positivo, se mueve el extremo f(x) del mismo signo, es decir a
2da Iteración
a=20 , b = 24, c = (20+24)/2 = 22
f(a) = 0.119, f(b) = -0.6, f(c) = -0.251
error = |24-20|= 4
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b
3ra Iteración
a=20 , b = 22, c = (20+22)/2 = 21
f(a) = 0.119, f(b) = -0.251, f(c) = -0.068
error = |22-20| = 2
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b
y así sucesivamente hasta que error< que 10-3
Usando el algoritmo en python se obtendrá la raiz en 20.632 con la tolerancia requerida.
Revisión de resultados
Usando como base los algoritmos desarrollados en clase:
# 1Eva_IT2019_T1 Oxígeno y temperatura en marimport numpy as np
import matplotlib.pyplot as plt
import sympy as sym
definterpola_dfinitas(xi,fi):
'''
Interpolación de diferencias finitas avanzadas
resultado: polinomio en forma simbólica
'''# Tabla de diferencias finitas
titulo = ['i','xi','fi']
n = len(xi)
# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
xi = np.transpose([xi])
fi = np.transpose([fi])
# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
# Aumenta el título para cada columna
titulo.append('df'+str(c-2))
# calcula cada diferencia por fila
f = 0
while (f < diagonal):
tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
f = f+1
diagonal = diagonal - 1
c = c+1
# POLINOMIO con diferencias finitas# caso: puntos en eje x equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
x = sym.Symbol('x')
h = xi[1,0]-xi[0,0]
polinomio = fi[0,0]
for c inrange(1,n,1):
denominador = np.math.factorial(c)*(h**c)
factor = dfinita[0,c-1]/denominador
termino=1
for f inrange(0,c,1):
termino = termino*(x-xi[f])
polinomio = polinomio + termino*factor
# simplifica polinomio, multiplica los (x-xi)
polinomio = polinomio.expand()
return(tabla, polinomio)
# INGRESO
tm = [0.,8,16,24,32,40]
ox = [14.6,11.5,9.9,8.4,7.3,6.4]
xi = [8,16,24,32]
fi = [11.5,9.9,8.4,7.3]
# PROCEDIMIENTO
x = sym.Symbol('x')
# literal a
tabla, polinomio = interpola_dfinitas(xi,fi)
p15 = polinomio.subs(x,15)
# literal b
deriva = polinomio.diff(x,1)
dp16 = deriva.subs(x,16)
px = sym.lambdify(x,polinomio)
xk = np.linspace(np.min(xi),np.max(xi))
pk = px(xk)
# SALIDAprint('Literal a')
print(tabla)
print(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(deriva)
print('dp(16) =', dp16)
# gráfica
plt.plot(tm,ox,'ro')
plt.plot(xk,pk)
plt.axhline(9,color="green")
plt.xlabel('temperatura')
plt.ylabel('concentracion de oxigeno')
plt.grid()
plt.show()
# --------literal c ------------# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = tolera# se convierte forma de símbolos a numéricos
buscar = polinomio-9
fx = sym.lambdify(x,buscar)
# INGRESO
a = 16
b = 24
tolera = 0.001
# PROCEDIMIENTO
tabla = []
tramo = b-a
fa = fx(a)
fb = fx(b)
i = 1
while (tramo > tolera):
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
i = i+1
cambia = np.sign(fa)*np.sign(fc)
if (cambia<0):
b = c
fb = fc
else:
a=c
fa = fc
tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)
raiz = c
# SALIDAprint('\n literal c')
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), paso]')
# print(tabla)# Tabla con formato
n=len(tabla)
for i inrange(0,n,1):
unafila = tabla[i]
formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
unafila = formato.format(*unafila)
print(unafila)
print('raiz: ',raiz)
Siguiendo el desarrollo analítico tradicional, para adecuar la ecuación para los algoritmo de búsquda de raíces de ecuaciones, se reemplazan los valores en la fórmula.
P = A\Big(\frac{1-(1+i)^{-n}}{i} \Big) 70000 = 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)
Como ambos lados de la ecuación deben ser iguales, si se restan ambos se obtiene una ecuación que tiene como resultado cero, que es la forma ideal para usar en el algoritmo que representa f(x) o en este caso f(i)
El intervalo de existencia correspondería a la tasa de interés mínimo y el interés máximo.
[izquierda, derecha] = [a,b]
Para el intervalo se deben tomar en cuenta algunas consideraciones descritas a continuación:
izquierda:
En el extremo izquierdo, las tasas no son negativas, lo que se interpreta en que un banco paga por que le presten dinero.
Tampoco tiene mucho sentido el valor cero, que son prestamos sin intereses. A menos que sean sus padres quienes le dan el dinero.
Un valor inicial para el interés puede ser por ejemplo 1% ó 0.01, siempre que se cumpla que existe cambio de signo en la función a usar.
derecha:
En el extremo derecho, si se propone por ejemplo i con 100%, o 1.00, no tendría mucho sentido un préstamo con intereses al 100% anual, que resulta en el doble del valor inicial en tan solo un periodo o año.
La tasa de interés de consumo que son de las más alto valor, se encuentran reguladas. En Ecuador es un valor alrededor del 16% anuales o 0.16.
Considerando las observaciones iniciales del problema, se propone empezar el análisis para la búsqueda de la raíz en el intervalo en un rango más amplio:
[ 0.01, 0.50]
Ser realiza la comprobación que existe cambio de signo en los extremos del intervalo.
fx(0.001) =- 43935.86
fx(0.50) = 67600.0
Para el ejercicio se hace notar que la es tasa nominal anual, pero los pagos son mensuales. Por lo que se debe unificar las tasas de interes a mensuales. Una aproximación es usar las tasas anuales divididas para los 12 meses del año.
Tolerancia al error
La tolerancia se considera en éste ejercicio como el valor de diferencias (tramo) entre iteraciones con precisión satisfactoria.
Por ejemplo si no negociaremos más con el banco por variaciones de tasas del 0.1% , entonces la tolerancia será de 0.001.
Las publicaciones de tasas en el mercado incluyen dos decimales, por lo que para el ejercicio aumentamos la precisión a : 0.0001
tolera = 1×10-4
Literal c
Se presentan dos formas se solución para el litera c:
– c.1 la requerida en el enunciado con Newton-Raphson
cuyo valr de error está casi dentro del valor de tolerancia,
que permite tomar el último valor como respuesta de tasa mensual
raiz = tasa mensual = 0.01703
Convirtiendo a la tasa tasa anual que es la publicada por las instituciones financieras se tiene que:
tasa anual nominal = 0.01703*12 = 0.2043
Teniendo como resultado una tasa anual de 20.43%
Algoritmo en Python
# 1ra Evaluación II Término 2018# Tema 4. Tasa de interes para hipotecaimport numpy as np
import matplotlib.pyplot as plt
defnewton_raphson(fx, dfx, xi, tolera):
'''
funciónx y fxderiva son de forma numérica
xi es el punto inicial de búsqueda
'''
tramo = abs(2*tolera)
while (tramo>=tolera):
xnuevo = xi - fx(xi)/dfx(xi)
tramo = abs(xnuevo-xi)
print(xi,xnuevo, tramo)
xi = xnuevo
return(xi)
# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fx = lambda i: P*i - A*(1-(1+i)**(-n))
dfx = lambda i: P + A*(-n)*(i+1)**(-n-1)
a = 0.01/12
b = 0.25/12
muestras = 51
tolera = 0.0001
# PROCEDIMIENTO
tasa = np.linspace(a,b,muestras)
fi = fx(tasa)
raiz = newton_raphson(fx, dfx, b, tolera)
tanual = 12*raiz
# SALIDAprint('raiz encontrada en: ', raiz)
print('tasa anual: ',tanual)
# Gráfica
plt.plot(tasa*12,fi)
plt.title('tasa anual de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.axhline(0, color='green')
plt.show()
c.2. Desarrollo con el método de la Bisección
Desarrollo Analítico con Bisección
Como parte del desarrollo del ejercicio se presenta las iteraciones para el algoritmo, tradicionalmente realizadas con una calculadora.
y se continuaría con las iteraciones, hasta cumplir que tramo<=tolera
Tabla de datos obtenidos
tabla para Bisección
i
a
c
b
f(a)
f(c)
f(b)
tramo
1
0.01
0.255
0.5
-43935.86
65294.11
67600.0
0.49
2
0.01
0.1325
0.255
-43935.86
60943.39
65294.11
0.215
3
0.01
0.07125
0.1325
-43935.86
53157.89
60943.39
0.1225
hasta lo calculado la raiz se encontraría en el intervalo [0.01,0.07125] con error esitmado de 0.06125, aún por mejorar con más iteraciones.
Algoritmo en Python para Bisección
El algoritmo bisección usa las variables a y b, por lo que los limites en el intervalo usados son [La,Lb]
para el problema la variable ‘i’ se usa en el eje x.
La selección de cambio de rango [a,b] se hace usando solo el signo del valor.
El algoritmo presentado es tal como se explica en la parte conceptual
Se deja como tarea convertir el algoritmo a funcion def-return de Python.
# 1Eva_IIT2018_T4 Tasa de interés en hipotecaimport numpy as np
import matplotlib.pyplot as plt
# INGRESO
P = 70000.00
A = 1200.00
n = 25*12
fi = lambda i: P - A*(1-((1+i)**-n))/i
# Intervalo de observación# e inicio de Bisección
La = 0.01
Lb = 0.50
tolera = 0.0001 #grafica
muestras = 21
# PROCEDIMIENTO# Método de Bisección
a = La
b = Lb
c = (a+b)/2
tramo = np.abs(b-a)
while (tramo>tolera):
fa = fi(a)
fb = fi(b)
fc = fi(c)
cambio = np.sign(fc)*np.sign(fa)
if (cambio>0):
a = c
b = b
else:
b = c
a = a
c = (a+b)/2
tramo = np.abs(b-a)
# Para la gráfica
tasa = np.linspace(La,Lb,muestras)
fr = fi(tasa)
# SALIDAprint('a, f(a):', a,fa)
print('c, f(c):', c,fc)
print('b, f(b):', b,fb)
print('la raiz esta entre: \n',a,b)
print('con un error de: ', tramo)
print('raiz es tasa buscada: ', c)
print('tasas anual buscada: ',c*12)
# Gráfica
plt.plot(tasa,fr)
plt.axhline(0, color='green')
plt.title('tasa de interes mensual')
plt.show()
la ejecución del algoritmo da como resultado
>>>
RESTART: D:/MATG1052Ejemplos/HipotecaInteres.py
a, f(a): 0.016998291015625 -385.52828922150366
c, f(c): 0.0170281982421875 -145.85350695741363
b, f(b): 0.01705810546875 92.28034212642524
la raiz esta entre:
0.016998291015625 0.01705810546875
con un error de: 5.981445312500111e-05
raiz es tasa buscada: 0.0170281982421875
tasas anual buscada: 0.20433837890625
El tema de integración para primera evaluación se realiza de forma analítica.
Como introducción a la Unidad 7 de la 2da Evaluación, se usa el método del trapecio para el polinomio encontrado en el literal anterior:
v(t) = -0.45125*t**2 + 8.7475*t
distancia recorrida: 202.8912640000001
Realice la integración analítica del polinomio encontrado y determine el error entre los métodos.
Desarrollo con Python
las instrucciones en Python para el ejercicio son:
# 1ra Evaluación II Término 2018# Tema 1. Interpolar velocidad del paracaidistaimport numpy as np
import matplotlib.pyplot as plt
import sympy as sym
# Literal a)definterpola_lagrange(xi,yi):
'''
Interpolación con método de Lagrange
resultado: polinomio en forma simbólica
'''# PROCEDIMIENTO
n = len(xi)
x = sym.Symbol('x')
# Polinomio
polinomio = 0
for i inrange(0,n,1):
# Termino de Lagrange
termino = 1
for j inrange(0,n,1):
if (j!=i):
termino = termino*(x-xi[j])/(xi[i]-xi[j])
polinomio = polinomio + termino*yi[i]
# Expande el polinomio
polinomio = polinomio.expand()
return(polinomio)
# INGRESO
t = [0.0, 2, 4, 6, 8]
v = [0.0, 16.40, 27.77, 35.64, 41.10]
xi = [t[0],t[2],t[4]]
yi = [v[0],v[2],v[4]]
muestras = 51
# PROCEDIMIENTO
polinomio = interpola_lagrange(xi,yi)
velocidad = polinomio.subs('x','t')
# Para graficar
vt = sym.lambdify('t',velocidad)
a = t[0]
b = t[-1]
ti = np.linspace(a, b, muestras)
vi = vt(ti)
# SALIDAprint('v(t) = ', velocidad)
# Grafica
plt.plot(t,v,'ro')
plt.plot(ti,vi)
plt.title('Interpolar velocidad de paracidista')
plt.xlabel('t')
plt.ylabel('v')
plt.show()
# Literal bdefintegratrapecio(funcionx,a,b,tramos):
h = (b-a)/tramos
x = a
suma = funcionx(x)
for i inrange(0,tramos-1,1):
x = x+h
suma = suma + 2*funcionx(x)
suma = suma + funcionx(b)
area = h*(suma/2)
return(area)
# INGRESO# El ingreso es el polinomio en forma lambda# se mantienen las muestras
tramos = muestras-1
# PROCEDIMIENTO
distancia = integratrapecio(vt,a,b,tramos)
# SALIDAprint('distancia recorrida: ', distancia)
Se requiere calcular una norma de suma de filas. es suficiente para demostrar el conocimiento del concepto el usar A.
Se adjunta el cálculo del número de condición y la solución al sistema de ecuaciones:
suma de columnas: [3. 4.75 7.51]
norma A: 7.51
numero de condicion: 97.03737354737129
solucion:
[ 2.03272727 -0.90787879 0.71515152]
El comentario importante corresponde al número de condición, que es un número muy alto para usar un método iterativo, por lo que la solución debe ser un método directo.
Se puede estimar será un número mucho mayor que 1, pues la matriz no es diagonal dominante.
Instrucciones en Python
# 1ra Evaluación II Término 2018# Tema 3. Interpolar con sistema de ecuacionesimport numpy as np
import matplotlib.pyplot as plt
# --------------------------# forma matricial para interpolar
A = np.array([[1, 1. , 1. ],
[1, 1.5, 2.25],
[1, 2.1, 4.41]])
B = np.array([1.84, 2.28, 3.28])
# literal b
sumacolumnas = np.sum(A, axis =1)
norma = np.max(sumacolumnas)
print('suma de columnas: ', sumacolumnas)
print('norma A: ', norma)
numerocondicion = np.linalg.cond(A)
print('numero de condicion: ', numerocondicion)
solucion = np.linalg.solve(A,B)
print('solucion: ')
print(solucion)
Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]
Al analizar las distancias de ex y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el rango de x = [a,b] para distancias se encuentra en:
a > 0, a = 0.1
b < 1, b = 0.7
El ejercicio usa la fórmula de distancia entre dos puntos:
d = \sqrt{(x_2-x_1)^2+(y_2- y_1)^2}
en los cuales:
[x1,y1] = [1,1]
[x2,y2] = [x, ex]
que al sustituir en la fórmula se convierte en:
d = \sqrt{(x-1)^2+(e^x- 1)^2}
que es lo requerido en el literal a
Literal b
Para encontrar el punto más cercano, se debe encontrar el mínimo de la distancia, se podría derivar la función y encontrar la raiz en cero.
Considere simplificar la función a un polinomio, donde tiene dos opciones:
b.1 Polinomio de Taylor, que también requiere derivadas (descartado)
b.2 evaluar la función en varios puntos, interpolar y obtener un polinomio.
Al reutilizar el algoritmo del tema 1 se obtiene lo planteado en b.2, usando un polinomio de grado 3, con muestras de 4 puntos equidistantes en el eje x.
En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .
c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático
Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.
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.0],
[40],
[70],
[50]]
reemplazando los valores iniciales en cada ecuación sin cambios.
