6.4 Métodos EDO con gráficos animados en Python

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

Edo con Taylor de 3 términos GIF animado

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 Runge-Kutta 2do orden primera derivada _animado

 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
.

Edo Presa Predador GIF animado

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 ]