Ejercicio: 1Eva_IIT2018_T4 Tasa de interés en hipoteca

literal a

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)

$70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) = 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:

$i \Big[70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big) \Big]= i (0)$ $70000 i - 1200 (1-(1+i)^{-300}) = 0$

La ecuación es la utilizada en el algoritmo de búsqueda de raíces pueden ser:

$fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)$ $fx(i) = 70000i - 1200(1-(1+i)^{-300})$

---

literal b

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

– c.2 una alterna con el método de la Bisección.

---

c.1. Desarrollo del ejercicio con el método del enunciado Newton-Raphson

---

Para el método de Newton-Raphson se tiene que:

$x_{i+1} = x_{i} - \frac{f(x_0i)}{f'(x_i)}$

Se requiere la derivada de la función planteada en el literal a:

$fx(i) = 70000i - 1200(1-(1+i)^{-300})$ $f'x(i) = 70000 + 1200(300)(1+i)^{-301})$

tomando como valor inicial xi = 0.16/12 ≈ 0.013

Se realizan las iteraciones suponiendo que tolera = 1×10^-4

iteración 1

$fx(0.013) = 70000(0.013) - 1200(1-(1+0.013)^{-300})$

 = -265.09

$f'x(0.013) = 70000 + 1200(300)(1+0.013)^{-301})$

= 62623.3

$x_{2} = 0.013 - \frac{-265.09}{62623.34} = 0.017233$

error = |0.013 – 0.01723| = 0.004331

iteración 2

$fx(0.01723) = 70000i - 1200(1-(1+0.0.01723)^{-300})$

= 13.446

$f'x(0.01723) = 70000 + 1200(300)(1+0.01723)^{-301}$

= 67897.5

$x_{3} = 0.017233 - \frac{13.446}{67897.5} = 0.017031$

error = |0.017233 – 0.017031| = 0.000198

cuyo valor 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

El resultado con el algoritmo es:

método de Newton-Raphson
i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 1.30000e-02 -2.65091e+02  6.26233e+04  1.72331e-02  4.23311e-03]
1 [1.72331e-02 1.34468e+01 6.78975e+04 1.70351e-02 1.98045e-04]
2 [1.70351e-02 1.24433e-02 6.77706e+04 1.70349e-02 1.83609e-07]
raiz encontrada en:  0.017034880749732726
tasa anual:  0.20441856899679273

Instrucciones en Python

# 1ra Evaluación II Término 2018
# Tema 4. Tasa de interes para hipoteca
import numpy as np

def newton_raphson(fx,dfx,xi, tolera, iteramax=100,
                   vertabla=False, precision=4):
    '''fx y dfx en forma numérica lambda
    xi es el punto inicial de búsqueda
    '''
    itera=0
    tramo = abs(2*tolera)
    if vertabla==True:
        print('método de Newton-Raphson')
        print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
        np.set_printoptions(precision)
    while (tramo>=tolera):
        fi = fx(xi)
        dfi = dfx(xi)
        xnuevo = xi - fi/dfi
        tramo = abs(xnuevo-xi)
        if vertabla==True:
            print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
        xi = xnuevo
        itera = itera + 1

    if itera>=iteramax:
        xi = np.nan
        print('itera: ',itera,
              'No converge,se alcanzó el máximo de iteraciones')

    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)

x0 = 0.013 # 0.16/12
tolera = 0.0001

# PROCEDIMIENTO
raiz   = newton_raphson(fx, dfx, x0, tolera, vertabla=True, precision=5)
tanual = 12*raiz

# SALIDA
print('raiz encontrada en: ', raiz)
print('tasa anual: ',tanual)

# GRAFICA
import matplotlib.pyplot as plt
a = 0.01/12
b = 0.25/12
muestras = 21

tasa = np.linspace(a,b,muestras)
fi   = fx(tasa)

plt.plot(tasa*12,fi, label="tasa anual")
plt.axhline(0, color='green')
plt.title('tasa anual de interes para Hipoteca')
plt.xlabel('tasa')
plt.ylabel('fx(tasa)')
plt.grid()
plt.legend()
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.

$fx(i) = 70000 - 1200\Big(\frac{1-(1+i)^{-300}}{i} \Big)$

iteración 1

$a = 0.01, b = 0.5$ $c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255$ $fx(0.01) = 70000 - 1200\Big(\frac{1-(1+(0.01))^{-300}}{0.01} \Big) = -43935.86$ $fx(0.255) = 70000 - 1200\Big(\frac{1-(1+(0.255))^{-300}}{0.255} \Big) = 65294.11$ $fx(0.5) = 70000 - 1200\Big(\frac{1-(1+(0.5))^{-300}}{0.5} \Big) = 67600.0$ $tramo = 0.5-0.01 =0.49$

cambio de signo a la izquierda

$a = 0.01, b=0.255$

iteración 2

$c = \frac{a+b}{2} = \frac{0.01+0.225}{2} = 0.1325$ $fx(0.1325) = 70000 - 1200\Big(\frac{1-(1+(0.1325))^{-300}}{0.1325} \Big) = 60943.39$ $tramo = 0.225-0.01 =0.215$

cambio de signo a la izquierda

$a = 0.01, b=0.1325$

iteración 3

$c = \frac{a+b}{2} = \frac{0.01+0.1325}{2} = 0.07125$ $fx(0.07125) = 70000 - 1200\Big(\frac{1-(1+(0.07125))^{-300}}{0.07125} \Big) = 53157.89$ $tramo = 0.1325-0.01 =0.1225$

cambio de signo a la izquierda

$a = 0.01, b=0.07125$

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 hipoteca
import 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)

# SALIDA
print('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

y la gráfica obtenida es:

Ejercicio: 1Eva_IIT2018_T1 Interpolar velocidad del paracaidista

El ejercicio tiene dos partes: la interpolación y el integral.

Literal a

t	[s]	0	2	4	6	8
v(t)	[m/s]	0.0	16.40	27.77	35.64	41.10

No se especifica el método a seguir, por lo que se puede seleccionar el de mayor preferencia.

Por ejemplo. usando el método de Lagrange, con los puntos primero, medio y último, para cubrir todo el intervalo:

$p_2(t) = 0\frac{(t-4)(t-8)}{(0-4)(0-8)} +$ $+ 27.77\frac{(t-0)(t-8)}{(4-0)(4-8)} +$ $+ 41.10\frac{(t-0)(t-4)}{(8-0)(8-4)}$ $p_2(t) = 0 + 27.77\frac{t(t-8)}{-16}) +$ $+ 41.10\frac{t(t-4)}{32}$ $p_2(t) = -1.73(t^2-8t) + 1.28(t^2-4t)$ $p_2(t) = -0.45 t^2 + 8.74t$

---

Literal b

El tema de integración para primera evaluación se realiza de forma analítica.

Una de las formas, que es independiente si se resolvió el literal a, es usar los datos proporcionados en la tabla el ejercicio:

t	[s]	0	2	4	6	8
v(t)	[m/s]	0.0	16.40	27.77	35.64	41.10

Se podría usar el método de Simpson de 1/3, puesto que los tamaños de paso en t son equidistantes se puede aplicar: h=2-0=2

$\int_0^8 v(t)dt = \frac{2}{3}\Big( 0+ 4(16.40)+27.77\Big)$ $+ \frac{2}{3}\Big( 27.77+ 4(35.64)+41.10\Big)$ $=203.2$

con error del orden de O(h⁵) 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 paracaidista
import 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 in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(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)

# SALIDA
print('    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

# PROCEDIMIENTO
def integratrapecio_fi(xi,fi):
    ''' sobre muestras de fi para cada xi
        integral con método de trapecio
    '''
    n = len(xi)
    suma = 0
    for i in range(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)

# SALIDA
print('Método del trapecio')
print('distancia recorrida: ', distancia)

Ejercicio: 1Eva_IIT2018_T3 Interpolar con sistema de ecuaciones

Los datos del ejercicio proporcionados son:

i	0	1	2	3	4	5
x	1.0	1.1	1.3	1.5	1.9	2.1
y(x)	1.84	1.90	2.10	2.28	2.91	3.28

Literal a

El tema es semejante al tema 1, cambiando el método de interpolación.
Se usan los puntos de las posiciones 0, 3 y 5.

$p_2(x) = b_0 + b_1x + b_2 x^2$

en la fórmula:

punto x[0] = 1, y[0]= 1.84

$1.84 = b_0 + b_1(1) + b_2 (1)^2$ $1.84 = b_0 + b_1 + b_2$

punto x[3] = 1.5, y[3]= 2.28

$2.28 = b_0 + b_1(1.5) + b_2 (1.5)^2$ $2.28 = b_0 + 1.5 b_1 + 2.25 b_2$

punto x[5] = 2.1, y[5]= 3.28

$3.28= b_0 + b_1(2.1) + b_2 (2.1)^2$ $3.28= b_0 + 2.1 b_1 + 4.41 b_2$

se obtiene el sistema de ecuaciones:

$b_0 + b_1 + b_2 = 1.84$ $b_0 + 1.5 b_1 + 2.25 b_2 = 2.28$ $b_0 + 2.1 b_1 + 4.41 b_2 = 3.28$

Con lo que se plantea la forma Ax=B:

$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & 1.5 & 2.25 \\1 & 2.1 & 4.41 \end{bmatrix}$ $B = \begin{bmatrix} 1.84\\ 2.28 \\ 3.28 \end{bmatrix}$

Matriz Aumentada

$AB = \begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 1 & 1.5 & 2.25 & 2.28 \\1 & 2.1 & 4.41 &3.28 \end{bmatrix}$

Pivoteo parcial por filas

Para el primer pivote no se requieren cambio de filas.
para el segundo pivote de la diagonal se deben intercambiar la fila segunda con la tercera

$\begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 1 & 2.1 & 4.41 &3.28 \\ 1 & 1.5 & 2.25 & 2.28 \end{bmatrix}$

Se aplica eliminación hacia adelante:

fila = 0, columna=0  pivote = AB[0,0]=1

factor entre las filas es 1/1=1.

$\begin{bmatrix}1 & 1 & 1 & 1.84 \\ 1-1 & 2.1-1 & 4.41 -1 &3.28 -1.84 \\ 1-1 & 1.5 -1 & 2.25 -1 & 2.28 - 1.84 \end{bmatrix}$ $\begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 & 1.25 & 0.44 \end{bmatrix}$

fila =1,  columna=1, pivote=AB[1,1] =1.1

factor entre filas es 0.5/1.1 = 1/2.2

$\begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0.5 -\frac{0.5}{1.1}(1.1)& 1.25 -\frac{0.5}{1.1}(3.41)& 0.44-\frac{0.5}{1.1}(1.44) \end{bmatrix}$ $\begin{bmatrix} 1 & 1 & 1 & 1.84 \\ 0 & 1.1 & 3.41 &1.44 \\ 0 & 0 & -0.3 & -0.214545 \end{bmatrix}$

aplicando sustitución hacia atrás

$b2 = -0.21/(-0.3) = 0.71515$ $b1= \frac{1.44-3.41 b_2}{1.1} = \frac{1.44-3.41( 0.71515)}{1.1}=-0.9078$ $b3= \frac{1.84-b_1-b_2}{1} = \frac{1.84-(-0.9078)-(0.71515)}{1} =2.0327$

con lo que el polinomio buscado es:

$p_2(x) = 2.0327 -0.9078 x + 0.71515 x^2$

y se obtiene el resultado de la interpolación.

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ón
import 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 in range(0,n,1):
    for j in range(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 in range(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])
    
# SALIDA
print('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)

.

Ejercicio: 1Eva_IIT2018_T2 Distancia mínima a un punto

Literal a

Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]

Al analizar las distancias de e^x y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el intervalo 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:

[x₁,y₁] = [1,1]
[x₂,y₂] = [x, e^x]

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 usar un método de búsqueda de raíces, se requiere encontrar el valor cuando f(x) = d’ = 0.

Un método como el de Newton Raphson requiere también f'(x) = d''

$f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}}$ $f'(x)= \frac{(e^x - 1)e^x + e^{2x} + 1 - \frac{(x + (e^x - 1)e^x - 1)^2}{(x - 1)^2 + (e^x - 1)^2}} {\sqrt{(x - 1)^2 + (e^x - 1)^2}}$

expresiones obtenidas usando Sympy

f(x) :
(x + (exp(x) - 1)*exp(x) - 1)/sqrt((x - 1)**2 + (exp(x) - 1)**2)
f'(x) :
((exp(x) - 1)*exp(x) + exp(2*x) + 1 - (x + (exp(x) - 1)*exp(x) - 1)**2/((x - 1)**2 + (exp(x) - 1)**2))/sqrt((x - 1)**2 + (exp(x) - 1)**2)

f(x) :
       / x    \  x        
   x + \e  - 1/*e  - 1    
--------------------------
    ______________________
   /                    2 
  /         2   / x    \  
\/   (x - 1)  + \e  - 1/  
f'(x) :
                                              2
                         /    / x    \  x    \ 
/ x    \  x    2*x       \x + \e  - 1/*e  - 1/ 
\e  - 1/*e  + e    + 1 - ----------------------
                                             2 
                                 2   / x    \  
                          (x - 1)  + \e  - 1/  
-----------------------------------------------
               ______________________          
              /                    2           
             /         2   / x    \            
           \/   (x - 1)  + \e  - 1/            

lo que permite observar la raíz de f(x) en una gráfica:

con las siguientes instrucciones:

# Eva_IIT2018_T2 Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
x = sym.Symbol('x')
fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)

a = 0
b = 1
muestras = 21

# PROCEDIMIENTO
dfx = sym.diff(fx,x,1)
d2fx = sym.diff(fx,x,2)

f = sym.lambdify(x,dfx)
xi = np.linspace(a,b,muestras)
fi = f(xi)


# SALIDA
print('f(x) :')
print(dfx)
print("f'(x) :")
print(d2fx)
print()
print('f(x) :')
sym.pprint(dfx)
print("f'(x) :")
sym.pprint(d2fx)

# GRAFICA
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid()
plt.show()

Usando el método de la bisección para el intervalo dado, se tiene:

$f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}}$

itera = 0 , a = 0, b=1

$c= \frac{0+1}{2} = 0.5$ $f(0) = \frac{0 + (e^0 - 1)e^0 - 1}{\sqrt{(0 - 1)^2 + (e^0 - 1)^2}} = -1$ $f(1) = \frac{1 + (e^1 - 1)e^1 - 1}{\sqrt{(1 - 1)^2 + (e^1 - 1)^2}} 2.7183$ $f(0.5) = \frac{(0.5) + (e^(0.5) - 1)e^(0.5) - 1}{\sqrt{((0.5) - 1)^2 + (e^(0.5) - 1)^2}} = 0.6954$

cambio de signo a la izquierda,

a= 0, b=c=0.5

tramo = |0.5-0| = 0.5

itera = 1

$c= \frac{0+0.5}{2} = 0.25$ $f(0.25) = \frac{(0.25) + (e^(0.25) - 1)e^(0.25) - 1}{\sqrt{((0.25) - 1)^2 + (e^(0.25) - 1)^2}} = -0.4804$

cambio de signo a la derecha,

a=c= 0.25, b=0.5

itera = 2

$c= \frac{0.25+0.5}{2} = 0.375$ $f(0.375) = \frac{(0.375) + (e^(0.375) - 1)e^(0.375) - 1}{\sqrt{((0.375) - 1)^2 + (e^(0.375) - 1)^2}} = 0.0479$

cambio de signo a la izquierda,

a= 0.25, b=c=0.375

se continúan las iteraciones con el algoritmo, para encontrar la raíz en 0.364:

método de Bisección
i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [0, 0.5, 1] [-1.      0.6954  2.7183]
   0.5
1 [0, 0.25, 0.5] [-1.     -0.4804  0.6954]
   0.25
2 [0.25, 0.375, 0.5] [-0.4804  0.0479  0.6954]
   0.125
3 [0.25, 0.3125, 0.375] [-0.4804 -0.2388  0.0479]
   0.0625
4 [0.3125, 0.34375, 0.375] [-0.2388 -0.1004  0.0479]
   0.03125
5 [0.34375, 0.359375, 0.375] [-0.1004 -0.0274  0.0479]
   0.015625
6 [0.359375, 0.3671875, 0.375] [-0.0274  0.01    0.0479]
   0.0078125
7 [0.359375, 0.36328125, 0.3671875] [-0.0274 -0.0088  0.01  ]
   0.00390625
8 [0.36328125, 0.365234375, 0.3671875] [-0.0088  0.0006  0.01  ]
   0.001953125
9 [0.36328125, 0.3642578125, 0.365234375] [-0.0088 -0.0041  0.0006]
   0.0009765625
raíz en:  0.3642578125

Al algoritmo anterior se complementa con las instrucciones de la función para la bisección.

# Eva_IIT2018_T2 Distancia mínima a un punto
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
x = sym.Symbol('x')
fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)

a = 0
b = 1
muestras = 21

# PROCEDIMIENTO
dfx = sym.diff(fx,x,1)
d2fx = sym.diff(fx,x,2)

f = sym.lambdify(x,dfx)
xi = np.linspace(a,b,muestras)
fi = f(xi)


# SALIDA
print('f(x) :')
print(dfx)
print("f'(x) :")
print(d2fx)
print()
print('f(x) :')
sym.pprint(dfx)
print("f'(x) :")
sym.pprint(d2fx)

# GRAFICA
plt.plot(xi,fi, label='f(x)')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid()
plt.show()

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera
import numpy as np

def biseccion(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 respuesta
        if (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)

# INGRESO
tolera = 0.001

# PROCEDIMIENTO
respuesta = biseccion(f,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

.

Ejercicio: 1Eva_IT2018_T1 Tanque esférico canchas deportivas

a) Para seleccionar el rango para h=[a,b], se observa que el tanque puede estar vacío, a=0 o lleno al máximo, b=2R = 2(3)=6, con lo que se obtiene:

h =[0.0, 6.0]

conociendo la proporción con el valor máximo, se tiene un valor inicial para h₀ para las iteraciones.

$Vmax = \frac{\pi}{3} (2R)^2 (3R-2R)$ $Vmax = \frac{4\pi }{3}R^{3}= 113.09$ $h_0 = (6)*30/113.09 = 1.59$

b) Usar Newton-Raphson con tolerancia 1e-6

$f(h) = \frac{\pi }{3}h^2 (3(3)-h)-30$ $f(h) = \frac{\pi }{3}(9h^2 -h^3-28.647)$ $f'(h) = \frac{\pi }{3}(18h-3h^2)$

el punto siguiente de iteración es:

$h_{i+1} = h_{i} -\frac{f(h)}{f'(h)}$ $= h_{i}-\frac{ \frac{\pi }{3}(9h^2 -h^3-28.647)}{ \frac{\pi }{3}(18h-3h^2)}$ $h_{i+1} = h_{i} -\frac{(9h^2 -h^3-28.647)}{(18h-3h^2)}$

con lo que se realiza la tabla de iteraciones:

hi	hi+1	error	orden
1.590	2.061	0.47	10-1
2.061	2.027	-0.034	10-2
2.027	2.02686	-0.00014	10-4
2.02686	2.0268689	-2.32E-09	10-9

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.

Ejercicio: 1Eva_IT2018_T3 Temperatura en nodos de placa

a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:

$a=\frac{50+c+100+b}{4}$ $b=\frac{a+30+50+d}{4}$ $c=\frac{a+60+100+d}{4}$ $d=\frac{b+60+c+30}{4}$

que reordenando se convierte en:

$4a=150+c+b$ $4b=a+80+d$ $4c=a+160+d$ $4d=b+c+90$

simplificando:

$4a-b-c= 150$ $a-4b+d = -80$ $a-4c+d = -160$ $b+c-4d = -90$

que a forma matricial se convierte en:

A = [[ 4, -1, -1, 0.0],
     [ 1, -4,  0, 1.0],
     [ 1,  0, -4, 1.0],
     [ 0,  1,  1,-4.0]]
B = [[ 150.0],
     [ -80.0],
     [-160.0],
     [ -90.0]]

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.

iteración 1
$a=\frac{50+70+100+40}{4} = 65$

$b=\frac{60+30+50+50}{4} = 47.5$ $c=\frac{60+60+100+50}{4} = 67.5$ $d=\frac{40+60+70+30}{4} = 50$

X(1) = [[65],
        [47.5],
        [67.5],
        [50]]

vector de error = 
     [|65-60|,
      |47.5-40|,
      |67.5-70|,
      |50-50|]
  = [|5|,
     |7.5|,
     |-2.5|,
     |0|]
errormax = 7.5

iteración 2
$a=\frac{50+67.5+100+47.5}{4} = 66.25$

$b=\frac{65+30+50+50}{4} = 48.75$ $c=\frac{65+60+100+50}{4} = 68.75$ $d=\frac{47.5+60+67.7+30}{4} = 51.3$

X(2) = [[66.25],
        [48.75],
        [68.75],
        [51.3]]

vector de error = 
       [|66.25-65|,
        |48.75-47.5|,
        |68.75-67.5|,
           |51.3-50|] 
  = [|1.25|,
     |1.25|,
     |1.25|,
     |1.3|]
errormax = 1.3

iteración 3
$a=\frac{50+68.75+100+48.75}{4} = 66.875$

$b=\frac{66.25+30+50+51.3}{4} = 49.38$ $c=\frac{66.25+60+100+51.3}{4} = 69.3875$ $d=\frac{48.75+60+68.75+30}{4} = 51.875$

X(2) = [[66.875],
        [49.38],
        [69.3875],
        [51.875]]

vector de error = 
      [|66.875-66.25|,
       |49.38-48.75|,
       |69.3875-68.75|,
       |51.875-51.3|]
    = [|0.655|,
       |0,63|,
       |0.6375|,
        |0.575|]
errormax = 0.655
con error relativo de:
100*0.655/66.875 = 0.97%

siguiendo las iteraciones se debería llegar a:

>>> np.linalg.solve(A,B)
array([[ 67.5],
       [ 50. ],
       [ 70. ],
       [ 52.5]])

Ejercicio: 1Eva_IT2018_T4 El gol imposible

Tabla de datos:

ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]

Observe que, un gol simplificado con física básica es un tiro parabólico cuya trayectoria se compone de movimientos en los ejes, Y y Z. Sin embargo, lo «imposible» del gol mostrado es añadir el movimiento en X. (Referencia de la imagen en el enunciado)

El movimiento de «profundidad» o dirección hacia el arco y(t) es semejante a un polinomio de primer grado, y el movimiento de «altura» z(t) es un polinomio de segundo grado. El movimiento «imposible» en el eje X, podría ser descrito por un polinomio de segundo o mayor grado.

a) Encontrar t para altura máxima, que se encuentra al igualar la derivada dz/dt a cero. Para interpolar el polinomio z(t), de segundo grado, se puede usar tres puntos de los sugeridos: 0, 0.3 y 0.6, que además son equidistantes en t (facilita usar cualquier método de interpolación).

Por ejemplo, con diferencias finitas avanzadas:

t	z(t)	d1	d2	d3
0.00	0.00	1.40	-0.89
0.30	1.40	0.51
0.60	1.91

$z(t) = 0 + 1.40\frac{(t-0)}{1!(0.3)} - 0.89 \frac{(t-0)(t-0.3)}{2!(0.3)^2}$ $= 0 + 4.66 t - 4.94(t^2-0.3t)$ $= 4.66 t + 1.48 t - 4.94 t^2$ $z(t) = 6.42 t - 4.94 t^2$

para encontrar el valor máximo se encuentra $\frac{dz}{dt} = 0$

$\frac{dz}{dt} = 6.42 - 2(4.94) t$ $6.42 - 2(4.94) t = 0$ $t = \frac{6.42}{2(4.94)}$ $t = 0.649$

Observe que el resultado tiene sentido, pues según la tabla, el máximo debería estar entre 0.60 y 0.75

b) El cruce por la «barrera», corresponde al desplazamiento del balón en el eje Y a 9 metros: y(t) = 9.
El polinomio modelo de primer grado usa dos puntos para la interpolación, de los sugeridos se escogen dos, por ejemplo: 0.0 y 0.3.

Usando diferencias finitas avanzadas :

$d1 = (8.88-0) = 8.88$ $y(t) = 0 + 8.88\frac{(t-0)}{1!(0.3)}$ $y(t) = 29.6 t$

usando y(t) = 9

29.6 t = 9
t = 0.30
z(0.30) = 1.40
(de la tabla de datos)

cuya respuesta con solo dos dígitos decimales es coherente, al estar cercano el valor a una distancia y=8.88 o aproximado a 9.
La respuesta puede ser mejorada usando más digitos significativos en las operaciones.

c)  La desviación máxima en eje X.
Determine un polinomio para la trayectoria en el eje X y obtenga el máximo igualando la derivada del polinomio x(t) a cero.

Por simplicidad, se usa un polinomio de segundo grado, alrededor de los valores máximos en el eje X

t	x(t)	d1	d2	d3
0.60	1.80	0.20	-0.30
0.75	2.00	-0.10
0.90	1.90

$x(t) = 1.80 + 0.20 \frac{(t-0.60)}{1!(0.15)} -0.30 \frac{(t-0.60)(t-0.75)}{2!(0.15)^2}$ $x(t) = 1.80 + 1.33 (t-0.60) - 6.67(t-0.60)(t-0.75)$

como se busca el máximo, se usa la derivada $\frac{dx}{dt} = 0$

$\frac{dx}{dt} = 1.33 - 6.67(2t +(-0.60-0.75))$ $1.33 - 13.34t + 9.00 = 0$ $-13.34t + 10.33 = 0$

t = 0.77

$x(0.77) = 1.80 + 1.33(0.77-0.60) - 6.67(0.77-0.60)(0.77-0.75)$ $x(0.77) = 2.003$

lo que es coherente con la tabla para el eje x, pues el máximo registrado es 2, y el próximo valor es menor, la curva será decreciente.

---

Desarrollo extra para observar y verificar resultados:

Usando los puntos y un graficador 3D se puede observar la trayectoria:

Tarea: Realice el ejercicio usando los algoritmos en Python, muestre los polinomios obtenidos y compare.

Nota: La gráfica 3D presentada usa interpolación de Lagrange con todos los puntos. Realice las observaciones y recomendaciones necesarias y presente su propuesta como tarea.

Ejercicio: 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

Se dispone de tres puntos para la gráfica.

x	f(x)
0	1
0.2	1.6
0.4	2.0

Si el polinomio de Taylor fuera de grado 0, sería una constante, que si se evalúa en x₀ = 0 para eliminar los otros términos, se encuentra que sería igual a 1

Como se pide el polinomio de grado 2, se tiene la forma:

$p(x) = a + bx + c x ^2$ $p(x) = 1 + bx + c x^2$

Se disponen de dos puntos adicionales que se pueden usar para determinar b y c:

$p(0.2) = 1 + 0.2 b + (0.2)^2 c = 1.6$ $p(0.4) = 1 + 0.4 b + (0.4)^2 c = 2.0$

simplificando:

$0.2 b + (0.2)^2 c = 1.6 - 1 = 0.6$ $0.4 b + (0.4)^2 c = 2.0 - 1 = 1$

multiplicando la primera ecuación por 2 y restando la segunda ecuación:

$0 - 0.08 c = 1.2-1 = 0.2$ $c = - 0.2/0.08 = -2.5$

sustituyendo el valor de c obtenido en la primera ecuación

$0.2 b + 0.04(-2.5) = 0.6$ $0.2 b = 0.6 - 0.04(-2.5) = 0.6 + 0.1 = 0.7$ $b = 0.7/0.2 = 3.5$

con lo que el polinomio queda:
$p(x) = 1 + 3.5 x - 2.5 x^2$

validando con python:
tomando los puntos de prueba:

xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]
>>>

se obtiene la gráfica:

se adjunta las instrucciones usadas para validar que el polinomio pase por los puntos requeridos.

# 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]

px = lambda x: 1 + 3.5*x - 2.5*(x**2)
a = -0.5
b = 1
muestras = 21

# PROCEDIMIENTO
xj = np.linspace(a,b,muestras)
pxj = px(xj)

# SALIDA
print(xj)
print(pxj)

# Gráfica
plt.plot(xj,pxj,label='p(x)')
plt.plot(xi,fi,'o', label='datos')

plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()

Nota: Se puede intentar realizar los polinomios aumentando el grado, sin embargo cada término agrega un componente adicional a los términos anteriores por la forma (x – x₀)^k

literal b

se requiere el integral aproximado de f(x) en el intervalo limitado por los 3 puntos de la tabla:

$\int_{0}^{0.4}f(x) dx$

Esta aproximación con un polinomio es el concepto de integración numérica con la regla de Simpson de 1/3, tema desarrollado en la unidad 5

$I_{S13} = \frac{0.2}{3} \Big(1+4(1.6)+ 2 \Big) = 0.62666$

Ejercicio: 1Eva_IT2017_T4 Componentes eléctricos

Desarrollo Analítico

Solo puede usar las cantidades disponibles de material indicadas, por lo que las cantidades desconocidas de producción por componente se convierten en las incógnitas x₀, x₁, x₂. Se usa el índice cero por compatibilidad con las instrucciones Python.

	Material 1	Material 2	Material 3
Componente 1	5 x0	9 x0	3 x0
Componente 2	7 x1	7 x1	16 x1
Componente 3	9 x2	3 x2	4 x2
Total	945	987	1049

Se plantean las fórmulas a resolver:

5 x0 +  7 x1 + 9 x2 = 945
9 x0 +  7 x1 + 3 x2 = 987
3 x0 + 16 x1 + 4 x2 = 1049

Se reescriben en la forma matricial AX=B

$\begin{bmatrix} 5 & 7 & 9\\ 9 & 7 & 3 \\ 3& 16 & 4\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 945 \\ 987 \\ 1049 \end{bmatrix}$

Se reordena, pivoteando por filas para tener la matriz diagonalmente dominante:

$\begin{bmatrix} 9 & 7 & 3\\ 3 & 16 & 4 \\ 5& 7 & 9\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 987 \\ 1049 \\ 945 \end{bmatrix}$

Se determina el número de condición de la matriz, por facilidad con Python:

numero de condicion: 4.396316324708121

Obtenido con:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]])

B = np.array([987.,1049.,945.])
# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('numero de condicion:', ncond)

Recordando que la matriz debe ser tipo real, se añade un punto a los dígitos.

El número de condición es cercano a 1, por lo que el sistema si debería converger a la solución.

Desarrollo con Python

La forma AX=B permite usar los algoritmos desarrollados, obteniendo la solución. Se verifica el resultado al realizar la multiplicación de A con el vector respuesta, debe ser el vector B con un error menor al tolerado.

respuesta con Jacobi
[[62.99996585]
 [44.99998037]
 [34.9999594 ]]
verificando:
[[ 986.99943346]
 [1048.99942111]
 [ 944.99932646]]

Si interpreta el resultado, se debe obtener solo la parte entera [63,45,35] pues las unidades producidas son números enteros.

instrucciones:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np

def jacobi(A,B,tolera,X0,iteramax=100, vertabla=False):
    ''' Método de Jacobi, tolerancia, vector inicial X0
        para mostrar iteraciones: tabla=1
    '''
    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)

    itera = 0
    if vertabla==True:
        print('itera,[X0],errado')
        print(itera, xnuevo, errado)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(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)
        itera = itera + 1
        if vertabla==True:
            print(itera, xnuevo, errado)
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)


# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]],dtype=float)

B = np.array([987.,1049.,945.],dtype=float)
tolera = 1e-4
X0 = [0.,0.,0.]

# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

respuesta = jacobi(A,B,tolera,X0,vertabla=True)

verifica = np.dot(A,respuesta)
verificaErrado = np.max(np.abs(B-np.transpose(verifica)))

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificar A.X:')
print(verifica)
print('max|A.X - B|')
print(verificaErrado)

al ejecutar el algoritmo se determina que se requieren 83 iteraciones para cumplir con con el valor de error tolerado.

Ejercicio: 1Eva_IT2017_T2 Tanque esférico-volumen

a. Planteamiento del problema

$V = \frac{\pi h^{2} (3r-h)}{3}$

Si r=1 m y V=0.75 m³,

$0.75 = \frac{\pi h^{2} (3(1)-h)}{3}$ $0.75 -\frac{\pi h^{2} (3(1)-h)}{3} = 0$ $f(h) = 0.75 -\frac{\pi h^{2} (3-h)}{3}$ $= 0.75 -\frac{\pi}{3}(h^{2} (3)-h^{2}h)$ $= 0.75 -\frac{\pi}{3}(3 h^{2} - h^{3})$ $f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}$

---

b. Intervalo de búsqueda de raíz

El tanque vacio tiene h=0 y completamente lleno h= 2r = 2(1) = 2, por lo que el intevalo tiene como extremos:

[0,2]

Verificando que exista cambio de signo en la función:

$f(0) = 0.75 -\pi (0)^{2} + \frac{\pi}{3}(0)^{3} = 0.75$ $f(2) = 0.75 -\pi (2)^{2} + \frac{\pi}{3}(2)^{3}= -3.4387$

y verificando al usar la gráfica dentro del intervalo:

Tolerancia

Se indica en el enunciado como 0.01 que es la medición mínima a observar con un flexómetro.

tolera = 0.01

---

c. Método de Newton-Raphson
d. Método de Punto Fijo

---

c. Método de Newton-Raphson

El método de Newton-Raphson requiere la derivada de la función:

$x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)}$ $f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}$ $f'(h) = -2\pi h + \pi h^{2}$

Tomando como punto inicial de búsqueda el extremo izquierdo del intervalo, genera una división para cero. Por lo que se mueve un poco a la derecha, algo más cercano a la raiz, viendo la gráfica por ejemplo 0.1

x₀ = 0.1

iteración 1

i =0

$f(0.1) = 0.75 -\pi (0.1)^{2} + \frac{\pi}{3}(0.1)^{3} =0.7196$ $f'(0.1) = -2\pi (0.1) + \pi (0.1)^{2} = -0.5969$ $x_{1} = x_0 -\frac{0.7496 }{-0.0625} = 1.3056$ $tramo = |x_{1}-x_{0}| = |0.1-1.3056 | = 1.2056$

iteración 2

i =1

$f(1.3056) = 0.75 -\pi (1.3056)^{2} + \frac{\pi}{3}(1.3056)^{3} = -2.2746$ $f'(1.3056) = -2\pi (1.3056) + \pi (1.3056)^{2} =-2.8481$ $x_{1} = x_0 -\frac{-2.2746}{-2.8481} = 0.5069$ $tramo = |x_{2}-x_{1}|=|0.5069-1.3056|=0.7987$

iteración 3

i =2

$f(0.5069) = 0.75 -\pi (0.5069)^{2} + \frac{\pi}{3}(0.5069)^{3} = 0.0789$ $f'(0.5069) = -2\pi (0.5069) + \pi (0.5069)^{2} =-2.3780$ $x_{1} = x_0 -\frac{0.0789}{-2.3780} = 0.5401$ $tramo = |x_{3}-x_{2}| =|0.5401 - 0.5069| = 0.0332$

Observe que el tramo disminuye en cada iteración , por lo que el método converge, si se siguen haciendo las operaciones se tiene que:

 [ xi, xnuevo, tramo]
[[1.00000000e-01 1.30560920e+00 1.20560920e+00]
 [1.30560920e+00 5.06991599e-01 7.98617601e-01]
 [5.06991599e-01 5.40192334e-01 3.32007350e-02]
 [5.40192334e-01 5.39518667e-01 6.73667593e-04]]
raiz 0.5395186666699257

Instrucciones en Python

para Método de Newton-Raphson

# 1Eva_IT2017_T2 Tanque esférico-volumen
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
dfx = lambda h: -2*np.pi*h+np.pi*(h**2)

# Parámetros de método
a = 0
b = 2
tolera = 0.01
x0 = 0.1
iteramax = 100

# parametros de gráfica
La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Newton-Raphson
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo
tabla = np.array(tabla)

# calcula para grafica
hi = np.linspace(La,Lb,muestras)
fi = fx(hi)
gi = dfx(hi)

# SALIDA
print(' [ xi, xnuevo, tramo]')
print(tabla)
print('raiz', xnuevo)
plt.plot(hi,fi)
plt.plot(xi,fx(xi),'ro')
plt.axhline(0, color='green')
plt.xlabel('h')
plt.ylabel('V')
plt.title('Método Newton-Raphson')
plt.show()

---

Planteo con Punto Fijo

---

d. Método de Punto Fijo

Del planteamiento del problema en el literal a, se tiene que:

$0.75 = \frac{\pi h^{2} (3(1)-h)}{3}$

de donde se despeja una h:

$\frac{3(0.75)}{\pi (3(1)-h) } = h^{2}$ $h = \sqrt{\frac{3*0.75}{\pi (3-h) }}$ $h = \sqrt{\frac{2.25}{\pi (3-h) }}$

con lo que se obtienen las expresiones a usar en el método

$identidad = h$ $g(h) = \sqrt{\frac{2.25}{\pi (3-h) }}$

El punto inicial de búsqueda debe encontrarse en el intervalo, se toma el mismo valor que x0 en el método de Newton-Raphson

x₀ = 0.10

Iteracion 1

$x_0= 0.10$ $g(0.10) = \sqrt{\frac{2.25}{\pi (3-(0.10) }}= 0.4969$ $tramo= |0.4969-0.10| = 0.3869$

Iteracion 2

$x_1= 0.4969$ $g(0.4969) = \sqrt{\frac{2.25}{\pi (3-(0.4969 ) }}= 0.5349$ $tramo= |0.5349- 0.4969| = 0.038$

Iteracion 3

$x_2 =0.5349$ $g(0.5349) = \sqrt{\frac{2.25}{\pi (3-(0.5349) }}= 0.5390$ $tramo= |0.5390 - 0.5349| = 0.0041$

con lo que se cumple el criterio de tolerancia, y se obtiene la raiz de:

raiz = 0.5390

Tabla de resultados, donde se observa que el tramo o error en cada iteración disminuye, por lo que el método converge.

 [i,xi,xi+1,tramo]
[[1.         0.1        0.4969553  0.3969553 ]
 [2.         0.4969553  0.5349116  0.03795631]
 [3.         0.5349116  0.53901404 0.00410243]]
raiz 0.5390140355891347
>>> 

Instrucciones en Python

para Método de Punto-Fijo, recordamos que el método puede diverger, por lo que se añade el parámetro iteramax

# 1Eva_IT2017_T2 Tanque esférico-volumen
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda h: h
gx = lambda h: np.sqrt(2.25/(np.pi*(3-h)))

a = 0.1
b = 2
tolera = 0.01
iteramax = 100

La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Punto Fijo
tabla = []
i = 1 # iteración
b = gx(a)
tramo = abs(b-a)
while not(tramo<=tolera or i>=iteramax):
    tabla.append([i,a,b,tramo])
    a = b
    b = gx(a)
    tramo = abs(b-a)
    i = i+1
tabla.append([i,a,b,tramo])
respuesta = b

# Validar respuesta
if (i>=iteramax):
    respuesta = np.nan
tabla = np.array(tabla)

# calcula para grafica
hi = np.linspace(La,Lb,muestras)
fi = fx(hi)
gi = gx(hi)

# SALIDA
print(' [i, xi, xi+1, tramo]')
print(tabla)
print('raiz', respuesta)
plt.plot(hi,fi, label = 'identidad', color='green')
plt.plot(hi,gi, label = 'g(h)', color = 'orange')
plt.plot(b,gx(b),'ro')
plt.axhline(0, color='green')
plt.xlabel('h')
plt.ylabel('V')
plt.title('Método Punto Fijo')
plt.legend()
plt.axhline(0, color='green')
plt.show()

Otra forma de probar la convergencia es que |g'(x)|<1 que se observa en la una gráfica adicional, lo que limita aún más el intervalo de búsqueda.

Desarrollo en la siguiente clase.

---