s2Eva_IIT2019_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

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

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

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

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 usando 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())