[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
Referencia: Sympy ODE solver. Stewart James. Cálculo de varias variables. 17.2 p1148.
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 ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
..
1. Ejemplo. Ecuación diferencial de un circuito RLC
Referencia: Lathi Ejemplo 2.1.a p155, Ejemplo 2.9 p175, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida
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.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}
1.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')
>>>
1.3 Ecuación homogénea
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 homogénea.
\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 ecuación homogénea 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
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
..
1.4. Ecuación auxiliar o característica.
En la ecuación homogénea , 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 así 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
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
..
1.5 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
>>>
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
..
2. Algoritmo 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)
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
..
3. 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:
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()
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
4. 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()
[ ejemplo ] ecuación: [ auxiliar ] [ general, complementaria ] [ algoritmo ] [ gráfica ]
