3.2.1 LTI CT – Respuesta a entrada cero con Sympy-Python

Referencia: Lathi 1.8-1 p111, Ejemplo 2.1.a p155. Oppenheim problema 2.61c p164, Ejemplo 9.24 p700

Para encontrar la  Respuesta a entrada cero, del sistema en el ejemplo 1 del Modelo entrada-salida,
representado por el circuito y  la ecuación diferencial expresada desde el termino de mayor orden, y haciendo x(t)=0:

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

Adicionalmente se indican las condiciones iniciales y'(0)=-5, y(0)=0. La expresión en operadores D queda:

(D^2 + 3D +2)y(t) = 0

..


Ejemplo 1. Respuesta entrada cero – Desarrollo analítico con Sympy-Python

Referencia: Lathi Ejemplo 2.1.a p155

La librería Sympy-Python permite usar expresiones de derivadas en forma simbólica y resolverlas siguiendo los pasos seguidos en la solución analítica.

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

Se define la variable independiente t como un símbolo y la función matemática y(t) como función resolver. La ecuación diferencial se escribe siguiendo el orden  de los lados de expresión como la parte izquierda, LHS, y la parte derecha, RHS.

import sympy as sym

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

# ecuacion EDO: 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,evaluate=False)
x0  = RHS.subs(x(t),0).doit()
ecuacion = sym.Eq(LHS,x0)

# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [-5,0]

Las condiciones iniciales y’(0)=-5, y(0)=0 , se agrupan en una lista cond_inicio=[y'(t0),y(t0)] siguiendo el orden descendente de derivada , semejante al orden seguido al escribir la ecuación diferencial. El orden usado facilita la compatibilidad con las soluciones numéricas a realizar con Scipy-Python.

Solución general de la ecuación diferencial

La solución general se obtiene aplicando la instrucción sym.dsolve() a la ecuación, resolviendo para y(t). La expresión se trabaja mejor en Sympy si se usan suma de  términos mas simples, para lo que se añade la instrucción sym.expand().

general = sym.dsolve(ecuacion, y(t))
general = general.expand()

El resultado de la solución general se muestra con sym.pprint():

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

el mismo que en notación matemática es:

y(t) = C_1 e^{-t} + C_2 e^{-2t}

Solución particular de la ecuación diferencial

Las condiciones iniciales se aplican a la solución general para encontrar la solución particular. Para aplicar cada condición inicial se usa el lado derecho (.rhs )de la ecuación general, por ejemplo:

y'_0 (0) = -5 y'(t)\Big|_{t=0} = \frac{d}{dt}\Big[ c_1 e^{-t} + c_2 e^{-2t})\Big]\Big|_{t=0} =\Big[ -c_1 e^{-t} -2c_2 e^{-2}\Big] \Big|_{t=0} = -c_1 e^{-0} -2c_2 e^{-2(0)} -5 = -c_1 -2c_2

(realizar lo mismo para y(0)=0)

Cada condición inicial genera una ecuación que se agrupa en un sistema de ecuaciones, eq_condicion=[]. El sistema de ecuaciones se resuelva mediante la instrucción sym.solve(), obteniendo los valores para las incógnitas C1 y C2 en un dicionario.

# EDO Solucion particular
grado = len(cond_inicio)
eq_condicion = []
for k in range(0,grado,1):
    ck  = cond_inicio[(grado-1)-k]
    dyk = general.rhs.diff(t,k)
    dyk = dyk.subs(t,t0)
    eq_k = sym.Eq(ck,dyk)
    eq_condicion.append(eq_k)
    
constante = sym.solve(eq_condicion)

teniendo como resultados, las ecuaciones de condiciones iniciales y los resultados de las constantes,

 0 =  C1 +   C2
-5 = -C1 - 2*C2
{C1: -5, C2: 5}

La solución particular de la ecuación diferencial se encuentra al sustituir las constantes encontradas en la ecuación general.

particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

obteniendo como solución Particular mostrada con sym.pprint():

            -t      -2*t
y(t) = - 5*e   + 5*e
y(t) = -5 e^{-t} + 5 e^{-2t}

Instrucciones en Python

Agrupando las instrucciones descritas se presenta la solución a la ecuación diferencial ordinaria usando las condiciones iniciales.

# Respuesta a entrada cero con Sympy
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
               'numpy']

# INGRESO
t = sym.Symbol('t')
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,evaluate=False)
x0  = RHS.subs(x(t),0).doit()
ecuacion = sym.Eq(LHS,x0)

# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [-5,0]

# PROCEDIMIENTO
# EDO solución general
general = sym.dsolve(ecuacion, y(t))
general = general.expand()

# EDO Solucion particular
grado = len(cond_inicio)
eq_condicion = []
for k in range(0,grado,1):
    ck  = cond_inicio[(grado-1)-k]
    dyk = general.rhs.diff(t,k)
    dyk = dyk.subs(t,t0)
    eq_k = sym.Eq(ck,dyk)
    eq_condicion.append(eq_k)
    
constante = sym.solve(eq_condicion)

# reemplaza valores constantes en general
particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

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

print('\nSolución general: ')
sym.pprint(general)

print('\nEcuaciones de condiciones iniciales:')
for eq_k in eq_condicion:
    sym.pprint(eq_k)               
print('constantes encontradas: ', constante)
print('\nSolución particular:')
sym.pprint(particular)

resultado del algoritmo:

La ecuación diferencial a resolver: 
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           

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

Ecuaciones de condiciones iniciales:
0 = C1 + C2
-5 = -C1 - 2*C2
constantes encontradas:  {C1: -5, C2: 5}

Solución particular:
            -t      -2*t
y(t) = - 5*e   + 5*e    
>>> 

gráfica de la solución particular con Sympy

Adicionalmente se puede realizar la gráfica, convirtiendo la ecuación de forma simbólica a forma numérica numpy con sym.lambdify(t,particular,’numpy’).

Se establece el intervalo a=0 y b=5 de observación para evaluar la función y se envia a graficar, con suficiente número de muestras para que la gráfica sea suave

Por lo que se añaden las siguientes instrucciones al algoritmo anterior:

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

# INGRESO
a = 0 ; b = 5 # intervalo tiempo [a,b]
muestras = 51

yt = sym.lambdify(t, particular.rhs,
                  modules=equivalentes)

# PROCEDIMIENTO
ti = np.linspace(a,b,muestras)
yi = yt(ti)

# grafica
plt.plot(ti,yi, color='orange', label='y0(t)')
untitulo = 'Respuesta entrada cero: y0(t)='
plt.title(untitulo + str(particular.rhs))
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend()
plt.grid()
plt.show()

Respuesta entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]


Ejemplo 2. EDO con raíces repetidas

Referencia: Lathi Ejemplo 2.1.b p155

Encontrar y0(t), la respuesta a entrada cero para un sistema LTI CT descrito  por:

(D^2 + 6D +9)y(t) = (3D+5)x(t)

con condiciones iniciales y(0) = 3, y'(0)=-7

El sistema tiene un polinomio característico de raíces repetidas. El resultado usando el algoritmo se muestra a continuación:

La ecuación diferencial a resolver: 
                        2          
           d           d           
9*y(t) + 6*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           

Solución general: 
           -3*t         -3*t
y(t) = C1*e     + C2*t*e    

Ecuaciones de condiciones iniciales:
3 = C1
-7 = -3*C1 + C2
constantes encontradas:  {C1: 3, C2: 2}

Solución particular:
            -3*t      -3*t
y(t) = 2*t*e     + 3*e    

Instrucciones en Python

# Respuesta a entrada cero con Sympy
# Lathi 2.1.b pdf 155, (D^2+ 6D + 9)y = (3D+5)x
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
               'numpy']
# INGRESO
t = sym.Symbol('t')
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) + 6*sym.diff(y(t),t,1) + 9*y(t)
RHS = 3*sym.diff(x(t),t,1,evaluate=False)+5*x(t)
x0  = RHS.subs(x(t),0).doit()
ecuacion = sym.Eq(LHS,x0)

# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [-7,3]

# Grafica: intervalo tiempo [a,b]
a = 0 ; b = 5
muestras = 51

# PROCEDIMIENTO
# EDO solución general
general = sym.dsolve(ecuacion, y(t))
general = general.expand()

# EDO Solucion particular
grado = len(cond_inicio)
eq_condicion = []
for k in range(0,grado,1):
    ck  = cond_inicio[(grado-1)-k]
    dyk = general.rhs.diff(t,k)
    dyk = dyk.subs(t,t0)
    eq_k = sym.Eq(ck,dyk)
    eq_condicion.append(eq_k)
    
constante = sym.solve(eq_condicion)

# reemplaza valores constantes en general
particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

# Para graficar la Salida
yt = sym.lambdify(t, particular.rhs,
                  modules=equivalentes)
ti = np.linspace(a,b,muestras)
yi = yt(ti)

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

print('\nSolución general: ')
sym.pprint(general)

print('\nEcuaciones de condiciones iniciales:')
for eq_k in eq_condicion:
    sym.pprint(eq_k)               
print('constantes encontradas: ', constante)
print('\nSolución particular:')
sym.pprint(particular)

# GRAFICA ------------------
plt.plot(ti,yi, color='orange', label='y0(t)')
untitulo = 'Respuesta entrada cero: y0(t)='
plt.title(untitulo + str(particular.rhs))
plt.xlabel('t')
plt.ylabel('y(t)')
plt.legend()
plt.grid()
plt.show()

Ejemplo 3. EDO con raíces complejas

Referencia: Lathi Ejemplo 2.1.c p155

Encontrar y0(t), la respuesta a entrada cero para un sistema LTI CT descrito  por:

(D^2 + 4D +40)y(t) = (D+2)x(t)

con condiciones iniciales y(0) = 2, y'(0)=16.78

Tarea: Use el algoritmo del ejercicio anterior y modifique lo necesario para realizar el ejercicio. Realice el desarrollo analítico en papel y lápiz y compare. En caso de ser necesario, proponga mejoras al algoritmo presentado.

Los resultados con el algoritmo son:

La ecuación diferencial a resolver: 
                         2          
            d           d           
40*y(t) + 4*--(y(t)) + ---(y(t)) = 0
            dt           2          
                       dt           

Solución general: 
           -2*t                -2*t         
y(t) = C1*e    *sin(6*t) + C2*e    *cos(6*t)

Ecuaciones de condiciones iniciales:
2 = C2
16.78 = 6*C1 - 2*C2
constantes encontradas:  {C1: 3.46333333333333, C2: 2.00000}

Solución particular:
                         -2*t                 -2*t         
y(t) = 3.46333333333333*e    *sin(6*t) + 2.0*e    *cos(6*t)