Categoría: Solución 2da Eva

Ejercicio: 2Eva_IT2008_T1_MN Producción petroleo

Literal a

Para el cálculo de las derivadas se hace uso de las fórmulas de diferenciación presentadas en la unidad 6, y basadas en el polinomio de Taylor:

f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h) f''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+f(x_i)}{h^2} + O(h)

[ dia, prod, dprod, d2prod]
[[ 1.000e+00  3.345e+03 -1.000e+02  6.600e+01]
 [ 2.000e+00  3.245e+03 -3.400e+01  1.320e+02]
 [ 3.000e+00  3.211e+03  9.800e+01 -5.600e+01]
 [ 4.000e+00  3.309e+03  4.200e+01  1.900e+01]
 [ 5.000e+00  3.351e+03  6.100e+01 -2.430e+02]
 [ 6.000e+00  3.412e+03 -1.820e+02  8.700e+01]
 [ 7.000e+00  3.230e+03 -9.500e+01  9.200e+01]
 [ 8.000e+00  3.135e+03 -3.000e+00  0.000e+00]
 [ 9.000e+00  3.132e+03 -3.000e+00  0.000e+00]
 [ 1.000e+01  3.129e+03  0.000e+00  0.000e+00]]

representados en las siguientes gráfica:

literal b

Dado que las fórmlas de error usadas tienen error del orden h: O(h), el error de las fórmulas es del orden de:

h= dia[1]-dia[0] = 2-1 = 1

literal c

Para el día dos se observa un decrecimiento en la producción, tal como lo refleja el valor negativo de la primera derivada.
Sin embargo para el día siguiente, la producción no mantiene la tasa de decrecimiento, se observa la segunda derivada positiva, Empieza a «acelerar».

---

Las instrucciones en Python para la tabla presentada son:

# 2Eva_IT2008_T1_MN Producción petroleo
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
dia = np.array( [1., 2, 3, 4, 5, 6, 7, 8, 9, 10])
produccion = [3345., 3245, 3211, 3309, 3351, 3412, 3230, 3135, 3132, 3129]
produccion = np.array(produccion, dtype = float)

# PROCEDIMIENTO
n = len(dia)

# primera derivada
dp = np.zeros(n,dtype=float)
for i in range(0,n-1,1):
    dp[i] = (produccion[i+1]-produccion[i])/(dia[i+1]-dia[i])

# segunda derivada
d2p = np.zeros(n,dtype=float)
h = dia[1]-dia[0]
for i in range(0,n-2,1):
    d2p[i] = (produccion[i]-2*produccion[i+1]+ produccion[i+2])/(h**2)

tabla = np.concatenate(([dia],[produccion],[dp],[d2p]),axis=0)
tabla = np.transpose(tabla)

# SALIDA
print(" [ dia, prod, dprod, d2prod]")
print(tabla)

# gráfica
plt.subplot(121)
plt.plot(dia,produccion)
plt.xlabel('dia')
plt.ylabel('producción')
plt.grid()
plt.subplot(122)
plt.plot(dia,dp,color='green',label='dp')
plt.xlabel('dia')
plt.plot(dia,d2p,color='orange',label='d2p')
plt.axhline(0)
plt.legend()
plt.grid()
plt.show()

Ejercicio: 2Eva_IIT2007_T2_AN Lanzamiento vertical proyectil

la ecuación del problema:

m \frac{\delta v}{\delta t} = -mg - kv|v|

se despeja:

\frac{\delta v}{\delta t} = -g - \frac{k}{m}v|v|

y usando los valores indicados en el enunciado:

\frac{\delta v}{\delta t} = -9,8 - \frac{0.002}{0.11}v|v|

con valores iniciales de:

t₀ = 0 , v₀ = 8 , h=0.2

Como muestra inicial, se usa Runge-Kutta de 2do Orden

iteración 1

K_1 = h\frac{\delta v}{\delta t}(0, 8) = 0.2[-9,8 - \frac{0.002}{0.11}8|8|] = -2.1927 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 8 -2.1927 ) = 0.2[-9,8 - \frac{0.002}{0.11}(8 -2.1927)|8 -2.1927|] =-2.0826 v_1 = -9,8 +\frac{-2.1927-2.0826 }{2} = 5.8623 t_1 = t_0 + h = 0 + 0.2 = 0.2 error = O(h^3) = O(0.2^3) = O(0.008)

iteración 2

K_1 = h\frac{\delta v}{\delta t}(0.2, 5.8623) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623)|5.8623|] = -2.085 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 5.8623 -2.085) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623 -2.085)|5.8623 -2.085|] =-2.0119 v_2 = -9,8 +\frac{-2.085-2.0119}{2} = 3.8139 t_2 = t_1 + h = 0.2 + 0.2 = 0.4

iteración 3

K_1 = h\frac{\delta v}{\delta t}(0.4, 3.8139) = 0.2[-9,8 - \frac{0.002}{0.11}( 3.8139)| 3.8139|] = -2.0129 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 3.8139 -2.0129) = 0.2[-9,8 - \frac{0.002}{0.11}(3.8139 -2.0129)|3.8139 -2.0129|] =-1.9718 v_3 = -9,8 +\frac{-2.0129-1.9718}{2} = 1.8215 t_3 = t_2 + h = 0.4 + 0.2 = 0.6

---

Tabla y gráfica del ejercicio para todo el intervalo:

 [xi,      yi,     K1,     K2    ]
[[ 0.      8.      0.      0.    ]
 [ 0.2     5.8623 -2.1927 -2.0826]
 [ 0.4     3.8139 -2.085  -2.0119]
 [ 0.6     1.8215 -2.0129 -1.9718]
 [ 0.8    -0.1444 -1.9721 -1.9599]
 [ 1.     -2.0964 -1.9599 -1.9439]]

Algoritmo en Python

# 3Eva_IT2009_T2 EDO Taylor Seno(x)
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda t,v: -9.8-(0.002/0.11)*v*np.abs(v)
x0 = 0
y0 = 8
h = 0.2
a = 0
b = 1

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h*d1y(xi,yi)
    K2 = h*d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)
# Gráfica
plt.plot(xg,yg)
plt.xlabel('ti')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar iteraciones para Runge-Kutta de 4to Orden

Ejercicio: 2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

La ecuación a resolver con Taylor es:

s'- \frac{26s}{200-t} - \frac{5}{2} = 0

Para lo que se plantea usar la primera derivada:

s'= \frac{26s}{200-t}+\frac{5}{2}

con valores iniciales de s(0) = 0, h=0.1

La fórmula de Taylor para tres términos es:

s_{i+1}= s_{i}+s'_{i}h + \frac{s''_{i}}{2}h^2 + error

Para el desarrollo se compara la solución con dos términos, tres términos y Runge Kutta.

1. Solución con dos términos de Taylor

Iteraciones

i = 0, t₀ = 0, s(0)=0

s'_{0}= \frac{26s_{0}}{200-t_{0}}+\frac{5}{2} = \frac{26(0)}{200-0}+\frac{5}{2} = \frac{5}{2} s_{1}= s_{0}+s'_{0}h = 0+ \frac{5}{2}*0.1= 0.25

t₁ =  t₀+h = 0+0.1 = 0.1

i=1

---

s'_{1}= \frac{26s_{1}}{200-t_{1}}+\frac{5}{2} = \frac{26(0.25)}{200-0.1}+\frac{5}{2} = 2.5325 s_{2}= s_{1}+s'_{1}h = 0.25 + (2.5325)*0.1 = 0.5032

t₂ =  t₁+h = 0.1+0.1 = 0.2

i=2,

resolver como tarea

---

2. Resolviendo con Python

estimado
 [xi,yi Taylor,yi Runge-Kutta, diferencias]
[[ 0.0  0.0000e+00  0.0000e+00  0.0000e+00]
 [ 0.1  2.5000e-01  2.5163e-01 -1.6258e-03]
 [ 0.2  5.0325e-01  5.0655e-01 -3.2957e-03]
 [ 0.3  7.5980e-01  7.6481e-01 -5.0106e-03]
 [ 0.4  1.0197e+00  1.0265e+00 -6.7714e-03]
 [ 0.5  1.2830e+00  1.2916e+00 -8.5792e-03]
 [ 0.6  1.5497e+00  1.5601e+00 -1.0435e-02]
 [ 0.7  1.8199e+00  1.8322e+00 -1.2339e-02]
 [ 0.8  2.0936e+00  2.1079e+00 -1.4294e-02]
 [ 0.9  2.3710e+00  2.3873e+00 -1.6299e-02]
 [ 1.0  2.6519e+00  2.6703e+00 -1.8357e-02]
 [ 1.1  2.9366e+00  2.9570e+00 -2.0467e-02]
 [ 1.2  3.2250e+00  3.2476e+00 -2.2632e-02]
 [ 1.3  3.5171e+00  3.5420e+00 -2.4853e-02]
 [ 1.4  3.8132e+00  3.8403e+00 -2.7129e-02]
 [ 1.5  4.1131e+00  4.1426e+00 -2.9464e-02]
 [ 1.6  4.4170e+00  4.4488e+00 -3.1857e-02]
 [ 1.7  4.7248e+00  4.7592e+00 -3.4310e-02]
 [ 1.8  5.0368e+00  5.0736e+00 -3.6825e-02]
 [ 1.9  5.3529e+00  5.3923e+00 -3.9402e-02]
 [ 2.0  5.6731e+00  5.7152e+00 -4.2043e-02]]
error en rango:  0.04204310894163932

---

2. Algoritmo en Python

# EDO. Método de Taylor 3 términos 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor2t(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) # + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    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] = [xi,yi]
    return(estimado)

# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

# INGRESO.
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: 26*y/(200-x)+5/2
x0 = 0
y0 = 0
h = 0.1
muestras = 20

# PROCEDIMIENTO
puntos = edo_taylor2t(d1y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# Con Runge Kutta
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xiRK2 = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# diferencias
diferencias = yi-yiRK2
error = np.max(np.abs(diferencias))
tabla = np.copy(puntos)
tabla = np.concatenate((puntos,np.transpose([yiRK2]),
                        np.transpose([diferencias])),
                       axis = 1)

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi Taylor,yi Runge-Kutta, diferencias]')
print(tabla)
print('error en rango: ', error)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o',
         color='r', label ='[x0,y0]')
plt.plot(xi[0:],yi[0:],'-',
         color='g',
         label ='y Taylor 2 términos')
plt.plot(xiRK2[0:],yiRK2[0:],'-',
         color='blue',
         label ='y Runge-Kutta 2Orden')
plt.axhline(y0/2)
plt.title('EDO: Taylor 2T vs Runge=Kutta 2Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

---

Usando Taylor con 3 términos

estimado
 [xi,        yi,        d1yi,      d2yi      ]
[[0.         0.         2.5        0.325     ]
 [0.1        0.251625   2.53272761 0.32958302]
 [0.2        0.50654568 2.56591685 0.33423301]
 [0.3        0.76480853 2.59957447 0.33895098]
 [0.4        1.02646073 2.63370731 0.34373796]
 [0.5        1.29155015 2.66832233 0.348595  ]
 [0.6        1.56012536 2.70342658 0.35352316]
 [0.7        1.83223563 2.73902723 0.35852351]
 [0.8        2.10793097 2.77513155 0.36359715]
 [0.9        2.38726211 2.81174694 0.36874519]
 [1.         2.67028053 2.84888087 0.37396876]
 [1.1        2.95703846 2.88654098 0.37926901]
 [1.2        3.24758891 2.92473497 0.3846471 ]
 [1.3        3.54198564 2.96347069 0.39010422]
 [1.4        3.84028323 3.00275611 0.39564157]
 [1.5        4.14253705 3.04259931 0.40126036]
 [1.6        4.44880328 3.08300849 0.40696184]
 [1.7        4.75913894 3.12399199 0.41274727]
 [1.8        5.07360187 3.16555827 0.41861793]
 [1.9        5.39225079 3.2077159  0.42457511]
 [2.         5.71514526 0.         0.        ]]