La respuesta total del sistema integra las respuestas obtenidas para entrada cero y estado cero. Se resume el desarrollo de un ejercicio con todos los componentes para el ejemplo 1 Modelo entrada-salida para circuitos RLC realizada con Sympy-Python y las funciones del curso en telg1001.py.
Para la respuesta se considera que la respuesta total del sistema se puede describir en forma gráfica como el resultado de dos componentes:
Se plantea el ejercicio para encontrar la respuesta a entrada cero ZIR y la respuesta a estado cero ZSR
Es de considerar que para las gráficas de las respuestas se considera que el sistema es observable, aplicable, usable desde t=0, dado que la solución para entrada cero no se ha acotado en los resultados. Imagine sistema con un capacitor entre sus componentes que a la salida que se descarga desde un voltaje de 1 voltio desde t=0, esta situación no necesariamente implica que podríamos conocer si viene descargandose desde 1.5, 3, 9, 12, o 1000 voltios. Observación a considerar con en la interpretación de los ejercicios para valores obtenidos antes de t=0 que es el punto inicial de observacion.
# Respuesta total del sistema# y(t) = ZIR(t) + ZSR(t)# https://blog.espol.edu.ec/telg1001/lti-ct-yszirzsr-respuesta-del-sistema-con-sympy-python/# Revisar causalidad de x(t) y h(t)import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)
# 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)
ecuacion = sym.Eq(LHS,RHS)
# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [-5,0]
# entrada x(t)
x = 10*sym.exp(-3*t)*u
# grafica intervalo [t_a,t_b]
t_a = 0; t_b = 5
muestras = 201
# PROCEDIMIENTO# Respuesta entrada cero ZIR
sol_ZIR = fcnm.respuesta_ZIR(ecuacion,cond_inicio,t0)
ZIR = sol_ZIR['ZIR']
# Respuesta al impulso h(t)
sol_h = fcnm.respuesta_impulso_h(ecuacion)
h = sol_h['h']
# respuesta a estado cero ZSR
sol_ZSR = fcnm.respuesta_ZSR(x,h)
ZSR = sol_ZSR['ZSR']
xh = sol_ZSR['xh']
# respuesta a y(t) = ZIR(t)+ZSR(t)
y_total = ZIR+ZSR
# revisa si grafica ZSRifnot(sol_ZSR['cond_graf']):
print('revisar acortar x(t) o h(t) para BIBO')
# SALIDAprint('\n ZIR(t):')
sym.pprint(ZIR)
print('\n h(t):')
sym.pprint(h)
print('\n ZSR(t):')
sym.pprint(ZSR)
print('xcausal: ',sol_ZSR['xcausal'])
print('hcausal: ',sol_ZSR['hcausal'])
print('limites de integral:', sol_ZSR['[tau_a,tau_b]'])
print('\n y(t) = ZIR(t)+ZSR(t):')
sym.pprint(y_total)
# GRAFICA
figura_ZIR = fcnm.graficar_ft(ZIR,t_a,t_b,
muestras,'ZIR')
figura_h = fcnm.graficar_ft(h,t_a,t_b,
muestras,'h')
# grafica animada de convolución
n_archivo = ''# sin crear archivo gif animado # n_archivo = 'LTIC_YTotal_ZIR_Ej01' # requiere 'imagemagick'if sol_ZSR['cond_graf']:
fig_ZSR = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
muestras,y_nombre='ZSR')
fig_ytotal = fcnm.graficar_xhy(ZIR,ZSR,y_total,t_a,t_b,
muestras,x_nombre='ZIR',
h_nombre='ZSR',y_nombre='y')
figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,-t_b,t_b,
muestras, reprod_x = 4,y_nombre='ZSR',
archivo_nombre = n_archivo)
plt.show()
El algoritmo en Sympy incorpora los pasos realizados con el desarrollo analítico de solución a ecuaciones diferenciales lineales.
Existen otros métodos como el de Transformada de Laplace que simplifica el paso de realizar el integral de convolución, pues al cambiar al dominio ‘s’ en lugar de tiempo las operaciones son mayoritariamete sobre polinomios. La siguiente unidad desarrolla el tema.
El ejemplo es continuación del presentado para respuesta a entrada cero, que tiene la ecuación:
(D^2 + 3D +2)y(t) = Dx(t)
Al diagrama desarrollado en la sección Respuesta entrada cero – Diagrama Bloques, se añaden los bloques para formar la señal de entrada x(t). La variable tiempo se indica con el Bloque Reloj, el escalón μ(t) se configura para empezar en cero y el valor -3 del exponente se los indica con un bloque de ganancia antes de eu.
Para diferenciación de los bloques por sección, la entrada se usa el color azul, el sistema en color naranja y la salida en color morado.
Para continuar con la simulación, asegúrese que no existan condiciones iniciales en el sistema, es decir los valores de estado inicial para los bloques 1/s se encuentren en cero.
Los osciloscopios de observación, tienen su propia referencia de tiempo con otro reloj. El osciloscopio de la señal de entrada x(t) se muestra primero, y la salida en la gráfica posterior.
La señal de salida se obtiene de forma semejante a la respuesta a impulso, es decir en el punto donde está la primera derivada. Revisar teoría para aplicar la forma: P(D)[y(t)μ(t)], que en éste ejercicio se aplica una D.
El resultado de y(t), solo como referencia, es el que se obtiene en el osciloscopio ubicado a la derecha abajo del diagrama.
Con lo que se muestra que es posible observar resultados usando algunos bloques en un simulador.
2. La integral de convolución con Numpy.convolve()
El algoritmo requiere la expresión de las funciones x(t) y h(t) para determinar y(t). Se realiza el cálculo de la convolución, que para aproximar a la integral requiere multiplicar por dt.
Note que debe mantener el valor dt pequeño, en el ejercicio se ha definido como 0.01 que corresponden a casi 1000 muestras en el intervalo [-5,5].
Para usar la función np.convolve(), debido a la naturaleza numérica del algoritmo, se requieren muestras en tiempos, simétricos alrededor de cero. Al mostrar la gráfica, se puede prescindir de la sección negativa, sin embargo, para recordarlo se ha dejado la parte negativa del intervalo de tiempo.
# Señales contínuas en convolución# Compare con el algoritmo en forma discreta# Lathi ejemplo 2.6 p166import numpy as np
import matplotlib.pyplot as plt
# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)
# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01
# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)
# Integral de Convolucion x[t]*h[t]# corrección de magnitud por dt para en integral
yi = np.convolve(xi,hi,'same')*dt
# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')
plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()
3. La integral de convolución con algoritmo numérico
Esta sección desarrolla el algoritmo para la integral de convolución, tomando paso a paso las operaciones para obtener la gráfica del resultado.
El algoritmo se presenta con propósitos didácticos, pues parece ser un poco más lento debido a que se realizan las operaciones de forma básica (interpretada) de Python. La función np.convolve() usa elementos compilados a lenguaje de maquina que se ejecutan de forma más eficiente.
Note que la única sección que cambia es la de convolución, por lo que la gráfica de salida es la misma que la sección anterior:
# Señales contínuas en convolución. Algoritmo numérico.# compare con el algoritmo en forma discreta# Lathi ejemplo 2.6 p166import numpy as np
import matplotlib.pyplot as plt
# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)
# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01
# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)
# Integral de Convolución x[t]*h[t]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i inrange(0,muestras):
suma=0
for k inrange(0,muestras):
suma = suma + x(ti[k])*h(ti[i]-ti[k])
yi[i]= suma*dt
# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')
plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()
Se desarrolla el integral de convolución con Sympy con los criterios indicados en el desarrollo analítico. Para determinar los límites del integral se requiere revisar la Causalidad de h(t) descrita con la función es_causal(ft). Los pasos para los resultados se completan con la función desarrollada en la sección anterior para Integral de convolución. ..
Ejemplo 1. Convolución con x(t) causal y h(t) causal
El límite inferior del integral se usa como 0–, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.
Para iniciar el algoritmo se definen las variables t y tau, simplificando la expresión del escalón unitario o Heaviside con u.
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
Se definen las señales de entrada x(t) y respuesta al impulso h(t) usando las expresiones en Sympy. Para seleccionar los límites del integral de convolución se indica si x(t) y h(t) son causales.
# entrada x(t), respuesta impulso h(t)
x = 10*sym.exp(-3*t)*u
h = (2*sym.exp(-2*t)-sym.exp(-t))*u
Se realiza la comprobación de causalidad revisando cada término de la función, que contenga un escalón o un impulso, tomando como referencia la forma en que se construye h(t) con el impulso y los modos caraterísticos.
En el procedimiento se definen los límites del integral de convolución, la expresión dentro del integral xh(t) y con el resultado se aplica la instrucción expand() para obtener términos suma mas sencillos.
En el caso que la respuesta al impulso h(t) no sea causal y la entrada x(t) es causal, se intercambian las funciones para mantener la simplicidad del integral con límite superior t en lugar de infinito. Se aplica la propiedad conmutativa de la convolución.
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
Se añaden las instrucciones al algoritmo del ejercicio para la convolución de la sección anterior y convierte a una función respuesta_ZSR(x,h) .
Con los resultados del algoritmo, queda añadir las instrucciones para mostrar las ecuaciones y las gráficas como los siguientes:
El integral de convolución de dos funciones x(t) y h(t) usa como operador el símbolo ⊗ que representa el integral de la ecuación mostrada.
La función h(t) representa la repuesta del sistema, función de transferencia o respuesta a impulso, que por la forma de obtenerla se considera causal por tener componente μ(t) y δ(t) que aseguraría que la función se desarrolle en el lado derecho del plano.
La función x(t) representa la entrada del sistema y el integral de convolución se utiliza para determinar la respuesta de estado cero de un sistema ZSR.
Si las funciones x(t) y el sistema h(t) son causales, la respuesta también será causal. Teniendo como límites del integral τa = 0 y τb = t
El límite inferior del integral se usa como 0–, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.
Ejemplo 1. Algoritmo para el integral de convolución con x(t) y h(t) causales con Sympy
Realizar la convolución entre la respuesta al impulso h(t) y la entrada x(t) mostradas:
x(t) = e^{-t} \mu(t) h(t) = e^{-2t} \mu(t)
Para el algoritmo, se define las variables a usar, escalón e impulso unitarios, definiendo las funciones x(t) y y(t) como:
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-t)*u
h = sym.exp(-2*t)*u
Se revisa la causalidad de las señales para actualizar los límites del integral de convolución.
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
El integral de convolución se aplica a la multiplicación entre x(t) y h(t), realizando el cambio de variable por τ y (t-τ). Para facilitar la operación del integral, se simplifica la expresión xh como términos de suma
# integral de convolucion
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
y = sym.integrate(xh,(tau,tau_a,tau_b))
y = sym.expand(y,t)
y = y.subs(u**2,u) # Heaviside**2=Heaviside
El resultado del integral se muestra como,
xcausal: True
hcausal: True
limites de integral: [ 0 , t ]
x(t)*h(t-tau):
-2*t tau
e *e *Heaviside(tau)*Heaviside(t - tau)
y(t):
-t -2*t
e *Heaviside(t) - e *Heaviside(t)
Instrucciones en Python
# Integral de convolucion con Sympy# Revisar causalidad de x(t) y h(t)import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-t)*u
h = sym.exp(-2*t)*u
# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 201
# PROCEDIMIENTO
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
# integral de convolucion x(t)*h(t)
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
y = sym.integrate(xh,(tau,tau_a,tau_b))
y = sym.expand(y,t)
ifnot(y.has(sym.Integral)):
y = fcnm.simplifica_escalon(y)
#ZSR = ZSR.subs(u**2,u) # Heaviside**2=Heaviside# SALIDAprint('xcausal: ',xcausal)
print('hcausal: ',hcausal)
print('limites de integral: [',tau_a,',',tau_b,']')
print('x(t)*h(t-tau):')
sym.pprint(xh)
print('\n y(t):')
sym.pprint(y)
if hcausal==False:
print('revisar causalidad de h(t)')
if xcausal==False:
print('revisar causalidad de x(t)')
Grafica para el integral de convolución entre x(t) y h(t)
Con los resultados se puede construir la gráfica, para observar en un intervalo simétrico que permita visualizar la operación de convolución.
Instrucciones complementarias para las gráficas, se realiza como una función a ser añadida a telg1001.py por ser reutilizada en los ejercicios.
# GRAFICA -----defgraficar_xh_y(xt,ht,yt,t_a,t_b,
muestras=101,x_nombre='x',
h_nombre='h',y_nombre='y'):
'''dos subgraficas, x(t) y h(t) en superior
h(t) nn inferior
'''# grafica evaluacion numerica
x_t = sym.lambdify(t,xt,modules=equivalentes)
h_t = sym.lambdify(t,ht,modules=equivalentes)
ti = np.linspace(t_a,t_b,muestras)
xi = x_t(ti)
hi = h_t(ti)
y_t = sym.lambdify(t,yt,modules=equivalentes)
yi = y_t(ti)
colorlinea_y = 'green'if y_nombre =='ZSR':
colorlinea_y = 'dodgerblue'
fig_xh_y, graf2 = plt.subplots(2,1)
untitulo = y_nombre+'(t) = $'+ str(sym.latex(yt))+'$'
graf2[0].set_title(untitulo)
graf2[0].plot(ti,xi, color='blue', label='x(t)')
graf2[0].plot(ti,hi, color='magenta', label='h(t)')
graf2[0].legend()
graf2[0].grid()
graf2[1].plot(ti,yi,colorlinea_y,
label = y_nombre+'(t)')
graf2[1].set_xlabel('t')
graf2[1].legend()
graf2[1].grid()
#plt.show()return(fig_xh_y)
figura_xh_y = graficar_xh_y(x,h,y,t_a,t_b,muestras)
plt.show()
Ejercicio 2. Convolución entre h(t) causal y x(t) no causal
Referencia: Opennheim 2.8 p101
Realizar la convolución entre h(t) y x(t) mostradas:
x(t) = e^{2t} \mu (-t) h(t) = \mu (t-3)
La señal x(t) es creciente desde el lado izquierdo del plano, es decir es no causal, también es acotada en cero por el escalón unitario invertido en el tiempo.
Mientras que h(t) presenta una respuesta con retraso de 3 unidades de tiempo con el escalón unitario, se encuentra en el lado derecho del plano, por lo que se la considera causal.
Para el análisis, el gráfico se presenta la animación de la operación de convolución
Para el ejercicio, es importante considerar que en un sistema lineal causal, solo se desarrolla una salida y(t) ante una entrada x(t). La señal de entrada x(t) no se registra a partir de tiempo con valor cero. Observe en la gráfica que la señal de salida y(t) se desarrolla con retardo de tres unidades respecto a la entrada.
Desde luego el sistema responderá en cualquier momento que se presenta la señal, por lo que x(t) podría ser no causal, La referencia de t=0 se considera desde el punto de vista del sistema h(t).
# entrada x(t), respuesta impulso h(t)
x = sym.exp(-2*t)*u
h = u
# grafica intervalo [t_a,t_b] plano simétrico
t_b = 5 ; t_a = -t_b
muestras = int(t_b)*10+1
Esto implica que de darse el caso que exista un h(t) no causal y x(t) causal, por la propiedad conmutativa de la convolución, para el desarrollo del integral se pueden intercambiar las entradas y hacerlo mas simple.
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True
El resultado con el algoritmo es:
xcausal: False
hcausal: True
limites de integral: [ -oo , t ]
x(t)*h(t-tau):
2*tau
e *Heaviside(-tau)*Heaviside(t - tau - 3)
ZSR(t):
/ -6 2*t \
|e *e 1| 1
|-------- - -|*Heaviside(3 - t) + -
\ 2 2/ 2
revisar causalidad de x(t)
En la respuesta se observa que el témino u(3-t)/2 en -∞ se convierte en 1/2 que anula la constante de 1/2 del tercer término. Al revisar la convergencia del primer término se encuentra que tiende a cero, en consecuencia en la señal hacia el lado izquierdo tiende a cero.
gráfica con el algoritmo
Instrucciones en Python
Para los ejercicios se agrupan las instrucciones en una función respuesta_ZSR(x,h) que recibe las dos señales x(t) y h(t), como se desarrolla en la siguiente sección. En la función se analiza con los algoritmos anteriores si son de tipo causal incluida en telg1001.py.
# Integral de convolucion con Sympy# como Respuesta a estado cero ZSR con x(t) y h(t)# https://blog.espol.edu.ec/telg1001/lti-ct-respuesta-a-estado-cero-con-sympy-python/import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
tau = sym.Symbol('tau',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t), respuesta impulso h(t)
x = sym.exp(2*t)*u.subs(t,-t)
h = u.subs(t,t-3)
# grafica intervalo [t_a,t_b] plano simétrico
t_b = 7 ; t_a = -t_b
muestras = int(t_b)*10+1
# PROCEDIMIENTO# revisa causalidad de señales
xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)
# intercambia si h(t) no es_causal# con x(t) es_causal por propiedad conmutativa
intercambia = Falseif hcausal==Falseand xcausal==True:
temporal = h
h = x
x = temporal
xcausal = False
hcausal = True
intercambia = True# limites de integral de convolución
tau_a = -sym.oo ; tau_b = sym.oo
if hcausal==True:
tau_b = t
if (xcausal and hcausal)==True:
tau_a = 0
# integral de convolución x(t)*h(t)
xh = x.subs(t,tau)*h.subs(t,t-tau)
xh = sym.expand(xh,t)
ZSR = sym.integrate(xh,(tau,tau_a,tau_b))
ZSR = sym.expand(ZSR,t)
ifnot(ZSR.has(sym.Integral)):
ZSR = fcnm.simplifica_escalon(ZSR)
#ZSR = ZSR.subs(u**2,u) # Heaviside**2=Heaviside
lista_escalon = ZSR.atoms(sym.Heaviside)
ZSR = sym.expand(ZSR,t) # terminos suma
ZSR = sym.collect(ZSR,lista_escalon)
# restaura si hubo intercambio de x(t) y h(t)if intercambia == True:
xcausal = True
hcausal = False# graficar si no tiene Integral o error
cond_graf = not(ZSR.has(sym.Integral))
cond_graf = cond_graf andnot(ZSR.has(sym.oo))
cond_graf = cond_graf andnot(ZSR.has(sym.nan))
# SALIDAprint('xcausal: ',xcausal)
print('hcausal: ',hcausal)
print('limites de integral: [',tau_a,',',tau_b,']')
print('x(t)*h(t-tau):')
if xh.is_Add:
i=0
for term_suma in xh.args:
if i>0:
print('+')
sym.pprint(term_suma)
i=i+1
else:
sym.pprint(xh)
print('\n ZSR(t):')
sym.pprint(ZSR)
if hcausal==False:
print('revisar causalidad de h(t)')
if xcausal==False:
print('revisar causalidad de x(t)')
if intercambia:
print('considere intercambiar h(t)con x(t)')
ifnot(cond_graf):
print('revisar acortar x(t) o h(t) para BIBO')
# GRAFICAif cond_graf:
figura_xh_y = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
muestras,y_nombre='ZSR')
else: # terminos de integrales sin evaluar, infinito o nan
figura_x = fcnm.graficar_ft(x,t_a,t_b,muestras,'x')
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,'h')
#plt.show()# grafica animada de convolución# n_archivo = '' # sin crear archivo gif animado
n_archivo = 'convolucion02'# requiere 'imagemagick'if cond_graf:
figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,t_a,t_b,
muestras, reprod_x = 4,y_nombre='y',
archivo_nombre = n_archivo)
plt.show()
Para el desarrollo analítico, se propone considerar usar las funciones por intervalos separados:
y(t) = \begin{cases} 0 , & t \lt 0 \\ \frac{1}{2} t^2 , & 0 \lt t \lt T \\ Tt-\frac{1}{2}T^2 , & T \lt t \lt 2T \\ - \frac{1}{2}t^2 + Tt + \frac{3}{2} T^2 , & 2T \lt t \lt 3T \\ 0 . & 3T \lt t \end{cases}
Con las indicaciones se propone como tarea realizar el desarrollo analítico. Observando la gráfica animada y considerando los intervalos propuestos.
El resultado con el algoritmo del ejercicio anterior es:
En el resultado del algoritmo se encuentra que es necesario simplificar la expresión resultante al menos para los términos con escalón unitario, tales como u*u o u*u(t-1). Sympy ‘1.11.1’ aún no incorpora esta simplificación, por lo que se realiza una función que revise término a término la expresión.
La función h(t) o respuesta al impulso se compone de un término de impulso y un escalón que multiplica a los modos característicos de la ecuación diferencial.
h(t)=b_0 \delta (t)+ [P(D)y_n (t)] \mu (t)
La comprobación de causalidad considera tomar como referencia la búsqueda de un δ(t) o un μ(t) del lado derecho del plano, t≥0. Se considera que el sistema responde solo cuando le llega una señal a partir de t≥0, por lo que los componentes de h(t) deben encontrarse en el lado derecho del plano.
La causalidad de h(t) se usa para determinar los límites del integral de convolución que facilitan las operaciones para obtenerlo. El algoritmo de comprobación de causalidad se compone de las partes para analizar la expresión de h(t) descritas en la imagen.
Se propone desarrollar algunos ejercicios de forma progresiva, empezando de lo pequeño hacia lo grande.
Para generalizar el algoritmo, en adelante h(t) se denominará f(t). ..
Ejemplo 1. busca el punto t donde se produce el impulso δ(t)
Referencia: h(t) de Lathi Ej. 2.13 p200
El primer objetivo es determinar si tan solo un término simple contiene un impulso δ(t), no existen otros terminos sumados. El ´término puede estar desplazado o multiplicado por otros factores.
# entrada x(t)# h = d
h = 2*sym.pi*d.subs(t,t-2)
# h = 2*d.subs(t,t-1)# h = 3*d.subs(t,t+1)# h = 5*sym.exp(-t)*d
La instrucción para revisar que aún no existan de términos de suma es not(ft.is_Add).
Si f(t) tiene un impulso se puede revisar con ft.has(sym.DiracDelta), que inicialmente sería suficiente para los objetivos del ejercicio. Sin embargo en ejercicios posteriores y de evaluación se encontrará que estas funciones se pueden desplazar hacia un lado del plano, tal como un retraso en el tiempo de procesamiento. Por lo que se considera necesario determinar el punto t donde se aplica el impulso y observar si es en el lado derecho del plano, t≥0.
En el caso de disponer de varios valores, es de interés el valor que se encuentre más hacia la izquierda, es decir los valores encontrados deben ser ordenados previamente. Si los impulsos ocupan solamente el lado derecho del plano, el sistema es causal.
donde_d = fcnm.busca_impulso(h)
# Causal en RHS
hcausal = Falseiflen(donde_d)>0:
donde=donde_d[0] # extremo izquierdoif donde>=0: # lado derecho del plano
hcausal=True
Para el siguiente ejemplo se debería obtener:
h(t):
2*pi*DiracDelta(t - 2)
hcausal: True
con gráfica mostrando la causalidad de h(t) como un retraso de la entrada.
Instrucciones en Python
# Revisar causalidad de ft(t) en sympy# considera causal h(t)=b0*d + (modos caracteristicos)*u# Parte 1 comprueba solo impulso d(t) import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t)# h = d
h = 2*sym.pi*d.subs(t,t-2)
# h = d.subs(t,t-1)# h = 3*d.subs(t,t+1)# h = 5*sym.exp(-t)*d# h = 4# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTO
donde_d = fcnm.busca_impulso(h)
# Causal en RHS
hcausal = Falseiflen(donde_d)>0:
donde=donde_d[0] # extremo izquierdoif donde>=0: # lado derecho del plano
hcausal=True# SALIDAprint('h(t):')
sym.pprint(h)
print('hcausal: ',hcausal)
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()
..
Ejemplo 2. f(t) es causal si contiene escalón unitario μ(t) en un termino simple, sin sumas, N>M
Referencia: Lathi 2.6-3 p206, Ej.2.13 p188
Ahora considere un termino de los modos característicos multiplicado por un escalón unitario μ(t) según el modelo de respuesta a impulso h(t).
h(t)=0 \delta (t)+ [P(D)y_n (t)] \mu (t)
la causalidad del escalón require primero encontrar desde dónde se aplican y luego la dirección o sentido. De forma semejante a la búsqueda de impulsos, se crea una función que busque el escalón e indique [donde, sentido]
Para un escalón unitario hacia la derecha se tiene como punto de partida que causal=False, cambiando a True en el caso que el impulso se encuentre a la derecha del plano y el sentido sea positivo.
causal = False
h = fcnm.simplifica_escalon(h) # h(t-a)*h(t-b)if h.has(sym.Heaviside): #sin factores
donde_u = busca_escalon(h)
iflen(donde_u)>0:
donde = donde_u[0][0]
sentido = donde_u[0][1]
if donde>=0 and sentido>0:
causal = True
El resultado esperado del algoritmo, para un caso es:
en la gráfica se observa el desarrollo hacia el lado derecho,
debiendo realizar otras pruebas para verificar o mejorar el algoritmo.
# Revisar causalidad de ft(t) en sympy# considera causal h(t)=b0*d + (modos caracteristicos)*u# Parte 1 comprueba solo impulso d(t) import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
k = sym.Symbol('k',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t)# h = u
h = u.subs(t,t-1)
# h = u.subs(t,t+1)# h = 5# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTO
donde_u = fcnm.busca_escalon(h)
# Causal en RHS
causal = Falseiflen(donde_u)>0: # existe un escalon unitario
donde = donde_u[0][0]
sentido = donde_u[0][1]
if donde>=0 and sentido>0: # lado derecho del plano
causal = True# SALIDAprint('h(t):')
sym.pprint(h)
print('busca_impulso:', fcnm.busca_impulso(h))
print('busca_escalon(ft):', donde_u)
print('xcausal: ',causal)
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()
..
Ejemplo 3. Revisa si f(t) es causal un término de factores que se multiplican
Referencia: h(t) de Lathi Ej. 2.10 p181
Se amplia el algoritmo anterior para incorporar la revisión de f(t) con factores que se multiplican.
Se muestras algunos ejemplos a usar para analizar esta parte:
# entrada x(t)# h = u# h = 3*sym.pi*u.subs(t,t-1)
h = sym.exp(-2*t)*u.subs(t,t+1)
# h = sym.exp(-2*t)*sym.cos(t)# h = 5# h = u.subs(t,-t+2)*u.subs(t,-t+1)# h = u.subs(t,-t+1)*u.subs(t,-t+1)# h = 3*u.subs(t,-t+1)*u.subs(t,-t+1)# h = u.subs(t,t+1)*u.subs(t,-t+1)# h = u.subs(t,t-1)*u.subs(t,-t+2)# h = u.subs(t,t+2)*u.subs(t,-t-1)# h = u.subs(t,t-2)*u.subs(t,-t-1)
El resultado esperado del algoritmo para uno de los ejemplos presentados en la entrada del algoritmo es:
h(t):
-2*t
e *Heaviside(t + 1)
hcausal: False
El resultado indica que no es causal, que es concordante con la gráfica pues este sistema tiene respuesta antes de t=0, en el lado izquierdo del plano. Es decir el sistema se ‘adelanta’ a la llegada de la señal en el tiempo.
La revisión de factores de cada término de suma, se realiza descomponiendo la expresión en sus factores con ft.args o se obtienen los término suma con term_suma = sym.Add.make_args(ft)
En el caso que el termino sea un escalon o impulso elevado al cuadrado, se extrae solo la base del término para revisar.
if ft.is_Pow:
ft = ft.as_base_exp()[0]
La multiplicación de funciones escalón μ(t)μ(t-1) se da cuando se aplica el integral de convolución entre x(t) y h(t), siendo x(t) un escalón. Sin embargo Sympy mantiene las expresiones sin simplificar, asunto a considerar en los siguientes ejercicios. por ahora se lo realiza con fcnm.simplifica_escalon(ft).
# Revisar causalidad de ft(t) en sympy# considera causal h(t)=b0*d + (modos caracteristicos)*u# Parte 1 comprueba solo impulso d(t) import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t)# h = u# h = 3*sym.pi*u.subs(t,t-1)
h = sym.exp(-2*t)*u.subs(t,t+1)
# h = sym.exp(-2*t)*sym.cos(t)# h = 5# h = u.subs(t,-t+2)*u.subs(t,-t+1)# h = u.subs(t,-t+1)*u.subs(t,-t+1)# h = 3*u.subs(t,-t+1)*u.subs(t,-t+1)# h = u.subs(t,t+1)*u.subs(t,-t+1)# h = u.subs(t,t-1)*u.subs(t,-t+2)# h = u.subs(t,t+2)*u.subs(t,-t-1)# h = u.subs(t,t-2)*u.subs(t,-t-1)# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTOdefes_causal(ft):
''' h(t) es causal si tiene
b0*d(t)+u(t)*(modos caracterisicos)
'''defes_causal_impulso(ft):
''' un termino en lado derecho del plano
'''
causal = False
donde_d = fcnm.busca_impulso(ft)
iflen(donde_d)>0: # tiene impulsoif donde_d[0]>=0: # derecha del plano
causal = Truereturn(causal)
causal = True
term_suma = sym.Add.make_args(ft)
for term_k in term_suma:
if term_k.has(sym.Heaviside):
term_k = fcnm.simplifica_escalon(term_k) # h(t-a)*h(t-b)
causal_k = True
donde_u = fcnm.busca_escalon(term_k)
iflen(donde_u)==0: # sin escalon?
causal_k = Falseiflen(donde_u)==1:
sentido1 = donde_u[0][0]
donde1 = donde_u[0][1]
# en lado izquierdo del planoif donde1<0 or sentido1<0:
causal_k = Falseiflen(donde_u)==2:
donde1 = donde_u[0][0]
sentido1 = donde_u[0][1]
donde2 = donde_u[1][0]
sentido2 = donde_u[1][1]
# rectangulo lado derecho del planoif (donde1<donde2):
if donde1<0: # u(t+1)*u(-t+1)
causal_k = Falseif (donde2<donde1):
if donde2<0: # u(-t+1)*u(t+1)
causal_k = Falseelse: # un termino, sin escalon unitario# causal depende si tiene un impulso
causal_k = es_causal_impulso(term_k)
causal = causal and causal_k
return(causal)
hcausal = es_causal(h)
# SALIDAprint('h(t):')
sym.pprint(h)
print('busca_impulso(ft):', fcnm.busca_impulso(h))
print('busca_escalon(ft):', fcnm.busca_escalon(h))
print('xcausal: ',hcausal)
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()
..
Ejemplo 4. Revisa si f(t) es causal de varios término de suma
Finalmente la revisión de la expresión f(t) revisa cada uno de los terminos de una suma. Dentro de cada termino de la suma se revisa cada factor que multiplica. Dentro de cada factor que multiplica de analiza si se encuentra un escalón o impulso unitario que opere solo en el lado derecho del plano.
Las expresiones usadas para revisar el ejercicio son
# entrada x(t)# h = u + d# h = 3*sym.pi*u.subs(t,t-1)+ 2*d# h = sym.exp(-2*t)*u.subs(t,t+1) + sym.exp(-4*t)*u.subs(t,t-3)# h = sym.exp(-2*t)*sym.cos(t) + 4# h = 5# h = u.subs(t,-t+2)*u.subs(t,-t+1)# h = u.subs(t,-t+1)*u.subs(t,-t+1)# h = 3*u.subs(t,-t+1)*u.subs(t,-t+1)
h = u.subs(t,t-1)-u.subs(t,t-2)
el resultado del algoritmo para el ejercicio corresponde a:
La gráfica correspondiente que permite comprobar visualmente el resultado es
Instrucciones en Python
El algoritmo básicamente realiza el seguimiento de los resultados que se obtienen con las funciones de las secciones anteriore para cada término de suma.
# Revisar causalidad de ft(t) en sympy# considera causal h(t)=b0*d + (modos caracteristicos)*u# Parte 1 comprueba solo impulso d(t) import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# entrada x(t)# h = u + d# h = 3*sym.pi*u.subs(t,t-1)+ 2*d# h = sym.exp(-2*t)*u.subs(t,t+1) + sym.exp(-4*t)*u.subs(t,t-3)# h = sym.exp(-2*t)*sym.cos(t) + 4# h = 5
h = u.subs(t,t-1)-u.subs(t,t-2)
# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51
# PROCEDIMIENTO
hcausal = fcnm.es_causal(h)
# SALIDAprint('h(t):')
sym.pprint(h)
print('busca_impulso(ft):', fcnm.busca_impulso(h))
print('busca_escalon(ft):', fcnm.busca_escalon(h))
print('xcausal: ',hcausal)
# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()
El estado cero del sistema, «Zero-State», supone no hay energía almacenada, que los capacitores están descargados, que recien sale el equipo de la caja. Para éste caso, la respuesta del sistema se conoce como respuesta a estado cero, «Zero-State Response» ZSR.
Respuesta total
=
respuesta a entrada cero
+
respuesta a estado cero
ZSR = h(t) ⊗ x(t)
Para los problemas presentados se asume que el sistema es lineal, causal e invariante en el tiempo. En la práctica, muchos de los sistemas son causales, pues su respuesta no inicia antes de aplicar una entrada, es decir, todas las entradas a evaluar empiezan en t=0.
La respuesta del sistema y(t) para un LTIC se determina con la convolución entre x(t) y h(t), definida mediante el operador ⊗ como :
El límite inferior del integral se usa como 0–, implica aunque se escriba solo 0 se pretende evitar la dificultad cuando x(t) tiene un impulso en el origen.
Recordar que, basado en la condición de causalidad, si x(t)=0 para t<0, esta señal es causal. En el caso contrario, si x(t)=0 para t>0 la señal es No causal o anticausal.
Para éste ejercicio se usa la tabla de integrales convolución, línea 4, así el enfoque del desarrollo se mantiene sobre la forma de la señal resultante. El siguiente ejemplo se desarrolla con el integral de convolución.
De las gráficas se observa que la entrada es semejante a conectar en la entrada un capacitor con carga, que la pierde en el tiempo.
En la salida se observa el efecto, la parte inicial corresponde a la corriente en el circuito mientras el capacitor de la entrada entrega energía al sistema. Note que en el sistema o circuito se debe ir cargando el capacitor del sistema. Luego, un poco más del segundo 1, la corriente invierte el sentido volviéndose negativa por la carga almacenada en el capacitor del sistema.
La explicación breve realizada debería ser comprobada en los experimentos de laboratorio, preferiblemente a escala menor con componentes tipo electrónico.
Para usar el simulador con diagrama de bloques para respuesta a impulso unitario, es necesario usar la expresión:
h(t) = b_0 \delta (t) +P(D)[y_n(t) \mu(t)]
Donde si el orden de P(D) es menor al orden de Q(D) entonces b0=0.
Para el ejercicio desarrollado en el sistema modelo, P(D)=D, por lo que hay una sola derivada, que se puede obtener de la entrada del primer integrador de y(t).
( D^2 + 3D +2 ) y(t) = Dx(t)
Para obtener h(t), el impulso unitario x(t), es semejante a tener una condicion inicial sin señal de entrada, siguiendo la teoría descrita de sobre condiciones iniciales.
El impulso unitario representa que y'(t) inicia con 1, por ser el término de más alto orden en P(D) y se configura el diagrama de bloques:
Con lo que la respuesta del sistema al impulso unitario que se obtiene es:
El parámetros inicial se configuró en:
Como referencia para observación se muestra la salida y(t):
2. Respuesta a impulso con sistema LTI definido en Scipy.signal
La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación P, Q y las condiciones de inicio. La respuesta a un impulso unitario δ(t) o función de transferencia h(t) se calcula con una versión muestreada xi a partir del tiempo muestreado en el intervalo [a,b].
sistema = senal.lti(P,Q)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)
– scipy.signal.impulse() donde solo se da el vector ti
– aplicar la teoria de los coeficientes donde solo la derivada de mayor orden de P(D) evaluada en cero, tiene el valor de uno. Todos los demás términos tienen el valor de cero. Revisar la parte conceptual de la sección para respuesta a impulso.
Para el ejercicio y como comprobación se calculan las dos formas, con entrada cero y condición inicial [1,0] y con la función de scipy, que son las usadas en la gráfica.
Instrucciones en Python
# Sistema lineal usando scipy.signal# Q(D)y(t)=P(D)x(t)# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)import numpy as np
import scipy.signal as senal
import matplotlib.pyplot as plt
# INGRESO# Señal de entrada x(t)
x = lambda t: t*0
# Sistema LTI# coeficients Q P de la ecuación diferencial
Q = [1., 3., 2.]
P = [1., 0.]
# condiciones iniciales [y'(0),y(0)]
cond_inicio = [1,0]
# grafica, parametros
a = 0 ; b = 5 # intervalo observacion [a,b]
muestras = 101
# PROCEDIMIENTO# intervalo observado
ti = np.linspace(a, b, muestras)
xi = x(ti)
sistema = senal.lti(P,Q)
# respuesta entrada cero
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, cond_inicio)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)
# SALIDA - GRAFICA
plt.subplot(211)
plt.suptitle('Respuesta a Impulso Unitario: h(t) = dy/dt')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')
plt.legend()
plt.grid()
plt.subplot(212)
# plt.plot(t_y, yi, color='magenta', label='y(t)')# plt.plot(t_y, yc[:,0], 'b-.', label='h(t)')
plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')
plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()
3. Algoritmo de Runge-Kutta d2y/dx2 de Análisis numérico
Para el caso del algoritmo de Runge-Kutta, se aplica las condiciones iniciales de y(0)= 0 y y'(0)=1,
siguiendo lo descrito en la parte teórica y se toman los valores de z=y’ para obtener el resultado.
Las instrucciones aplicadas al algoritmo de entrada cero se muestran a continuación:
Instrucciones en Python
# Respuesta a entrada cero# solucion para (D^2+ D + 1)y = 0import numpy as np
import matplotlib.pyplot as plt
defrungekutta2_fg(f,g,x0,y0,z0,h,muestras):
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,3),dtype=float)
# incluye el punto [x0,y0]
estimado[0] = [x0,y0,z0]
xi = x0
yi = y0
zi = z0
for i inrange(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
estimado[i] = [xi,yi,zi]
return(estimado)
# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y
t0 = 0
y0 = 0
z0 = 1
h = 0.1
tn = 5
muestras = int((tn-t0)/h)
tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)
# GRAFICA#plt.plot(ti,yi, label='y(t)')
plt.plot(ti,zi, color = 'Red', label='h(t)=dy/dt')
plt.ylabel('dy/dt')
plt.xlabel('t')
plt.title('Respuesta impulso: Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()
Recuerde que se aplican las modificaciones de acuerdo al planteamiento del problema, por lo que revise los conceptos antes de aplicar el algoritmo.
La respuesta al impulso reutiliza el algoritmo de para encontrar la solución homogénea de la ecuación diferencial lineal. Se reutiliza la función respuesta_ZIR() de la sección anterior.
Siguiendo el método simplificado al emparejar impulsos, las condiciones iniciales, dado que el orden de las derivadas de la izquierda es 2, se establece como y'(0)=1, y(0)=0.
N = 1
yn(0) = 1
N = 2
yn(0) = 0, y’n(0) = 1
N = 3
yn(0) = 0, y’n(0) = 0, y"n(0) = 1
N = 4
…
Siendo N el grado mayor de las derivadas de y(t) o lado izquierdo LHS y M el grado de mayor orden de las derivadas para x(t) o lado derecho RHS.
Si N>M, se tiene que b0=0. Si N=M el valor de b0 es el coeficiente de la derivada de mayor grado para x(t).
Desarrollo del algoritmo en Python
Se empieza buscando el orden N, para y(t) de las derivadas de la ecuación diferencial lineal. Con N se crea un vector ceros para condiciones de inicio y se escribe el valor de 1 a la primera casilla qu representa la posición de mayor orden de derivada .
# Método simplificado al emparejar impulsos
N = sym.ode_order(ecuacion.lhs,y)
M = sym.ode_order(ecuacion.rhs,x)
# Condiciones iniciales para respuesta a impulso
cond_inicio = [0]*N # lista de ceros tamano N
cond_inicio[0] = 1 # condicion de mayor orden
Se debe buscar el coeficiente b0, que es el del término de mayor orden de la derivada para el lado derecho de la ecuación.
En la ecuacion parte derecha eq_RHS, se busca cada término hasta encontrar el de orden M. Se extrae el coeficiente del término encontrado. es decir todas las partes que no contienen el término de la derivada, ejemplo 3π.
# coeficiente de derivada de x(t) de mayor orden
b0 = sym.nan
if N>M: # orden de derivada diferente
b0 = 0
if N==M: # busca coeficiente de orden mayor
eq_RHS = sym.expand(ecuacion.rhs)
term_suma = sym.Add.make_args(eq_RHS)
for term_k in term_suma:
# coeficiente derivada mayorif (M == sym.ode_order(term_k,x)):
b0 = 1 # para separar coeficiente
factor_mul = sym.Mul.make_args(term_k)
for factor_k in factor_mul:
ifnot(factor_k.has(sym.Derivative)):
b0 = b0*factor_k
Con el valor de b0 y las cond_inicio se usa el algoritmo respuesta_ZIR() realizado para encontrar la respuesta a entrada cero a partir de la ecuación homogenea.
A partir de la solución homogénea, se crea la función h(t) aplicando la expresión:
h(t)=b_0 \delta (t)+ [P(D)y_n (t)] \mu (t)
y dado que para el ejercicio N>M, el orden de las derivadas de la izquierda es mayor que el orden de las derivadas de la derecha, se tiene que b0=0
h(t)=0 \delta (t) + [D y_n (t)] \mu (t)
# ecuacion homogenea x(t)=0, entrada cero y# condiciones de impulso unitario
sol_ht = fcnm.respuesta_ZIR(ecuacion,cond_inicio)
# Respuesta a impulso h(t)
P_y = ecuacion.rhs.subs(x(t),sol_ht['ZIR']).doit()
h = P_y*u + b0*sym.DiracDelta(t)
# h = sym.expand(h)
con lo que se llega a la respuesta de h(t),
h(t)= (-e^{-t} + 2e^{-2t})\mu (t)
La respuesta del algoritmo en Python es:
clasifica EDO:
factorable
nth_linear_constant_coeff_variation_of_parameters
nth_linear_constant_coeff_variation_of_parameters_Integral
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
1 = -C1 - 2*C2
constante : {C1: 1, C2: -1}
ZIR :
-t -2*t
e - e
h :
/ -t -2*t\
\- e + 2*e /*Heaviside(t)
>>>
y con gráfica:
Instrucciones en Python
# Respuesta a impulso h(t) con Sympy-Python# https://blog.espol.edu.ec/telg1001/lti-ct-respuesta-al-impulso-con-sympy-python/# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dximport numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t', real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# 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)
ecuacion = sym.Eq(LHS,RHS)
# cond_inicial con Método simplificado# de emparejar impulsos# Grafica: intervalo tiempo [t_a,t_b]
t_a = 0 ; t_b = 5
muestras = 51
# PROCEDIMIENTO# Método simplificado al emparejar términos
N = sym.ode_order(ecuacion,y)
M = sym.ode_order(ecuacion,x)
# coeficiente de derivada de x(t) de mayor orden
b0 = sym.nan
if N>M: # orden de derivada diferente
b0 = 0
if N==M: # busca coeficiente de orden mayor
eq_RHS = sym.expand(ecuacion.rhs)
term_suma = sym.Add.make_args(eq_RHS)
for term_k in term_suma:
# coeficiente derivada mayorif (M == sym.ode_order(term_k,x)):
b0 = 1 # para separar coeficiente
factor_mul = sym.Mul.make_args(term_k)
for factor_k in factor_mul:
ifnot(factor_k.has(sym.Derivative)):
b0 = b0*factor_k
# Condiciones iniciales para respuesta a impulso
cond_inicio = [0]*N # lista de ceros tamano N
cond_inicio[0] = 1 # condicion de mayor orden# ecuacion homogenea x(t)=0, entrada cero y# condiciones de impulso unitario
sol_yc = fcnm.respuesta_ZIR(ecuacion,cond_inicio)
yc = sol_yc['ZIR']
# Respuesta a impulso h(t)
P_y = ecuacion.rhs.subs(x(t),yc).doit()
h = P_y*u + b0*d
sol_yc['h'] = h
edo_tipo = sym.classify_ode(ecuacion, y(t))
# SALIDAprint('clasifica EDO:')
for elemento in edo_tipo:
print(' ',elemento)
fcnm.print_resultado_dict(sol_yc)
# GRAFICA ------------------# Para graficar la Salida
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,'h')
plt.show()
Ejemplo 2. Ecuación diferencial lineal con Orden N=M
Referencia: Lathi. Ejercicio 2.4.a p167
Determine la respuesta al impulso de un sistema LTI C descrito por la siguiente ecuación diferencial ordinaria:
(D+2)y(t) = (3D+5) x(t)
El orden N=M=1, por lo que aplican las condiciones iniciales de t0=0, y(0)=1, y'(0)=0 para resolver usando el algoritmo de respuesta a entrada cero para obtener y(t) y aplicar luego la expresión:
h(t)=b_0 \delta (t)+ [P(D)y_n (t)] \mu (t)
siendo b0=3, que es el coeficiente de la derivada de mayor grado para el lado derecho de la expresión.
clasifica EDO:
1st_linear
almost_linear
nth_linear_constant_coeff_variation_of_parameters
1st_linear_Integral
almost_linear_Integral
nth_linear_constant_coeff_variation_of_parameters_Integral
homogenea :
d
2*y(t) + --(y(t)) = 0
dt
general :
-2*t
y(t) = C1*e
eq_condicion :
1 = C1
constante : {C1: 1}
ZIR :
-2*t
e
N : 1
M : 1
b0 : 3
h :
-2*t
3*DiracDelta(t) - e *Heaviside(t)
>>>
Instrucciones en Python
El ejercicio se desarrolla creando la función edo_resp_impulso() para ser incluida en telg1001.py y así simplificar el algoritmo para el próximo ejercicio.
# Respuesta a impulso h(t) con Sympy-Python# https://blog.espol.edu.ec/telg1001/lti-ct-respuesta-al-impulso-con-sympy-python/# Lathi 2.1.a pdf 155, (D+2)y = (3D+5)ximport numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t', real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,1) + 2*y(t)
RHS = 3*sym.diff(x(t),t,1,evaluate=False)+5*x(t)
ecuacion = sym.Eq(LHS,RHS)
# cond_inicial con Método simplificado# de emparejar impulsos# Grafica: intervalo tiempo [t_a,t_b]
t_a = 0 ; t_b = 5
muestras = 51
# PROCEDIMIENTOdefrespuesta_impulso_h(ecuacion,t0=0,
y = sym.Function('y'),
x = sym.Function('x')):
''' respuesta a impulso h(t) de un
sistema con Ecuacion Diferencial lineal
'''# Método simplificado al emparejar términos
N = sym.ode_order(ecuacion,y)
M = sym.ode_order(ecuacion,x)
# coeficiente de derivada de x(t) de mayor orden
b0 = sym.nan
if N>M: # orden de derivada diferente
b0 = 0
if N==M: # busca coeficiente de orden mayor
eq_RHS = sym.expand(ecuacion.rhs)
term_suma = sym.Add.make_args(eq_RHS)
for term_k in term_suma:
# coeficiente derivada mayorif (M == sym.ode_order(term_k,x)):
b0 = 1 # para separar coeficiente
factor_mul = sym.Mul.make_args(term_k)
for factor_k in factor_mul:
ifnot(factor_k.has(sym.Derivative)):
b0 = b0*factor_k
# Condiciones iniciales para respuesta a impulso
cond_inicio = [0]*N # lista de ceros tamano N
cond_inicio[0] = 1 # condicion de mayor orden# ecuacion homogenea x(t)=0, entrada cero y# condiciones de impulso unitario
sol_yc = fcnm.respuesta_ZIR(ecuacion,cond_inicio)
yc = sol_yc['ZIR']
# Respuesta a impulso h(t)
P_y = ecuacion.rhs.subs(x(t),yc).doit()
h = P_y*u + b0*d
sol_yc['N'] = N
sol_yc['M'] = M
sol_yc['cond_inicio'] = cond_inicio
sol_yc['b0'] = b0
sol_yc['h'] = h
return(sol_yc)
edo_tipo = sym.classify_ode(ecuacion, y(t))
# Respuesta a impulso h(t)
sol_h = respuesta_impulso_h(ecuacion)
h = sol_h['h']
# SALIDAprint('clasifica EDO:')
for elemento in edo_tipo:
if'linear'in elemento.split('_'):
print(' ',elemento)
fcnm.print_resultado_dict(sol_h)
# GRAFICA ------------------# Para graficar la Salida
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,'h')
plt.show()
Ejemplo 3. Ecuación diferencial lineal con Orden de N>M
Referencia: Lathi. Ejercicio 2.4.b p167
Determine la respuesta al impulso de un sistema LTI C descrito por la siguiente ecuación:
D(D+2)y(t) = (D+4) x(t)
El orden N>M, por lo que aplican las condiciones iniciales de t0=0, y(0)=0, y'(0)=1 para resolver usando el algoritmo de respuesta a entrada cero para obtener y(t) y aplicar luego la expresión:
h(t)=b_0 \delta (t)+ [P(D)y_n (t)] \mu (t)
siendo b0=0, que es el coeficiente de la derivada de mayor grado para el lado derecho de la expresión.
clasifica EDO:
nth_linear_constant_coeff_variation_of_parameters
nth_linear_constant_coeff_variation_of_parameters_Integral
homogenea :
2
d d
2*--(y(t)) + ---(y(t)) = 0
dt 2
dt
general :
-2*t
y(t) = C1 + C2*e
eq_condicion :
0 = C1 + C2
1 = -2*C2
constante : {C1: 1/2, C2: -1/2}
ZIR :
-2*t
1 e
- - -----
2 2
N : 2
M : 1
cond_inicio : [ 1 0 ]
b0 : 0
h :
/ -2*t\
\2 - e /*Heaviside(t)
>>>
Instrucciones en Python
# Respuesta a impulso h(t) con Sympy-Python# https://blog.espol.edu.ec/telg1001/lti-ct-respuesta-al-impulso-con-sympy-python/# Lathi. Ejercicio 2.4.b p167 D(D+2)y(t) = (D+4)x(t)import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import telg1001 as fcnm
# INGRESO
t = sym.Symbol('t', real=True)
y = sym.Function('y')
x = sym.Function('x')
h = sym.Function('h')
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 2*sym.diff(y(t),t,1)
RHS = sym.diff(x(t),t,1,evaluate=False)+4*x(t)
ecuacion = sym.Eq(LHS,RHS)
# cond_inicial con Método simplificado de emparejar impulsos# Grafica: intervalo tiempo [t_a,t_b]
t_a = 0 ; t_b = 5
muestras = 51
# PROCEDIMIENTO
edo_tipo = sym.classify_ode(ecuacion, y(t))
# Respuesta a impulso h(t)
sol_h = fcnm.respuesta_impulso_h(ecuacion)
h = sol_h['h']
# SALIDAprint('clasifica EDO:')
for elemento in edo_tipo:
if'linear'in elemento.split('_'):
print(' ',elemento)
fcnm.print_resultado_dict(sol_h)
# GRAFICA ------------------
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,'h')
plt.show()
