g(x) es creciente en todo el intervalo, con valor minimo en g(1) = 2, y máximo en g(3) =2.449. Por observación de la gráfica, la pendiente g(x) es menor que la recta identidad en todo el intervalo
Verifique la cota de g'(x)
g(x)=3+x=(3+x)1/2g′(x)=21(3+x)−1/2g′(x)=23+x1
Tarea: verificar que g'(x) es menor que 1 en todo el intervalo.
Literal c
Usando el algoritmo del punto fijo, iniciando con el punto x0=2
y tolerancia de 10-4, se tiene que:
Iteración 1: x0=2
g(x0)=3+2=2.2361
error = |2.2361 – 2| = 0.2361
Iteración 2: x1 = 2.2361
g(x1)=3+2.2361=2.2882
error = |2.2882 – 2.2361| = 0.0522
Iteración 3: x2 = 2.2882
g(x2)=3+2.2882=2.2996
error = |2.2996 – 2.28821| = 0.0114
Iteración 4: x3 = 2.2996
g(x3)=3+2.2996=2.3021
error = |2.3021- 2.2996| = 0.0025
Iteración 5: x4 = 2.3021
g(x4)=3+2.3021=2.3026
error = |2.3021- 2.2996| = 5.3672e-04
con lo que determina que el error en la 5ta iteración es de 5.3672e-04 y el error se reduce en casi 1/4 entre iteraciones. El punto fijo converge a 2.3028
El ejercicio se puede simplificar con una matriz de 3×3 dado que una de las corrientes I2 es conocida con valor -10, queda resolver el problema para
[I1 ,I3 ,I4 ]
A = [[ 55.0, 0, -25],
[ 0 , -37, -4],
[-25 , -4, 29]]
B = [-200,-250,100]
sin embargo, para verificar la respuesta se aplica A.X=B
verificar que A.X = B, obteniendo nuevamente el vector B.
[-200. -250. 100.]]
literal b
La norma de la matriz infinito se determina como:
∣∣x∣∣=max[∣xi∣]
considere que en el problema el término en A de magnitud mayor es 55.
El vector suma de filas es:
[[| 55|+| 0|+|-25|], [[80],
[| 0|+|-37|+| -4|], = [41],
[[-25|+| -4|+| 29|]] [58]]
por lo que la norma ∞ ejemplo ||A||∞
es el maximo de suma de filas: 80
para revisar la estabilidad de la solución, se observa el número de condición
>>> np.linalg.cond(A)
4.997509004325602
En éste caso no está muy alejado de 1. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.
Tarea: Matriz de transición de Jacobi
Literal c
En el método de Gauss-Seidel acorde a lo indicado, se inicia con el vector cero. Como no se indica el valor de tolerancia para el error, se considera tolera = 0.0001
el vector solución X es:
[[ -2.7277672 ]
[-10. ]
[ 6.54065815]
[ 1.99891216]]
sin embargo, para verificar la respuesta se aplica A.X=B.
[[-200.]
[-250.]
[ 100.]
[ -10.]]
Se revisa el número de condición de la matriz:
>>> np.linalg.cond(A)
70.21827416891405
Y para éste caso, el número de condición se encuentra alejado del valor 1, contrario a la respuesta del la primera forma de solución con la matriz 3×3. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.
Algoritmo en Python
Presentado por partes para revisión:
Para el método de Gauss, los resultados del algoritmo se muestran como:
Instrucciones en Python usando las funciones creadas en la unidad:
# 1Eva_IIT2019_T3 Circuito eléctrico# Método de Gauss# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport numpy as np
defpivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
if vertabla==True:
print('Matriz aumentada')
print(AB)
print('Pivoteo parcial:')
# Pivoteo por filas AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print(AB)
return(AB)
defgauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
''' Gauss elimina hacia adelante
tarea: verificar términos cero
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
if vertabla==True:
print('Elimina hacia adelante:')
for i inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
if vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(adelante,n,1):
if (np.abs(pivote)>=casicero):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - factor*AB[i,:]
if vertabla==True:
print(' factor: ',factor,' para fila: ',k)
else:
print(' pivote:', pivote,'en fila:',i,
'genera division para cero')
if vertabla==True:
print(AB)
return(AB)
defgauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
''' Gauss sustituye hacia atras
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Sustitución hacia atras
X = np.zeros(n,dtype=float)
ultfila = n-1
ultcolumna = m-1
for i inrange(ultfila,0-1,-1):
suma = 0
for j inrange(i+1,ultcolumna,1):
suma = suma + AB[i,j]*X[j]
X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
return(X)
# INGRESO
A = [[ 55.0, 0, -25],
[ 0 , -37, -4],
[-25 , -4, 29]]
B = [-200,-250,100]
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
AB = gauss_eliminaAdelante(AB,vertabla=True)
X = gauss_sustituyeAtras(AB,vertabla=True)
# SALIDAprint('solución: ')
print(X)
# 1Eva_IIT2019_T3 Circuito eléctrico# 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)
# INGRESO
A = [[ 55.0, 0, -25],
[ 0 , -37, -4],
[-25 , -4, 29]]
B = [-200,-250,100]
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)
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.
p3(x)=12.9−0.15x−0.00390625x2+0.00009765625x3
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_dfinitasAvz(xi,fi, vertabla=False,
precision=6, casicero = 1e-15):
'''Interpolación de diferencias finitas
resultado: polinomio en forma simbólica,
redondear a cero si es menor que casicero
'''
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
n = len(xi)
# revisa tamaños de paso equidistantes
h_iguales = pasosEquidistantes(xi, casicero)
if vertabla==True:
np.set_printoptions(precision)
# POLINOMIO con diferencias Finitas avanzadas
x = sym.Symbol('x')
polisimple = sym.S.Zero # expresión del polinomio con Sympyif h_iguales==True:
tabla,titulo = dif_finitas(xi,fi,vertabla)
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
if vertabla==True:
print('dfinita: ',dfinita)
print(fi[0],'+')
n = len(dfinita)
polinomio = fi[0]
for j inrange(1,n,1):
denominador = np.math.factorial(j)*(h**j)
factor = np.around(dfinita[j-1]/denominador,precision)
termino = 1
for k inrange(0,j,1):
termino = termino*(x-xi[k])
if vertabla==True:
txt1='';txt2=''if n<=2 or j<=1:
txt1 = '('; txt2 = ')'print('+(',np.around(dfinita[j-1],precision),
'/',np.around(denominador,precision),
')*',txt1,termino,txt2)
polinomio = polinomio + termino*factor
# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()
if vertabla==True:
print('polinomio simplificado')
print(polisimple)
return(polisimple)
defdif_finitas(xi,fi, vertabla=False):
'''Genera la tabla de diferencias finitas
resultado en: [título,tabla]
Tarea: verificar tamaño de vectores
'''
xi = np.array(xi,dtype=float)
fi = np.array(fi,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
if vertabla==True:
print('tabla de diferencias finitas')
print(titulo)
print(tabla)
return(tabla, titulo)
defpasosEquidistantes(xi, casicero = 1e-15):
''' Revisa tamaños de paso h en vector xi.
True: h son equidistantes,
False: h tiene tamaño de paso diferentes y dónde.
'''
xi = np.array(xi,dtype=float)
n = len(xi)
# revisa tamaños de paso equidistantes
h_iguales = Trueif n>3:
dx = np.zeros(n,dtype=float)
for i inrange(0,n-1,1): # calcula hi como dx
dx[i] = xi[i+1]-xi[i]
for i inrange(0,n-2,1): # revisa diferencias
dx[i] = dx[i+1]-dx[i]
if dx[i]<=casicero: # redondea cero
dx[i]=0
ifabs(dx[i])>0:
h_iguales=Falseprint('tamaños de paso diferentes en i:',i+1,',',i+2)
dx[n-2]=0
return(h_iguales)
# PROGRAMA ----------------# 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
polinomio = interpola_dfinitasAvz(xi,fi, vertabla=True)
p15 = polinomio.subs(x,15)
# literal b
derivap = polinomio.diff(x,1)
dp16 = derivap.subs(x,16)
px = sym.lambdify(x,polinomio)
xk = np.linspace(np.min(xi),np.max(xi))
pk = px(xk)
# SALIDAprint('Literal a')
print(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(derivap)
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 ------------defbiseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
'''
Algoritmo de Bisección
Los valores de [a,b] son seleccionados
desde la gráfica de la función
error = tolera
'''
fa = fx(a)
fb = fx(b)
tramo = np.abs(b-a)
itera = 0
cambia = np.sign(fa)*np.sign(fb)
if cambia<0: # existe cambio de signo f(a) vs f(b)if vertabla==True:
print('método de Bisección')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
np.set_printoptions(precision)
while (tramo>=tolera and itera<=iteramax):
c = (a+b)/2
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if vertabla==True:
print(itera,[a,c,b],np.array([fa,fc,fb]))
if (cambia<0):
b = c
fb = fc
else:
a = c
fa = fc
tramo = np.abs(b-a)
if vertabla==True:
print(' ',tramo)
itera = itera + 1
respuesta = c
# Valida respuestaif (itera>=iteramax):
respuesta = np.nan
else:
print(' No existe cambio de signo entre f(a) y f(b)')
print(' f(a) =',fa,', f(b) =',fb)
respuesta=np.nan
return(respuesta)
# 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
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDAprint('raíz en: ', respuesta)
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(i1−(1+i)−n)70000=1200(i1−(1+i)−300)
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)
70000−1200(i1−(1+i)−300)=0
Para evitar inconvenientes con la división para cero en caso que i tome el valor de cero, dado se multiplica toda la ecuación por 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
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 estImado 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
con error del orden de O(h5) que al considerar h=2 no permite hacer una buena estimación del error. Sin embargo la respuesta es bastante cercana si se usa el método el trapecio con el algoritmo:
valores de fi: [ 0. 27.77 41.1 ]
divisores en L(i): [ 32. -16. 32.]
Polinomio de Lagrange, expresiones
-1.735625*x*(x - 8.0) + 1.284375*x*(x - 4.0)
Polinomio de Lagrange:
-0.45125*x**2 + 8.7475*x
Método del trapecio
distancia recorrida: 193.28
>>>
El error entre los métodos es |203.2-193.28|= 9.92
Revisar el resultado usando un método con mayor precisión que el trapecio.
Algoritmo 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 sympy as sym
import matplotlib.pyplot as plt
# Literal a)# Interpolacion de Lagrange# divisoresL solo para mostrar valores# INGRESO
t = [0.0, 2, 4, 6, 8]
v = [0.0, 16.40, 27.77, 35.64, 41.10]
cuales = [0,2,4]
# PROCEDIMIENTO
xi = np.array(t,dtype=float)
fi = np.array(v,dtype=float)
xi = xi[cuales]
fi = fi[cuales]
# 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 = 51
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(t,v,'o', label = 'Puntos')
plt.plot(xi,fi,'o', label = 'Puntos en polinomio')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.grid()
plt.show()
# Literal b# INGRESO# El ingreso es el polinomio en forma lambda# se mantienen las muestras# intervalo de integración# a, b seleccionados para la gráfica anterior
tramos = muestras -1
# PROCEDIMIENTOdefintegratrapecio_fi(xi,fi):
''' sobre muestras de fi para cada xi
integral con método de trapecio
'''
n = len(xi)
suma = 0
for i inrange(0,n-1,1):
dx = xi[i+1]-xi[i]
untrapecio = dx*(fi[i+1]+fi[i])/2
suma = suma + untrapecio
return(suma)
tramos = muestras-1
# PROCEDIMIENTO
distancia = integratrapecio_fi(xi,fi)
# SALIDAprint('Método del trapecio')
print('distancia recorrida: ', distancia)
Observación: En la gráfica se muestra que el polinomio pasa por los puntos seleccionados de la tabla. En los otros puntos hay un error que se puede calcular como la resta del punto y su valor con p(x). Queda como tarea.
Usando el algoritmo del polinomio de interpolación con la matriz de Vandermonde se obtiene:
Matriz Vandermonde:
[[1. 1. 1. ]
[2.25 1.5 1. ]
[4.41 2.1 1. ]]
los coeficientes del polinomio:
[ 0.71515152 -0.90787879 2.03272727]
Polinomio de interpolación:
0.715151515151516*x**2 - 0.907878787878792*x + 2.03272727272728
formato pprint
2
0.715151515151516*x - 0.907878787878792*x + 2.03272727272728
suma de columnas: [3. 4.75 7.51]
norma D: 7.51
numero de condicion: 97.03737354737122
solucion:
[ 0.71515152 -0.90787879 2.03272727]
>>>
Literal b
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 condición: 97.03737354737129
solución:
[ 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
# 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones# El polinomio de interpolaciónimport numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO
xj = [1.0, 1.1, 1.3, 1.5, 1.9, 2.1 ]
yj = [1.84, 1.90, 2.10, 2.28, 2.91, 3.28]
cuales = [0, 3, 5]
# muestras = tramos+1
muestras = 51
# PROCEDIMIENTO# Convierte a arreglos numpy
xi = np.array(xj,dtype=float)
fi = np.array(yj,dtype=float)
xi = xi[cuales]
fi = fi[cuales]
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
# Polinomio a forma Lambda# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)
# Para graficar el polinomio en [a,b]
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])# 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)
# Grafica
plt.plot(xj,yj,'o', label='Puntos')
plt.plot(xi,B,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()
# literal b
sumacolumnas = np.sum(D, axis =1)
norma = np.max(sumacolumnas)
print('suma de columnas: ', sumacolumnas)
print('norma D: ', norma)
numerocondicion = np.linalg.cond(D)
print('numero de condicion: ', numerocondicion)
solucion = np.linalg.solve(D,B)
print('solucion: ')
print(solucion)
