Ejercicio: 2Eva2019TII_T2 EDO, problema de valor inicial
la ecuación del problema planteado es:
y'(t) = f(t,y) = \frac{y}{2t^3}0 ≤ t ≤ 1
y(0.5) = 1.5
literal a
La solución empieza usando la Serie de Taylor por ejemplo para tres términos:
y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)Se observa que se tiene el valor inicial y la primera derivada, si usamos tres términos se puede usar la segunda derivada.
y'(t) = f(t,y) = \frac{y}{2t^3} y''(t) = f'(t,y) = \frac{y'}{2t^3} + y \Big(\frac{-3}{2t^4}\Big) y''(t) = \frac{1}{2t^3}\Big( \frac{y}{2t^3} \Big) - \frac{3y}{2t^4} y''(t) = \frac{y}{4t^6} -\frac{3y}{2t^4}por lo que la ecuación de Taylor a usar queda de la siguiente forma:
y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i y_{i+1} = y_{i} + h\frac{y}{2t^3}+\frac{h^2}{2} \Big( \frac{y}{4t^6} -\frac{3y}{2t^4}\Big)que es la ecuacion que se usará con un error de O(h3)
Reemplazando los valores en la fórmula se obtiene la siguiente tabla:
estimado
[ti, yi, d1yi, d2yi]
[[ 0.5 1.5 6. -12. ]
[ 0.6 2.04 4.7222 -12.68004115]
[ 0.7 2.4488 3.5697 -10.09510219]
[ 0.8 2.7553 2.6907 -7.46259891]
[ 0.9 2.9870 2.0487 -5.42399041]
[ 1. 3.1648 0. 0. ]]cuya gráfica es:
literal b
Para desarrollar Runge-Kutta de 2do orden se dispone de los siguientes datos:
y'(t) = f(t,y) = \frac{y}{2t^3}t0 = 0.5, y0 = 1.5, h = 0.1
pasos del algoritmo,
K_1 = h * y'(t_i) K_2 = h * y'(t_i+h, y_i + K_1) y_{i+1} = y_i + \frac{K_1+K_2}{2} t_i = t_i + h</p> <p>iteración 1: i = 0
K_1 = 0.1 * y'(0.5) = 0.6 K_2 = 0.1 * y'(0.5+0.1, 1.5 + 0.6) = 0.4861 y_{1} = 1.5 + \frac{0.6+0.4861}{2} = 2.0430 t_1 = 0.5 + 0.1 = 0.6</p> <p>iteración 2: i = 1
K_1 = 0.1 * y'(0.6) = 0.4729 K_2 = 0.1 * y'(0.6+0.1, 2.0430 + 0.4729) = 0.3667 y_{1} = 2.0430 + \frac{0.4729+0.3667}{2} = 2.4629 t_1 = 0.6 + 0.1 = 0.7</p> <p>iteración 3: i = 2
K_1 = 0.1 * y'(0.7) = 0.3590 K_2 = 0.1 * y'(0.7+0.1, 2.4629 + 0.3590) = 0.3667 y_{1} = 2.4629 + \frac{0.3590+0.3667}{2} = 2.7802 t_1 = 0.7 + 0.1 = 0.8</p> <p>obteniendo la siguiente tabla:
estimado
[ti, yi, K1, K2]
[[0.5 1.5 0. 0. ]
[0.6 2.0430 0.6 0.4861]
[0.7 2.4629 0.4729 0.3667]
[0.8 2.7802 0.3590 0.2755]
[0.9 3.0206 0.2715 0.2093]
[1. 3.2048 0.2071 0.1613]]
Diferencias entre Taylor y Runge-Kutta2
[ 0. -0.0030 -0.0140 -0.0248 -0.0335 -0.0400]Algoritmo en Python
# EDO. Método de Taylor 3 términos
# Runge-Kutta de 2 Orden
# 2Eva_IIT2019_T2 EDO, problema de valor inicial
import numpy as np
# Funciones desarrolladas en clase
def edo_taylor3t(d1y,d2y,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]
x = x0
y = y0
for i in range(1,tamano,1):
estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
x = x+h
estimado[i,0:2] = [x,y]
return(estimado)
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] = [xi,yi,K1,K2]
return(estimado)
# PROGRAMA PRUEBA
# INGRESO
# d1y = y', d2y = y''
d1y = lambda t,y: y/(2*t**3)
d2y = lambda t,y: y/(4*t**6)-(3/2)*y/(t**4)
t0 = 0.5
y0 = 1.5
h = 0.1
muestras = 5
# PROCEDIMIENTO
# Edo con Taylor
puntos = edo_taylor3t(d1y,d2y,t0,y0,h,muestras)
ti = puntos[:,0]
yi = puntos[:,1]
# Runge-Kutta
puntosRK2 = rungekutta2(d1y,t0,y0,h,muestras)
# ti = puntosRK2[:,0] # lo mismo del anterior
yiRK2 = puntosRK2[:,1]
# diferencias
diferencia = yi-yiRK2
# SALIDA
print('estimado[ti, yi, d1yi, d2yi]')
print(puntos)
print('estimado[ti, yi, K1, K2]')
print(puntosRK2)
print('Diferencias entre Taylor y Runge-Kutta2')
print(diferencia)
# Gráfica
import matplotlib.pyplot as plt
plt.plot(ti[0],yi[0],'o', color='r',
label ='[t0,y0]')
plt.plot(ti[1:],yi[1:],'o', color='g',
label ='y con Taylor 3 términos')
plt.plot(ti[1:],yiRK2[1:],'o', color='m',
label ='y Runge-Kutta 2 Orden')
plt.title('EDO: Solución numérica')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()