Autor: Edison Del Rosario

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.

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 y_c

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

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

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 ordinara EDO requieren identificar el tipo de la ecuación para establecer el método de solución. Sympy incorpora librerias 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.

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 y_c(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 se x(t) = 0

Es necesario crear la ecuación homogenea 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 y_p(t)

En el caso de la solución particular, se simplifica la ecuación diferencial al sustituir x(t) con la entrada particular x_p(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 homogenea 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

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