3.7 Simplifica multiplicación entre impulso δ(t) o escalón unitario μ(t) con Sympy

La multiplicación entre impulso δ(t) o escalón unitario μ(t) aparece al desarrollar en integral de convolución cuando ambas señales son de tipo causal.

x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Las expresiones resultantes pueden incluir algunas de éstas operaciones como se muestra en el ejemplo 3 de una convolución entre señales tipo rectangular y rampa de caracter causal.

x(t) = \mu (t) - \mu (t-1) h(t) = t\mu (t) - t\mu (t-2)

Con respuestas al usar el algoritmo con Sympy que integran varias multiplicaciones de escalón unitario.

Como Sympy ‘1,11,1’ revisada en Enero del 2023, aún no incluye este tipo de operaciones, se realizan dos funciones:  simplifica_impulso() y simplifica_escalon(). Se desarrolla como un ejercicio de análisis de expresiones con Sympy.

>>> sym.__version__
'1.11.1'

Simplificar multiplicación de impulso unitario δ(t)

Las operaciones encontradas en los ejercicios del curso dan una idea básica de por dónde empezar. También se probó las instrucciones doit(), evalf() sin cambios.

>>> sym.DiracDelta(t)*sym.DiracDelta(t-1)
DiracDelta(t)*DiracDelta(t - 1)
>>> sym.simplify(sym.DiracDelta(t)*sym.DiracDelta(t-1))
DiracDelta(t)*DiracDelta(t - 1)
>>> sym.DiracDelta(t)**2
DiracDelta(t)**2
>>> sym.simplify(sym.DiracDelta(t)**2)
DiracDelta(t)**2
>>>

Como punto de partida de plantea encontrar todos los puntos de la expresión donde intervienen los impulsos. La función busca_impulso() permite revisar si existiendo impulsos en diferentes tiempos, el resultado deberá ser cero.

x = d
h = d.subs(t,t-1)

xh = x*h

Para facilitar el análisis se realizan las gráficas de las dos funcines que se multiplican y la operación resultante del algoritmo:

simplifica impulso Sympy 01

otro caso a considerar es:

x = sym.sin(t)
h = d.subs(t,t-2)

xh = x*h

simplifica impulso Sympy 02

Instrucciones con Python

El algoritmo se desarrolla como funciones, para ser incorporadas a telg1001.py

# Simplificar multiplicacion de impulso unitario d(t)*d(t-1)
# http://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/
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 = d
# h = d.subs(t,t-1)
# h = d # h= d**2

# x = 2
# h = d
# h = sym.pi*d.subs(t,t-2)

x = sym.sin(t)
h = d.subs(t,t-2)
# h = 5*d.subs(t,t-2)**2

# x = d.subs(t,t-1)
# h = u.subs(t,t-1)
# h = d*u.subs(t,t-2)+ u*u.subs(t,t-2)
# h = d*d.subs(t,t-1)*u.subs(t,t+2)+3

# grafica intervalo [t_a,t_b] plano simétrico
t_b = 4 ; t_a = -t_b
muestras = 81

# PROCEDIMIENTO
def simplifica_impulso(ft):
    ''' simplifica d**2, d(t-1)*d
    '''
    def simplifica_d_d(ft):
        '''un termino de suma  d**2, d(t+1)*d,
        '''
        respuesta = ft
        if ft.has(sym.DiracDelta):# tiene impulsos
            impulso_en = fcnm.busca_impulso(ft)
            if len(impulso_en)==0: # anulado por d(t-a)*d(t-b)
               respuesta = 0*t
            elif len(impulso_en)>0: # tiene impulsos
                respuesta = 1
                factor_mul = sym.Mul.make_args(ft)
                for factor_k in factor_mul:
                    if not(factor_k.has(sym.DiracDelta)):
                        if not(factor_k.has(sym.Heaviside)):
                            termino = factor_k.subs(t,impulso_en[0])
                        else: # tiene escalón
                            termino = factor_k.subs(t,impulso_en[0])
                            if termino == 1/2: #coinciden d,u
                                termino = 1
                        respuesta = respuesta*termino
                    else:  # factor con impulso
                        if factor_k.is_Pow: # tiene exponente
                            respuesta = respuesta*factor_k.args[0]
                        else: # termino sin exponente
                            respuesta = respuesta*factor_k
        return(respuesta)
    
    # revisa terminos suma
    respuesta = 0*t
    ft = sym.expand(ft,t)
    term_suma = sym.Add.make_args(ft)
    for term_k in term_suma:
        respuesta = respuesta + simplifica_d_d(term_k)
        
    return(respuesta)

xh = x*h
ft = simplifica_impulso(xh)

# SALIDA
print('xh inicial:')
sym.pprint(xh)
print('h simplificado:')
sym.pprint(ft)

# grafica
figura = fcnm.graficar_xh_y(x,h,ft,t_a,t_b,muestras,y_nombre='xh')
plt.show()

 


Simplificar multiplicación de escalón unitario μ(t)

Las operaciones encontradas en los ejercicios del curso dan una idea básica de por dónde empezar

>>> sym.Heaviside(t)*sym.Heaviside(t-1)
Heaviside(t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(t)*sym.Heaviside(t-1))
Heaviside(t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(-t+1)*sym.Heaviside(t-1))
Heaviside(1 - t)*Heaviside(t - 1)
>>> sym.simplify(sym.Heaviside(t)**2)
Heaviside(t)**2
>>> sym.simplify(sym.DiracDelta(t-1)*sym.Heaviside(t))
DiracDelta(t - 1)*Heaviside(t)
>>>

Como punto de partida de plantea encontrar todos los puntos de la expresión donde intervienen el escalón unitario. La función busca_escalon() permite revisar la ubicación y el sentido de desarrollo de cada escalón unitario, con estos datos se puede determinar si la función se anula con la otra o la parte que permanece.

x = u
h = u.subs(t,t-1)

xh = x*h

Con una grafica se puede comprobar el escalón que se superpone en la multiplicación de x(t)h(t).

simplifica escalon Sympy 01

otro ejercicio a considerar es cuando los escalones unitarios se pueden complementar o anular. En el siguiente ejemplo x(t) h(t) se anulan entre si.

x = u.subs(t,t-1)
h = u.subs(t,-t-1)

xh = x*h

simplifica escalon Sympy 02

mientras en en otro caso, las señales tienen una región donde pueden coexistir

x = u.subs(t,-t+1)
h = u.subs(t,t-0)

Instrucciones con Python

# Simplificar multiplicacion de escalon unitario u(t)*u(t-1)
# http://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/
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), respuesta impulso h(t)
# x = u
# h = u.subs(t,t-1)

# x = u.subs(t,t-1)
# h = u.subs(t,-t-1)

x = u.subs(t,-t+1)
h = u.subs(t,t-0)

#x = 3*u.subs(t,-t+2)
#h = 2*u.subs(t,-t+3)

# x = u.subs(t,t+1)*u.subs(t,-t+2)
# h = 1
# h = d.subs(t,-t+2)#*u.subs(t,t-3)


# grafica intervalo [t_a,t_b] plano simétrico
t_b = 4 ; t_a = -t_b
muestras = 301

# PROCEDIMIENTO
def simplifica_escalon(ft):
    ''' simplifica multiplicaciones
        Heaviside(t-a)*Heaviside(t-b) en f(t)
    '''
    def simplifica_u_u(ft):
        '''solo dos pares de u(t-1)*u(t-2),
           sin coeficientes
        '''
        donde_u = fcnm.busca_escalon(ft)
        donde_u = np.array(donde_u)
        # direccion donde_[:,1],
        # lugar donde[:,0]
        # analiza multiplicación
        resultado = ft   # todo igual
        if donde_u[0,1]*donde_u[1,1] > 0: # direccion igual
            k = 0
            if donde_u[0,1]>0: # hacia derecha
                k = np.argmax(donde_u[:,0])
                k_signo = 1
            else: # hacia izquierda
                k = np.argmin(donde_u[:,0])
                k_signo = -1
            ubica = donde_u[k,1]*t-k_signo*donde_u[k][0]
            resultado = sym.Heaviside(ubica)
        else: # direccion diferente
            if donde_u[0][1]>0 and (donde_u[0,0]>donde_u[1,0]):
                    resultado = 0
            if donde_u[0][1]<0 and (donde_u[0,0]<=donde_u[1,0]):
                    resultado = 0
        return(resultado)

    def simplifica_u_term(ft):
        ''' simplifica un termino de varios
            factores que multiplican 2*pi*u*u(t-1)
        '''
        respuesta = ft
        if ft.has(sym.Heaviside): # tiene escalon
            escalon_en = fcnm.busca_escalon(ft)
            revisa = 1 ; otros = 1 ; cuenta = 0
            factor_mul = sym.Mul.make_args(ft)
            for factor_k in factor_mul:
                if factor_k.has(sym.Heaviside):
                    if factor_k.is_Pow: # con exponente
                        revisa = revisa*factor_k.args[0]
                        cuenta = cuenta + 1
                    else: # sin exponente
                        revisa = revisa*factor_k
                        cuenta = cuenta + 1
                    if cuenta>1: # simplificar
                        revisa = simplifica_u_u(revisa)
                        cuenta = len(fcnm.busca_escalon(revisa))
                else: # factor sin Heaviside
                    otros = otros*factor_k
            respuesta = otros*revisa
        return(respuesta)

    # revisa terminos suma
    respuesta = 0*t
    ft = sym.expand(ft,t)
    if ft.has(sym.DiracDelta): # tiene impulsos
        ft = fcnm.simplifica_impulso(ft)
    term_suma = sym.Add.make_args(ft)
    for term_k in term_suma:
        respuesta = respuesta + simplifica_u_term(term_k)

    return(respuesta)

xh = x*h
ft = simplifica_escalon(xh)

# SALIDA
print('busca_impulso(ft)', fcnm.busca_impulso(h))
print('busca_escalon(ft)', fcnm.busca_escalon(h))
print('x*h inicial:')
sym.pprint(xh)
print('simplificado:')
sym.pprint(ft)

# grafica
figura = fcnm.graficar_xh_y(x,h,ft,t_a,t_b,y_nombre='xh')
plt.show()

..

3.6 Busca impulso unitario δ(t) o escalón unitario μ(t) con Sympy-Python

En los sistemas lineales, la causalidad se muestra con respuestas para t≥0, por lo que se usan términos multiplicados por escalón unitario μ(t). Mientras que la respuesta al impulso genera al menos un término con impulso unitario δ(t) con un coeficiente constante.

Para el desarrollo de los ejercicios, es de gran utilidad disponer de una función que permita buscar el impulso unitario dentro de una expresión de varios términos. Listar las ubicaciones de los impulsos es un resultado que se usa en la generación de gráficas, pues no siempre la muestra en el intervalo de tiempo usado se encuentra exactamente sobre la aplicación del impulso unitario.

Se propone usar dos funciones. una para encontrar dónde ocurre el impulso y otra para buscar donde y el sentido del escalón unitario.


1. Busca impulso unitario, DiracDelta(t) o δ(t) con Sympy

Si f(t) fuese solo δ(t), δ(t-1) o δ(t+1), el valor de desplazamiento de obtiene con el primer argumento de la función mediante ft.args[0] que al evaluarlo en cero y multiplicarlo por (-1) se tendría donde se aplica.

ecuacion = sym.Eq(ft.args[0],0)
donde = sym.solve(ecuacion,t)

Para los casos como δ(-t-1), δ(-t+1) es necesario trabajar con el signo de la variable t. Una forma de evaluar la expresión interiores de los impulsos es obtener la ecuacion de la forma at+b=0 con sym.Eq(ft.args[0],0) y luego despejar el valor de t con la instrucción sym.solve(ecuacion,t).

En expresiones simples de un solo término con coeficientes, se tiene:

h = 2*sym.pi*d.subs(t,t-2)

con resultados como:

h(t):
2*pi*DiracDelta(t - 2)
busca_impulso(h):  [2]

busca impulso sympy 01

En el caso que se usen términos de varios impulsos, el valor de interés se encuentra usando el mínimo de la lista donde y observando que aún sea ≥0.

h = d.subs(t,t-1) + 2*sym.pi*d.subs(t,t-2)

con resuldado esperado de algoritmo

h(t):
2*pi*DiracDelta(t - 2) + DiracDelta(t - 1)
busca_impulso(h):  [1, 2]

busca impulso Sympy 02

Instrucciones en Python

# Busca impulso unitario en h(t)
# http://blog.espol.edu.ec/telg1001/busca-impulso-unitario-o-escalon-unitario-con-sympy/
import sympy as sym
import matplotlib.pyplot as plt
import telg1001 as fcnm

# INGRESO
t = sym.Symbol('t',real=True)
u = sym.Heaviside(t)
d = sym.DiracDelta(t)
a = sym.Symbol('a',real=True)

# busca impulso en h(t)
# h = d
# h = 2*sym.pi*d.subs(t,t-2)
h = d.subs(t,t-1) + 2*sym.pi*d.subs(t,t-2)
# 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
def busca_impulso(ft):
    ''' busca en f(t) impulsos sym.DiracDelta
        entrega una lista ordenada de resultados
    '''
    def impulso_donde(ft):
        ''' busca posicion de impulso en ft simple d, d**2
            un solo termino,sin factores
        '''
        ft = sym.sympify(ft,t) # convierte a sympy si es constante
        donde = [] # revisar f(t) sin impulsos
        if ft.has(sym.DiracDelta):
            if not(ft.is_Pow): # sin exponentes
                ecuacion = sym.Eq(ft.args[0],0)
            else: # con exponentes d**2
                ecuacion = sym.Eq(ft.args[0].args[0],0)
            donde = sym.solve(ecuacion,t)
        return(donde)
    
    def impulso_donde_unterm(ft):
        ''' revisa impulso en un termino suma de ft
        '''
        donde = [] # revisar f(t) sin impulsos
        factor_mul = sym.Mul.make_args(ft)
        for factor_k in factor_mul:
            if factor_k.has(sym.DiracDelta):
                donde_k = impulso_donde(factor_k)
                if len(donde) == 0: # vacio
                    donde.extend(donde_k)
                if len(donde)>0: # sin repetir
                    if not(donde_k[0] in donde): 
                        donde = [] # anulado por d(t-a)*d(t-b)
        return(donde)

    # revisa terminos suma
    ft = sym.sympify(ft,t) # convierte a sympy si es constante
    ft = sym.expand(ft)
    respuesta = []
    term_suma = sym.Add.make_args(ft)
    for term_k in term_suma:
        donde = impulso_donde_unterm(term_k)
        if not(donde in respuesta) and len(donde)>0:
            respuesta.extend(donde)
    if len(respuesta)>1: # ordena ascendente
        respuesta.sort()
    respuesta = list(respuesta)
    return(respuesta)

donde_d = busca_impulso(h)

# SALIDA
print('h(t):')
sym.pprint(h)
print('busca_impulso(h): ',donde_d)

# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,muestras,f_nombre='h')
plt.show()

1. Busca escalón unitario, Heaviside(t) o μ(t) con Sympy

Para pruebas se tienen algunas expresiones de h(t), recordando que se analiza expresiones de un solo término, sin sumas o mutiplicaciones.

# entrada x(t)
# h = u
h = u.subs(t,t-1)
# h = u.subs(t,t+1)
# h = 5

Semejante al ejemplo anterior, se busca dónde se produce la transición del escalón unitario y el sentido si es de subida o de bajada como el signo de ‘t’.

\frac{d}{dt}\mu (t) = \delta(t)

Con la función busca_impulso() obtiene donde se produce el impulso.

El sentido del escalón se obtiene al derivar la expresión del argumento. Si el escalón se desarrolla hacia de donde se produce, será positivo, y hacia la izquierda será negativo. Al aplicar escalon_donde(ft) permite determinar el punto de subida y bajada.

Para determinar el sentido del escalón unitario, se procede a tomar la expresion del argumento ft.args[0] y luego derivar con la  instrucción sym.diff(ft.args[0]) para obtener el signo_t. con esto se observa si el desarrollo es hacia el lado derecho del plano (RHS). Con los parametros de donde y el sentido del escalón se puede aplicar un criterio de causalidad.

    def escalon_donde(ft):
        ''' ft sin factores o terminos suma
        '''
        ft = sym.sympify(ft,t) # convierte a sympy si es constante
        respuesta = []
        if ft.has(sym.Heaviside):
            eq_u      = sym.Eq(ft.args[0],0)
            sentido   = sym.diff(eq_u.lhs,t,1)
            donde     = sym.solve(eq_u,t)[0]
            respuesta = [donde,sentido]
        return(respuesta)

Instrucciones en Python

# Busca escalon unitario en h(t)
# http://blog.espol.edu.ec/telg1001/busca-impulso-unitario-o-escalon-unitario-con-sympy/
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 h(t)
# h = u
h = u.subs(t,t-1)
# h = u.subs(t,t+1)
# h = 5
# 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 = 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)

# grafica intervalo [t_a,t_b] simetrico a 0
t_b = 4 ; t_a = -t_b
muestras = 51

# PROCEDIMIENTO
def busca_escalon(ft):
    ''' busca en f(t) el donde,sentido de escalon unitario
        en un termino simple  sin sumas, para ubicar
        lado del plano de f(t)
    '''
    
    def escalon_donde(ft):
        ''' ft sin factores o terminos suma
        '''
        ft = sym.sympify(ft,t) # convierte a sympy si es constante
        respuesta = []
        if ft.has(sym.Heaviside):
            eq_u      = sym.Eq(ft.args[0],0)
            sentido   = sym.diff(eq_u.lhs,t,1)
            donde     = sym.solve(eq_u,t)[0]
            respuesta = [donde,sentido]
        return(respuesta)

    def escalon_donde_unterm(ft):
        '''revisa termino de factores multiplica
        '''
        respuesta = []
        factor_mul = sym.Mul.make_args(ft)
        for factor_k in factor_mul:
            if factor_k.has(sym.Heaviside):
                ubicado = []
                if not(factor_k.is_Pow): # sin exponente
                    ubicado = escalon_donde(factor_k)
                else:  # con exponentes d**k
                    ubicado = escalon_donde(factor_k.args[0])
                if len(ubicado)>0:
                    respuesta.append(ubicado)
        return(respuesta)
        
    # revisa terminos suma
    ft = sym.sympify(ft,t) # convierte a sympy si es constante
    ft = sym.expand(ft)
    respuesta = []
    term_suma = sym.Add.make_args(ft)
    for term_k in term_suma:
        donde = escalon_donde_unterm(term_k)
        if not(donde in respuesta) and len(donde)>0:
            respuesta.extend(donde)
    if len(respuesta)>1: # ordena ascendente
        respuesta = np.array(respuesta)
        columna   = respuesta[:,0]
        ordena    = np.argsort(columna)
        respuesta = respuesta[ordena]

    return(respuesta)

donde_u = busca_escalon(h)

# SALIDA
print('h(t):')
sym.pprint(h)
print('busca_impulso:',fcnm.busca_impulso(h))
if len(donde_u)<=1:
    print('busca_escalon:',donde_u)
else:
    print('busca_escalon:')
    print(busca_escalon(h))

# Grafica
figura_h = fcnm.graficar_ft(h,t_a,t_b,f_nombre='h')
plt.show()

Las funciones se las incorpora a los algoritmos telg1001.py

..

 

3.5 LTI CT Respuesta del Sistema Y(s)=ZIR+ZSR con Sympy-Python

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.

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero
ZSR = h(t) ⊗ x(t)

Referencia: Lathi 1.8-1 p111. Oppenheim problema 2.61c p164 Ejemplo 9.24 p700FIEC05058_RLC

Para el ejemplo, se plantea determinar la corriente de lazo y(t) del circuito mostrado en la imagen.

\frac{dy(t)}{dt} +3 y(t) + 2\int_{-\infty}^t y(\tau)d\tau = x(t)

Como se había indicado en los desarrollos preliminares, para tener todo expresado con un solo operador, se derivan ambos lados de la ecuación:

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

En la entrada de sistema se aplica:

x(t) = 10 e^{-3t} \mu (t)

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.

Resultados con el Algoritmo en Sympy-Python

 ZIR(t):
     -t      -2*t
- 5*e   + 5*e    

 h(t):
/   -t      -2*t\             
\- e   + 2*e    /*Heaviside(t)

 ZSR(t):
     -t                    -2*t                    -3*t             
- 5*e  *Heaviside(t) + 20*e    *Heaviside(t) - 15*e    *Heaviside(t)
xcausal:  True
hcausal:  True
limites de integral: [ 0 , t ]

 y(t) = ZIR(t)+ZSR(t):
     -t                   -t       -2*t                   -2*t       -3*t     
- 5*e  *Heaviside(t) - 5*e   + 20*e    *Heaviside(t) + 5*e     - 15*e    *Heaviside(t)
>>> 

Graficas de y(t) = ZIR(t)+ZSR(t)

LTIC Ejercicio01 Y_total Sympy

LTIC Ejercicio 01 ZIR Sympy

LTIC Ejercicio 01 ht Sympy

LTIC Ejercicio 01 ZSR Sympy

LTIC YTotal ZIR Ej01 animado

Instrucciones en Python

# Respuesta total del sistema
# y(t) = ZIR(t) + ZSR(t)
# http://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 ZSR
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')
    
# SALIDA
print('\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.

3.4.3 LTI CT – Respuesta a estado cero ZSR con Sympy-Python

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

Referencia Lathi ejemplo 2.9 p175, Lathi ejemplo 2.6 p166

En la entrada de sistema se aplica:

x(t) = 10 e^{-3t} \mu (t)

El sistema tiene respuesta a impulso:

h(t) = \big( 2e^{-2t} -e^{-t}\big)\mu (t)

El ejemplo es la continuación de lo realizado en respuesta a entrada cero, que tiene la ecuación:

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

La respuesta del sistema y(t) para un LTIC se determina como:

y(t) = x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Si la entrada 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

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau , & t\ge 0\\ 0, & t<0 \end{cases}

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.

xcausal = fcnm.es_causal(x)
hcausal = fcnm.es_causal(h)

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 = False
    if hcausal==False and 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:

Resultados con Python

xh :
      -t  -2*tau                                     
- 10*e  *e      *Heaviside(tau)*Heaviside(t - tau) + 

    -2*t  -tau                                  
+ 20*e    *e    *Heaviside(tau)*Heaviside(t - tau)
xcausal : True
hcausal : True
[tau_a,tau_b] : [0, t]
intercambia : False
cond_graf : True
ZSR :
/     -t       -2*t       -3*t\             
\- 5*e   + 20*e     - 15*e    /*Heaviside(t)

respuesta Estado Cero 01 Sympy

Grafica animada de la convolución

LTIC ZSR Ej01 animado

Instrucciones en Python

# Respuesta estado cero ZSR con x(t) y h(t)
# Integral de convolucion con Sympy
# http://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 = 10*sym.exp(-3*t)*u
h = (2*sym.exp(-2*t)-sym.exp(-t))*u

# grafica intervalo [t_a,t_b] plano simétrico
t_b = 5 ; t_a = -t_b
muestras = 101

# PROCEDIMIENTO
def respuesta_ZSR(x,h):
    '''Respuesta a estado cero x(t) y h(t)
    '''
    # 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 = False
    if hcausal==False and 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)
    if not(ZSR.has(sym.Integral)):
        ZSR = fcnm.simplifica_escalon(ZSR)

    lista_escalon = ZSR.atoms(sym.Heaviside)
    ZSR = sym.expand(ZSR,t) # terminos suma
    ZSR = sym.collect(ZSR,lista_escalon)

    if intercambia == True:
        xcausal = True
        hcausal = False

    # graficar si no tiene Integral o error
    cond_graf = ZSR.has(sym.Integral)
    cond_graf = cond_graf or ZSR.has(sym.oo)
    cond_graf = cond_graf or ZSR.has(sym.nan)
    cond_graf = not(cond_graf)
    
    sol_ZSR = {'xh'      : xh,
               'xcausal' : xcausal,
               'hcausal' : hcausal,
               '[tau_a,tau_b]': [tau_a,tau_b],
               'intercambia'  : intercambia,
               'cond_graf'    : cond_graf,
               'ZSR' : ZSR,}
    return(sol_ZSR)

sol_ZSR = respuesta_ZSR(x,h)
ZSR = sol_ZSR['ZSR']

# SALIDA
fcnm.print_resultado_dict(sol_ZSR)
if sol_ZSR['hcausal']==False:
    print('revisar causalidad de h(t)')
if sol_ZSR['xcausal']==False:
    print('revisar causalidad de x(t)')
if sol_ZSR['intercambia']:
    print('considere intercambiar h(t)con x(t)')
if not(sol_ZSR['cond_graf']):
    print('revisar acortar x(t) o h(t) para BIBO')

# GRAFICA
if sol_ZSR['cond_graf']:
    fig_xh_y = fcnm.graficar_xh_y(x,h,ZSR,t_a,t_b,
                            muestras,y_nombre='ZSR')
#plt.show()

# grafica animada de convolución
n_archivo = '' # sin crear archivo gif animado 
# n_archivo = 'LTIC_ZSR_Ej01' # requiere 'imagemagick'
if sol_ZSR['cond_graf']:
    figura_animada = fcnm.graf_animada_xh_y(x,h,ZSR,t_a,t_b,
                      muestras, reprod_x = 4,y_nombre='ZSR',
                      archivo_nombre = n_archivo)
plt.show()

[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]
..


Ejemplo 2. Respuesta entrada cero ZSR entre exponencial y escalón unitario

Referencia: Oppenheim Ej 2.6 p98

Sea la y(t) la respuesta a entrada cero entre las siguientes señales:

x(t) = e^{-2t} \mu (t) h(t) = \mu (t)

Para interpretar el integral de convolución en el ejercicio aa desarrollar se adjunta la animación siguiente:

LTIC ZSR Ej02 animado

Las funciones para el algoritmo se describen como:

# entrada x(t), respuesta impulso h(t)
x = sym.exp(-2*t)*u
h = u

Teniendo como resultado:

xh :
 -2*tau                                  
e      *Heaviside(tau)*Heaviside(t - tau)
xcausal : True
hcausal : True
[tau_a,tau_b] : [0, t]
intercambia : False
cond_graf : True
ZSR :
/     -2*t\             
|1   e    |             
|- - -----|*Heaviside(t)
\2     2  /                

Gráfica de algoritmo

LTIC ZSR Ej02 Sympy

[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]
..


Ejemplo 3. Ejercicio con Filtros

Tarea:  Lathi 2.6-4 filtros p207

Para pruebas de algoritmo


[ ejemplo 1. x(t) y h(t) causal ]  [ ejemplo 2. h(t) causal ]  [ ejemplo 3. x(t) y h(t) causal]

3.4.2 Integral de Convolución entre x(t) y h(t) causal con Sympy-Python

Referencia: Lathi 2.4-1 p170, Ej 2.8 p173.

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.

x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Integral de Convolucion 01 animado

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

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau , & t\ge 0\\ 0, & t<0 \end{cases}

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.

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]

..


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)
if not(y.has(sym.Integral)):
    y = fcnm.simplifica_escalon(y)
#ZSR  = ZSR.subs(u**2,u) # Heaviside**2=Heaviside

# SALIDA
print('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 -----
def graficar_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()

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectángulo y rampa]

..


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

convolucion 02 animado

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 = False
    if hcausal==False and 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

LTI C ZSR Ej03 Sympy

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)
# http://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 = False
if hcausal==False and 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)
if not(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 and not(ZSR.has(sym.oo))
cond_graf = cond_graf and not(ZSR.has(sym.nan))

# SALIDA
print('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)')
if not(cond_graf):
    print('revisar acortar x(t) o h(t) para BIBO')

# GRAFICA
if 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()

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]

..


Ejemplo 3. Convolución entre  rectangular y rampa causales

Referencia: Oppenheim Ejemplo 2.7 p99

Realizar la convolución entre las siguientes dos señales:

x(t) = \mu (t) - \mu (t-1) h(t) = t\mu (t) - t\mu (t-2)

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.
convolucion entre rectangulo y rampa truncada
El resultado con el algoritmo del ejercicio anterior es:

xcausal:  True
hcausal:  True
limites de integral: [ 0 , t ]
x(t)*h(t-tau):
t*Heaviside(tau)*Heaviside(t - tau)
+
t*Heaviside(tau - 1)*Heaviside(t - tau - 2)
+
tau*Heaviside(tau)*Heaviside(t - tau - 2)
+
tau*Heaviside(t - tau)*Heaviside(tau - 1)
+
-t*Heaviside(tau)*Heaviside(t - tau - 2)
+
-t*Heaviside(t - tau)*Heaviside(tau - 1)
+
-tau*Heaviside(tau)*Heaviside(t - tau)
+
-tau*Heaviside(tau - 1)*Heaviside(t - tau - 2)

 ZSR(t):
   2                                  2                
  t *Heaviside(t)*Heaviside(t - 1)   t *Heaviside(t)   
- -------------------------------- + --------------- + 
                 2                          2          

   2                                      2                 2   
  t *Heaviside(t - 3)*Heaviside(t - 2)   t *Heaviside(t - 2)    
+ ------------------------------------ - -------------------- + 
                   2                              2            

+ t*Heaviside(t)*Heaviside(t - 1) - t*Heaviside(t - 3)*Heaviside(t - 2) 
                                                                                           
   Heaviside(t)*Heaviside(t - 1)   3*Heaviside(t- 3)*Heaviside(t - 2) 
- ----------------------------- - ----------------------------------- + 
                 2                                  2                

                    2
+ 2*Heaviside(t - 2) 

>>>

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 para simplificar escalón unitario multiplicados se desarrolla en Simplificar multiplicación entre impulso δ(t) o escalón unitario μ(t) con Sympy
De donde se obtiene la instrucción para el algoritmo,

ZSR = fcnm.simplifica_escalon(ZSR)

El resultado luego de aplicar la función sobre ZSR se muestra un poco mas simple:

 ZSR(t):
 2                 2                     2                     2              
t *Heaviside(t)   t *Heaviside(t - 3)   t *Heaviside(t - 2)   t *Heaviside(t - 1)
--------------- + ------------------- - ------------------- - -------------------
       2                   2                     2                     2      

                                                                              
                                            3*Heaviside(t - 3)   Heaviside(t - 1)
- t*Heaviside(t - 3) + t*Heaviside(t - 1) - ------------------ - ----------------
                                                    2                   2        

+ 2*Heaviside(t - 2)

y si reagrupan las expresiones por cada escalón desplazado con sym.collect() se puede obtener:

lista_escalon = ZSR.atoms(sym.Heaviside)
ZSR = sym.expand(ZSR,t) # terminos suma
ZSR = sym.collect(ZSR,lista_escalon)

con el siguiente resultado:

 ZSR(t):
 2                /     2\                    /   2        \                  
t *Heaviside(t)   |    t |                    |  t        1|                  
--------------- + |2 - --|*Heaviside(t - 2) + |- -- + t - -|*Heaviside(t - 1) 
       2          \    2 /                    \  2        2/                  

  / 2        \                 
  |t        3|                 
+ |-- - t - -|*Heaviside(t - 3)
  \2        2/                 

convolucion ejercicio 03

Convolución: [ ejemplo 1: x(t) y h(t) causales ] [ ejemplo 2 : x(t) no causal y h(t) causal ] [ejemplo 3: rectangulo y rampa]


But what is a convolution? 3Blue1Brown. 18 nov 2022

3.4.1 LTI CT Causalidad de h(t) para Integral del Convolución con Sympy-Python

Referencia: Lathi Ej 2.8 p173.

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.

causal h(t) funciones

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.

La función a usar se desarrollada para buscar impulsos unitarios en una expresión es:

donde_d = fcnm.busca_impulso(h)

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 = False
if len(donde_d)>0:
    donde=donde_d[0] # extremo izquierdo
    if 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.
causalidad h(t) 01 Sympy

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 = False
if len(donde_d)>0:
    donde=donde_d[0] # extremo izquierdo
    if donde>=0:  # lado derecho del plano
        hcausal=True

# SALIDA
print('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)
    if len(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:

h(t):
Heaviside(t - 1)
busca_impulso: []
busca_escalon(ft): [[1, 1]]
xcausal:  True

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 = False
if len(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

# SALIDA
print('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.

causalidad h(t) 03 Sympy

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

# PROCEDIMIENTO
def es_causal(ft):
    ''' h(t) es causal si tiene
        b0*d(t)+u(t)*(modos caracterisicos)
    '''
    def es_causal_impulso(ft):
        ''' un termino en lado derecho del plano
        '''
        causal  = False
        donde_d = fcnm.busca_impulso(ft)
        if len(donde_d)>0:    # tiene impulso
            if donde_d[0]>=0: # derecha del plano
                causal = True 
        return(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)
            if len(donde_u)==0: # sin escalon?
                causal_k = False
            if len(donde_u)==1: 
                sentido1 = donde_u[0][0]
                donde1   = donde_u[0][1]
                # en lado izquierdo del plano
                if donde1<0 or sentido1<0:
                    causal_k = False
            if len(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 plano
                if (donde1<donde2): 
                    if donde1<0:  # u(t+1)*u(-t+1)
                        causal_k = False
                if (donde2<donde1):
                    if donde2<0: # u(-t+1)*u(t+1)
                        causal_k = False

        else: # 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)

# SALIDA
print('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

Referencia: Ejercicio Lathi 2.6-3 p206, Lathi Ej. 2.12 p185,

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:

h(t):
-Heaviside(t - 2) + Heaviside(t - 1)
busca_impulso(ft): []
busca_escalon(ft):
[[1 1]
 [2 1]]
hcausal:  True

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)

# SALIDA
print('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()

 

3.4 LTI CT – Respuesta a estado cero ZSR – Desarrollo analítico

Referencia: Lathi 2.4 p168, Oppenheim 2.2.2 p94 , Hsu 2.5.B p60

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 :

y(t) = x(t) \circledast h(t) = \int_{-\infty}^{+\infty} x(\tau)h(t-\tau) \delta \tau

Es importante observar que el integral de convolución se realiza con respecto a τ en lugar de t.

Si h(t) es causal, lineal, continua

Es decir, multiplicada por μ(t), se tiene que, h(t)=0 para t<0 y el integral puede expresarse como:

y(t) = \int_{-\infty}^{+\infty} h(\tau)\mu (\tau) x(t-\tau) \delta \tau = \int_{0}^{\infty} h(\tau) x(t-\tau) \delta \tau

O de forma alterna, aplicando la condición de causalidad (Oppenheim 2.3.6 p112)

y(t) = \int_{0}^{\infty} h(\tau) x(t-\tau) \delta \tau = \int_{-\infty}^{t} x(\tau) h(t-\tau) \delta \tau

Si la entrada x(t) y el sistema h(t) son causales

la respuesta también será causal.

y(t)=\begin{cases}\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau , & t\ge 0\\ 0, & t<0 \end{cases} y(t)=\int_{0^{-}}^{t} x(\tau)h(t-\tau) \delta \tau = \int_{0^{-}}^{t} h(\tau)x(t-\tau) \delta \tau

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.

Respuesta estado cero:  [Desarrollo Analítico]  [Scipy-Python]  [Numpy-Convolución]  [algoritmo Python-Convolución]  [Simulador]

..


Ejemplo 1. Respuesta Estado Cero ZSR con h(t) causal y x(t) causal

Referencia:  Lathi ejemplo 2.6 p166

Encuentre la corriente y(t) del circuito RLC, cuando todas las condiciones iniciales son cero y en la entrada se tiene la señal x(t) descrita por:

x(t) = 10 e^{-3t} \mu (t)

Además el sistema tiene respuesta a impulso:

h(t) = \big( 2e^{-2t} -e^{-t}\big)\mu (t)

El ejemplo es la continuación del presentado para respuesta a entrada cero ZIR, que tiene la ecuación:

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

Desarrollo Analítico

La respuesta se obtiene aplicando convolución entre la señal de entrada x(t) y la respuesta al impulso h(t) del sistema:

y(t) = x(t) \circledast h(t) = [ 10 e^{-3t} \mu (t)] \circledast [(2e^{-2t} - e^{-t}) \mu (t)]

usando la propiedad distributiva de la convolución:

y(t) = [10e^{-3t} \mu (t) \circledast 2e^{-2t} \mu (t)] - [10e^{-3t} \mu (t) \circledast e^{-t} \mu (t)] = 20[e^{-3t}\mu (t) \circledast e^{-2t} \mu (t)] - 10[e^{-3t} \mu(t) \circledast e^{-t} \mu (t)]

Para éste ejercicio se usa la tabla de integrales convolución, línea 4,  asi el enfoque del desarrrollo se mantiene sobre la forma de la señal resultante. El siguiente ejemplo se desarrolla con el integral de convolución.

y(t) = 20\frac{e^{-3t} - e^{-2t}}{-3-(-2)}\mu (t) - 10\frac{e^{-3t} - e^{-t}}{-3-(-1)}\mu (t) y(t) = 20\frac{e^{-3t} - e^{-2t}}{-3+2}\mu (t) - 10\frac{e^{-3t} - e^{-t}}{-3+1}\mu (t) = -20\big[e^{-3t} - e^{-2t}\big]\mu (t) + 5\big[e^{-3t} - e^{-t}\big]\mu (t) = \big[-5e^{-t} + 20e^{-2t} - 15e^{-3t}\big]\mu (t)

respuesta Estado Cero 01 Sympy ZSR

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.

Proponer como tarea.

Respuesta estado cero:  [Desarrollo Analítico]  [Scipy-Python]  [Numpy-Convolución]  [algoritmo Python-Convolución]  [Simulador]


Ejemplo 2. Respuesta Estado Cero ZSR con h(t) causal y x(t) causal

Referencia: Lathi Ejemplo 2.8 p173

Para un sistema LTIC, si la respuesta al impulso es

h(t) = e^{-2t} \mu (t)

determine la respuesta y(t) para la entrada

x(t) = e^{-t} \mu (t)

Desarrollo Analítico

Para éste ejercicio se desarrollará la integral de convolución. La entrada y respuesta al impulso se convierte a:

x(\tau) = e^{-\tau} \mu(\tau) h(t-\tau) = e^{-2(t-\tau)} \mu(t-\tau)

recuerde que la integración es respecto a τ enel intervalo 0≤τ≤t.

y(t) = \begin{cases} \int_{0}^{t} \Big[e^{-\tau}\mu(\tau) e^{-2(t-\tau)}\mu(t-\tau)\Big] \delta\tau , & t\ge 0 \\0, & t \lt 0 \end{cases}

Los valores de u(τ) =1 debido a se convierte a 0 para τ<0 y en el caso de u(t-τ)=1 se convierte a 0 cuando τ≥t.

y(t) = \int_{0}^{t}e^{-\tau} e^{-2(t-\tau)} \delta\tau = \int_{0}^{t} e^{-\tau} e^{-2t} e^{2\tau} \delta\tau = e^{-2t} \int_{0}^{t} e^{\tau} \delta\tau = e^{-2t} e^{\tau} \Big|_0^t = e^{-2t} (e^{t} - 1) = e^{-t} - e^{-2t}

para t≥0, y además como y(t) = 0 para t<0

y(t) = \big( e^{-t} - e^{-2t} \big) \mu(t)

Integral de Convolucion 01 animado


Respuesta estado cero:  [Desarrollo Analítico]  [Scipy-Python]  [Numpy-Convolución]  [algoritmo Python-Convolución]  [Simulador]


Ejemplo 3. Respuesta entrada cero ZSR entre exponencial y escalón unitario

Referencia: Oppenheim Ejemplo 2.6 p98

Sea x(t) la entrada a un sistema LTI con respuesta a impulso unitario h(t),

x(t) = e^{-at} \mu (t) \text{ , } a>0 h(t) = u(t)

dado que la señal de entrada tiene valores para t≥0 al tener un componente μ(t), se tiene que:

y(t) = \int_{0}^{t} x(\tau)h(t-\tau) \delta \tau x(\tau) = e^{-a\tau} \mu (\tau) h(t-\tau) = \mu (t-\tau)

observando que en la región de integración sobre τ se encuentra [0,t], entonces τ≥0 se tiene que μ(τ)=1 y para t-τ≥0 se tiene también que μ(t-τ) = 1.

y(t) = \int_{0}^{t} e^{-a\tau} \delta \tau = -\frac{1}{a} e^{-a\tau} \Big|_0^t = -\frac{1}{a}(e^{-at}-e^{0}) y(t) = \frac{1}{a}(1 - e^{-at})

recordando que esta definido para t≥0

y(t) = \frac{1}{a}(1 - e^{-at}) \mu (t)

Para la gráfica, se define a=2 y se obtiene,

LTIC ZSR_ Ej02 Sympy

El proceso de convolución se observa en la animación realizada al desplazar el varo de tau

 

Respuesta estado cero:  [Desarrollo Analítico]  [Scipy-Python]  [Numpy-Convolución]  [algoritmo Python-Convolución]  [Simulador]

3.3.1 LTI CT – Respuesta a impulso unitario h(t) con Sympy-Python

La respuesta al impulso reutiliza el algoritmo de para encontrar la solución homogenea de la ecuación diferencial lineal. Se reutiliza la función respuesta_ZIR() de la sección anterior.

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]


.. ..


Ejemplo 1. Respuesta a impulso de un sistema RLC

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

Encontrar la Respuesta al impulso h(t), del sistema en el ejemplo 1 Modelo entrada-salida,
representado por el circuito y la ecuación diferencial lineal expresada desde el termino de mayor orden:

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

La expresión en operadores D es:

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

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 mayor
        if (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:
                if not(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.

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

Con la ecuación homogenea  y con las condiciones iniciales, y'(0)=1, y(0)=0, de una entrada impulso, se obtiene como respuesta:

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

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:

respuesta impulso 01 Sympy h(t)

Instrucciones en Python

# Respuesta a impulso h(t) con Sympy-Python
# http://blog.espol.edu.ec/telg1001/lti-ct-respuesta-al-impulso-con-sympy-python/
# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dx
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) + 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 mayor
        if (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:
                if not(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))

# SALIDA
print('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()

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

..


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

respuesta impulso 02 Sympy h(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
# http://blog.espol.edu.ec/telg1001/lti-ct-respuesta-al-impulso-con-sympy-python/
# Lathi 2.1.a pdf 155, (D+2)y = (3D+5)x
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,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

# PROCEDIMIENTO
def respuesta_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 mayor
            if (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:
                    if not(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']

# SALIDA
print('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()

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

..


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

respuesta impulso 03 Sympy

Instrucciones en Python

# Respuesta a impulso h(t) con Sympy-Python
# http://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']

# SALIDA
print('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()

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

3.3 LTI CT – Respuesta a impulso unitario h(t) – Desarrollo analítico

Referencia: Lathi 2.3 p163, Ejemplo 2.4 p159, Oppenheim problema 2.2.2 p97, Hsu 2.5.E p61

La respuesta de un sistema al impulso h(t), se obtiene al aplicar un impulso unitario δ(t) en la entrada x(t), semejante a un destello durante un tiempo muy pequeño.

El impulso es semejante a tomar una foto con flash en una habitación oscura, con todo tranquilo, sin movimientos. Luego del flash, a pesar que ya no exista más la luz, todos las personas o seres vivos reaccionan al destello de luz durante un tiempo y visualizaron el contenido de la habitación, ésta es la respuesta al impulso del sistema.

La respuesta del sistema se puede conformar como la suma de una señal contínua generada mediante una secuencia de impulsos cuando el tiempo entre impulsos tiende a cero (dt→0).

Respuesta
total
= respuesta a
entrada cero
+ respuesta a
estado cero
ZSR = h(t) ⊗ x(t)

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

..


1. Respuesta a impulso h(t): Desarrollo Analítico

Un sistema de orden 2 se describe con la expresión de operadores D como:

Q(D)y(t) = P(D) x(t) (a_0 D^2+a_1 D + a_2)y(t)=(b_0 D^2 + b_1 D + b_2) x (t)

para disminuir la cantidad de variables se hace a0=1, al dividir la ecuación en ambos lados para a0. El resultado como aún tiene variables desconocidas, se mantiene la nomenclatura, que de forma simplificada se presenta como:

( D^2+a_1 D + a_2)y(t)=(b_0 D^2 + b_1 D + b_2) x (t)

La respuesta de un sistema al impulso unitario δ(t) aplicada en t=0, con todas las condiciones iniciales en cero en t=0-, se conoce como h(t).

(D^2+a_1 D + a_2)h(t)=(b_0 D^2 + b_1 D + b_2) \delta (t)

Para el análisis se consideran dos intervalos:

– Para t≥0+ donde x(t)=0, cuya respuesta es conocida como los términos de modos característicos o respuesta homogenea desarrollada en respuesta entrada cero ZIR

(D^2+a_1 D + a_2)h(t)=0, t \ge 0

h(t) = términos de modos característicos

– Para t=0 donde existe x(t)=δ(t), que suma a la respuesta anterior una constante A0 por un impulso.

h(t) = (términos de modos característicos) + A0 δ(t)

que al insertar esta ecuación en la expresión original tenemos:

(D^2+a_1 D + a_2)(A_0 \delta(t) + \text{terminos modos caracteristicos}) = =(b_0 D^2 + b_1 D + b_2) \delta (t)

que al multiplicar la parte de la izquierda, y comparando términos de ambos lados se determina que A0 = b0

que es como expresar

h(t) = b_0 \delta (t) + \text{modos caracteristicos}

y si el orden del lado derecho M es menor que el del izquierdo N, b0=0.

h(t)=b_0 \delta (t)+ [P(D)y_n (t)] \mu (t)

Revisar Lathi 2.3 pdf 163 para una descripción más detallada del método.


Método simplificado al emparejar impulsos

Referencia: Lathi 2.3 p165

El método permite reducir el procedimiento para determinar h(t) de un sistema LTIC.

h(t) = b_0 \delta (t) + [P(D) y_n(t)] \mu (t)

donde yn(t) es una combinación de modos característicos del sistema sujeto a las condiciones iniciales siguientes:

y_n(0) = y'_n(0) = y''_n(0) = y_n^{(N-2)}(0) = 0 y_n^{(N-1)}(0) = 1

Donde y(k)n(0) es el valor de la k-ésima derivada de yn(t) para t=0.

Si orden de P(D) es menor que el orden de Q(D) entonces b0=0 y el término del impulso b0δ(t) en h(t) es cero.

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

Ejemplo 1. Respuesta al impulso de una ecuación diferencial lineal de 2do orden

Referencia: Lathi Ejemplo 2.6 p166, Oppenheim 2.2.2 p94

Un sistema se describe con la ecuación mostrada y se pide determinar la respuesta al impulso unitario h(t).

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

Desarrollo Analítico

Semejante al análisis realizado para respuesta a entrada cero ZIR, el problema tiene un sistema se segundo orden con polinomio característico:

\lambda ^2 + 3\lambda +2 = (\lambda +1)(\lambda + 2)
Raíces características Modos característicos
λ1 = -1 e-t
λ2 = -2 e-2t

con lo que la solución general tienen la forma:

y_0 (t) = c_1 e^{-t} + c_2 e^{-2t}

se toman las condiciones iniciales de la tabla del método simplificado para emparejar impulsos para N=2, por ser un sistema de segundo orden

y(0) = 0 , y'(0) = 1

y se aplican en la solución general:

y_0 (t) = c_1 e^{-t} + c_2 e^{-2t} y'_0 (t)= -c_1 e^{-t} -2c_2 e^{-2t}

al sustituir los valores iniciales en la solución general con t=0, se obtienen los valores de las constantes,

0 =  c1 + c2 
1 = -c1 - 2c2 

c1 = 1 
c2 = -1

obteniendo como solución de la ecuación homogenea:

y_n (t) = 1 e^{-t} - 1e^{-2t}

De la ecuación del problema, se tiene que del lado derecho de la ecuación, P(D)=D, entonces:

P(D)y_n (t) = Dy_n (t) = y'_n (t) = -e^{-t} + 2e^{-2t}

Dado el caso que en la ecuacion, el orden de la derivadas de la derecha M es menor que el orden de derivadas de la izquierda N, (N>M), el término de segundo orden no se encuentra en P(D), entonces: b0=0

h(t) = b_0 \delta (t) + [P(D) y_n(t)] \mu (t) h(t)= 0\delta (t) + [-e^{-t} + 2e^{-2t}] \mu (t) h(t)= (-e^{-t} + 2e^{-2t})\mu (t)

LTIC Ejercicio 01 ht Sympy

Respuesta a Impulso: [Desarrollo Analítico]  [Sympy-Python]  [SciPy-Python]  [Runge-Kutta]  [Simulador]

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

Para encontrar la  Respuesta a entrada cero del sistema lineal del Modelo entrada-salida, circuito RLC ejemplo respuesta impulsorepresentado por un circuito, se plantea la ecuación diferencial lineal escrita desde el termino de mayor orden. Por ejemplo:

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

Luego se se simplifica haciendo x(t)=0, obteniendo ecuación complementaria relacionada a la forma homogénea de la ecuación diferencial.

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


Ejemplo 1. Respuesta a entrada cero ZIR para un sistema LTI CT – Desarrollo analítico con Sympy

Referencia: Lathi Ejemplo 2.1.a p155, Oppenheim 2.4.1 p117, ejemplo 1 de Modelo entrada-salida, Sympy ODE solver

El circuito que representa la respuesta a entrada cero tiene x(t)=0 como se muestra en la figura.

La ecuación diferencial con operadores D es:

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

Adicionalmente, para el desarrollo se indican las condiciones iniciales y'(0)=-5, y(0)=0.

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

Para iniciar el algoritmo, se define la variable independiente t como un símbolo y las función matemática y(t) para la salida y x(t) para la entrada.

La ecuación diferencial se escribe siguiendo el orden de los lados de expresión denominados para la parte izquierda, LHS, y la parte derecha, RHS.

import sympy as sym

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

# ecuacion EDO: lado izquierdo = lado derecho
#               Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,evaluate=False)
ecuacion_edo = sym.Eq(LHS,RHS)

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

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

Solución general de la ecuación diferencial lineal homogénea

Se obtiene la ecuación homogénea aplicando x(t)=0 al lado derecho de la ecuación diferencial (RHS).

La solución general de la EDO se obtiene mediante la instrucción sym.dsolve() aplicada a la ecuación, indicando que se busca resolver para y(t). Para facilitar el procesamiento con Sympy, las ecuaciones se expresan con  términos mas simples de sumas, aplicando la instrucción sym.expand().

# ecuacion homogenea 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 entrada cero
general = sym.dsolve(homogenea, y(t))
general = general.expand()

Con lo que se obtiene la solución general mostrada con sym.pprint():

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

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

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

Se realiza lo mismo para y(0)=0

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

# aplica condiciones iniciales 
N = sym.ode_order(LHS,y) # orden Q(D)
eq_condicion = []
for k in range(0,N,1):
    ck   = cond_inicio[(N-1)-k]
    dyk  = general.rhs.diff(t,k)
    dyk  = dyk.subs(t,t0)
    eq_k = sym.Eq(ck,dyk)
    eq_condicion.append(eq_k)
    
constante = sym.solve(eq_condicion)

El resultado de eq_condicion y las constantes obtenidas se muestra a continuación

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

La solución homogénea o  ZIR de la ecuación diferencial se obtiene al sustituir las constantes encontradas en la ecuación general.

# reemplaza las constantes en general
y_c = general
for Ci in constante:
    y_c = y_c.subs(Ci, constante[Ci])
# ecuacion complementaria o respuesta a entrada cero ZIR
ZIR = y_c.rhs

obteniendo solución a la ecuación homogena ‘y‘, o respuesta a entrada cero ZIR, mostrada con sym.pprint(y_homogenea). Para ZIR se toma el dado derecho de y_homogenea.

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

resultado del algoritmo:

ecuación diferencial a resolver: 
                        2                 
           d           d          d       
2*y(t) + 3*--(y(t)) + ---(y(t)) = --(x(t))
           dt           2         dt      
                      dt                  
clasifica EDO:
  factorable
  nth_linear_constant_coeff_variation_of_parameters
  nth_linear_constant_coeff_variation_of_parameters_Integral
ecuación diferencial homogenea: 
                        2          
           d           d           
2*y(t) + 3*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           

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

Ecuaciones de condiciones iniciales:
0 = C1 + C2
-5 = -C1 - 2*C2
constantes en solucion general:  {C1: -5, C2: 5}
ZIR(t):
     -t      -2*t
- 5*e   + 5*e       
>>> 

Instrucciones en Python

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

# 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,evaluate=False)
ecuacion = sym.Eq(LHS,RHS)

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

# 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
general = sym.dsolve(homogenea, y(t))
general = general.expand()

# aplica condiciones iniciales 
N = sym.ode_order(ecuacion,y(t)) # orden Q(D)
eq_condicion = []
for k in range(0,N,1):
    ck   = cond_inicio[(N-1)-k]
    dyk  = general.rhs.diff(t,k)
    dyk  = dyk.subs(t,t0)
    eq_k = sym.Eq(ck,dyk)
    eq_condicion.append(eq_k)
    
constante = sym.solve(eq_condicion)

# reemplaza las constantes en general
y_c = general
for Ci in constante:
    y_c = y_c.subs(Ci, constante[Ci])
# respuesta a entrada cero ZIR
ZIR = y_c.rhs

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

print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)

print('ecuación diferencial homogenea: ')
sym.pprint(homogenea)

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

print('\nEcuaciones de condiciones iniciales:')
for eq_k in eq_condicion:
    sym.pprint(eq_k)
print('constantes en solucion general: ', constante)
print('\nZIR(t):')
sym.pprint(ZIR)

Gráfica de respuesta a entrada cero ZIR

respuesta Entrada Cero 01 Sympy ZIR

Adicionalmente se realiza la gráfica, convirtiendo la ecuación de forma simbólica a forma numérica numpy con sym.lambdify(t,ZIR,'numpy').

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

Se añaden las siguientes instrucciones al algoritmo anterior:

# 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,ZIR, modules=equivalentes) 
ti = np.linspace(t_a,t_b,muestras)
yi = yt(ti)

# Grafica
plt.plot(ti,yi, color='orange', label='ZIR(t)')
untitulo = 'ZIR(t)=$'+ str(sym.latex(ZIR))+'$'
plt.title(untitulo)
plt.xlabel('t')
plt.ylabel('ZIR(t)')
plt.legend()
plt.grid()
plt.show()

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

..


Ejemplo 2. Respuesta a entrada cero ZIR con raíces características repetidas

Referencia: Lathi Ejemplo 2.1.b p155

Encontrar y0(t), la respuesta a entrada cero para un sistema LTI CT descrito  por la ecuación diferencial lineal de orden 2:

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

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

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

clasifica EDO:
  factorable
  nth_linear_constant_coeff_variation_of_parameters
  nth_linear_constant_coeff_variation_of_parameters_Integral
homogenea :
                        2          
           d           d           
9*y(t) + 6*--(y(t)) + ---(y(t)) = 0
           dt           2          
                      dt           
general :
           -3*t         -3*t
y(t) = C1*e     + C2*t*e    
eq_condicion :
3 = C1
-7 = -3*C1 + C2
constante : {C1: 3, C2: 2}
ZIR :
     -3*t      -3*t
2*t*e     + 3*e        

respuesta Entrada Cero 02 Sympy ZIR

Instrucciones en Python

Se reordena el algoritmo del ejercicio anterior para disponer de funciones que simplifiquen las instrucciones principales para obtener las respuestas de los ejercicios. Se crean las funciones respuesta_ZIR() y graficar_ft() que luego pasarán a ser parte de los algoritmos del curso en telg1001.py.

Para las definiciones de RHS y LHS dentro de la función se usa ecuación.rhs y ecuacion.lhs.

En el diccionario solución, se incorpora la respuesta en la entrada 'ZIR'.

# Respuesta a entrada cero ZIR con Sympy-Python
# http://blog.espol.edu.ec/telg1001/lti-ct-respuesta-entrada-cero-con-sympy-python/
# Lathi 2.1.b pdf 155, (D^2+ 6D + 9)y = (3D+5)x
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'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')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 6*sym.diff(y(t),t,1) + 9*y(t)
RHS = 3*sym.diff(x(t),t,1,evaluate=False)+5*x(t)
ecuacion = sym.Eq(LHS,RHS)

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

# Grafica: intervalo tiempo [t_a,t_b]
t_a = 0 ; t_b = 5
muestras = 51

# PROCEDIMIENTO
def respuesta_ZIR(ecuacion,cond_inicio=[],t0,
                y = sym.Function('y'),x = sym.Function('x')):
    ''' Sympy: ecuacion: diferencial en forma(LHS,RHS),
        condiciones de inicio t0 y [y'(t0),y(t0)]
        cond_inicio en orden descendente derivada
    '''
    # ecuacion homogenea 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 entrada cero
    general = sym.dsolve(homogenea,y(t))
    general = general.expand()

    # aplica condiciones iniciales 
    N = sym.ode_order(ecuacion,y(t)) # orden Q(D)
    eq_condicion = []
    for k in range(0,N,1):
        ck   = cond_inicio[(N-1)-k]
        dyk  = general.rhs.diff(t,k)
        dyk  = dyk.subs(t,t0)
        eq_k = sym.Eq(ck,dyk)
        eq_condicion.append(eq_k)
        
    constante = sym.solve(eq_condicion)

    # reemplaza las constantes en general
    y_c = general
    for Ci in constante:
        y_c = y_c.subs(Ci, constante[Ci])
    # respuesta a entrada cero ZIR
    ZIR = y_c.rhs
    
    sol_ZIR = {'homogenea'   : sym.Eq(homogenea,0),
               'general'     : general,
               'eq_condicion': eq_condicion,
               'constante'   : constante,
               'ZIR' : ZIR,}
    return(sol_ZIR)

edo_tipo = sym.classify_ode(ecuacion, y(t))
# Respuesta entrada cero ZIR
sol_ZIR = respuesta_ZIR(ecuacion,cond_inicio,t0)
ZIR = sol_ZIR['ZIR']

# SALIDA
print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)
fcnm.print_resultado_dict(sol_ZIR)

# GRAFICA ------------------
figura_ZIR = fcnm.graficar_ft(ZIR,t_a,t_b,muestras,'ZIR')
plt.show()

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


Ejemplo 3. EDO y respuesta a entrada cero ZIR con raíces complejas

Referencia: Lathi Ejemplo 2.1.c p155

Encontrar y0(t), la respuesta a entrada cero para un sistema LTI CT descrito  por la ecuación diferencial ordinaria de orden 2:

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

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

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

Los resultados con el algoritmo son:

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

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

Ecuaciones de condiciones iniciales:
2 = C2
16.78 = 6*C1 - 2*C2
constantes en solución general:  {C1: 3.46333333333333, C2: 2.00000000000000}

Solución ZIR:
                  -2*t                 -2*t         
3.46333333333333*e    *sin(6*t) + 2.0*e    *cos(6*t)
>>>

respuesta Entrada Cero 03 Sympy ZIR
Instrucciones con Python

# Respuesta a entrada cero ZIR con Sympy-Python
# http://blog.espol.edu.ec/telg1001/lti-ct-respuesta-entrada-cero-con-sympy-python/
# Lathi 2.1.c pdf 155, (D^2+ 4D + 40)y = (D+2)x
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')

# ecuacion: lado izquierdo = lado derecho
#           Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 4*sym.diff(y(t),t,1) + 40*y(t)
RHS = sym.diff(x(t),t,1,evaluate=False)+2*x(t)
ecuacion = sym.Eq(LHS,RHS)

# condiciones iniciales [y'(t0),y(t0)]
t0 = 0
cond_inicio = [16.78,2]

# 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 entrada cero ZIR
sol_ZIR = fcnm.respuesta_ZIR(ecuacion,cond_inicio,t0)
ZIR = sol_ZIR['ZIR']

# SALIDA
print('clasifica EDO:')
for elemento in edo_tipo:
    print(' ',elemento)
fcnm.print_resultado_dict(sol_ZIR)

# GRAFICA ------------------
figura_ZIR = fcnm.graficar_ft(ZIR,t_a,t_b,muestras,'ZIR')
plt.show()

Resumen de funciones de Señales y Sistemas telg1001.py

Las funciones desarrolladas para uso posterior se agrupan en un archivo resumen telg1001.py.

Use una copia del archivo como telg1001.py en el directorio de trabajo, e incorpore las funciones  con import y use en el algoritmo llamando a la funcion con fcnm.funcion().

iimport telg1001 as fcnm

respuesta = fcnm.funcion(parametros)

Al desarrollar más funciones en cada unidad, se incorporan al archivo del curso de Señales y Sistemas.

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() para los demás elementos del resultado

def print_resultado_dict(resultado):
    ''' print de diccionario resultado
        formato de pantalla
    '''
    eq_sistema = ['ZIR','h','ZSR','xh']
    for entrada in resultado:
        tipo = type(resultado[entrada])
        cond = (tipo == sym.core.relational.Equality)
        cond = cond or (entrada in eq_sistema)
        if cond:
            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()