6.2 EDO Runge-Kutta 2do Orden dy/dx

[ Runge Kutta dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

1. EDO Runge-Kutta 2do Orden para primera derivada dy/dx

Referencia: Burden 5.4 p209, Chapra 25.3 p740, Rodríguez 9.1.7 p354, Boyce DiPrima 4Ed 8.4 p450

Los métodos de Runge-Kutta mejoran la aproximación a la respuesta sin requerir determinar las expresiones de las derivadas de orden superior, como fué necesario en EDO con Taylor.

Para una ecuación diferencial de primera derivada (primer orden) con una condición de inicio: $\frac{\delta y}{\delta x} + etc =0$

y'(x) = f(x_i,y_i)

y(x_0) = y_0

El método Runge-Kutta de 2do Orden usa una corrección a la derivada tomando valores de puntos xi y xi+h, un tamaño de paso h hacia adelante, calculados como términos K₁ y K₂.

K_1 = h f(x_i,y_i)

K_2 = h f(x_i+h, y_i + K_1)

y_{i+1} = y_i + \frac{K_1+K_2}{2}

x_{i+1} = x_i + h

Runge-Kutta 2do Orden tiene error de truncamiento O(h³)

Las iteraciones se repiten para encontrar el siguiente punto en X[i+1] como se muestra en el gráfico animado siguiente:

Como tema de introducción observar dos minutos del video sugerido a partir de donde se encuentra marcado el enlace. En este caso, en combate aéreo, las armas se encuentran fijas en las alas.

Video Revisar:

Luego de observar el video de introducción conteste las siguientes preguntas:
¿ Que trayectoria siguen los proyectiles al salir del cañón?
¿ Que trayectoria siguen los aviones, el perseguido y el que caza?
¿ Cuál es la relación entre las trayectorias de los dos aviones?

[ Runge Kutta dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

2. Ejercicio

Referencia: Rodríguez 9.1.1 ejemplo p335. Chapra 25.1.3 p731

Se pide encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x₀=0, y₀=1

\frac{dy}{dx}-y -x +x^2 -1 = 0

y'-y -x +x^2 -1 = 0

[ Runge Kutta dy/dx] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

3. Desarrollo Analítico

Se reordena la expresión de forma que la derivada se encuentre en el lado izquierdo:

f(x,y) = y' = y -x^2 +x +1

Se usa las expresiones de Runge-Kutta en orden, K1 corresponde a una corrección de EDO con Taylor de dos términos (método de Euler). K2 considera el cálculo a un tamaño de paso más adelante. iteración:

K_1 = h f(x_i,y_i)

K_2 = h f(x_i+h, y_i + K_1)

y_{i+1} = y_i + \frac{K_1+K_2}{2}

x_{i+1} = x_i + h

Runge-Kutta 2do Orden tiene error de truncamiento O(h³)

para el ejercicio, el tamaño de paso h=0.1

itera = 0 , x0 = 0, y0 = 1

K_1 = 0.1 f(0,1) = 0.1 \Big( 1 -0^2 +0 +1 \Big) = 0.2

K_2 = 0.1 f(0+0.1, 1+ 0.2)

K_2 = 0.1 \Big( (1+ 0.2) - (0+0.1) ^2 +(0+0.1) +1\Big) = 0.229

y_1 = 1 + \frac{0.2+0.229}{2} = 1.2145

x_1 = 0 + 0.1 = 0.1

itera = 1 , x0 = 0.1, y0 = 1.2145

K_1 = 0.1 f(0.1,1.2145) = 0.1( 1.2145 -0.1^2 +0.1 +1) = 0.23045

K_2 = 0.1 f(0.1+0.1, 1.2145 + 0.23045)

K_2 = 0.1 \Big((1.2145 + 0.23045) -(0.1+0.1)^2 +(0.1+0.1) +1\Big) = 0.260495

y_2 = 1.2145 + \frac{0.23045+0.260495}{2} = 1.4599725

x_2 = 0.1 +0.1 = 0.2

itera = 2 , x0 = 0.2, y0 = 1.4599725

completar como tarea …

completando los puntos con el algoritmo y realizando la gráfica se obtiene

EDO con Runge-Kutta 2 Orden
 [xi, yi, K1, K2]
[[0.         1.         0.         0.        ]
 [0.1        1.2145     0.2        0.229     ]
 [0.2        1.4599725  0.23045    0.260495  ]
 [0.3        1.73756961 0.26199725 0.29319698]
 [0.4        2.04856442 0.29475696 0.32723266]
 [0.5        2.39436369 0.32885644 0.36274209]]
>>>

Compare los resultados con Taylor de 2 y 3 términos.

Runge-Kutta 2do Orden tiene error de truncamiento O(h³)

[ Runge Kutta dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

4. Algoritmo en Python

Se adjunta el programa de prueba que usa la función rungekutta2(d1y,x0,y0,h,muestras) :

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

def rungekutta2(d1y,x0,y0,h,muestras):
    # Runge Kutta de 2do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    
    # incluye el punto [x0,y0]
    tabla[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 + (1/2)*(K1+K2)
        xi = xi + h
        
        tabla[i] = [xi,yi,K1,K2]
    return(tabla)
# 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'
d1y = lambda x,y: y -x**2 + x + 1
x0 = 0
y0 = 1
h  = 0.1
muestras = 5

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

# SALIDA
# np.set_printoptions(precision=4)
print( 'EDO con Runge-Kutta 2 Orden')
print(' [xi, yi, K1, K2]')
print(tabla)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi,yiRK2)
plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='[xi,yi] 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() #comentar para la siguiente gráfica

[ Runge Kutta dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

5. Cálculo de Error con la solución conocida

La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.

y = e^x + x + x^2

Note que el error crece al distanciarse del punto inicial.

Error máximo estimado:  0.004357584597315167
entre puntos: 
[0.         0.00067092 0.00143026 
 0.0022892  0.00326028 0.00435758]
>>>

Para las siguientes instrucciones, comente la última línea #plt.show() antes de continuar con:

# ERROR vs solución conocida -----------------
y_sol = lambda x: ((np.e)**x) + x + x**2

yi_psol  = y_sol(xi)
errores  = yi_psol - yiRK2
errormax = np.max(np.abs(errores))

# SALIDA
print('Error máximo estimado: ',errormax)
print('entre puntos: ')
print(errores)

# GRAFICA [a,b+2*h]
a = x0
b = h*muestras+2*h
muestreo = 10*muestras+2
xis = np.linspace(a,b,muestreo)
yis = y_sol(xis)

plt.plot(xis,yis, label='y solución conocida',
         linestyle='dashed')
plt.legend()
plt.show()

[ Runge Kutta dy/dx ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

6. Ejercicio de Evaluación:

2Eva_IT2018_T1 Paracaidista wingsuit

Solución Propuesta: s2Eva_IT2018_T1 Paracaidista wingsuit , Runge-Rutta para primera derivada.