[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
..
1. EDO \frac{\delta^2 y}{\delta x^2} con Runge-Kutta
Para una ecuación diferencial de segunda derivada (segundo orden) con condiciones de inicio en x0, y0, y'0
\frac{\delta ^2 y}{\delta x^2} = \frac{\delta y}{\delta x} + etcForma estandarizada de la ecuación:
y'' = y' + etcSe puede sustituir la variable y' por z, lo que se convierte a dos expresiones que forman un sistema de ecuaciones:
\begin{cases} z= y' = f_x(x,y,z) \\ z' = (y')' = z + etc = g_x(x,y,z) \end{cases}y se pueden reutilizar los métodos para primeras derivadas, por ejemplo Runge-Kutta de 2do y 4to orden para las variables x,y,z de forma simultanea.
Runge-Kutta 2do Orden tiene error de truncamiento O(h3)
Runge-Kutta 4do Orden tiene error de truncamiento O(h5)
Al aplicar Runge-Kutta para las variables dependientes, se tiene:
y'' = y' + etc \begin{cases} f_x(x,y,z) = z \\ g_x(x,y,z) = z + etc \end{cases} K_{1y} = h f(x_i, y_i, z_i) K_{1z} = hg(x_i, y_i, z_i) K_{2y} = h f(x_i +h, y_i + K_{1y} , z_i + K_{1z}) K_{2z} = h g(x_i +h, y_i + K_{1y}, z_i + K_{1z}) y_{i+1}=y_i+\frac{K_{1y}+K_{2y}}{2} z_{i+1}=z_i+\frac{K_{1z}+K_{2z}}{2} x_{i+1} = x_i +h[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
2. Ejercicio
Referencia: Chapra Ejercicio 25.23 p265, 2Eva_IT2018_T1 Paracaidista wingsuit 
Si suponemos que el arrastre es proporcional al cuadrado de la velocidad, se puede modelar la altura (y) de un objeto que cae, como un paracaidista, por medio de la ecuación diferencial ordinaria:
\frac{d^2y}{dt^2} = g - \frac{Cd}{m} \Big( \frac{dy}{dt} \Big) ^2Que es una EDO de 2do orden o como 2da derivada.
Resuelva para la altura que recorre un objeto de 90 Kg con coeficiente de arrastre Cd =0.225 kg/m.
Si la velocidad inicial es 0 y la altura inicial es 1 Km, determine la velocidad y posición en cada tiempo, usando un tamaño de paso de 2s.
[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
3. Desarrollo analítico
La solución propone resolver de forma simultanea para t,y,v con Runge Kutta para segunda derivada donde:
Al sustituir los valores de las constantes en la ecuación como gravedad, masa e índice de arrastre se tiene:
f(t,y,v) = -v g(t,y,v) = 9.8 - \frac{0.225}{90} v^2con las condiciones iniciales del ejercicio t0 = 0 , y0 = 1000, v0 = 0
la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.
Para las iteraciones, recuerde que
t se corrige con t+h (en el algoritmo era la posición para x)
y se corrige con y+K1y
v se corrige con v+K1v (en el algoritmo era la posición para z)
itera = 0
K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6..
K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2
..
y_1 = y_0 + \frac{K1y+K2y}{2} = 1000 + \frac{0-39.2}{2}= 980.4
| ti | yi | vi | K1y | K1v | K2y | K2v |
| 0 | 1000 | 0 | 0 | 0 | 0 | 0 |
| 2 | 980.4 | 18.63 | 0 | 19.6 | -39.2 | 17.6792 |
itera = 1
K1y = 2(-(18.63)) = -37.2792 K1v = 2(9.8 - \frac{0.225}{90} (18.63)^2) = 17.8628..
K2y =2(-(18.6396+17.8628)) =-73.00
..
y_2 =980.4 + \frac{ -37.2792+(-73.00)}{2}= 925.25
| ti | yi | vi | K1y | K1v | K2y | K2v |
| 0 | 1000 | 0 | 0 | 0 | 0 | 0 |
| 2 | 980.4 | 18.63 | 0 | 19.6 | -39.2 | 17.6792 |
| 4 | 925.25 | 34.0399 | -37.2792 | 17.8628 | -73.00 | 12.9378 |
itera = 2
K1y = h f(t,y,v) = 2(-(34.0399)) = -68.0798 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (34.0399)^2) = 13.8064..
K2y = h f(t+h, y+K1y, v + K1v) = 2(-(34.0399+13.8064)) =-95.6927
..
y_2 = 925.25 + \frac{ -68.0798+(-95.6927)}{2}= 843.3716
| ti | yi | vi | K1y | K1v | K2y | K2v |
| 0 | 1000 | 0 | 0 | 0 | 0 | 0 |
| 2 | 980.4 | 18.63 | 0 | 19.6 | -39.2 | 17.6792 |
| 4 | 925.25 | 34.0399 | -37.2792 | 17.8628 | -73.00 | 12.9378 |
| 6 | 843.3716 | 45.0199 | -68.0798 | 13.8064 | -95.6927 | 8.1536 |
Usando el algoritmo se obtiene el siguiente resultado:
EDO f,g con Runge-Kutta 2 Orden i [ xi, yi, zi ] [ K1y, K1z, K2y, K2z ] 0 [ 0. 1000. 0.] [0. 0. 0. 0.] 1 [ 2. 980.4 18.6396] [ 0. 19.6 -39.2 17.6792] 2 [ 4. 925.258 34.0399] [-37.2792 17.8628 -73.0049 12.9379] 3 [ 6. 843.3717 45.02 ] [-68.0799 13.8064 -95.6927 8.1536] 4 [ 8. 743.8657 52.1312] [ -90.0399 9.466 -108.972 4.7564] 5 [ 10. 633.5917 56.4855] [-104.2623 6.0117 -116.2857 2.697 ] 6 [ 12. 516.9737 59.0692]
con la siguiente gráfica
[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
4. Algoritmo en Python para \frac{\delta^2 y}{\delta x^2} Runge-Kutta 2do Orden
Se presenta las instrucciones en Python con una función.
# EDO d2y/dx2. Método de RungeKutta 2do Orden # estima la solucion para muestras espaciadas h en eje x # valores iniciales x0,y0,z0 entrega tabla[xi,yi,zi,K1y,K1z,K2y,K2z] import numpy as np def rungekutta2_fg(f,g,x0,y0,z0,h,muestras, vertabla=False, precision=6): ''' solucion a EDO d2y/dx2 con Runge-Kutta 2do Orden, f(x,y,z) = z #= y' g(x,y,z) = expresion d2y/dx2 con z=y' tambien es solucion a sistemas edo f() y g() x0,y0,z0 son valores iniciales, h es tamano de paso, muestras es la cantidad de puntos a calcular. ''' tamano = muestras + 1 tabla = np.zeros(shape=(tamano,3+4),dtype=float) # incluye el punto [x0,y0,z0,K1y,K1z,K2y,K2z] tabla[0] = [x0,y0,z0,0,0,0,0] xi = x0 # valores iniciales yi = y0 zi = z0 for i in range(1,tamano,1): K1y = h * f(xi,yi,zi) K1z = h * g(xi,yi,zi) K2y = h * f(xi+h, yi + K1y, zi + K1z) K2z = h * g(xi+h, yi + K1y, zi + K1z) yi = yi + (K1y+K2y)/2 zi = zi + (K1z+K2z)/2 xi = xi + h tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z] if vertabla==True: np.set_printoptions(precision) print('EDO f,g con Runge-Kutta 2 Orden') print('i ','[ xi, yi, zi',']') print(' [ K1y, K1z, K2y, K2z ]') for i in range(0,tamano,1): txt = ' ' if i>=10: txt = ' ' print(str(i),tabla[i,0:3]) print(txt,tabla[i,3:]) return(tabla) # PROGRAMA PRUEBA ------------------- # 2Eva_IT2018_T1 Paracaidista wingsuit # INGRESO f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre g = lambda t,y,v: 9.8 - (0.225/90)*(v**2) t0 = 0 y0 = 1000 v0 = 0 h = 2 muestras = 15+1 # PROCEDIMIENTO tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras, vertabla=True, precision=4) ti = tabla[:,0] yi = tabla[:,1] vi = tabla[:,2] # SALIDA # print('tabla de resultados') # print(tabla) # GRAFICA import matplotlib.pyplot as plt plt.subplot(211) plt.plot(ti,vi,label='velocidad v(t)', color='green') plt.plot(ti,vi,'o') plt.ylabel('velocidad (m/s)') plt.title('paracaidista Wingsuit con Runge-Kutta') plt.legend() plt.grid() plt.subplot(212) plt.plot(ti,yi,label='Altura y(t)',) plt.plot(ti,yi,'o',) plt.axhline(0, color='red') plt.xlabel('tiempo (s)') plt.ylabel('Altura (m)') plt.legend() plt.grid() plt.show()
[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]
..
5. Algoritmo en Python para \frac{\delta^2 y}{\delta x^2} Runge-Kutta 4to Orden
Se adjunta la función para 4to orden que se puede probar con el mismo ejercicio anterior.
# EDO d2y/dx2. Método de RungeKutta 4to Orden # estima la solucion para muestras espaciadas h en eje x # valores iniciales x0,y0,z0 entrega # tabla[xi,yi,zi,K1y,K1z,K2y,K2z,K3y,K3z,K4y,K4z] import numpy as np def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras, vertabla=False, precision=6): ''' solucion a EDO d2y/dx2 con Runge-Kutta 4to Orden, f(x,y,z) = z #= y' g(x,y,z) = expresion d2y/dx2 con z=y' tambien es solucion a sistemas edo f() y g() x0,y0,z0 son valores iniciales, h es tamano de paso, muestras es la cantidad de puntos a calcular. ''' tamano = muestras + 1 tabla = np.zeros(shape=(tamano,3+8),dtype=float) # incluye el punto [x0,y0] tabla[0] = [x0,y0,z0,0,0,0,0,0,0,0,0] xi = x0 # valores iniciales yi = y0 zi = z0 for i in range(1,tamano,1): K1y = h * fx(xi,yi,zi) K1z = h * gx(xi,yi,zi) K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2) K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2) K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2) K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2) K4y = h * fx(xi+h, yi + K3y, zi + K3z) K4z = h * gx(xi+h, yi + K3y, zi + K3z) yi = yi + (K1y+2*K2y+2*K3y+K4y)/6 zi = zi + (K1z+2*K2z+2*K3z+K4z)/6 xi = xi + h tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z,K3y,K3z,K4y,K4z] if vertabla==True: np.set_printoptions(precision) print('EDO f,g con Runge-Kutta 4 Orden') print('i ','[ xi, yi, zi',']') print(' [ K1y, K1z, K2y, K2z ]') print(' [ K3y, K3z, K4y, K4z ]') for i in range(0,tamano,1): txt = ' ' if i>=10: txt = ' ' print(str(i),tabla[i,0:3]) print(txt,tabla[i,3:7]) print(txt,tabla[i,7:]) return(tabla)
[ Runge Kutta \frac{\delta^2 y}{\delta x^2} ] [ Ejercicio ] [ analítico ]
Algoritmo: [ RK 2do Orden ] [ RK 4to Orden]