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.
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
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
Instrucciones en Python
