s3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

Ejercicio: 3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

la ecuación a resolver es:

\frac{\delta^2 y(t)}{\delta t^2}+3 \frac{\delta y(t)}{ \delta t}+2 y(t) =0

con valores iniciales: y0(t)=0 , y’0(t) =-5

se puede escribir como:

y"+3 y'+2y = 0 y" = -3y'-2y

sutituyendo las expresiones de las derivadas como las funciones f(x) por z y g(x) por z’:

y’ = z = f(x)

y» = z’= -3z-2y = g(x)

Los valores iniciales de t0=0, y0=0, z0=-5 se usan en el algoritmo.

En este caso también se requiere conocer un intervalo de tiempo a observar [0,tn=6] y definir el tamaño de paso o resolución del resultado

h = \delta t = \frac{b-a}{n} = \frac{6-0}{60} = 0.1

t0 = 0, y0 = 0,  y’0 = z0 = -5

iteración 1

K1y = h * f(ti,yi,zi) = 0.1 (-5) = -0.5

K1z = h * g(ti,yi,zi) ) = 0.1*(-3(-5)-2(0)) = 1.5

K2y = h * f(ti+h, yi + K1y, zi + K1z) = 0.1(-5+1.5) = -0.35

K2z = h * g(ti+h, yi + K1y, zi + K1z) = 0.1 ( -3(-5+1.5) – 2(0-0.5)) = 1.15

yi = yi + (K1y+K2y)/2 =0+(-0.5-0.35)/2=-0.425

zi = zi + (K1z+K2z)/2 = -5+(1.5+1.15)/2 = -3.675

ti = ti + h = 0+0.1 = 0.1

iteración 2

K1y = 0.1 (-3.675) = -0.3675

K1z = 0.1*(-3(-3.675)-2(-0.425)) = 1.1875

K2y = 0.1(-3.675+ 1.1875) = -0.24875

K2z = 0.1 ( -3(-3.675+ 1.1875) – 2(-0.425-0.3675)) = 0.90475

yi = -0.425+(-0.3675-0.24875)/2=-0.7331

zi = -3.675+( 1.1875+0.90475)/2 = -2.6288

ti =0.1+0.1 = 0.2

iteración 3

continuar como ejercicio

El algoritmo permite obtener la gráfica y la tabla de datos.

los valores de las iteraciones son:

t, y, z
[[ 0.        0.       -5.      ]
 [ 0.1      -0.425    -3.675   ]
 [ 0.2      -0.733125 -2.628875]
 [ 0.3      -0.949248 -1.807592]
 [ 0.4      -1.093401 -1.167208]
 [ 0.5      -1.18168  -0.67202 ]
 [ 0.6      -1.226984 -0.293049]
 [ 0.7      -1.239624 -0.006804]
 [ 0.8      -1.227806  0.205735]
 [ 0.9      -1.19804   0.359943]
 [ 1.       -1.155465  0.468225]
 [ 1.1      -1.104111  0.540574]
 [ 1.2      -1.047121  0.585021]
 [ 1.3      -0.986923  0.608001]
 [ 1.4      -0.925374  0.614658]
 [ 1.5      -0.863874  0.609087]
 [ 1.6      -0.803463  0.594537]
 [ 1.7      -0.744893  0.573574]
 [ 1.8      -0.68869   0.548208]
 [ 1.9      -0.635205  0.520011]
 [ 2.       -0.584652  0.490193]
 [ 2.1      -0.53714   0.459683]
 [ 2.2      -0.492695  0.42918 ]
 [ 2.3      -0.451288  0.399206]
 [ 2.4      -0.412843  0.370135]
 [ 2.5      -0.377253  0.342233]
 [ 2.6      -0.34439   0.315674]
 [ 2.7      -0.314114  0.290567]
 [ 2.8      -0.286275  0.266966]
 [ 2.9      -0.26072   0.244887]
 [ 3.       -0.237297  0.224314]
 [ 3.1      -0.215858  0.205211]
 [ 3.2      -0.196256  0.187526]
 [ 3.3      -0.178354  0.171195]
 [ 3.4      -0.162019  0.156149]
 [ 3.5      -0.147126  0.142312]
 [ 3.6      -0.133558  0.129611]
 [ 3.7      -0.121206  0.117969]
 [ 3.8      -0.109966  0.107312]
 [ 3.9      -0.099745  0.097569]
 [ 4.       -0.090454  0.08867 ]
 [ 4.1      -0.082013  0.080549]
 [ 4.2      -0.074346  0.073146]
 [ 4.3      -0.067385  0.066401]
 [ 4.4      -0.061067  0.06026 ]
 [ 4.5      -0.055334  0.054673]
 [ 4.6      -0.050134  0.049591]
 [ 4.7      -0.045417  0.044972]
 [ 4.8      -0.04114   0.040776]
 [ 4.9      -0.037263  0.036964]
 [ 5.       -0.033748  0.033503]
 [ 5.1      -0.030563  0.030362]
 [ 5.2      -0.027677  0.027512]
 [ 5.3      -0.025062  0.024926]
 [ 5.4      -0.022692  0.022581]
 [ 5.5      -0.020546  0.020455]
 [ 5.6      -0.018602  0.018527]
 [ 5.7      -0.016841  0.01678 ]
 [ 5.8      -0.015246  0.015196]
 [ 5.9      -0.013802  0.013761]
 [ 6.       -0.012494  0.012461]]

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = -5

h = 0.1
tn = 6
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, label='y(t)')

plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

s3Eva_2021PAOI_T2 Tensiones mínimas en cables por carga variable

Ejercicio: 3Eva_2021PAOI_T2 Tensiones mínimas en cables por carga variable

El ejercicio usa el resultado del tema anterior, planteando una función de Python como la solución para valores dados. Se requiere una función, para disponer de los valores solución en cada llamada para el intervalo de análisis.

Por lo que básicamente lo que se pide es usar algún algoritmo de búsqueda de raíces. Para simplicidad en la explicación se toma el algoritmo de la bisección.

Los resultados se grafican como theta vs Tensiones, y el punto a buscar es cuando las tensiones en los cables son de igual magnitud, es decir:

TCA = TCB

Resultando en :

Resultado: [TCA, TCB], diferencia
[array([3.46965006e-14, 4.00000000e+02]), -399.99999999999994]
tetha, TCA, TCB
[[-2.61799388e-01  3.46965006e-14  4.00000000e+02]
 [-1.74532925e-01  3.70996817e+01  3.85789041e+02]
 [-8.72664626e-02  7.39170124e+01  3.68641994e+02]
 [ 8.32667268e-17  1.10171790e+02  3.48689359e+02]
 [ 8.72664626e-02  1.45588094e+02  3.26082988e+02]
 [ 1.74532925e-01  1.79896384e+02  3.00994928e+02]
 [ 2.61799388e-01  2.12835554e+02  2.73616115e+02]
 [ 3.49065850e-01  2.44154918e+02  2.44154918e+02]
 [ 4.36332313e-01  2.73616115e+02  2.12835554e+02]
 [ 5.23598776e-01  3.00994928e+02  1.79896384e+02]
 [ 6.10865238e-01  3.26082988e+02  1.45588094e+02]
 [ 6.98131701e-01  3.48689359e+02  1.10171790e+02]
 [ 7.85398163e-01  3.68641994e+02  7.39170124e+01]
 [ 8.72664626e-01  3.85789041e+02  3.70996817e+01]]
       raiz en:  0.34898062924398343 radianes
       raiz en:  19.995117187500004 grados
error en tramo:  8.522115488257542e-05
>>> 

Instrucciones en Python

se añaden las instrucciones de la bisección al algoritmo del tema anterior, para encontrar el punto de intersección,

import numpy as np
import matplotlib.pyplot as plt

# Tema 1
def funcion(P,theta,alfa,beta):
    # ecuaciones
    A = np.array([[np.cos(alfa), -np.cos(beta)],
                  [np.sin(alfa),  np.sin(beta)]])
    B = np.array([P*np.sin(theta), P*np.cos(theta)])

    # usar un algoritmo directo
    X = np.linalg.solve(A,B)
    
    diferencia = X[0]-X[1]
    return([X,diferencia])    

# INGRESO
alfa = np.radians(35)
beta = np.radians(75)
P = 400

# PROCEDIMIENTO
theta = beta-np.radians(90)
resultado = funcion(P,theta,alfa, beta)

# SALIDA
print("Resultado: [TCA, TCB], diferencia")
print(resultado)

# Tema 1b --------------
# PROCEDIMIENTO
dtheta = np.radians(5)
theta1 = beta-np.radians(90)
theta2 = np.radians(90)-alfa

tabla = []
theta = theta1
while not(theta>=theta2):
    X = funcion(P,theta,alfa,beta)[0] # usa vector X
    tabla.append([theta,X[0],X[1]])
    theta = theta + dtheta
    
tabla = np.array(tabla)
thetai = np.degrees(tabla[:,0])
Tca = tabla[:,1]
Tcb = tabla[:,2]

# SALIDA
print('tetha, TCA, TCB')
print(tabla)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
# plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()


# Tema 2 -------------------------
# busca intersección con Bisección
diferencia = Tca-Tcb
donde_min  = np.argmin(np.abs(diferencia))
a = tabla[donde_min-1,0]
b = tabla[donde_min+1,0]
tolera = 0.0001

tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = funcion(P,a,alfa,beta)[1] # usa delta
    fb = funcion(P,b,alfa,beta)[1]
    fc = funcion(P,c,alfa,beta)[1]
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: 
        a = a
        b = c
    if cambia > 0:
        a = c
        b = b
    tramo = b-a

# SALIDA
print('       raiz en: ', c,"radianes")
print('       raiz en: ', np.degrees(c),"grados")
print('error en tramo: ', tramo)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()

s3Eva_2021PAOI_T1 Tensiones en cables por carga variable

Ejercicio: 3Eva_2021PAOI_T1 Tensiones en cables por carga variable

Planteamiento del problema

Las ecuaciones de equilibrio del sistema corresponden a:

-T_{CA} \cos (\alpha) + T_{CB} \cos (\beta) + P \sin (\theta) = 0 T_{CA} \sin (\alpha) + T_{CB} \sin (\beta) - P \cos (\theta) = 0

se reordenan considerando que P y θ son valores constantes para cualquier caso

T_{CA} \cos (\alpha) - T_{CB} \cos (\beta) = P \sin (\theta) T_{CA} \sin (\alpha) + T_{CB} \sin (\beta) = P \cos (\theta)

se convierte a la forma matricial

\begin{bmatrix} \cos (\alpha) & -\cos (\beta) \\ \sin (\alpha) & \sin (\beta) \end{bmatrix} \begin{bmatrix} T_{CA} \\ T_{CB} \end{bmatrix} = \begin{bmatrix} P \sin (\theta) \\ P \cos (\theta) \end{bmatrix}

tomando valores por ejemplo:

α = 35°, β = 75°, P = 400 lb, Δθ = 5°

θ = 75-90 = -15

\begin{bmatrix} \cos (35) & -\cos (75) \\ \sin (35) & \sin (75) \end{bmatrix} \begin{bmatrix}T_{CA} \\ T_{CB} \end{bmatrix} = \begin{bmatrix} 400 \sin (-15) \\ 400 \cos (15) \end{bmatrix}

Desarrollo analítico

matriz aumentada

\begin{bmatrix} \cos (35) & - \cos (75) & 400 \sin (-15) \\ \sin (35) & \sin (75 ) & 400 \cos (15 ) \end{bmatrix}
A = 
[[ 0.81915204 -0.25881905]
 [ 0.57357644  0.96592583]]
B = 
[-103.52761804  386.37033052]

AB = 
[[ 0.81915204 -0.25881905 -103.52761804]
 [ 0.57357644  0.96592583 386.37033052]]

Pivoteo parcial por filas

cos(-15°) tendría mayor magnitud que sin(-15°), por lo que la matriz A se encuentra pivoteada

Eliminación hacia adelante

pivote =  0.81915204/0.57357644
[[ 0.81915204 -0.25881905 -103.52761804]
 [ 0.0         1.63830408  655.32162903]]

Sustitución hacia atras

usando la última fila:

1.63830408 TCB = 655.32162903
TCB = 400

luego la primera fila:

0.81915204TCA -0.25881905TCB = -103.52761804

0.81915204TCA = 0.25881905TCB  -103.52761804

TCA = 2,392 x10-6 ≈ 0

Se deja como tarea realizar el cálculo para:  θ+Δθ

Instrucciones en Python

Resultado:

Resultado: [TCA, TCB], diferencia 
[array([3.46965006e-14, 4.00000000e+02]), -399.99999999999994]

usando el intervalo para θ1 y θ2:

con las instrucciones:

import numpy as np
import matplotlib.pyplot as plt

# Tema 1
def funcion(P,theta,alfa,beta):
    # ecuaciones
    A = np.array([[np.cos(alfa), -np.cos(beta)],
                  [np.sin(alfa),  np.sin(beta)]])
    B = np.array([P*np.sin(theta), P*np.cos(theta)])

    # usar un algoritmo directo
    X = np.linalg.solve(A,B)
    
    diferencia = X[0]-X[1]
    return([X,diferencia])    

# INGRESO
alfa = np.radians(35)
beta = np.radians(75)
P = 400

# PROCEDIMIENTO
theta = beta-np.radians(90)
resultado = funcion(P,theta,alfa, beta)

# SALIDA
print("Resultado: [TCA, TCB], diferencia")
print(resultado)

# Tema 1b --------------
# PROCEDIMIENTO
dtheta = np.radians(5)
theta1 = beta-np.radians(90)
theta2 = np.radians(90)-alfa

tabla = []
theta = theta1
while not(theta>=theta2):
    X = funcion(P,theta,alfa,beta)[0] # usa vector X
    tabla.append([theta,X[0],X[1]])
    theta = theta + dtheta
    
tabla = np.array(tabla)
thetai = np.degrees(tabla[:,0])
Tca = tabla[:,1]
Tcb = tabla[:,2]

# SALIDA
print('tetha, TCA, TCB')
print(tabla)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
# plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()

s2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Ejercicio: 2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Dada la EDP elíptica

\frac{\partial ^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2} = 0 0 \lt x \lt \frac{1}{2}, 0 \lt y\lt \frac{1}{2}

Se convierte a la versión discreta usando diferencias divididas centradas:

 


\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} =

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2

\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = 0

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1
u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = 0
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = 0

que permite plantear las ecuaciones para cada punto en posición [i,j]

En cada iteración se requiere el uso de los valores en la frontera

u(x,0)=0, 0 \leq x \leq \frac{1}{2} u(0,y)=0 , 0\leq y \leq \frac{1}{2} u\Big(x,\frac{1}{2} \Big) = 200 x , 0 \leq x \leq \frac{1}{2} u\Big(\frac{1}{2} ,y \Big) = 200 y , 0 \leq y \leq \frac{1}{2}

Iteraciones

i=1, j=1

u[0,1]-4u[1,1]+u[2,1] + + u[1,0]+u[1,2] = 0

usando los valores en la frontera,

0-4u[1,1]+u[2,1] + 0+u[1,2] = 0 -4u[1,1]+u[2,1] + u[1,2] = 0

i=2, j=1

u[1,1]-4u[2,1]+u[3,1] + + u[2,0]+u[2,2] = 0

usando los valores en la frontera,

u[1,1]-4u[2,1]+200(1/6) + 0+u[2,2] = 0 u[1,1]-4u[2,1] + u[2,2] = -200(1/6)

i=1, j=2

u[0,2]-4u[1,2]+u[2,2] + + u[1,1]+u[1,3] = 0 0 - 4u[1,2]+u[2,2] + u[1,1]+200\frac{1}{6} = 0 - 4u[1,2] + u[2,2]+u[1,1] = -200\frac{1}{6}

i=2, j=2

u[1,2]-4u[2,2]+u[3,2] + + u[2,1]+u[2,3] = 0 u[1,2]-4u[2,2]+200\frac{2}{6} + u[2,1]+200\frac{2}{6} = 0 u[1,2]-4u[2,2]+ u[2,1] = -(2)200\frac{2}{6}

Sistema de ecuaciones a resolver:

\begin{bmatrix} -4 & 1 &1&0\\1 &-4&0&1\\1&0&-4&1 \\0&1&1&-4\end{bmatrix} \begin{bmatrix} u[1,1]\\u[2,1] \\u[1,2]\\u[2,2] \end{bmatrix} = \begin{bmatrix} 0\\-200(1/6)\\-200(1/6)\\-200(4/6) \end{bmatrix}

Resolviendo el sistema se tiene:

[11.11111111 22.22222222 22.22222222 44.44444444]

Instrucciones Python

import numpy as np

A = np.array([[-4, 1, 1, 0],
              [ 1,-4, 0, 1],
              [ 1, 0,-4, 1],
              [ 0, 1, 1,-4]])

B = np.array([0,-200*(1/6),-200*(1/6),-200*(4/6)])

x= np.linalg.solve(A,B)

print(x)

s2Eva_2021PAOI_T1 Masa transportada por tubo

Ejercicio: 2Eva_2021PAOI_T1 Masa transportada por tubo

Las expresiones siguientes se usan dentro de la expresión del integral

Q(t)=9+4 \cos ^2 (0.4t) c(t)=5e^{-0.5t}+2 e^{-0.15 t} M = \int_{t_1}^{t_2} Q(t)c(t) dt

literal a

Usando los valores dados para el intervalo [2,8] con 6 tramos h = (8-2)/6 =1

Se usa los valores de cada ti en Se puede obtener una tabla de valores muestreados para integrar f(t) = Q(i)c(t)

[ti,	 Qi,	 Ci,	 fi]
[ 2.     10.9416  3.321  36.3374]
[ 3.      9.5252  2.3909 22.7739]
[ 4.      9.0034  1.7743 15.9747]
[ 5.      9.6927  1.3552 13.1352]
[ 6.     11.175   1.0621 11.8687]
[ 7.     12.5511  0.8509 10.6793]
[ 8.     12.9864  0.694   9.0121]

Para el integral se usan los valores por cada dos tramos

I\cong \frac{1}{3}[36.3374+4(22.7739) + 15.9747] + \frac{1}{3}[15.9747+4(13.1352) + 11.8687] + \frac{1}{3}[11.8687+4(10.6793) + 9.0121] I = 95.7965

literal b

L acota de error de truncamiento por cada fórmula usada, se estima como O(h5),

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

para un valor de z entre [a,b]

por lo que al usar 3 veces la formula de Simpson se podría estimar en:

error_{trunca} = 3(1^5/90) = 0.033

literal c

El resultado se puede mejorar de dos formas:

1. Dado que el número de tramos es múltiplo de 3, se puede cambiar la fórmula a Simpon de 3/8, que tendría una cota de error menor

2. Aumentar el número de tramos disminuyendo el valor de h para que el error disminuya. Por ejemplo si se reduce a 0.5, el error disminuye en el orden de 0.55

Podría recomendar la segunda opión, pues a pesar que se aumenta la cota de error por cada vez que se usa la fórmula, el error de cada una disminuye en ordenes de magnitud 0,03125


La gráfica del ejercicio es:

Instrucciones en Python

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
Q = lambda t: 9+4*(np.cos(0.4*t)**2)
C = lambda t: 5*np.exp(-0.5*t)+2*np.exp(-0.15*t)

t1 = 2
t2 = 8
n  = 6

# PROCEDIMIENTO
muestras = n+1 
dt = (t2-t1)/n
ti = np.arange(t1,t2+dt,dt)
Qi = Q(ti)
Ci = C(ti)
fi = Qi*Ci

# integración con Simpson 1/3
h= dt
I13 = 0
for i in range(0,6,2):
    S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
    I13 = I13 + S13

# SALIDA
np.set_printoptions(precision=4)
print("[ti,\t Qi,\t Ci,\t fi]")
for i in range(0,muestras,1):
    print(np.array([ti[i],Qi[i],Ci[i],fi[i]]))
print("Integral S13: ",I13)
# grafica
plt.plot(ti,Qi, label = "Q(t)")
plt.plot(ti,Ci, label = "c(t)")
plt.plot(ti,fi, label = "f(t)")
plt.plot(ti,Qi,'.b')
plt.plot(ti,Ci,'.r')
plt.plot(ti,fi,'.g')
plt.xlabel("t")
plt.ylabel("f(t)=Q(t)*c(t)")
plt.show()

s2Eva_2021PAOI_T2 EDO para cultivo de peces

Ejercicio: 2Eva_2021PAOI_T2 EDO para cultivo de peces

Siendo la captura una constante mas una función periódica,

h(t) = a + b \sin (2 \pi t)

La ecuación EDO del ejercicio, junto a las constantes a=0.9 y b=0.75, r=1

\frac{\delta y(t)}{\delta t} = r y(t)-h(t)

se convierte en:

\frac{\delta y(t)}{\delta t} = (1) y(t)- \Big( 0.9 + .75 \sin (2 \pi t)\Big) \frac{\delta y(t)}{\delta t} = y(t)- 0.9 - .75 \sin (2 \pi t)

Considerando que la población inicial de peces es 1 o 100%, y(0)=1

literal a

h=1/12
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,2),dtype=float)
estimado[0] = [0,1]
ti = 0
yi = 1
for i in range(1,tamano,1):
    K1 = 1/12 * d1y(ti,yi)
    K2 = 1/12 * d1y(ti+1/24, yi + K1/2)
    K3 = 1/12 * d1y(ti+1/24, yi + K2/2)
    K4 = 1/12 * d1y(ti+1/12, yi + K3)

    yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
    ti = ti + 1/12
        
    estimado[i] = [ti,yi]

literal b

iteración i=0

t(0) = 0

y(0) = 1

K1 = \frac{1}{12} \Big(1- 0.9 - .75 \sin (2 \pi 0)\Big) = 0,008333 K2 = \frac{1}{12} \Big(1- 0.9 - .75 \sin \Big(2 \pi (0+\frac{1}{12})\Big)\Big) = -0.02222 y(1) = 0 + \frac{0.008333+(-0.02222)}{2} = 0.9930 t(1) = 0 + \frac{1}{12} = \frac{1}{12}

iteración i=1

t(1) = \frac{1}{12}

y(1) = 0.9930

K1 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big( 2 \pi\frac{1}{12}\Big)\Big) = -0.02349 K2 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big(2 \pi (\frac{1}{12}+\frac{1}{12})\Big)\Big) = -0.04832 y(1) = 0.9930 + \frac{-0.02349+(-0.04832)}{2} = 0.9571 t(1) = \frac{1}{12} + \frac{1}{12} = \frac{2}{12}

iteración i=2

t(2) = \frac{2}{12}

y(1) = 0.9571

K1 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big( 2 \pi\frac{2}{12}\Big)\Big) = -0.04936 K2 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big(2 \pi (\frac{2}{12}+\frac{1}{12})\Big)\Big) = -0.06185 y(1) = 0.9571 + \frac{-0.04936+(-0.06185)}{2} = 0.9015 t(3) = \frac{2}{12} + \frac{1}{12} = \frac{3}{12}

literal c

Resultado del algoritmo, muestra que la estragegia de cosecha, en el tiempo no es sostenible, dado que la población de peces en el tiempo decrece.

estimado[xi,yi,K1,K2]
[[ 0.0000e+00  1.0000e+00  8.3333e-03 -2.2222e-02]
 [ 8.3333e-02  9.9306e-01 -2.3495e-02 -4.8330e-02]
 [ 1.6667e-01  9.5714e-01 -4.9365e-02 -6.1852e-02]
 [ 2.5000e-01  9.0153e-01 -6.2372e-02 -5.9196e-02]
 [ 3.3333e-01  8.4075e-01 -5.9064e-02 -4.1109e-02]
 [ 4.1667e-01  7.9066e-01 -4.0361e-02 -1.2475e-02]
 [ 5.0000e-01  7.6425e-01 -1.1313e-02  1.8994e-02]
 [ 5.8333e-01  7.6809e-01  2.0257e-02  4.4822e-02]
 [ 6.6667e-01  8.0063e-01  4.5845e-02  5.8039e-02]
 [ 7.5000e-01  8.5257e-01  5.8547e-02  5.5053e-02]
 [ 8.3333e-01  9.0937e-01  5.4907e-02  3.6606e-02]
 [ 9.1667e-01  9.5513e-01  3.5844e-02  7.5807e-03]
 [ 1.0000e+00  9.7684e-01  6.4031e-03 -2.4313e-02]
 [ 1.0833e+00  9.6788e-01 -2.5593e-02 -5.0602e-02]
 [ 1.1667e+00  9.2978e-01 -5.1645e-02 -6.4322e-02]
 [ 1.2500e+00  8.7180e-01 -6.4850e-02 -6.1881e-02]
 [ 1.3333e+00  8.0844e-01 -6.1757e-02 -4.4027e-02]
 [ 1.4167e+00  7.5554e-01 -4.3288e-02 -1.5645e-02]
 [ 1.5000e+00  7.2608e-01 -1.4494e-02  1.5549e-02]
 [ 1.5833e+00  7.2661e-01  1.6800e-02  4.1077e-02]
 [ 1.6667e+00  7.5554e-01  4.2089e-02  5.3969e-02]
 [ 1.7500e+00  8.0357e-01  5.4464e-02  5.0630e-02]
 [ 1.8333e+00  8.5612e-01  5.0470e-02  3.1799e-02]
 [ 1.9167e+00  8.9725e-01  3.1021e-02  2.3563e-03]
 [ 2.0000e+00  9.1394e-01  1.1619e-03 -2.9991e-02]
 [ 2.0833e+00  8.9953e-01 -3.1289e-02 -5.6773e-02]
 [ 2.1667e+00  8.5550e-01 -5.7835e-02 -7.1028e-02]
 [ 2.2500e+00  7.9107e-01 -7.1578e-02 -6.9169e-02]
 [ 2.3333e+00  7.2069e-01 -6.9069e-02 -5.1948e-02]
 [ 2.4167e+00  6.6018e-01 -5.1235e-02 -2.4254e-02]
 [ 2.5000e+00  6.2244e-01 -2.3130e-02  6.1924e-03]
 [ 2.5833e+00  6.1397e-01  7.4142e-03  3.0909e-02]
 [ 2.6667e+00  6.3313e-01  3.1888e-02  4.2918e-02]
 [ 2.7500e+00  6.7053e-01  4.3378e-02  3.8619e-02]
 [ 2.8333e+00  7.1153e-01  3.8421e-02  1.8746e-02]
 [ 2.9167e+00  7.4012e-01  1.7926e-02 -1.1830e-02]
 [ 3.0000e+00  7.4317e-01  0.0000e+00  0.0000e+00]]

Instrucciones en Python

# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        estimado[i-1,2:]=[K1,K2]
        estimado[i] = [xi,yi,0,0]
    return(estimado)

# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
# d1y = y' = f, d2y = y'' = f'
a =0.9; b=0.75; r=1
d1y = lambda t,y: r*y-(a+b*np.sin(2*np.pi*t))
x0 = 0
y0 = 1
h  = 1/12
muestras = 12*3

# PROCEDIMIENTO
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi,K1,K2]')
print(puntosRK2)


# Gráfica
import matplotlib.pyplot as plt


plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='y Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

s1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

Ejercicio: 1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

literal a

Se plantea el sistema de ecuaciones de acuerdo a la cantidad de medicamentos que se administra a cada tipo de paciente, siendo las variables X=[a,b,c,d] de acuerdo a la tabla:

Niños Adolescentes Adultos Adultos Mayores Medicamentos /semana
Medicamento_A 0.3 0.4 1.1 4.7 3500
Medicamento_B 1 3.9 0.15 0.25 3450
Medicamento_C 0 2.1 5.6 1.0 6500
0.3 a + 0.4 b + 1.1 c + 4.7 d = 3500 a + 3.9 b + 0.15 c + 0.25 d = 3450 0 a +2.1 b + 5.6 c + 1.0 d = 6500

Mostrando que el número de incógnitas no es igual al número de ecuaciones, por lo que sería necesario de disponer de otra ecuación, o  usar una variable libre.


literal b

Siguiendo las indicaciones sobre la variable libre que sea para el grupo de pacientes niños, te tiene que

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K 3.9 b + 0.15 c + 0.25 d = 3450 - K 2.1 b + 5.6 c + 1.0 d = 6500 - 0 k

y haciendo K = 100, se convierte en :

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3(100) = 3470 3.9 b + 0.15 c + 0.25 d = 3450 - 100 = 3350 2.1 b + 5.6 c + 1.0 d = 6500 - 0 = 6500

literal c

Antes de desarrollar el ejercicio para el algoritmo se requiere convertir el sistema de ecuaciones a su forma matricial. Por lo que se plantea la matriz aumentada:

\begin{pmatrix} 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

Se realiza el pivoteo parcial por filas, que al aplicar en la primera columa, se intercambia la primera y segunda fila, de tal forma que el mayor valor quede en la primera casilla de la diagonal.

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

se continua el mismo proceso para la segunda casilla de la diagonal hacia abajo, se aplica también pivoteo parcial,

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \end{pmatrix}

Con lo que el sistema queda listo para resolver por cualquier método.

Si seleccionamos Gauss-Seidel, el vector inicial de acuerdo al enunciado será X0=[100,100,100]

y las ecuaciones a resolver son:

b = \frac{3350 - 0.15 c - 0.25 d}{3.9} c = \frac{6500 - 2.1 b - 1.0 d}{5.6} d = \frac{3470- 0.4 b - 1.1 c}{4.7}

iteración 1

X0=[100,100,100]

b = \frac{3350 - 0.15(100) - 0.25 (100)}{3.9} = 848.71 c = \frac{6500 - 2.1 (848.71 ) - 1.0 (100)}{5.6} = 824.58 d = \frac{3470- 0.4 (848.71 ) - 1.1 (824.58)}{4.7} = 473.07

X1=[848.71, 824.58, 473.07]

diferencia = [848.71, 824.58, 473.07] – [100,100,100]

diferencia = [748.71, 724.58, 373.07]

error = max|diferencia| = 748.71

resultado del error es mayor que tolerancia de 0.01, se continua con la siguiente iteración.

iteración 2

X1=[848.71, 824.58, 473.07]

b = \frac{3350 - 0.15 (824.58) - 0.25 (473.07)}{3.9} = 796.93 c = \frac{6500 - 2.1 (796.93) - 1.0 (473.07)}{5.6} = 777.38 d = \frac{3470- 0.4 (796.93) - 1.1 (777.38)}{4.7} = 488.53

X2 = [796.93, 777.38, 488.53]

diferencia = [796.93, 777.38, 488.53]  – [848.71, 824.58, 473.07]

diferencia = [-51.78 ,  -47.2, 15.52]

error = max|diferencia| = 51.78

iteración 3

X2 = [796.93, 777.38, 488.53]

b = \frac{3350 - 0.15 (777.38) - 0.25 (488.53)}{3.9} = 797.75 c = \frac{6500 - 2.1 (797.75) - 1.0 (488.53)}{5.6} = 774.31 d = \frac{3470- 0.4 (797.75) - 1.1 (774.31)}{4.7} = 489.18

x3 = [ 797.75, 774.31, 489.18]

diferencias = [ 797.75, 774.31, 489.18] – [796.93, 777.38, 488.53]

diferencias = [0.82, -3.07, 0.65]

error = max|diferencia| = 3.07

Observación, el error disminuye en cada iteración, por lo que el sistema converge.

solución con el algoritmo, solo se toma la parte entera de la respuesta, pues los pacientes son números enteros.

respuesta X: [797.83, 774.16, 489.20]

con lo que el primer resultado es:

 X = [797, 774, 489]

literal d

Si la cantidad de pacientes presentados en una seman es la que se indica en el enunciado [350,1400,1500,1040], se determina que el sistema hospitalario estatal no podrá atender a todos los pacientes al compara con la capacidad encontrada en el literal c.

Hay un exceso de 2129 pacientes encontrados por grupo:

exceso = [350, 1400, 1500, 1040] – [100, 797, 774, 489] = [250, 603, 726, 551]

con lo que se recomienda una medida de confinamiento para disminuir los contagios y poder atender a los pacientes que se presenten.

literal e

Siendo K = 0, al estar vacunados todos los niños, y suponiendo que la vacuna es una cura, se tiene:

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K = 3500 3.9 b + 0.15 c + 0.25 d = 3450 - K = 3450 2.1 b + 5.6 c + 1.0 d = 6500 - 0 k = 6500

pivoteado parcialmente por filas se tiene:

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3450 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3500 \end{pmatrix}

Resolviendo por un método directo:

se obtiene:

respuesta X: [823.46, 763.35, 495.94]

que al compararse con la capacidad anterior en números enteros se encuentra una diferencia de un incremento de 21 pacientes neto ante la condición de usar todos los medicamentos disponibles.

diferencia = [823, 763, 495] – [797, 774, 489] = [26, -11, 6]


 Instrucciones en Python

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# INGRESO
A = np.array([[0.4, 1.10, 4.7 ],
              [3.9, 0.15, 0.25],
              [2.1, 5.60, 1.0  ]])
k = 100
B = np.array([3500.-0.3*k,3450-1*k,6500-0*k])

X0  = np.array([100.0,100,100])

tolera = 0.001
iteramax = 100

# PROCEDIMIENTO
# Matriz aumentada
Bcolumna = np.transpose([B])
AB  = np.concatenate((A,Bcolumna),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)
    
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

A = np.copy(AB[:,:n])
B = np.copy(AB[:,n])

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
    print(X)
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)

s1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

Ejercicio: 1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

Literal a

Evaluando las sucesión de la forma recursiva:

xi
[-0.45       -0.4383     -0.4458     -0.441      -0.4441 
 -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
 -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
errores
[ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
 -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
  1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]

literal b

Se puede afirmar que converge, observe la diferencia entre cada dos valores consecutivos de la sucesión … (continuar de ser necesario)

literal c

Para el algoritmo se requiere la función f(x) y su derivada f'(x)

f(x) = x +ln(x+2) f'(x) = 1 + \frac{1}{x+2}

x0 = -0.45

itera = 1

x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} x_{1} = x_0 -\frac{f(x_0)}{f'(x_0)} = -0.45 -\frac{-0.45+ln(-0.45+2)}{1 + \frac{1}{-0.45+2}}

x1 = -0.44286

error = |x1-x0| = |-0.44286 -(-0.45)| = 0.007139

itera = 2

x_{2} = -0.4428 -\frac{-0.4428 + ln(-0.45+2)}{1 + \frac{1}{-0.4428+2}}

x1 = -0.44286

error = |x1-x0| = |-0.44285 -(-4.4286)| = 6.4394e-06

con lo que se cumple el valor de tolerancia y no se requiere otra iteración

la raiz se encuentra en x = -0.44286

Solución con algoritmo

xi
[-0.45       -0.4383     -0.4458     -0.441      -0.4441 
 -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
 -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
errores
[ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
 -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
  1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]
['xi', 'xnuevo', 'tramo']
[[-4.5000e-01 -4.4286e-01  7.1392e-03]
 [-4.4286e-01 -4.4285e-01  6.4394e-06]]
raiz en:  -0.44285440100759543
con error de:  6.439362322474551e-06
>>> 

Instrucciones en Python

import numpy as np
import matplotlib.pyplot as plt

# literal a sucesión en forma recursiva
def secuenciaL(n,x0):
    if n == 0:
        xn = x0
    if n>0:
        xn = np.log(1/(2+secuenciaL(n-1,x0)))
    return(xn)
x0 = -0.45
n = 15
xi = np.zeros(n,dtype=float)
xi[0] = x0
errado = np.zeros(n,dtype=float)
for i in range(1,n,1):
    xi[i] = secuenciaL(i,x0)
    errado[i-1] = xi[i] - xi[i-1]
   
np.set_printoptions(precision=4)
print('xi: ')
print(xi)
print('errado: ')
print(errado)

#Grafica literal a y b
plt.plot(xi)
plt.plot(xi,'o')
plt.xlabel('n')
plt.ylabel('xn')
plt.show()

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np

# INGRESO
fx  = lambda x: x+np.log(x+2)
dfx = lambda x: 1+ 1/(x+2)

x0 = -0.45
tolera = 0.0001

a = -0.5
b = -0.2
muestras = 21
# PROCEDIMIENTO
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

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# para la gráfica
xj = np.linspace(a,b,muestras)
fj = fx(xj)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

plt.plot(xj,fj)
plt.axhline(0, color='grey')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

litera d: tarea

s1eva_2021PAOI_T3 Interpolar, modelo de contagios 2020

Ejercicio: 1eva_2021PAOI_T3 Interpolar, modelo de contagios 2020

literal a

Los datos de los pacientes casos graves entre las semanas 11 a la 20, que son el intervalo donde será válido el polinomio de interpolación son:

semana 11 12 13 14 15 16 17 18 19 20
casos graves 1503 3728 7154 6344 4417 3439 2791 2576 2290 2123

de los cuales solo se usarán los indicados en el literal a : 11,13,16,18,20.

xi0 = [    9,   10,   11,   12,   13,   14,
          15,   16,   17,   18,   19,   20,
          21,   22,   23,   24,   25,   26 ])
fi0 = [ 1435, 1645, 1503, 3728, 7154, 6344,
        4417, 3439, 2791, 2576, 2290, 2123,
        2023, 2067, 2163, 2120, 2125, 2224 ])
xi = [  11,   13,   16,   18,   20]
fi = [1503, 7154, 3439, 2576, 2123]

Se observa que los datos estan ordenados en forma ascendente respecto a la variable independiente, tambien se determina que no se encuentran equidistantes entre si (13-11=2, 16-13=3). Por lo que se descarta usar el método de diferencias finitas avanzadas.

Los métodos que se podrían usar con puntos no equidistantes en el eje semanas serían el método de diferencias divididas de Newton o el  método de Lagrange.

Seleccionando por ejemplo, Diferencias divididas de Newton, donde primero se realiza la tabla:

xi fi f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0] f[x4,x3,x2,x1,x0]
11 1503 =(7154-1503) /(13-11) = 2825.5 =(-1238.33-2835.5) /(16-11) = -812.76 =(161.36-(-812.76)) /(18-11) = 139.16 =(-15.73-139.16) /(20-11) = -17.21
13 7154 (3439-7154) /(16-13) = -1238.33 (-431.5-(-1238.33)) /(18-13) = 161.36 (51.25-161.36) /(20-13)= -15.73 —-
16 3439 (2576-3439) /(18-16) = -431.5 (-226.5-(-431.5)) /(20-16) = 51.25 —-
18 2576 (2123-2576) /(20-18) = -226.5 —-
20 2123 —-

con lo que se puede contruir el polinomio usando las diferencias divididas para el intervalo dado:

[2825.5    -812.76  139.16  -17.21]
p_4(x) = 1503 + 2825.5(x-11) - 812.76(x - 13)(x - 11) + 139.16(x - 16)(x - 13)(x - 11) - 17.21(x - 18)(x - 16) (x - 13) (x - 11)

Simplificando el algoritmo se tiene:

p_4(x) = - 1172995.28 + 298304.50 x - 27840.50x^2 + 1137.36x^3 - 17.21 x^4


literal b

El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.

xi fi p4(x) |error|
11 1503 1503 0
12 3728 6110.96 2382.96
13 7154 7154 0
14 6344 6293.18 50.81
15 4417 4776.52 359.52
16 3439 3439 0
17 2791 2702.53 88.46
18 2576 2576 0
19 2290 2655.22 365.22
20 2123 2123 0

literal c

Podría aplicarse uno de varios criterios, lo importante por lo limitado del tiempo en la evaluación son las conclusiones y recomendaciones expresadas en el literal e, basadas en lo realizado en los literales c y d. Teniendo como opciones:

– cambiar uno de los puntos selecionados, mateniendo así el grado del polinomio
– aumentar el número de puntos usados para armar el polinomio con grado mayor
– dividir el intervalo en uno o mas segmentos, con el correspondiente número de polinomios.

Se desarrolla la opción de cambiar uno de los puntos seleccionados, usando para esta ocasión como repaso la interpolación de Lagrange. Para los puntos se usa el punto con mayor error de la tabla del literal anterior y se elimina el punto penúltimo, es decir se usa la semana 12 en lugar de la semana 18 de la siguiente forma:

xi = [  11,   12,   13,   16,   20]
fi = [1503, 3728, 7154, 3439, 2123]
p_4(x) = 1503 \frac{(x-12)(x-13)(x-16)(x-20)}{(11-12)(11-13)(11-16)(11-20)} + 3728\frac{(x-11)(x-13)(x-16)(x-20)}{(12-11)(12-13)(12-16)(12-20)} + 7154\frac{(x-11)(x-12)(x-16)(x-20)}{(13-11)(13-12)(13-16)(13-20)} + 3439\frac{(x-11)(x-12)(x-13)(x-20)}{(16-11)(16-12)(16-13)(16-20)} + 2123\frac{(x-11)(x-12)(x-13)(x-16)}{(20-11)(20-12)(20-13)(20-16)}

Simplificando el polinomio:

p_4(x) = \frac{46927445}{21} - \frac{1655552687}{2520} x + \frac{715457663}{10080}x^2 - \frac{8393347}{2520} x^3 + \frac{577153}{10080} x^4

literal d

El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.

xi fi p4(x) error
11 1503 1503 0
12 3728 3728 0
13 7154 7154 0
14 6344 8974.01 2630.01
15 4417 7755.22 3338.22
16 3439 3439 0
17 2791 -2659.13 -5450.13
18 2576 -7849.45 -10425.45
19 2290 -8068.1 -10358.1
20 2123 2123 0

literal e

El cambio aplicado a los puntos usados en el modelo del polinomio disminuyó el error entre las semanas 11 a 13. Sin embargo la magnitud del error aumentó  para las semanas posteriores a la 13, es decir aumentó la distorsión de la estimación y se recomienda realizar otras pruebas para mejorar el modelo aplicando los otros criterios para determinar el que tenga mejor desempeño respecto a la medida de error.



Intrucciones en Python

Literal a y b. Desarrollado a partir del algoritmo desarrollado en clases:

# Polinomio interpolación
# Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi0 = np.array([    9,   10,   11,   12,   13,   14,
                   15,   16,   17,   18,   19,   20,
                   21,   22,   23,   24,   25,   26 ])
fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344,
                 4417, 3439, 2791, 2576, 2290, 2123,
                 2023, 2067, 2163, 2120, 2125, 2224 ])

xi1 = np.array([   11,   12,   13,   14,   15,   16,
                   17,   18,   19,   20 ])
fi1 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439,
                 2791, 2576, 2290, 2123 ])

xi = np.array([  11,   13,   16,   18,   20])
fi = np.array([1503, 7154, 3439, 2576, 2123])

# PROCEDIMIENTO

# Tabla de Diferencias Divididas Avanzadas
titulo = ['i   ','xi  ','fi  ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('F['+str(j-2)+']')

    # cada fila de columna
    i = 0
    paso = j-2 # inicia en 1
    while (i < diagonal):
        denominador = (xi[i+paso]-xi[i])
        numerador = tabla[i+1,j-1]-tabla[i,j-1]
        tabla[i,j] = numerador/denominador
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Divididas
# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    factor = dDividida[j-1]
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# calcula errores en intervalo usado
pfi1 = px(xi1)
errado1 = np.abs(fi1-pfi1)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)
print('errores en intervalo:')
print(xi1)
print(errado1)

# Gráfica
plt.plot(xi0,fi0,'o', label = 'Puntos')
plt.plot(xi,fi,'ro', label = 'Puntos')
for i in range(0,n,1):
    etiqueta = '('+str(xi[i])+','+str(fi[i])+')'
    plt.annotate(etiqueta,(xi[i],fi[i]))
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.grid()
plt.show()

Literal c y d. Se puede continuar con el algoritmo anterior. Como repaso se adjunta un método diferente al anterior.

# Interpolacion de Lagrange
# divisores L solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi0 = np.array([    9,   10,   11,   12,   13,   14,
                   15,   16,   17,   18,   19,   20,
                   21,   22,   23,   24,   25,   26 ])
fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344,
                 4417, 3439, 2791, 2576, 2290, 2123,
                 2023, 2067, 2163, 2120, 2125, 2224 ])

xi2 = np.array([   11,   12,   13,   14,   15,   16,
                   17,   18,   19,   20 ])
fi2 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439,
                 2791, 2576, 2290, 2123 ])

xi = np.array([  11,   12,   13,   16,   20])
fi = np.array([1503, 3728, 7154, 3439, 2123])

# PROCEDIMIENTO
# 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)

# calcula errores en intervalo usado
pfi2 = px(xi2)
errado2 = np.abs(fi2-pfi2)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
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)
print('errores en intervalo:')
print(xi2)
print(errado2)

# Gráfica
plt.plot(xi0,fi0,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

s1Eva_IT2016_T3_MN Tasa interés anual

Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

Propuesta de Solución  empieza con el planteamiento del problema, luego se desarrolla con el método de Bisección y método del Punto Fijo solo con el objetivo de comparar resultados.


Planteamiento del problema

La fórmula del enunciado para el problema es:

A = P \frac{i(1+i)^{n}}{(1+i)^{n} -1}

que con los datos dados se convierte a:

5800 = 35000 \frac{i(1+i)^8}{(1+i)^8 -1} 35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0

que es la forma de f(x) = 0

f(i)=35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800

Intervalo de búsqueda

Como el problema plantea la búsqueda de una tasa de interés, consideramos que:

  • Las tasas de interés no son negativas. ⌉(i<0)
  • Las tasas de interés no son cero en las instituciones bancarias (i≠0)
  • Las tasas muy grandes 1 = 100/100 = 100% tampoco tienen mucho sentido

permite acotar la búsqueda a un intervalo (0,1].
Sin embargo tasas demasiado altas tampoco se consideran en el problema pues el asunto es regulado (superintendencia de bancos), por lo que se podría intentar entre un 1% = 0.01l y la mitad del intervalo 50%= 0.5 quedando

[0.01,0.5]

Tolerancia

si la tolerancia es de tan solo menor a un orden de magnitud que el valor buscado, se tiene que las tasas de interés se representan por dos dígitos después del punto decimal, por lo que la tolerancia debe ser menor a eso.

Por ejemplo: tolerancia < 0.001 o aumentando la precisión

tolera = 0.0001


Método de la Bisección

itera = 1

a = 0.01, b = 0.5

c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 f(0.01) = 35000 \frac{0.01(1+0.01)^8}{(1+0.01)^8-1}-5800 = -1225.83 f(0.5) = 35000 \frac{0.5(1+0.5i)^8}{(1+0.5)^8-1}-5800 = 12410.54

con lo que se verifica que existe cambio de signo al evaluar f(x) en el intervalo y puede existir una raíz.

f(0.255) = 35000 \frac{0.255(1+0.255)^8}{(1+0.255)^8-1}-5800 = 4856.70

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.255]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.255-0.01 = 0.245

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 2

a = 0.01, b = 0.255

c = \frac{0.01+0.255}{2} = 0.1325 f(0.01) = -1225.83 f(0.255) = 4856.70 f(0.1325 ) = 35000 \frac{0.1325(1+0.1325)^8}{(1+0.1325)^8-1}-5800 = 1556.06

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.1325]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.1325-0.01 = 0.1225

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 3

a = 0.01, b = 0.1325

c = \frac{0.01+0.1225}{2} = 0.071 f(0.01) = -1225.83 f(0.1325) = 1556.06 f(0.071 ) = 35000 \frac{0.071(1+0.071)^8}{(1+0.071)^8-1}-5800 = 89.79

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.071]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.071-0.01 = 0.061

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

Para una evaluación del tema en forma escrita es suficiente para mostrar el objetivo de aprendizaje, el valor final se lo encuentra usando el algoritmo.

       raiz en:  0.06724243164062502
error en tramo:  5.981445312500111e-05
iteraciones:  13
>>>

Algoritmo en Python

# 1Eva_IT2016_T3_MN Tasa interés anual
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
a = 0.01
b = 0.5
tolera = 0.0001

# PROCEDIMIENTO
cuenta = 0
np.set_printoptions(precision=3)
tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = fx(a)
    fb = fx(b)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: a = a b = c if cambia > 0:
        a = c
        b = b
    tramo = b-a
    cuenta = cuenta+1

# SALIDA
print('       raiz en: ', c)
print('error en tramo: ', tramo)
print('iteraciones: ',cuenta)


Método del Punto Fijo

El planteamiendo del punto fijo se realiza con x= g(x), por lo que se reordena la ecuación a:

35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0 35000 \frac{i(1+i)^8}{(1+i)^8 -1}=5800 \frac{i(1+i)^8}{(1+i)^8 -1}=\frac{5800}{35000} i=\frac{58}{350} \frac{(1+i)^8 -1} {i(1+i)^8}

con lo que g(x) es:

g(i)=\frac{58}{350} \frac{(1+i)^8 -1} {(1+i)^8}

valor inicial de búsqueda

Para el punto inicial i0 se puede usar uno de los extremos del intervalo propuesto en la sección de planteamiento. Para reducir aun más la búsqueda se pude seleccionar el punto intermedio

 i0 = 0.255

Itera = 1

 i0 = 0.255

g(0.255i)=\frac{58}{350} \frac{(1+0.255)^8 -1} {(1+0.255)^8} = 0.1387

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1387 - 0.255| = 0.1163

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 2

i1 = g(i0 ) = 0.1387

g(0.1387)=\frac{58}{350} \frac{(1+0.1387)^8 -1} {(1+0.1387)^8} = 0.1071

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1071 - 0.1387| = 0,0316

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 3

i2 = g(i1 ) = 0.1071

g(0.1071)=\frac{58}{350} \frac{(1+0.1071)^8 -1} {(1+0.1071)^8} = 0.0922

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.0922 - 0.1071| = 0,0149

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Observación: Como el error disminuye entre cada iteración, se considera que el método converge, si se realizan suficientes iteraciones se cumplierá con el valor de tolerancia y se habrá llegado a la precisión requerida.

Usando el algoritmo se tiene:

iteraciones: 17
    raiz en: 0.06754389199556779
>>>

Algoritmo en Python

# Algoritmo de Punto Fijo
# x0 valor inicial de búsqueda
# error = tolera

import numpy as np
def puntofijo(gx,antes,tolera, iteramax=50):
    itera = 1 # iteración
    nuevo = gx(antes)
    tramo = abs(nuevo-antes)
    while(tramo>=tolera and itera<=iteramax ):
        antes = nuevo
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        itera = itera+1
    respuesta = nuevo
    # Validar respuesta
    if (itera>=iteramax ):
        respuesta = np.nan
    print('iteraciones:',itera)
    return(respuesta)

# PROGRAMA ---------

# INGRESO
gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)

a = 0       # intervalo
b = 0.5

x0 = 0.255
tolera = 0.0001
iteramax  = 50      # itera máximo

# PROCEDIMIENTO
respuesta = puntofijo(gx,0.255,tolera,iteramax)

# SALIDA
print('    raiz en:',respuesta)