animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para la solución de Ecuaciones Diferenciales Ordinarias, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

..

---

EDO con Taylor de 3 términos

Tabla de resultados:

 EDO con Taylor 3 términos
 [xi,     yi,     d1yi,    d2yi,   término 1,   término 2 ]
[[0.         1.         0.         0.         0.         0.        ]
 [0.1        1.215      2.         3.         0.2        0.015     ]
 [0.2        1.461025   2.305      3.105      0.2305     0.015525  ]
 [0.3        1.73923262 2.621025   3.221025   0.2621025  0.01610513]
 [0.4        2.05090205 2.94923262 3.34923262 0.29492326 0.01674616]
 [0.5        2.39744677 3.29090205 3.49090205 0.32909021 0.01745451]]
>>> 

Instrucciones en Python

# EDO. Método de Taylor con3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras, vertabla=False, precision=6):
    ''' solucion a EDO usando tres términos de Taylor,
    x0,y0 son valores iniciales, h es el tamaño de paso,
    muestras es la cantidad de puntos a calcular.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,6),dtype=float)
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        d1yi = d1y(x,y)
        d2yi = d2y(x,y)
        y = y + h*d1yi + ((h**2)/2)*d2yi
        x = x + h
        
        term1 = h*d1yi
        term2 = ((h**2)/2)*d2yi
        
        tabla[i] = [x,y,d1yi,d2yi,term1,term2]
    if vertabla==True:
        titulo = ' [xi,     yi,     d1yi,   d2yi,'
        titulo = titulo + '   término 1,   término 2 ]'
        np.set_printoptions(precision)
        print(' EDO con Taylor 3 términos')
        print(titulo)
        print(tabla)
        
    return(tabla)

# PROGRAMA -----------------
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1
# INGRESO.
# d1y = y', d2y = y''
d1y = lambda x,y: y - x**2 + x + 1
d2y = lambda x,y: y - x**2 - x + 2
x0 = 0
y0 = 1
h  = 0.1
muestras = 5

# PROCEDIMIENTO
tabla = edo_taylor3t(d1y,d2y,x0,y0,h,muestras,
                     vertabla=True)

# SALIDA
##print(' EDO con Taylor 3 términos')
##print(' [xi,     yi,     d1yi,',
##      '   d2yi,   término 1,   término 2 ]')
##print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
xi = tabla[:,0]
yi = tabla[:,1]
plt.plot(xi,yi)
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
#plt.show() #comentar para la siguiente gráfica


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation

unmetodo = 'Edo con Taylor 3 términos'
narchivo = 'EdoTaylor3t' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
term1 = tabla[:,4]
term2 = tabla[:,5]
dfi = tabla[:,2]
n = len(xi)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])
# Lineas y puntos base
linea_fx, = graf_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf_ani.plot(xi[0], yi[0],'o',
                        color='green',
                        label='xi,yi')
puntoa, = graf_ani.plot(xi[0], yi[0],'o',
                        color='Blue')
puntob, = graf_ani.plot(xi[1], yi[1],'o',
                        color='orange')

linea_h, = graf_ani.plot(xi, xi, color='orange',
                         label='h',
                         linestyle='dashed')
linea_term1, = graf_ani.plot(xi, xi,
                             color='green',label="h*y'[i]",
                             linestyle='dashed')
linea_term2, = graf_ani.plot(xi, yi, linewidth=4,
                             color='magenta',
                             label="((h**2)/2!)*y''[i]")
# Aproximacion con tangente
b0 = yi[0] - dfi[1]*xi[0]
tangentei = dfi[1]*xi + b0
linea_tang, = graf_ani.plot(xi, tangentei, color='dodgerblue',
                             label="tangente",
                             linestyle='dotted')

# Cuadros de texto en gráfico
txt_i  = graf_ani.text(xi[0], yi[0]+0.03*deltay,'[x[i],y[i]]',
                       horizontalalignment='center')
txt_i1 = graf_ani.text(xi[1], xi[1]+0.03*deltay,'[x[i+1],y[i+1]]',
                       horizontalalignment='center')
# Configura gráfica
graf_ani.axhline(0, color='black')  # Linea horizontal en cero
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,dfi,term1,term2):
    
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    else:
        ymax = np.max([yi[0],yi[2]])
        ymin = np.min([yi[0],yi[2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    # actualiza cada punto
    puntoa.set_xdata(xi[i]) 
    puntoa.set_ydata(yi[i])
    puntob.set_xdata(xi[i+1]) 
    puntob.set_ydata(yi[i+1])
    puntof.set_xdata(xi[0:i]) 
    puntof.set_ydata(yi[0:i])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_term1.set_xdata([xi[i+1],xi[i+1]])
    linea_term1.set_ydata([yi[i],yi[i]+term1[i+1]])
    linea_term2.set_xdata([xi[i+1],xi[i+1]])
    linea_term2.set_ydata([yi[i]+term1[i+1],
                           yi[i]+term1[i+1]+term2[i+1]])
    
    b0 = yi[i] - dfi[i+1]*xi[i]
    tangentei = dfi[i+1]*xi + b0
    linea_tang.set_ydata(tangentei)
    # actualiza texto
    txt_i.set_position([xi[i], yi[i]+0.03*deltay])
    txt_i1.set_position([xi[i+1], yi[i+1]+0.03*deltay])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_tang,
            linea_term1,linea_term2,
            txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_term1.set_ydata(np.ma.array(xi, mask=True))
    linea_term2.set_ydata(np.ma.array(xi, mask=True))
    linea_tang.set_ydata(np.ma.array(xi, mask=True))
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_tang,
            linea_term1,linea_term2,
            txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,dfi,term1,term2),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

..

---

Runge Kutta de 2do Orden para primera derivada

 EDO con Runge-Kutta 2do Orden primera derivada
 [xi,     yi,     K1,    K2 ]
[[0.       1.       0.       0.      ]
 [0.1      1.2145   0.2      0.229   ]
 [0.2      1.459973 0.23045  0.260495]
 [0.3      1.73757  0.261997 0.293197]
 [0.4      2.048564 0.294757 0.327233]
 [0.5      2.394364 0.328856 0.362742]]

Instrucciones en Python

# EDO. Método de Runge-Kutta 2do Orden primera derivada 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras, vertabla=False, precision=6):
    ''' solucion a EDO con Runge-Kutta 2do Orden primera derivada,
        x0,y0 son valores iniciales
        muestras es la cantidad de puntos a calcular con tamaño de paso h.
    '''
    # 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]
    if vertabla==True:
        np.set_printoptions(precision)
        titulo = ' [xi,     yi,     K1,    K2 ]'
        print(' EDO con Runge-Kutta 2do Orden primera derivada')
        print(titulo)
        print(tabla)
    return(tabla)

# PROGRAMA -----------------
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1
# INGRESO.
# d1y = y', d2y = y''
d1y = lambda x,y: y - x**2 + x + 1
d2y = lambda x,y: y - x**2 - x + 2
x0 = 0
y0 = 1
h  = 0.1
muestras = 5

# PROCEDIMIENTO
tabla = rungekutta2(d1y,x0,y0,h,muestras,
                     vertabla=True)

# SALIDA
##print(' EDO con Runge-Kutta 2do Orden primera derivada')
##print(' [xi,     yi,     d1yi,',', K1,   K2 ]')
##print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
xi = tabla[:,0]
yi = tabla[:,1]
plt.plot(xi,yi)
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución Runge-Kutta 2do Orden primera derivada')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
#plt.show() #comentar para la siguiente gráfica


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation

unmetodo = 'EDO: Runge-Kutta 2do Orden primera derivada'
narchivo = 'EdoRK2df' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
K1 = tabla[:,2]
K2 = tabla[:,3]
n = len(xi)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
# Lineas y puntos base
linea_fx, = graf_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf_ani.plot(xi[0], yi[0],'o',
                        color='green',
                        label='xi,yi')
puntoa, = graf_ani.plot(xi[0], yi[0],'o',
                        color='Blue')
puntob, = graf_ani.plot(xi[1], yi[1],'o',
                        color='orange')

linea_h, = graf_ani.plot(xi, xi, color='orange',
                         label='h',
                         linestyle='dashed')
linea_K1, = graf_ani.plot(xi-0.02*deltax, xi-0.02*deltax,
                          color='green',label="K1",
                          linestyle='dashed')
linea_K2, = graf_ani.plot(xi+0.02*deltax, yi,
                          color='magenta',
                          label="K2",
                          linestyle='dashed')
linea_K12, = graf_ani.plot(xi, yi,
                          color='magenta')

# Cuadros de texto en gráfico
txt_i  = graf_ani.text(xi[0], yi[0]+0.05*deltay,'[x[i],y[i]]',
                       horizontalalignment='center')
txt_i1 = graf_ani.text(xi[1], xi[1]+0.05*deltay,'[x[i+1],y[i+1]]',
                       horizontalalignment='center')

# Configura gráfica
graf_ani.axhline(0, color='black')  # Linea horizontal en cero
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('f(x)')
graf_ani.legend()
graf_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,term1,term2):   
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    else:
        ymax = np.max([yi[0],yi[2]])
        ymin = np.min([yi[0],yi[2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])
    # actualiza cada punto
    puntoa.set_xdata(xi[i]) 
    puntoa.set_ydata(yi[i])
    puntob.set_xdata(xi[i+1]) 
    puntob.set_ydata(yi[i+1])
    puntof.set_xdata(xi[0:i]) 
    puntof.set_ydata(yi[0:i])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_K1.set_xdata([xi[i+1]-0.02*deltax,xi[i+1]-0.02*deltax])
    linea_K1.set_ydata([yi[i],yi[i]+K1[i+1]])
    linea_K2.set_xdata([xi[i+1]+0.02*deltax,xi[i+1]+0.02*deltax])
    linea_K2.set_ydata([yi[i],yi[i]+K2[i+1]])
    linea_K12.set_xdata([xi[i+1]-0.02*deltax,xi[i+1]+0.02*deltax])
    linea_K12.set_ydata([yi[i]+K1[i+1],yi[i]+K2[i+1]])
    # actualiza texto
    txt_i.set_position([xi[i], yi[i]+0.05*deltay])
    txt_i1.set_position([xi[i+1], yi[i+1]+0.05*deltay])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_K1,linea_K2,linea_K12,
            txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_K1.set_ydata(np.ma.array(xi, mask=True))
    linea_K2.set_ydata(np.ma.array(xi, mask=True))
    linea_K12.set_ydata(np.ma.array(xi, mask=True))
    
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,linea_K1,linea_K2,linea_K12,
            txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,K1,K2),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]
..

---

Sistemas EDO. modelo depredador-presa con Runge-Kutta 2do Orden
.

Instrucciones en Python

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 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]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    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:
            txt = ' '
            if i>=10:
                txt='  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35

# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y

# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1

# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras,vertabla=True)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
print('Sistemas EDO: Modelo presa-predador')
##np.set_printoptions(precision=6)
##print(' [ ti, xi, yi]')
##print(tabla[:,0:4])

# GRAFICA tiempos vs población
import matplotlib.pyplot as plt

fig_t, (graf1,graf2) = plt.subplots(2)
fig_t.suptitle('Modelo predador-presa')
graf1.plot(ti,xi, color='blue',label='xi presa')

#graf1.set_xlabel('t tiempo')
graf1.set_ylabel('población x')
graf1.legend()
graf1.grid()

graf2.plot(ti,yi, color='orange',label='yi predador')
graf2.set_xlabel('t tiempo')
graf2.set_ylabel('población y')
graf2.legend()
graf2.grid()

# gráfica xi vs yi
fig_xy, graf3 = plt.subplots()
graf3.plot(xi,yi)

graf3.set_title('Modelo presa-predador [xi,yi]')
graf3.set_xlabel('x presa')
graf3.set_ylabel('y predador')
graf3.grid()
#plt.show()


# GRAFICA CON ANIMACION --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

unmetodo = 'Sistemas EDO Presa-Predador con Runge-Kutta'
narchivo = 'EdoPresaPredador' # nombre archivo GIF
muestras = 51 

# Puntos para la gráfica
a = xi[0]
b = xi[1]
n = len(ti)

# Parametros de trama/foto
retardo = 1000   # milisegundos entre tramas
tramas = len(xi)

# GRAFICA animada en fig_ani
fig_ani, (graf1_ani,graf2_ani) = plt.subplots(2)
ymax = np.max([yi[0],yi[2]])
ymin = np.min([yi[0],yi[2]])
deltax = np.abs(xi[2]-xi[0])
deltay = np.abs(yi[2]-yi[0])
graf1_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.05*deltay])

zmax = np.max([zi[0],zi[2]])
zmin = np.min([zi[0],zi[2]])
deltax = np.abs(xi[2]-xi[0])
deltaz = np.abs(zi[2]-zi[0])
graf2_ani.set_xlim([xi[0],xi[2]+0.05*deltax])
graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.05*deltaz])
# Lineas y puntos base
linea_fx, = graf1_ani.plot(xi, yi,color='blue',
                          linestyle='dashed')
puntof, = graf1_ani.plot(xi[0], yi[0],'o',
                        color='blue',
                        label='xi,yi')
puntoa, = graf1_ani.plot(xi[0], yi[0],'o',
                        color='green')
puntob, = graf1_ani.plot(xi[1], yi[1],'o',
                        color='dodgerblue')
linea_h, = graf1_ani.plot(xi, xi, color='green',
                         label='h',
                         linestyle='dashed')
linea_term1, = graf1_ani.plot(xi, xi,
                             color='dodgerblue',label="(K1y+K2y)/2",
                             linestyle='dashed')
# Cuadros de texto en gráfico
#txt_i  = graf1_ani.text(xi[0], yi[0]+0.03*deltay,'[x[i],y[i]]',
#                       horizontalalignment='center')
#txt_i1 = graf1_ani.text(xi[1], xi[1]+0.03*deltay,'[x[i+1],y[i+1]]',
#                       horizontalalignment='center')

linea_gx, = graf2_ani.plot(xi, zi,color='orange',
                          linestyle='dashed')
puntog, = graf2_ani.plot(xi[0], zi[0],'o',
                        color='orange',
                        label='xi,zi')
puntog_a, = graf2_ani.plot(xi[0], zi[0],'o',
                        color='green')
puntog_b, = graf2_ani.plot(xi[1], zi[1],'o',
                        color='red')
lineag_h, = graf2_ani.plot(xi, xi, color='green',
                         label='h',
                         linestyle='dashed')
lineag_term1, = graf2_ani.plot(xi, xi,
                             color='red',label="(K1z+K2z)/2",
                             linestyle='dashed')

# Configura gráfica
graf1_ani.axhline(0, color='black')  # Linea horizontal en cero
graf1_ani.set_title(unmetodo)
graf1_ani.set_xlabel('x')
graf1_ani.set_ylabel('y(x)')
graf1_ani.legend()
graf1_ani.grid()

graf2_ani.axhline(0, color='black')  # Linea horizontal en cero
#graf2_ani.set_title(unmetodo)
graf2_ani.set_xlabel('x')
graf2_ani.set_ylabel('z(x)')
graf2_ani.legend()
graf2_ani.grid()

# Nueva Trama
def unatrama(i,xi,yi,zi):
    
    if i>1:
        ymax = np.max([yi[0:i+2]])
        ymin = np.min([yi[0:i+2]])
        deltax = np.abs(xi[i+1]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf1_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])

        zmax = np.max([zi[0:i+2]])
        zmin = np.min([zi[0:i+2]])
        deltaz = np.abs(zmax-zmin)
        graf2_ani.set_xlim([xi[0]-0.05*deltax,xi[i+1]+0.05*deltax])
        graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.1*deltaz])
    else:
        ymax = np.max([yi[0:2]])
        ymin = np.min([yi[0:2]])
        deltax = np.abs(xi[2]-xi[0])
        deltay = np.abs(ymax-ymin)
        graf1_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf1_ani.set_ylim([ymin-0.05*deltay,ymax+0.1*deltay])

        zmax = np.max([zi[0:2]])
        zmin = np.min([zi[0:2]])
        deltaz = np.abs(zmax-zmin)
        graf2_ani.set_xlim([xi[0]-0.05*deltax,xi[2]+0.05*deltax])
        graf2_ani.set_ylim([zmin-0.05*deltaz,zmax+0.1*deltaz])
    # actualiza cada punto
    puntoa.set_xdata(xi[i]) 
    puntoa.set_ydata(yi[i])
    puntob.set_xdata(xi[i+1]) 
    puntob.set_ydata(yi[i+1])
    puntof.set_xdata(xi[0:i]) 
    puntof.set_ydata(yi[0:i])
    # actualiza cada linea
    linea_fx.set_xdata(xi[0:i+1])
    linea_fx.set_ydata(yi[0:i+1])
    linea_h.set_xdata([xi[i],xi[i+1]])
    linea_h.set_ydata([yi[i],yi[i]])
    linea_term1.set_xdata([xi[i+1],xi[i+1]])
    linea_term1.set_ydata([yi[i],yi[i+1]])
    # actualiza texto
    #txt_i.set_position([xi[i], yi[i]+0.03*deltay])
    #txt_i1.set_position([xi[i+1], yi[i+1]+0.03*deltay])

    # actualiza cada punto
    puntog_a.set_xdata(xi[i]) 
    puntog_a.set_ydata(zi[i])
    puntog_b.set_xdata(xi[i+1]) 
    puntog_b.set_ydata(zi[i+1])
    puntog.set_xdata(xi[0:i]) 
    puntog.set_ydata(zi[0:i])
    # actualiza cada linea
    linea_gx.set_xdata(xi[0:i+1])
    linea_gx.set_ydata(zi[0:i+1])
    lineag_h.set_xdata([xi[i],xi[i+1]])
    lineag_h.set_ydata([zi[i],zi[i]])
    lineag_term1.set_xdata([xi[i+1],xi[i+1]])
    lineag_term1.set_ydata([zi[i],zi[i+1]])

    return (puntoa,puntob,puntof,
            linea_fx,linea_h,
            linea_term1,)
            #txt_i,txt_i1,)
# Limpia Trama anterior
def limpiatrama(): 
    puntoa.set_ydata(np.ma.array(xi, mask=True))
    puntob.set_ydata(np.ma.array(xi, mask=True))
    puntof.set_ydata(np.ma.array(xi, mask=True))
    linea_h.set_ydata(np.ma.array(xi, mask=True))
    linea_term1.set_ydata(np.ma.array(xi, mask=True))

    puntog_a.set_ydata(np.ma.array(xi, mask=True))
    puntog_b.set_ydata(np.ma.array(xi, mask=True))
    puntog.set_ydata(np.ma.array(xi, mask=True))
    lineag_h.set_ydata(np.ma.array(xi, mask=True))
    lineag_term1.set_ydata(np.ma.array(xi, mask=True))
    return (puntoa,puntob,puntof,
            linea_fx,linea_h,
            linea_term1,
            puntog_a,puntog_b,puntog,
            linea_gx,lineag_h,
            lineag_term1,)
            #txt_i,txt_i1,)

# Trama contador
i = np.arange(0,tramas-1,1)
ani = animation.FuncAnimation(fig_ani,
                              unatrama,
                              i ,
                              fargs = (xi,yi,zi),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

animación: [ EDO Taylor 3t ] [ Runge Kutta  dy/dx ] [ Sistemas EDO con RK ]

Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  d²y/dx² ]
..

---

1. Ejercicio

Referencia: Chapra 28.2 p831 pdf855, Rodriguez 9.2.1 p263
https://es.wikipedia.org/wiki/Ecuaciones_Lotka%E2%80%93Volterra
https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

Modelos depredador-presa y caos. Ecuaciones Lotka-Volterra. En el sistema de ecuaciones:

\frac{dx}{dt} = ax - bxy \frac{dy}{dt} = -cy + dxy

variables
x = número de presas
y = número de depredadores
t = tiempo de observación
coeficientes
a = razón de crecimiento de la presa, (0.5)
c = razón de muerte del depredador (0.35)
b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)

Considere como puntos iniciales en la observación de las especies:
t = 0, x = 2, y = 1, h = 0.5

Los términos que multiplican xy hacen que las ecuaciones sean no lineales.
Observe que la variable tiempo no se encuentra en las expresiones f y g, h se aplica a tiempo.

Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  d²y/dx² ]
..

---

2. Desarrollo analítico

Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:

f(t,x,y) = 0.5 x - 0.7 xy g(t,x,y) = -0.35y + 0.35xy

Las expresiones se adaptan al método de Runge-Kutta para primeras derivadas por cada variable de población. Se deben usar de forma simultánea para cada tiempo t.

K1x = h f(t,x,y) = 0.5 \Big( 0.5 x - 0.7 xy \Big) K1y = h g(t,x,y) = 0.5 \Big(-0.35y + 0.35xy \Big)

..

K2x = h f(t+h,x+K1x,y+K1y) = 0.5 \Big( 0.5 (x+K1x) - 0.7 (x+K1x)(y+K1y) \Big) K2y = h g(t+h,x+K1x,y+K1y) = 0.5 \Big(-0.35(y+K1y) + 0.35(x+K1x)(y+K1y) \Big)

..

x[i+1] = x[i] + \frac{K1x+K2x}{2} y[i+1] = y[i] + \frac{K1y+K2y}{2} t[i+1] = t[i] + h

con lo que se puede aplicar al ejercicio en cada iteración. dadas las condiciones iniciales.

Itera = 0

t = 0, x = 2, y = 1, h = 0.5

K1x = 0.5 f(0,2,1) = 0.5 \Big( 0.5 (2) - 0.7 (2)(1) \Big) = -0.2 K1y = 0.5 g(0,2,1) = 0.5 \Big(-0.35(1) + 0.35(2)(1) \Big) =0.175

..

K2x = 0.5 f(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big( 0.5 (2+(-0.2)) - 0.7 (2+(-0.2))(1+0.175) \Big) = -0.29025 K2y = 0.5 g(0+0.5, 2+(-0.2), 1+0.175) = 0.5 \Big(-0.35(1+0.175) + 0.35(2+(-0.2))(1+0.175) \Big) = 0.1645

..

x[1] = x[0] + \frac{K1x+K2x}{2} = 2 + \frac{-0.2+(-0.29025)}{2} = 1.7548 y[1] = y[0] + \frac{K1y+K2y}{2} = 1 + \frac{0.175+0.1645}{2}= 1.1697 t[1] = t[0] + h = 0 +0.5 = 0.5

itera = 1

t = 0.5, x = 1.7548, y = 1.1697, h = 0.5

K1x = 0.5 \Big( 0.5 (0,1.7548) - 0.7 (0,1.7548)(1.1697) \Big) = -0.2797 K1y = 0.5 \Big(-0.35(1.1697) + 0.35(1.7548)(1.1697) \Big) =0.1545

..

K2x = 0.5 \Big( 0.5 (1.7548+(-0.2797)) - 0.7 (1.7548+(-0.2797))(1.1697+0.1545) \Big) =-0.3149 K2y = 0.5 \Big(-0.35(1.1697+0.1545) + 0.35(1.7548+(-0.2797))(1.1697+0.1545) \Big) = 0.1645

..

x[2] = 1.7548 + \frac{-0.2797+(-0.3149)}{2} = 1.4575 y[2] = 1.1697 + \frac{0.1545+0.1645}{2} = 1.3020 t[2] = t[0] + h = 0.5 +0.5 = 1

itera=2

t = 1, x = 1.4575, y = 1.3020, h = 0.5

continuar como tarea …

Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  d²y/dx² ]

..

---

3. Algoritmo en Python

Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en

EDO con Runge-Kutta d2y/dx2

Al ejecutar el algoritmo se obtienen los siguientes resultados:

Runge-Kutta Segunda derivada
i  [ xi,  yi,  zi ]
   [ K1y,  K1z,  K2y,  K2z ]
0 [0. 2. 1.]
   [0. 0. 0. 0.]
1 [0.5      1.754875 1.16975 ]
  [-0.2      0.175   -0.29025  0.1645 ]
2 [1.       1.457533 1.302069]
  [-0.279749  0.154528 -0.314935  0.11011 ]
3 [1.5      1.167405 1.373599]
  [-0.29985   0.104254 -0.280406  0.038807]
4 [2.       0.922773 1.381103]
  [-0.26939   0.040241 -0.219874 -0.025233]
5 [2.5      0.734853 1.33689 ]
  [-0.215362 -0.018665 -0.160478 -0.069761]
6 [3.       0.598406 1.258434]
  [-0.160133 -0.062033 -0.11276  -0.09488 ]
... 

Los resultados de la tabla se muestran parcialmente, pues se usaron mas de 100 iteraciones.

Los resultados se pueden observar de diferentes formas:

a) Cada variable x_i, y_i versus t_i, es decir cantidad de animales de cada especie durante el tiempo de observación

b) Independiente de la unidad de tiempo, x_i vs y_i, muestra la relación entre la cantidad de presas y predadores. Relación que es cíclica y da la forma a la gráfica.

Las instrucciones del algoritmo en Python usadas en el problema son:

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 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]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    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:
            txt = ' '
            if i>=10:
                txt='  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35

# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y

# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1

# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras,vertabla=True)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
print('Sistemas EDO: Modelo presa-predador')
##np.set_printoptions(precision=6)
##print(' [ ti, xi, yi]')
##print(tabla[:,0:4])

Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:

# GRAFICA tiempos vs población
import matplotlib.pyplot as plt

fig_t, (graf1,graf2) = plt.subplots(2)
fig_t.suptitle('Modelo predador-presa')
graf1.plot(ti,xi, color='blue',label='xi presa')

#graf1.set_xlabel('t tiempo')
graf1.set_ylabel('población x')
graf1.legend()
graf1.grid()

graf2.plot(ti,yi, color='orange',label='yi predador')
graf2.set_xlabel('t tiempo')
graf2.set_ylabel('población y')
graf2.legend()
graf2.grid()

# gráfica xi vs yi
fig_xy, graf3 = plt.subplots()
graf3.plot(xi,yi)

graf3.set_title('Modelo presa-predador [xi,yi]')
graf3.set_xlabel('x presa')
graf3.set_ylabel('y predador')
graf3.grid()
plt.show()

Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  d²y/dx² ]

[ Runge Kutta  d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..

---

1. EDO Runge-Kutta para Segunda derivada d²y/dx²

Para una ecuación diferencial de segunda derivada (segundo orden) con condiciones de inicio en x₀, y₀, y’₀

\frac{\delta ^2 y}{\delta x^2} = \frac{\delta y}{\delta x} + etc

Forma estandarizada de la ecuación:

y'' = y' + etc

Se 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(h³)
Runge-Kutta 4do Orden tiene error de truncamiento O(h⁵)

[ Runge Kutta  d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]

..

---

2. Runge-Kutta 2do Orden para Segunda derivada d²y/dx² en Python

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

import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 son valores iniciales, h es tamano de paso,
        muestras es la cantidad de puntos a calcular.
        f(x,y,z) = z #= y'
        g(x,y,z) = expresion con z=y'
    '''
    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
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    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:
            txt = ' '
            if i>=10:
                txt = '  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)

[ Runge Kutta  d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..

---

3. Runge-Kutta 4do Orden para Segunda derivada d²y/dx² en Python

import numpy as np

def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras):
    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
    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]
    return(tabla)

[ Runge Kutta  d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..

---

4. Ejercicio

2Eva_IT2018_T1 Paracaidista wingsuit

Solución Propuesta: s2Eva_IT2018_T1 Paracaidista wingsuit

otro ejercicio, una aplicación del algoritmo en Señales y Sistemas:

LTI CT – Respuesta entrada cero – Desarrollo analítico, TELG1001-Señales y Sistemas

[ Runge Kutta  d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]

[ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

..

---

EDO Runge-Kutta 4to Orden de Primera derivada dy/dx

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:

\frac{\delta y}{\delta x} + etc =0 y'(x) = f(x_i,y_i) y(x_0) = y_0

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

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

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

y_{i+1} = y_i + h f(x_i,y_i) + + \frac{h^2}{2!} f'(x_i,y_i) + \frac{h^3}{3!} f''(x_i,y_i) + +\frac{h^4}{4!} f'''(x_i,y_i) + O(h^5) x_{i+1} = x_i + h

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

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:

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

iteración:

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

Las iteraciones se dejan como tarea

[ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

..

---

Algoritmo en Python como Función

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)
        titulo = ' [xi,     yi,     K1,    K2,     K3,     K4 ]'
        print(' EDO con Runge-Kutta 4do Orden primera derivada')
        print(titulo)
        print(tabla)
    return(tabla)

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.

[ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

..

---

Ejercicio

2Eva_IT2018_T1 Paracaidista wingsuit

Solución Propuesta: s2Eva_IT2018_T1 Paracaidista wingsuit

 

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

[ Runge Kutta 4to Orden ] [ Función ] [ Ejercicio en video ]

[ 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

Para una ecuación diferencial ordinaria de primera derivada, el método Runge-Kutta de 2do Orden usa una corrección sobre la derivada a partir de los puntos xi y xi+h,  es decir un tamaño de paso h hacia adelante, calculados como términos K₁ y K₂.

Considere una ecuación diferencial de primera derivada con una condición de inicio se reordena y se escribe como f(x,y) siguiendo los pasos:

\frac{\delta y}{\delta x} + etc =0 y'(x) = f(x_i,y_i) y(x_0) = y_0

Los términos K₁ y K₂ se calculan para predecir el próximo valor en y[i+1], observe que el término K₁ es el mismo que el método de Edo con Taylor de dos términos.

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.

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

Para observar al idea básica, considere observar un combate aéreo simulado en la década de 1940, donde las armas se encuentras fijas en las alas del avión. Observe dos minutos del video sugerido a partir de donde se encuentra marcado el enlace.

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) = 0.1 (y_i -x_i^2 +x_i +1) K_2 = h f(x_i+h, y_i + K_1) K_2 = 0.1 \Big((y_i+K_1) -(x_i+h)^2 +(x_i+h) +1 \Big) 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, se realizan tres iteraciones en las actividades del curso con lápiz y papel,

itera = 0 , x₀ = 0, y₀ = 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 , x₁ = 0.1, y₁ = 1.2145

K_1 = 0.1 f(0.1,1.2145) = 0.1( 1.2145 -0.1^2 +0.1 +1) K_1 = 0.2304 K_2 = 0.1 f(0.1+0.1, 1.2145 + 0.2304) =0.1 \Big((1.2145 + 0.2304) -(0.1+0.1)^2 +(0.1+0.1) +1\Big) K_2 = 0.2604 y_2 = 1.2145 + \frac{0.2304+0.2604}{2} = 1.4599 x_2 = 0.1 +0.1 = 0.2

itera = 2 , x₂ = 0.2, y₂ = 1.4599

K_1 = 0.1 f(0.2,1.4599) = 0.1( 1.4599 -0.2^2 +0.2 +1) K_1 = 0.2619 K_2 = 0.1 f(0.2+0.1, 1.4599 + 0.2619) =0.1 \Big((1.4599 + 0.2619) -(0.2+0.1)^2 +(0.2+0.1) +1\Big) K_2 = 0.2931 y_2 = 1.4599 + \frac{0.2619+0.2931}{2} = 1.7375 x_2 = 0.2 +0.1 = 0.3

luego de las 3 iteraciones en papel, se completan los demás puntos con el algoritmo obteniendo la gráfica resultante para y(x) correspondiente.

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.

[ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

---

1. Ecuaciones diferenciales ordinarias aproximadas con Taylor

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

En los métodos con Taylor para Ecuaciones Diferenciales Ordinarias (EDO) se aproxima el resultado a n términos de la serie, para lo cual se ajusta la expresión del problema a cada derivada correspondiente.

La solución empieza usando la Serie de Taylor para tres términos ajustada a la variable del ejercicio:

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)

A partir de la expresión de y'(x) y el punto inicial conocido en x[i],se busca obtener el próximo valor en x[i+1] al avanzar un tamaño de paso h. Se repite el proceso en el siguiente punto encontrado y se continua hasta alcanzar el intervalo objetivo.

En éstos métodos la solución siempre es una tabla de puntos xi,yi que se pueden usar para interpolar y obtener una función polinómica.

[ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]

..

---

2. Ejercicio

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

Se requiere 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

que con nomenclatura simplificada:

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

[ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

---

3. Desarrollo Analítico

Al despejar el valor de  y’ de expresión del ejercicio,

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

se puede obtener y" al derivar una vez,

y'' = y' -2x + 1

para luego combinar las expresiones en

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

simplificando:

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

Ecuaciones que permiten estimar nuevos valores y_i+1 para nuevos puntos  muestra distanciados en i*h desde el punto inicial siguiendo las siguientes expresiones de iteración:

y'_i = y_i -x_i^2 + x_i +1 y''_i = y_i -x_i^2 - x_i +2 y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h

Se empieza evaluando el nuevo punto a una distancia x₁= x₀+h del punto de origen con lo que se obtiene y₁ , repitiendo el proceso para el siguiente punto en forma sucesiva.

itera = 0 , x0 = 0, y0 = 1

y'_0 = 1 -0^2 +0 +1 = 2 y''_0 = 1 -0^2 -0 +2 = 3 y_1 = y_{0} + h y'_0 + \frac{h^2}{2!} y''_0 y_1 = 1 + 0.1 (2) + \frac{0.1^2}{2!} 3 = 1.215 x_1 = 0 + 0.1

itera = 1 , x = 0.1, y = 1.215

y'_1 = 1.215 - 0.1^2 + 0.1 +1 = 2.305 y''_1 = 1.215 - 0.1^2 - 0.1 +2 = 3.105 y_2 = 1.215 + 0.1 (2.305) + \frac{0.1^2}{2!} 3.105 = 1.461 x_2 = 0.1 + 0.1 = 0.2

itera = 2 , x = 0.2, y = 1.461

y'_2 = 1.461 - 0.2^2 + 0.2 +1 = 2.621 y''_2 = 1.461 - 0.2^2 - 0.2 +2 = 3.221 y_3 = 1.461 + 0.1 (2.621) + \frac{0.1^2}{2!} 3.221 = 1.7392 x_3 = 0.2 + 0.1 = 0.3

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

 EDO con Taylor 3 términos
[xi, yi, d1yi, d2yi]
[[0.         1.         0.         0.        ]
 [0.1        1.215      2.         3.        ]
 [0.2        1.461025   2.305      3.105     ]
 [0.3        1.73923262 2.621025   3.221025  ]
 [0.4        2.05090205 2.94923262 3.34923262]
 [0.5        2.39744677 3.29090205 3.49090205]]
>>>

Observación, note que los resultados de las derivadas, se encuentran desplazados una fila para cada iteración. Asunto a ser considerado en la gráfica de las derivadas en caso de incluirlas.

Nota: Compare luego los pasos del algoritmo con el método de Runge-Kutta de 2do orden.

[ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..

---

4. Algoritmo en Python

Para simplificar los cálculos se crea una función edo_taylor3t() para encontrar  los valores para una cantidad de muestras distanciadas entre si h veces del punto inicial [x₀,y₀]

# EDO. Método de Taylor con3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    ''' solucion a EDO usando tres términos de Taylor, x0,y0 son valores iniciales
        muestras es la cantidad de puntos a calcular con tamaño de paso h.
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        d1yi = d1y(x,y)
        d2yi = d2y(x,y)
        y = y + h*d1yi + ((h**2)/2)*d2yi
        x = x + h
        tabla[i] = [x,y,d1yi,d2yi]
    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', d2y = y''
d1y = lambda x,y: y - x**2 + x + 1
d2y = lambda x,y: y - x**2 - x + 2
x0 = 0
y0 = 1
h = 0.1
muestras = 5

# PROCEDIMIENTO
tabla = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = tabla[:,0]
yi = tabla[:,1]

# SALIDA
print(' EDO con Taylor 3 términos')
print('[xi, yi, d1yi, d2yi]')
print(tabla)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi,yi)
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show() # plt.show() #comentar para la siguiente gráfica

Tarea: Realizar el ejercicio con más puntos muestra, donde se visualice que el error aumenta al aumentar la distancia del punto inicial [x₀,y₀]

[ EDO Taylor ] [ 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

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 - yi
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,linestyle='dashed', label='y solución conocida')
plt.legend()
plt.show()

Se puede observar los siguientes resultados:

Error máximo estimado:  0.0012745047595
entre puntos: 
[ 0.  0.000170  0.000377  0.000626  0.000922  0.00127 ]

[ EDO Taylor ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]