6.3 EDO d2y/dx2 Runge-Kutta con Python

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) ^2

Que 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}$

Algoritmo:

3. Desarrollo analítico - paso a paso

La solución propone resolver de forma simultanea para t,y,v con Runge Kutta para segunda derivada, siendo las expresiones:

f(t,y,v) = -v

g(t,y,v) = g - \frac{Cd}{m} v^2

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^2

con las condiciones iniciales del ejercicio t₀ = 0 , y₀ = 1000, v₀ = 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

K2v = h g(t+h, y+K1y, v + K1v)

= 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792

y_1 = y_0 + \frac{K1y+K2y}{2}

= 1000 + \frac{0-39.2}{2}= 980.4

v_1 = v_0 + \frac{K1v+K2v}{2} = 0 + \frac{19.6-17.67}{2} = 18.63

t_1 =t_0 + h = 0+2 = 2

EDO Runge-Kutta 2Orden d2y/dx2 itera=0 gráfica

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	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

K2v = 2(9.8 - \frac{0.225}{90} (18.6396+17.8628)^2) =12.9378

y_2 =980.4 + \frac{ -37.2792+(-73.00)}{2}= 925.25

v_2 = 18.63 + \frac{17.8628+12.9378}{2} = 34.0399

t_2 =t_1 + h = 2+2 = 4

EDO Runge-Kutta 2Orden d2y/dx2 itera=1 gráfica

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	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

K2v = h g(t+h, y+K1y, v + K1v)

= 2(9.8 - \frac{0.225}{90} (34.0399+13.8064)^2) =8.1536

y_2 = 925.25 + \frac{ -68.0798+(-95.6927)}{2}= 843.3716

v_2 = 34.0399 + \frac{13.8064+8.1536}{2} = 45.0199

t_2 =t_1 + h = 4+2 = 6

EDO Runge-Kutta 2Orden d2y/dx2 itera=2 gráfica

ti	yi	vi	K1y	K1v	K2y	K2v
0	1000	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.37	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]
...
13 [  26.     -348.3079   62.5395]
   [-1.2497e+02  7.7276e-02 -1.2513e+02  2.8959e-02]
14 [  28.     -473.4309   62.5698]
   [-1.2508e+02  4.4071e-02 -1.2517e+02  1.6499e-02]
15 [  30.     -598.5956   62.587 ]
   [-1.2514e+02  2.5126e-02 -1.2519e+02  9.4015e-03]

con la siguiente gráfica

Runge Kutta $\frac{\delta^2 y}{\delta x^2}$

Algoritmo:

4. Algoritmo en Python para $\frac{\delta^2 y}{\delta x^2}$ Runge-Kutta 2do Orden

Se presenta las instrucciones en Python del algoritmo para segunda derivada 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

# INGRESO
# 2Eva_IT2018_T1 Paracaidista wingsuit
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 = 14+1

# Algoritmo como función
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)

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras,
                       vertabla=True, precision=4)
# SALIDA
# print('tabla de resultados')
# print(tabla)

Runge Kutta $\frac{\delta^2 y}{\delta x^2}$

Algoritmo:

5. Gráfica

Instrucciones adicionales al algoritmo anterior

# GRAFICA ---------------------
import matplotlib.pyplot as plt

titulo = 'EDO Runge-Kutta 2ord - Paracaidista Wingsuit'
i = muestras # iteración en gráfica

titulo = titulo+', i='+str(i)
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
K1y = tabla[:,3]
K1z = tabla[:,4]
K2y = tabla[:,5]
K2z = tabla[:,6]

plt.subplot(211)
plt.plot(xi[0],yi[0],'o',
         color='red', label ='[t0,y0]')
plt.plot(xi[1:i+2],yi[1:i+2],'o',
         color='green', label ='[t[i],y[i]]')
plt.plot(xi[0:i+2],yi[0:i+2],
         color='blue',label='y(t)')
if i<muestras: # gráfica para una iteración
    plt.plot(xi[i+1],yi[i+1],'o',color='orange',
             label ='[x[i+1],y[i+1]]')
    plt.plot(xi[i:i+3],yi[i:i+3],'.',color='gray')
    plt.plot(xi[i:i+2],[yi[i],yi[i]], color='orange',
             label='h',linestyle='dashed')
    plt.plot([xi[i+1]-0.02*h,xi[i+1]-0.02*h],
             [yi[i],yi[i]+K1y[i+1]],
             color='green',label='K1y',linestyle='dashed')
    plt.plot([xi[i+1]+0.02*h,xi[i+1]+0.02*h],
             [yi[i],yi[i]+K2y[i+1]],
             color='magenta',label='K2y',linestyle='dashed')
    plt.plot([xi[i+1]-0.02*h,xi[i+1]+0.02*h],
             [yi[i]+K1y[i+1],yi[i]+K2y[i+1]],
             color='magenta')
if np.min(yi[0:i+1])<0: # linea 0
    plt.axhline(0, color='red')
plt.ylabel('y = Altura')
plt.title(titulo)
plt.legend()
plt.grid()
plt.tight_layout()

plt.subplot(212)
plt.plot(xi[0],zi[0],'o',
         color='red', label ='[t0,v0]')
plt.plot(xi[1:i+2],zi[1:i+2],'o',
         color='green', label ='[t[i],v[i]]')
plt.plot(xi[0:i+2],zi[0:i+2],
         color='green',label='v(t)')
if i<muestras: # gráfica para una iteración
    plt.plot(xi[i+1],zi[i+1],'o',color='orange',
             label ='[t[i+1],v[i+1]]')
    plt.plot(xi[i:i+3],zi[i:i+3],'.',color='gray')
    plt.plot(xi[i:i+2],[zi[i],zi[i]], color='orange',
             label='h',linestyle='dashed')
    plt.plot([xi[i+1]-0.02*h,xi[i+1]-0.02*h],
             [zi[i],zi[i]+K1z[i+1]],
             color='green',label='K1z',linestyle='dashed')
    plt.plot([xi[i+1]+0.02*h,xi[i+1]+0.02*h],
             [zi[i],zi[i]+K2z[i+1]],
             color='magenta',label='K2z',linestyle='dashed')
    plt.plot([xi[i+1]-0.02*h,xi[i+1]+0.02*h],
             [zi[i]+K1z[i+1],zi[i]+K2z[i+1]],
             color='magenta')

plt.xlabel('x = tiempo')
plt.ylabel('z = velocidad')
plt.legend()
plt.grid()
plt.tight_layout()

plt.show() #comentar para la siguiente gráfica

Runge Kutta $\frac{\delta^2 y}{\delta x^2}$

Algoritmo:

6. 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}$

Algoritmo:

6.3 EDO d2y/dx2 Runge-Kutta con Python

1. EDO $\frac{\delta^2 y}{\delta x^2}$ con Runge-Kutta

2. Ejercicio

3. Desarrollo analítico - paso a paso

4. Algoritmo en Python para $\frac{\delta^2 y}{\delta x^2}$ Runge-Kutta 2do Orden

5. Gráfica

6. Algoritmo en Python para $\frac{\delta^2 y}{\delta x^2}$ Runge-Kutta 4to Orden

Métodos Numéricos

Análisis Numérico

6.3 EDO d2y/dx2 Runge-Kutta con Python

1. EDO \frac{\delta^2 y}{\delta x^2} con Runge-Kutta

2. Ejercicio

3. Desarrollo analítico - paso a paso

4. Algoritmo en Python para \frac{\delta^2 y}{\delta x^2} Runge-Kutta 2do Orden

5. Gráfica

6. Algoritmo en Python para \frac{\delta^2 y}{\delta x^2} Runge-Kutta 4to Orden

Métodos Numéricos

Análisis Numérico

1. EDO $\frac{\delta^2 y}{\delta x^2}$ con Runge-Kutta

4. Algoritmo en Python para $\frac{\delta^2 y}{\delta x^2}$ Runge-Kutta 2do Orden

6. Algoritmo en Python para $\frac{\delta^2 y}{\delta x^2}$ Runge-Kutta 4to Orden