Ejercicio: 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua
La solución obtenida se realiza con h=0.5 y usando dos métodos para comparar resultados.
\frac{dy}{dt} = -k \sqrt{y}
1. EDO con Taylor
Usando una aproximación con dos términos de Taylor:
y_{i+1}=y_{i}+ y'_{i} h+\frac{y"_{i}}{2}h^{2}
Por lo que se obtienen las derivadas necesarias:
y'_i= -k (y_i)^{1/2}
y"_i= \frac{-k}{2}(y_i)^{-1/2}
1.1 iteraciones
i=0, y0=3, t0=0
y'_0= -k(y_0)^{1/2} =-0.06(3)^{1/2} = -0.1039
y"_0= \frac{-0.06}{2}(3)^{-1/2} = -0.0173
y_{1}=y_{0}+ y'_{0} (1)+\frac{y"_{0}}{2}(1)^{2}
y_{1}=3+ (-0.1039) (0.5)+\frac{-0.0173}{2}(0.5)^{2}= 2.9458
t1=t0+h = 0+0.5= 0.5
i=1, y1=2.9458, t1=0.5
y'_1= -k(y_1)^{1/2} =-0.06(2.887)^{1/2} =-0.1029
y"_1= \frac{-0.06}{2}(2.887)^{-1/2} = -0.0174
y_{2}=y_{1}+ y'_{1} (1)+\frac{y"_{1}}{2}(1)^{2}
y_{1}=2.9458+ (-0.1029) (1)+\frac{-0.0174}{2}(1)^{2}= 2.8921
t2=t1+h = 0.5+0.5 = 1.0
i=2, y2=2.8921, t2=1.0
Resolver como Tarea
1.2 Resultados con Python
Realizando una tabla de valores usando Python y una gráfica, encuentra que el valor buscado del tanque a la mitad se obtiene en 16 minutos.
estimado[xi,yi]
[[ 0. 3. ]
[ 0.5 2.94587341]
[ 1. 2.89219791]
[ 1.5 2.83897347]
[ 2. 2.7862001 ]
...
[14. 1.65488507]
[14.5 1.61337731]
[15. 1.57231935]
[15.5 1.53171109]
[16. 1.49155239]
[16.5 1.45184313]
[17. 1.41258317]
[17.5 1.37377234]
[18. 1.33541049]
[18.5 1.29749744]
[19. 1.26003297]
[19.5 1.22301689]
[20. 1.18644897]]
Algoritmo en Python para Solución EDO con tres términos:
# 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_taylor3t(d1y,d2y,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)
# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua
# INGRESO.
k=0.06
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: -k*(y**0.5)
d2y = lambda x,y: -(k/2)*(y**(-0.5))
x0 = 0
y0 = 3
h = 1/2
muestras = 40
# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]
# SALIDA
print('estimado[xi,yi]')
print(puntos)
# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.axhline(y0/2)
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()
2. EDO con Runge-Kutta de 2do Orden dy/dx
Para éste método no se requiere desarrollar la segunda derivada, se usa el mismo h =0.5 con fines de comparación de resultados
2.1 ITeraciones
i = 1, y0=3, t0=0
K_1 = h y'(x_0,y_0) = (0.5)*(-0.06)(3)^{1/2} =-0.05196
K_2 = h y'(x_0+h,y_0+K_1) = (0.5)* y'(0.5,3-0.05196) = -0.05150
y_1 = y_0+\frac{K1+K2}{2} = 3+\frac{-0.05196-0.05150}{2} = 2.9482
i = 2, y1=2.9482, t1=0.5
K_1 = h y'(x_1,y_1) = (0.5)*(-0.06)(2.9482)^{1/2} =-0.05149
K_2 = h y'(x_1+h,y_1+K_1) = (0.5)* y'(0.5,2.9482-0.05149) = -0.05103
y_1 = y_0+\frac{K1+K2}{2} = 3+\frac{-0.05149-0.05103}{2} = -2.8946
i = 3, y1=2.8946, t1=1.0
Resolver como Tarea
2.2 Resultados con Python
Si comparamos con los resultados anteriores en una tabla, y obteniendo las diferencias entre cada iteración se tiene que:
estimado[xi,yi Taylor, yi Runge-Kutta, diferencias]
[[ 0.0 3.00000000 3.00000000 0.00000000e+00]
[ 0.5 2.94587341 2.94826446 -2.39104632e-03]
[ 1.0 2.89219791 2.89697892 -4.78100868e-03]
[ 1.5 2.83897347 2.84614338 -7.16990106e-03]
[ 2.0 2.78620010 2.79575783 -9.55773860e-03]
...
[ 14.0 1.65488507 1.72150488 -6.66198112e-02]
[ 14.5 1.61337731 1.68236934 -6.89920328e-02]
[ 15.0 1.57231935 1.64368380 -7.13644510e-02]
[ 15.5 1.53171109 1.60544826 -7.37371784e-02]
[ 16.0 1.49155239 1.56766273 -7.61103370e-02]
[ 16.5 1.45184313 1.53032719 -7.84840585e-02]
[ 17.0 1.41258317 1.49344165 -8.08584854e-02]
[ 17.5 1.37377234 1.45700611 -8.32337718e-02]
[ 18.0 1.33541049 1.42102058 -8.56100848e-02]
[ 18.5 1.29749744 1.38548504 -8.79876055e-02]
[ 19.0 1.26003297 1.35039950 -9.03665304e-02]
[ 19.5 1.22301689 1.31576397 -9.27470733e-02]
[ 20.0 1.18644897 1.28157843 -9.51294661e-02]]
error en rango: 0.09512946613018003
# 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_taylor3t(d1y,d2y,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.
k=0.06
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: -k*(y**0.5)
d2y = lambda x,y: -(k/2)*(y**(-0.5))
x0 = 0
y0 = 3
h = 1/2
muestras = 40
# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,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
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[1:],yi[1:],'o',
color='g',
label ='y Taylor 3 terminos')
plt.plot(xiRK2[1:],yiRK2[1:],'o',
color='blue',
label ='y Runge-Kutta 2Orden')
plt.axhline(y0/2)
plt.title('EDO: Taylor 3T vs Runge=Kutta 2Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()
