1. EDO  \frac{\delta y}{\delta x} Runge-Kutta 4to Orden

Referencia: Chapra 25.3.3 p746, Rodríguez 9.1.8 p358

Para una ecuación diferencial de primera derivada (primer orden) con una condición de inicio:

La fórmula de Runge-Kutta de 4to orden realiza una corrección con 4 valores de K:

K_1 = h f(x_i,y_i) K_2 = h f\Big(x_i+\frac{h}{2}, y_i + \frac{K_1}{2} \Big) K_3 = h f\Big(x_i+\frac{h}{2}, y_i + \frac{K_2}{2} \Big) K_4 = h f(x_i+h, y_i + K_3 ) y_{i+1} = y_i + \frac{K_1 + 2K_2 + 2K_3 + K_4}{6} x_{i+1} = x_i + h

debe ser equivalente a la serie de Taylor de 5 términos:

Runge-Kutta 4do Orden tiene error de truncamiento O(h⁵)

2. Ejercicio

Para el desarrollo analítico se tienen las siguientes expresiones para el ejercicio usado en Runge-Kutta de orden 2, que ahora será con orden 4:

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 medio tamaño de paso más adelante.

iteración:

Se inicia la tabla con las condiciones iniciales en primera fila:

itera=0

itera=1

Las iteraciones con los valores desde i=2 se dejan como tarea

Los resultados del algoritmo se presentan en la tabla:

EDO dy/dx con Runge-Kutta 4to Orden
i, [xi,     yi,     K1,    K2,     K3,     K4 ]
0 [0. 1. 0. 0. 0. 0.]
1 [0.1        1.21517063 0.2        0.21475    0.2154875  0.23054875]
2 [0.2        1.46140213 0.23051706 0.24579292 0.24655671 0.26217273]
3 [0.3        1.7398578  0.26214021 0.27799722 0.27879007 0.29501922]
4 [0.4        2.05182327 0.29498578 0.31148507 0.31231003 0.32921678]
5 [0.5        2.39871935 0.32918233 0.34639144 0.3472519  0.36490752]

3. Algoritmo en Python como Función

# EDO dy/dx. Metodo de RungeKutta 4to Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0, entrega tabla[xi,yi,K1,K2,K3,K4]
import numpy as np
   
# INGRESO
# d1y = y' = f
d1y = lambda x,y: y -x**2 + x + 1
x0 = 0 # condiciones iniciales
y0 = 1
h  = 0.1
muestras = 5
   
# algoritmos como funcion
def rungekutta4(d1y,x0,y0,h,muestras, vertabla=False, precision=6):
    ''' solucion a EDO con Runge-Kutta 4do Orden primera derivada,
        x0,y0 son valores iniciales, tamaño de paso h.
        muestras es la cantidad de puntos a calcular.
    '''
    # Runge Kutta de 4do orden
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+4),dtype=float)
      
    # incluye el punto [x0,y0,K1,K2,K3,K4]
    tabla[0] = [x0,y0,0,0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h/2, yi + K1/2)
        K3 = h * d1y(xi+h/2, yi + K2/2)
        K4 = h * d1y(xi+h, yi + K3)
  
        yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
        xi = xi + h
          
        tabla[i] = [xi,yi,K1,K2,K3,K4]
          
    if vertabla==True:
        np.set_printoptions(precision)
        print(' EDO con Runge-Kutta 4do Orden primera derivada')
        print('i, [xi,     yi,     K1,    K2,     K3,     K4 ]')
        for i in range(0,tamano,1):
            print(i,tabla[i])
    return(tabla)
   
# PROCEDIMIENTO
tabla = rungekutta4(d1y,x0,y0,h,muestras)
n = len(tabla)
   
# SALIDA
print('EDO dy/dx con Runge-Kutta 4to Orden')
print('i, [xi,     yi,     K1,    K2,     K3,     K4 ]')
for i in range(0,n,1):
    print(i,tabla[i])

Note que el método de Runge-Kutta de 4to orden es similar a la regla de Simpson 1/3. La ecuación representa un promedio ponderado para establecer la mejor pendiente.

4. 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.

Error máximo estimado: 1.917156652542218e-06
entre puntos:
[0.00000000e+00 2.93075648e-07 6.25886732e-07 1.00354959e-06 1.43181723e-06 1.91715665e-06]

5. Ejercicio en video

2Eva2018TI_T1 Paracaidista wingsuit

Solución Propuesta: s2Eva2018TI_T1 Paracaidista wingsuit

La segunda parte corresponde a Runge-Kutta de 4to Orden