6.5.1 EDO lineal – ecuaciones auxiliar, general y complementaria con Sympy-Python

Referencia: Sympy ODE solver. Stewart James. Cálculo de varias variables. 17.2 p1148 pdf545. Lathi Ejemplo 2.1.a p155, Ejemplo 2.9 p175, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida,

Continuando con el ejercicio de la sección anterior de una ecuación diferencial ordinaria, lineal y de coeficientes constantes, se obtuvo la solución homogénea o ecuación complementaria con dsolve().  Para el desarrollo analítico de la solución se requieren describir los pasos tales como la ecuación auxiliar, como se describe a continuación como un ejercicio de análisis de expresiones con Sympy.

Ejemplo 1. Ecuación diferencial de un circuito RLC, ecuación complementaria.circuito RLC

El circuito RLC con entrada x(t) de la figura tiene una corriente y(t) o salida del sistema que se representa  por medio de una ecuación diferencial.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = \frac{dx(t)}{dt}

Las condiciones iniciales del sistema a tiempo t=0 son y(0)=0 , y'(0)=-5

1. Ecuación diferencial y condiciones de frontera o iniciales

La ecuación diferencial del ejercicio con Sympy se escribe con el operador sym.diff(y,t,1). Indicando la variable independiente ‘t‘, que ‘y‘ es una función y el orden de la derivada con 1. esquema que se mantiene para la descripción de las condiciones iniciales.

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

2. Clasificar la ecuación diferencial ordinaria

Sympy permite revisar el tipo de EDO a resolver usando la instrucción:

>>> sym.classify_ode(ecuacion, y(t))
('factorable', 
 'nth_linear_constant_coeff_variation_of_parameters', 
 'nth_linear_constant_coeff_variation_of_parameters_Integral')
>>>

3. Ecuación homogenea

Para el análisis de la respuesta del circuito, se inicia haciendo la entrada x(t)=0, también conocida como ecuación para «respuesta a entrada cero» o ecuación homogenea.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = 0

En el algoritmo, La ecuación homogénea se obtiene al substituir x(t)=0 en cada lado de la ecuación, que también es un paso para encontrar la respuesta a entrada cero. La ecuacion homogenea se escribe de la forma f(t)=0.

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

4. Ecuación auxiliar o característica.

En la ecuación homogenea , se procede a sustituir en la ecuación el operador dy/dt de la derivada con una variable r elevada al orden de la derivada de cada término . El resultado es una ecuación algebraica que se analiza encontrando las raíces de r.

r ^2 + 3r +2 = 0

Los valores de ‘r’ que resuelven la ecuación permiten estimar la forma de la solución para y(t) conocida como la ecuación general

Se usan diferentes formas para mostrar la ecuación auxiliar, pues en algunos casos se requiere usar la forma de factores, en otros casos los valores de las raíces y las veces que se producen. Formas usadas para generar diagramas de polos y ceros, o expresiones de transformada de Laplace. Motivo por el que los resultados se los presenta en un diccionario, y asi usar la respuesta que sea de interés en cada caso.

homogenea :
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
auxiliar : r**2 + 3*r + 2
Q : r**2 + 3*r + 2
Q_factor : (r + 1)*(r + 2)
Q_raiz : {-1: 1, -2: 1}
>>> 

Instrucciones con Python

# Ecuación Diferencial Ordinaria EDO
# ecuación auxiliar o característica 
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# PROCEDIMIENTO
def edo_lineal_auxiliar(ecuacion,
                 t = sym.Symbol('t'),r = sym.Symbol('r'),
                 y = sym.Function('y'),x = sym.Function('x')):
    ''' ecuacion auxiliar o caracteristica de EDO
        t independiente
    '''
    # ecuación homogénea x(t)=0, entrada cero
    RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
    LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
    homogenea = LHSx0 - RHSx0
    homogenea = sym.expand(homogenea,t)

    # ecuación auxiliar o característica
    Q = 0*r
    term_suma = sym.Add.make_args(homogenea)
    for term_k in term_suma:
        orden_k = sym.ode_order(term_k,y)
        coef = 1 # coefientes del término suma
        factor_mul = sym.Mul.make_args(term_k)
        for factor_k in factor_mul:
            cond = factor_k.has(sym.Derivative)
            cond = cond or factor_k.has(y(t))
            if not(cond):
                coef = coef*factor_k
        Q = Q + coef*(r**orden_k)
               
    # Q factores y raices
    Q_factor = sym.factor(Q,r)
    Q_poly   = sym.poly(Q,r)
    Q_raiz   = sym.roots(Q_poly)
    
    auxiliar = {'homogenea' : sym.Eq(homogenea,0),
                'auxiliar'  : Q,
                'Q'         : Q,
                'Q_factor'  : Q_factor,
                'Q_raiz'    : Q_raiz }
    return(auxiliar)

# analiza la ecuación diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = edo_lineal_auxiliar(ecuacion,t,r)

# SALIDA
print('clasifica EDO:')
for untipo in edo_tipo:
    print(' ',untipo)
print('ecuacion auxiliar:')
for entrada in auxiliar:
    print(' ',entrada,':',auxiliar[entrada])

Otra forma de mostrar el resultado es usando un procedimiento creado para mostrar elementos del diccionario según corresponda a cada tipo. el procedimiento se adjunta al final.

Las funciones se incorporan a los algoritmos del curso en matg1052.py


6. Ecuación diferencial, ecuación general y complementaria con Sympy-Python

La ecuación general se encuentra con la instrucción sym.dsolve(), sin condiciones iniciales, por que el resultado presenta constantes Ci por determinar.

    # solucion general de ecuación homogénea
    general = sym.dsolve(homogenea, y(t))
    general = general.expand()

el resultado para el ejercicio es

general :
           -t       -2*t
y(t) = C1*e   + C2*e

Para encontrar los valores de las constantes, se aplica cada una de las condiciones iniciales a la ecuación general, obteniendo un sistema de ecuaciones.

     # Aplica condiciones iniciales o de frontera
    eq_condicion = []
    for cond_k in y_cond: # cada condición
        valor_k = y_cond[cond_k]
        orden_k = sym.ode_order(cond_k,y)
        if orden_k==0: # condicion frontera
            t_k = cond_k.args[0] # f(t_k)
            expr_k = general.rhs.subs(t,t_k)
        else: # orden_k>0
            # f.diff(t,orden_k).subs(t,t_k)
            subs_param = cond_k.args[2] # en valores
            t_k = subs_param.args[0]  # primer valor
            dyk = general.rhs.diff(t,orden_k)
            expr_k = dyk.subs(t,t_k)
        eq_condicion.append(sym.Eq(valor_k,expr_k))

con el siguiente resultado:

eq_condicion :
0 = C1 + C2
-5 = -C1 - 2*C2

El sistema de ecuaciones se resuelve con la instrucción

constante = sym.solve(eq_condicion)

que entrega un diccionario con cada valor de la constante

constante : {C1: -5, C2: 5}

La ecuación complementaria se obtiene al sustituir en la ecuación general los valores de las constantes.

El resultado del algoritmo

homogenea :
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
general :
           -t       -2*t
y(t) = C1*e   + C2*e    
eq_condicion :
0 = C1 + C2
-5 = -C1 - 2*C2
constante : {C1: -5, C2: 5}
complementaria :
            -t      -2*t
y(t) = - 5*e   + 5*e    
>>>

Instrucciones con Python

Las instrucciones detalladas se presentan en el algoritmo integrado.

# Solución complementaria a una Ecuación Diferencial Ordinaria EDO
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import matg1052 as fcnm

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# PROCEDIMIENTO
def edo_lineal_complemento(ecuacion,y_cond):
    # ecuación homogénea x(t)=0, entrada cero
    RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
    LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
    homogenea = LHSx0 - RHSx0

    # solucion general de ecuación homogénea
    general = sym.dsolve(homogenea, y(t))
    general = general.expand()

    # Aplica condiciones iniciales o de frontera
    eq_condicion = []
    for cond_k in y_cond: # cada condición
        valor_k = y_cond[cond_k]
        orden_k = sym.ode_order(cond_k,y)
        if orden_k==0: # condicion frontera
            t_k    = cond_k.args[0] # f(t_k)
            expr_k = general.rhs.subs(t,t_k)
        else: # orden_k>0
            subs_param = cond_k.args[2] # en valores
            t_k = subs_param.args[0]  # primer valor
            dyk = sym.diff(general.rhs,t,orden_k)
            expr_k = dyk.subs(t,t_k)
        eq_condicion.append(sym.Eq(valor_k,expr_k))

    constante = sym.solve(eq_condicion)

    # ecuacion complementaria
    # reemplaza las constantes en general
    yc = general
    for Ci in constante:
        yc = yc.subs(Ci, constante[Ci])
    
    complemento = {'homogenea'      : sym.Eq(homogenea,0),
                   'general'        : general,
                   'eq_condicion'   : eq_condicion,
                   'constante'      : constante,
                   'complementaria' : yc}
    return(complemento)

# analiza ecuacion diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = fcnm.edo_lineal_auxiliar(ecuacion,t,r)
complemento = edo_lineal_complemento(ecuacion,y_cond)
yc = complemento['complementaria'].rhs

# SALIDA
print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)

fcnm.print_resultado_dict(auxiliar)
fcnm.print_resultado_dict(complemento)

7. Gráfica de la ecuación complementaria yc

Para la gráfica se procede a convertir la ecuación yc en la forma numérica con sym.lambdify(). Se usa un intervalo de tiempo entre [0,5] y 51 muestras con el siguiente resultado:

solucion homogenea EDO

Instrucciones en Python

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

# INGRESO
t_a = 0 ; t_b = 5 # intervalo tiempo [t_a,t_b]
muestras = 51

# PROCEDIMIENTO
yt = sym.lambdify(t,yc, modules=equivalentes) 
ti = np.linspace(t_a,t_b,muestras)
yi = yt(ti)

# Grafica
plt.plot(ti,yi, color='orange', label='yc(t)')
untitulo = 'yc(t)=$'+ str(sym.latex(yc))+'$'
plt.title(untitulo)
plt.xlabel('t')
plt.ylabel('yc(t)')
plt.legend()
plt.grid()
plt.show()

8.  Mostrar el resultado del diccionario según el tipo de entrada

Para mostrar los resultados de una forma más fácil de leer, se usará sym.pprint() para las ecuaciones, y print() simple para los demás elementos del resultado

def print_resultado_dict(resultado):
    ''' print de diccionario resultado
        formato de pantalla
    '''
    for entrada in resultado:
        tipo = type(resultado[entrada])
        if tipo == sym.core.relational.Equality:
            print(entrada,':')
            sym.pprint(resultado[entrada])
        elif tipo==list or tipo==tuple:
            tipoelem = type(resultado[entrada][0])
            if tipoelem==sym.core.relational.Equality:
                print(entrada,':')
                for fila in resultado[entrada]:
                    sym.pprint(fila)
            else:
                print(entrada,':')
                for fila in resultado[entrada]:
                    print(' ',fila)
        else:
            print(entrada,':',resultado[entrada])
    return()

6.5 EDO lineal – solución complementaria y particular con Sympy-Python

Referencia: Sympy ODE solver. Stewart James. Cálculo de varias variables. 17.2 p1148 pdf545. Lathi Ejemplo 2.1.a p155, Ejemplo 2.9 p175, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida,

Los métodos analíticos para encontrar una solución y(t) a una ecuación diferencial ordinaria EDO requieren identificar el tipo de la ecuación para establecer el método de solución. Sympy incorpora librerías que pueden asistir en la solución con muchos de los métodos analíticos conocidos.

Ejemplo 1. Ecuación diferencial de un circuito RLC, ecuación complementaria.circuito RLC

El circuito RLC con entrada x(t) de la figura tiene una corriente y(t) o salida del sistema que se representa  por medio de una ecuación diferencial.

\frac{d^{2}y(t)}{dt^{2}} + 3\frac{dy(t)}{dt} + 2y(t) = \frac{dx(t)}{dt}

Las condiciones iniciales del sistema a tiempo t=0 son y(0)=0 , y'(0)=-5

Encuentre la solución considerando como entrada x(t)

x(t) = 10 e^{-3t} \mu(t)

1.1 Ecuación diferencial y condiciones de frontera o iniciales

La ecuación diferencial del ejercicio con Sympy se escribe con el operador sym.diff(y,t,k). Indicando la variable independiente, que ‘y‘ es una función y el orden de la derivada con k.

La ecuación puede escribirse como lado izquierdo (LHS) es igual al lado derecho (RHS). Una ecuación en Sympy se compone de las dos partes sym.Eq(LHS,RHS).

# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

Las condiciones de frontera o iniciales y(0)=0, y'(0)=-5 se ingresan en un diccionario en el siguiente formato

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

la entrada x(t) conocida se define como:

# ecuacion entrada x(t)
xp = 10*sym.exp(-3*t)*sym.Heaviside(t)

1.2 Clasificar la ecuación diferencial ordinaria

Para revisar el tipo de EDO a resolver se tiene la instrucción:

>>> sym.classify_ode(ecuacion, y(t))
('factorable', 
 'nth_linear_constant_coeff_variation_of_parameters', 
 'nth_linear_constant_coeff_variation_of_parameters_Integral')
>>>

Solución EDO como suma de solución complementaria y solución particular

La solución para una EDO puede presentarse también como la suma de una solución complementaria y una solución particular

y(t) = y_p(t) + y_c(t)

1.3 Solución complementaria yc(t)

La solución complementaria de la EDO se interpreta también como respuesta de entrada cero del circuito, donde la salida y(t) depende solo de las condiciones iniciales del circuito. Para el ejemplo, considera solo las energías internas almacenadas en el capacitor y el inductor. Para éste caso x(t) = 0

Es necesario crear la ecuación homogénea con x(t)=0 en la ecuación diferencial para encontrar la solución con las condiciones iniciales dadas en el enunciado del ejercicio.

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general de ecuación homogénea
yc = sym.dsolve(homogenea, y(t),ics=y_cond)
yc = yc.expand()

el resultado de la ecuación complementaria es

solucion complementaria y_c(t): 
            -t      -2*t
y(t) = - 5*e   + 5*e    

1.4 Solución particular yp(t)

En el caso de la solución particular, se simplifica la ecuación diferencial al sustituir x(t) con la entrada particular xp(t). Las condiciones iniciales en este caso consideran que el circuito no tenía energía almacenada en sus componentes (estado cero), por lo que las condiciones iniciales no se consideran.

# ecuación particular x(t)=0, estado cero
RHSxp = ecuacion.rhs.subs(x(t),x_p).doit()
LHSxp = ecuacion.lhs.subs(x(t),x_p).doit()
particular = LHSxp - RHSxp

# solucion particular de ecuación homogénea
yp = sym.dsolve(particular, y(t))

Se aplica nuevamente la búsqueda de la solución con sym.dsolve() y se obtiene la solución particular que incluye parte del modelo de la ecuación homogénea con los coeficientes Ci

>>> sym.pprint(yp.expand())
           -t       -2*t      -t                    -2*t                    -3*t
y(t) = C1*e   + C2*e     - 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)

>>> yp.free_symbols
{C2, C1, t}
>>>

En éste punto se incrementan los términos de las constantes C1 y C2 como símbolos, que para tratar exclusivamente la solución particular, se reemplazan con ceros. Para obtener las variables de la ecuación se usa la instrucción yp.free_symbols que entrega un conjunto. El conjunto se descarta el símbolo t y se sustituye todos los demás por cero en la ecuación.

    # particular sin terminos Ci
    y_Ci = yp.free_symbols
    y_Ci.remove(t) # solo Ci
    for Ci in y_Ci: 
        yp = yp.subs(Ci,0)

quedando la solución particular como:

            -t                    -2*t                    -3*t             
y(t) = - 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)

1.5 Solución EDO como suma de complementaria + particular

La solución de la ecuación diferencial ordinaria, ante la entrada x(t) y condiciones iniciales es la suma de los componentes complementaria y particular:

# solucion total = complementaria + particular
ytotal = yc.rhs + yp.rhs

El resultado de todo el algoritmo  permite interpretar la respuesta con los conceptos de las respuestas del circuito a entrada cero y estado cero, que facilitan el análisis de las soluciones en ejercicios de aplicación.

ecuación diferencial a resolver: 
                        2                 
           d           d          d       
2*y(t) + 3*--(y(t)) + ---(y(t)) = --(x(t))
           dt           2         dt      
                      dt                  
condiciones iniciales
  y(0)  =  0
  Subs(Derivative(y(t), t), t, 0)  =  -5
clasifica EDO:
  factorable
  nth_linear_constant_coeff_variation_of_parameters
  nth_linear_constant_coeff_variation_of_parameters_Integral
x(t): 
    -3*t             
10*e    *Heaviside(t)
solucion complementaria yc(t): 
            -t      -2*t
y(t) = - 5*e   + 5*e    
solucion particular yp(t): 
       /     -t       -2*t       -3*t\             
y(t) = \- 5*e   + 20*e     - 15*e    /*Heaviside(t)
solucion y(t) = yc(t) +yp(t): 
/     -t       -2*t       -3*t\                   -t      -2*t
\- 5*e   + 20*e     - 15*e    /*Heaviside(t) - 5*e   + 5*e    
>>> 

Instrucciones en Python

# Solución de unaEcuación Diferencial Ordinaria EDO
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
y = sym.Function('y')
x = sym.Function('x')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)

# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
          sym.diff(y(t),t,1).subs(t,0) : -5}

# ecuacion entrada x(t)
xp = 10*sym.exp(-3*t)*sym.Heaviside(t)

# PROCEDIMIENTO
edo_tipo = sym.classify_ode(ecuacion, y(t))

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general de ecuación homogénea
yc = sym.dsolve(homogenea, y(t),ics=y_cond)
yc = yc.expand()

def edo_lineal_particular(ecuacion,xp):
    ''' edo solucion particular con entrada x(t)
    '''
    # ecuación particular x(t)=0, estado cero
    RHSxp = ecuacion.rhs.subs(x(t),xp).doit()
    LHSxp = ecuacion.lhs.subs(x(t),xp).doit()
    particular = LHSxp - RHSxp

    # solucion particular de ecuación homogénea
    yp = sym.dsolve(particular, y(t))

    # particular sin terminos Ci
    y_Ci = yp.free_symbols
    y_Ci.remove(t) # solo Ci
    for Ci in y_Ci: 
        yp = yp.subs(Ci,0)
        
    # simplifica y(t) y agrupa por escalon unitario
    yp = sym.expand(yp.rhs,t)
    lista_escalon = yp.atoms(sym.Heaviside)
    yp = sym.collect(yp,lista_escalon)
    yp = sym.Eq(y(t),yp)
    
    return(yp)

yp = edo_lineal_particular(ecuacion,xp)
# solucion total
ytotal = yp.rhs + yc.rhs


# SALIDA
print('ecuación diferencial a resolver: ')
sym.pprint(ecuacion)

# condiciones iniciales
print('condiciones iniciales')
for caso in y_cond:
    print(' ',caso,' = ', y_cond[caso])

print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)

print('x(t): ')
sym.pprint(xp)

print('solucion complementaria yc(t): ')
sym.pprint(yc)

print('solucion particular yp(t): ')
sym.pprint(yp)

print('solucion y(t) = yc(t) +yp(t): ')
sym.pprint(ytotal)

1.6  Grafica de resultados

Se adjunta la gráfica de los resultados de las ecuaciones complementaria, particular y total

EDO lineal Complementaria Particular 01

Instrucciones adicionales en Python

# GRAFICA ------------------
import numpy as np
import matplotlib.pyplot as plt
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]

# INGRESO
t_a = 0 ; t_b = 5 # intervalo tiempo [t_a,t_b]
muestras = 51

# PROCEDIMIENTO
y_c = sym.lambdify(t,yc.rhs, modules=equivalentes)
y_p = sym.lambdify(t,yp.rhs, modules=equivalentes)
y_cp = sym.lambdify(t,ytotal, modules=equivalentes)
ti = np.linspace(t_a,t_b,muestras)
yci = y_c(ti)
ypi = y_p(ti)
ycpi = y_cp(ti)

# Grafica
plt.plot(ti,yci, color='orange', label='yc(t)')
plt.plot(ti,ypi, color='dodgerblue', label='yp(t)')
plt.plot(ti,ycpi, color='green', label='y(t)')
untitulo = 'y(t)=$'+ str(sym.latex(ytotal))+'$'
plt.title(untitulo)
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

Un desarrollo semejante del ejercicio aplicando conceptos del curso de señales y sistemas de proporciona en:

LTI CT Respuesta del Sistema Y(s)=ZIR+ZSR con Sympy-Python

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 ]

6.3 Sistemas EDO. modelo depredador-presa con Runge-Kutta y Python

Sistemas EDO [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Runge Kutta  d2y/dx2 ]
..


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  d2y/dx2 ]
..


2. Desarrollo analíticoEdo Presa Predador GIF animado

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  d2y/dx2 ]

..


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 xi, yi versus ti, es decir cantidad de animales de cada especie durante el tiempo de observación

b) Independiente de la unidad de tiempo, xi vs yi, 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  d2y/dx2 ]

6.2.2 EDO Runge-Kutta d2y/dx2 con Python

[ Runge Kutta  d2y/dx2 ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..


1. EDO Runge-Kutta para Segunda derivada d2y/dx2

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} + 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(h3)
Runge-Kutta 4do Orden tiene error de truncamiento O(h5)

[ Runge Kutta  d2y/dx2 ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]

..


2. Runge-Kutta 2do Orden para Segunda derivada d2y/dx2 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  d2y/dx2 ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..


3. Runge-Kutta 4do Orden para Segunda derivada d2y/dx2 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  d2y/dx2 ] 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  d2y/dx2 ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]

6.2.1 EDO Runge-Kutta 4to Orden dy/dx con Python

[ 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:
Runge Kutta 4to Orden

\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(h5)

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 ]

6.2 EDO Runge-Kutta 2do Orden dy/dx con Python

[ 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 K1 y K2.

EDO Runge-Kutta 2do orden primera derivada _animado

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 K1 y K2 se calculan para predecir el próximo valor en y[i+1], observe que el término K1 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 02

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

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 2do Orden primera derivada

[ 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 x0=0, y0=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(h3)

EDO Runge-Kutta 2do orden primera derivada _animado

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 , x0 = 0, y0 = 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 , x1 = 0.1, y1 = 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 , x2 = 0.2, y2 = 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]]
>>> 

ecuación diferencial ordinaria con Runge-Kutta de 2do orden

Compare los resultados con Taylor de 2 y 3 términos.

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

[ 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

EDO RungeKutta 2Orden 01

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.

6.1 EDO con Taylor de 3 términos con Python

[ 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)

Edo Taylor 3 términos GIF animado
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.

EDO Taylor 3 terminos

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 x0=0, y0=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 yi+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

Edo Taylor 3 términos GIF animado

Se empieza evaluando el nuevo punto a una distancia x1= x0+h del punto de origen con lo que se obtiene y1 , 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.

EDO_Taylor_3terminos01

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 [x0,y0]

# 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 [x0,y0]

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