4.2.3 Transformada Inversa de Laplace de F(s) con Sympy-Python

Referencia: Lathi 4.1. p330. Hsu 3.5.A p119, Oppenheim 9.3 p670

Se dice que la señal x(t) es la transformada inversa de Laplace X(s) se determina como:

x(t) = \frac{1}{2 \pi j} \int_{c - j \infty}^{c + j \infty} X(s) e^{st} ds

donde c es una constante seleccionada para asegurar la convergencia de la integral.

La operación para encontrar la transformada inversa de Laplace requiere un integral en el plano complejo. El camino de integración es a lo largo de c+jω, siendo que ω varía entre -∞ a ∞.

Para la señal e(-at)μ(t) wa posible si c>-a, por ejemplo para un punto c=1 con ω desde -∞ a ∞. Sin embargo esto requiere aplicar conocimientos en teoría de funciones de variable compleja. Es posible evitar estos ‘detalles’ usando la tabla de transformadas de Laplace, donde encontrar la transformada inversa consiste en buscar el modelo de expresión en el domino ‘s‘ y buscar la pareja en el dominio del tiempo ‘t‘.

La transformada inversa de Laplace con Sympy tiene la instrucción sym.inverse_laplace_transform(Fs,s,t), que para términos simples, facilita el proceso de desarrollar del integral hacia el dominio del tiempo. Para simplificar los términos de la expresión F(s) se usan las instrucciones como sym.expand(Fs,s) o para fracciones parciales sym.apart(Fs,s).

Transformada  Inversa de Laplace: [ ej1 H(s) función de transferencia ]  [ej2 con desplazamiento ] [ ej3 suma de términos s ]   [ ej4 Y(s)=H(s)*X(s) con escalones deplazados ]


Ejemplo 1. Transformada Inversa de Laplace de una función de transferencia en dominio ‘s’, H(s)

Referencia: Lathi 1.8-1 p111. Oppenheim Ejemplo 9.24 p700

El ejercicio tiene como referencia la función de transferencia del ejercicio desarrollado para el Ejemplo 1. Corriente en circuito RLC del modelo de entrada-salida. Se requiere obtener la transformada inversa de Laplace:

F(s)= \frac{s}{s^2+3s+2}
F(s): 
  2       1  
----- - -----
s + 2   s + 1

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

Instrucciones en Python

La expresión de Fs es una fracción que contiene polinomios de ‘s’ y se puede escribir en una sola línea para éste caso. Para disponer de términos más simples, se aplica sym.apart(Fs,s). Para mostrar el resultado de f(t) por cada término de fracción parcial se aplica sym.expand(Fs,s).

# Transformada Inversa de Laplace con Sympy
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True)
s = sym.Symbol('s')

Fs = s/(s**2+3*s+2)

# coeficientes como racional en dominio 'ZZ'
#k1 = sym.Rational(1/3).limit_denominator(100)
#k2 = sym.Rational(4/3).limit_denominator(100)
#Fs = 1 - k2*1/(s+1) + k1*1/(s-2)

#Fs = (1/s)*(1-sym.exp(-2*s))

#Hs = sym.exp(s)/s - 2*sym.exp(-s)/s + sym.exp(-3*s)/s
#Xs = 2*sym.exp(-s)/s - 2*sym.exp(-3*s)/s
#Fs = Hs*Xs

# PROCEDIMIENTO
# convierte a Sympy si es solo constante
Fs = sym.sympify(Fs)
# separa exponenciales constantes
Fs = sym.expand_power_exp(Fs)
Fs = sym.expand(Fs,s) # terminos suma

# Fracciones parciales
if not(Fs.has(sym.exp)):
    Fs = sym.apart(Fs,s)

ft  = sym.inverse_laplace_transform(Fs,s,t)

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

# SALIDA
print('\n F(s): ')
sym.pprint(Fs)

print('\n f(t): ')
sym.pprint(ft)

Transformada  Inversa de Laplace: [ ej1 H(s) función de transferencia ]  [ej2 con desplazamiento ] [ ej3 suma de términos s ]   [ ej4 Y(s)=H(s)*X(s) con escalones deplazados ]
..


Ejemplo 2. Transformada Inversa de Laplace con impulso δ(t) y suma de términos

Referencia: Oppenheim ejemplo 9.5 p661

Se requiere realizar el proceso contrario a lo desarrollado en el Ejemplo 4 de Transformadas de Laplace con Sympy. En el ejercicio se expone sobre uso de los coeficientes en forma de enteros o fracciones.

Obtener la función en el dominio del tiempo:

F(s) = 1 - \frac{4}{3}\frac{1}{s+1} + \frac{1}{3}\frac{1}{s-2}

La expresión de ingreso para F(s) en el algoritmo anterior se escribe como,

# coeficientes como racional en dominio 'ZZ'
k1 = sym.Rational(1/3).limit_denominator(100)
k2 = sym.Rational(4/3).limit_denominator(100)

Fs = 1 - k2*1/(s+1) + k1*1/(s-2)

con lo que se obteiene el resultado esperado y acorde al ejercicio original de la referencia:

F(s): 
        4           1    
1 - --------- + ---------
    3*(s + 1)   3*(s - 2)

 f(t): 
/ 2*t      -t\                             
|e      4*e  |                             
|---- - -----|*Heaviside(t) + DiracDelta(t)
\ 3       3  /                             
>>>  

Transformada  Inversa de Laplace: [ ej1 H(s) función de transferencia ]  [ej2 con desplazamiento ] [ ej3 suma de términos s ]   [ ej4 Y(s)=H(s)*X(s) con escalones deplazados ]


Ejemplo 3 Transformada Inversa de Laplace con términos de desplazamiento en tiempo

Referencia: Lathi Ejercicio 4.1.a p337

Se requiere la transformada inversa de Laplace del resultado ejercicio 2  desarrollado con Sympy para una función gate o compuerta:

F(s)= \frac{1}{s} \Big( 1 - e^{-2s} \Big)

En el caso de usar desplazamientos en tiempo, se recomienda también usar las expresiones simples de suma, para que por cada término de suma aplicar la instrucción de la transformada inversa. Recuerde usar sym.expand(Fs,s) y sym.apart(Fs,s) antes de aplicar la tranformada inversa.

La expresión en Sympy de la entrada es:

Fs = (1/s)*(1-sym.exp(-2*s))

con lo que el resultado a obtener con el algoritmo del ejemplo 4 es:

 F(s): 
     -2*s
1   e    
- - -----
s     s  

 f(t): 
Heaviside(t) - Heaviside(t - 2)
>>>   

La expresión de F(s) al aplicar directamente fracciones parciales con sym.apart() se muestra un error por tener un elemento exponencial sym.exp()en el numerador.

Como los ejercicios a resolver tienen varios términos que se multiplican o que se suman, se procede crear una función para los procesos de expansion en fracciones parciales y transformadas inversas con fracciones parciales.

Transformada  Inversa de Laplace: [ ej1 H(s) función de transferencia ]  [ej2 con desplazamiento ] [ ej3 suma de términos s ]   [ ej4 Y(s)=H(s)*X(s) con escalones deplazados ]
..


Ejemplo 4. y(t) desde h(t) y x(t) con términos escalón desplazados

Referencia: 1Eva2009TII_T3 LTI CT y(t) desde h(t) y x(t) con términos escalón desplazados

Se simplifica el enunciado del ejercicio, enfocandose en que la transformada se aplica a F(s)  y que es el producto de la señal de entrada X(s) y la respuesta al impulso H(s),

F(s) = H(s)*X(s) H(s) = \frac{1}{s}e^{s} - 2\frac{1}{s}e^{-s}+ \frac{1}{s} e^{-3s} X(s) = 2\frac{1}{s}e^{-s} - 2\frac{1}{s}e^{-3s}

La expresión a usar para la transformada inversa de Laplace se convierte en:

F(s) = \Bigg[ \frac{1}{s}e^{s} - 2\frac{1}{s}e^{-s}+ \frac{1}{s} e^{-3s} \Bigg ] \Bigg [ 2\frac{1}{s}e^{-s} - 2\frac{1}{s}e^{-3s} \Bigg]

se obtiene términos de la transformadas como:

 F(s): 
        -2*s      -4*s      -6*s
2    6*e       6*e       2*e    
-- - ------- + ------- - -------
 2       2         2         2  
s       s         s         s   

 f(t): 
2*t*Heaviside(t) + (12 - 6*t)*Heaviside(t - 2) + 
 (12 - 2*t)*Heaviside(t - 6) 
+ (6*t - 24)*Heaviside(t - 4)
>>> 

Transformada  Inversa de Laplace: [ ej1 H(s) función de transferencia ]  [ej2 con desplazamiento ] [ ej3 suma de términos s ]   [ ej4 Y(s)=H(s)*X(s) con escalones deplazados ]

4.2.2 Transformada de Laplace para f(t) con Sympy-Python

Sympy ofrece la instrucción sym.laplace_transform(ft,t,s) para expresiones de f(t) con términos simples. La instrucción desarrolla el integral unilateral con t entre [0,∞], es decir con entradas tipo causal con t>0 o con términos μ(t) y los desplazados hacia la derecha μ(t-1).

Partiendo de las variable ‘t‘ y ‘s‘ como símbolos , se establece la expresión correspondiente en f(t) para determinar F(s). La transformada se obtiene al usar la instrucción:

Fs_L = sym.laplace_transform(ft,t,s)
Fs = Fs_L[0] # solo expresion

El resultado contiene la expresión, el valor de un polo del plano de convergencia y una condición de convergencia auxiliar. Para los objetivos de los ejercicios el enfoque es sobre el primer componente Fs = Fs_L[0].

En los ejercicios desarrollados se describen las ventajas y restricciones al usar las instrucciones librería Sympy, versión 1.11.1 o superior al 2022.10.30. Se indica que el inconveniente está resuelto en la versión 1.12 22/11/2022 https://github.com/sympy/sympy/issues/24294

Por simplicidad, el análisis de polos y ceros se realiza en la sección de estabilidad del sistema con Python

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ] [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]

..


Ejemplo 1. Transformada de Laplace de una exponencial decreciente, un solo termino

Referencia: Lathi ejemplo 4.1. p331, Oppenheim Ejemplo 9.2 p656, Hsu Ejemplo 3.1 p111

Para una señal f(t),Transformada Laplace f(t) Ej01

f(t) = e^{-at} \mu (t)

se tiene la transformada F(s),

F(s) = \frac{1}{s+a}

realizar la transformada con la instrucción directa de Sympy:

Siendo las variables t, u de tipo símbolo, se definen las funciones como,

u = sym.Heaviside(t)
ft = sym.exp(-a*t)*u

el resultado de la operación será:

 f(t): 
 -a*t
e    *Heaviside(t)

 F(s): 
  1  
-----
a + s

Se puede verificar el resultado en la Tabla de Transformadas de Laplace

Instrucciones con Python  para Ejemplo 1

# Transformadas de Laplace Unilateral con Sympy
# supone f(t) causal, usa escalón u(t)
import sympy as sym

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

ft = sym.exp(-a*t)*u

# PROCEDIMIENTO
ft = sym.expand(ft,t) # términos de suma
Fs_L = sym.laplace_transform(ft,t,s)
Fs = Fs_L[0] # solo expresion

# SALIDA
print('\n f(t): ')
sym.pprint(ft)
print('\n F(s): ')
sym.pprint(Fs)

Realizado el primer ejemplo con las instrucciones Sympy, se obtiene una guia para continuar con otros casos de ejercicios.

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ] [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]


Ejemplo 2. Transformada de Laplace para suma de términos f(t) con desplazamiento, o función «gate» o compuerta

Referencia: Lathi práctica 4.1.a p337

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

Transformada Laplace f(t) Ej03El bloque de ingreso se expresa como:

u = sym.Heaviside(t)
ft = u - u.subs(t,t-2)

Al aplicar el algoritmo anterior, modificando la expresión ft, la Transformada de Laplace muestra que el término desplazamiento tiene un componente exponencial.

 f(t): 
Heaviside(t) - Heaviside(t - 2)

 F(s): 
     -2*s
1   e    
- - -----
s     s  
>>> 

Para considerar el término exponencial en el cálculo de polos,  se separa el ejercicio en partes con o sin  sym.exp(-a*s) cuando se realice el análisis de estabilidad del sistema.

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ] [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]

..


Ejemplo 3 Transformada de Laplace con cos(t)

Referencia: Oppenheim Ejemplo 9.4 p658

Encontrar la transformada de Laplace para:

x(t) = e^{-2t}\mu (t) + e^{-t} \cos (3t) \mu (t)

Transformada Laplace ft Ej05

Para el algoritmo, la expresión se escribe como

ft = sym.exp(-2*t)*u + sym.exp(-t)*sym.cos(3*t)*u

con lo que se obtiene como resultado,

 f(t): 
 -t                          -2*t             
e  *cos(3*t)*Heaviside(t) + e    *Heaviside(t)

 F(s): 
      2              
   2*s  + 5*s + 12   
---------------------
 3      2            
s  + 4*s  + 14*s + 20 
>>> 

Para disponer de expresiones mas simples de F(s) en fracciones parciales, se añade la instrucción

Fs = sym.apart(Fs) # separa en fracciones parciales

con lo que se obtiene el siguiente resultado para F(s),

 F(s): 
    s + 1         1  
------------- + -----
 2              s + 2
s  + 2*s + 10        
>>> 

El primer componente de la suma corresponde a la parte de sym.exp(-t)*sym.cos(3*t) y la segunda parte corresponde a sym.exp(-2*t)

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ] [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]


Ejemplo 4 Transformada de Laplace con Impulso unitario δ(t) y suma de exponenciales

Referencia: Oppenheim ejemplo 9.5 p661

x(t) = \delta(t) -\frac{4}{3} e^{-t} \mu (t) + \frac{1}{3} e^{2t} \mu (t)

La expresión de f(t) podría escribirse directamente como:

d = sym.DiracDelta(t)
u = sym.Heaviside(t)
ft = d - (4/3)*sym.exp(-t)*u + (1/3)*sym.exp(2*t)*u

Al usar la instrucción sym.laplace_transform(ft,t,s)  convierte las fracciones a números reales o con decimales.

Nota: Sympy hasta la versión 1.11.1, las operaciones en el dominio ‘s’ para la Transformadas Inversas de Laplace se encuentran implementadas para manejar principalmente números enteros y fracciones. Los resultados de expresiones combinadas con coeficientes enteros y coeficientes reales no necesariamente se simplifican entre si, pues se manejan diferentes dominios ‘ZZ’ o ‘QQ’. (Revisión 2022-Nov)

Para optimizar la simplificación de expresiones con coeficientes entre enteros y reales, los números reales se convierten a su aproximación racional con la instrucción sym.Rational(0.333333).limit_denominator(100).

Convirtiendo los coeficientes a racionales, se define ft como:

u = sym.Heaviside(t)
d = sym.DiracDelta(t)

# coeficientes como racional en dominio 'ZZ' enteros
k1 = sym.Rational(1/3).limit_denominator(100)
k2 = sym.Rational(4/3).limit_denominator(100)

ft = d - k2*sym.exp(-t)*u + k1*sym.exp(2*t)*u

Al observar los resultados con el algoritmo se puede observar que no se ha procesado el valor de la constante en la transformada.

 f(t): 
 2*t                                   -t             
e   *Heaviside(t)                   4*e  *Heaviside(t)
----------------- + DiracDelta(t) - ------------------
        3                                   3         

 F(s): 
      1       1    
1 - ----- + -----   #faltan las constantes...
    s + 1   s - 2
>>>

Se encuentra que la instrucción sym.laplace_transform(ft,t,s) no ha procesado la constante cuando f(t) tiene mas de dos componentes que se multiplican, (Sympy version 1.11.1 revisado hasta 2022-Nov). Asunto que se encuentra a la fecha bajo revisión segun el enlace:

https://github.com/sympy/sympy/issues/23360

Mientras tanto, para obtener resultados e identificado el asunto, se crea una función separa_constante() para un término, donde se separa la constante como el término multiplicador de las partes (args) que no contienen la variable ‘t’.

    def separa_constante(termino):
        ''' separa constante antes de usar
            sym.laplace_transform(term_suma,t,s)
            para incorporarla luego de la transformada
            inconveniente revisado en version 1.11.1
        '''
        constante = 1
        if termino.is_Mul:
            factor_mul = sym.Mul.make_args(termino)
            for factor_k in factor_mul:
                if not(factor_k.has(t)):
                    constante = constante*factor_k
            termino = termino/constante
        return([termino,constante])

usando el resultado previo del algoritmo, se prueba la función con el último termino de la suma. Luego de separar la constante, se aplica la transformada de Laplace de Sympy y se incorpora la constante al resultado.

>>> k1 = sym.Rational(1/3).limit_denominator(100)
>>> ft = k1*sym.exp(2*t)*u
>>> sym.laplace_transform(ft,t,s)
(1/(s - 2), 2, True)
>>> [termino,constante] = separa_constante(ft)
>>> termino
exp(2*t)*Heaviside(t)
>>> constante
1/3
>>> sym.laplace_transform(termino,t,s)[0]*constante
1/(3*(s - 2))

En el ejemplo se muestra que es necesario separar la constante en al menos dos términos de suma, por lo que  se debe considerar el caso de que f(t) sea de uno o varios términos suma. Para simplificar el proceso en los próximos ejercicios se crea la función laplace_term_suma(ft) que se encargará de realizar el proceso término a término.

    ft = sym.expand(ft)  # expresion de sumas
    ft = sym.powsimp(ft) # simplifica exponentes

    term_suma = sym.Add.make_args(ft)
    Fs = 0
    for term_k in term_suma:
        [term_k,constante] = separa_constante(term_k)
        Fsk = sym.laplace_transform(term_k,t,s)
        Fs  = Fs + Fsk[0]*constante

Al incorporar las funciones al algoritmo, se puede verificar que se obtienen los resultados obtenidos en la forma analítica.

 f(t): 
 2*t                                   -t             
e   *Heaviside(t)                   4*e  *Heaviside(t)
----------------- + DiracDelta(t) - ------------------
        3                                   3         

 F(s): 
        4           1    
1 - --------- + ---------
    3*(s + 1)   3*(s - 2)
>>>

Se prueban todos los ejercicios revisados con sus respuestas y se tiene como algoritmo:

# Transformadas de Laplace Unilateral con Sympy
# supone f(t) causal, usa escalón u(t)
import sympy as sym

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

# coeficientes como racional en dominio 'ZZ' enteros
k1 = sym.Rational(1/3).limit_denominator(100)
k2 = sym.Rational(4/3).limit_denominator(100)

ft = d - k2*sym.exp(-t)*u + k1*sym.exp(2*t)*u

#ft = d - 3*d.subs(t,t-2) + 2*d.subs(t,t-3)
#ft = k2*sym.exp(-2*t)*u - k1*sym.exp(-2*t)*u.subs(t,t-5)
#ft = k2*sym.exp(-2*t)*u.subs(t,t-5)
#ft = d.subs(t,t-2)
#ft = d
#ft = 3*sym.exp(-2*t)*u + sym.exp(-t)*sym.cos(3*t)*u
#ft = u
#ft = u - u.subs(t,t-2)
#ft = u.subs(t,t-2)

# PROCEDIMIENTO
def laplace_transform_suma(ft):
    '''transformada de Laplace de suma de terminos
       separa constantes para conservar en resultado
    '''
    def separa_constante(termino):
        ''' separa constante antes de usar
            sym.laplace_transform(term_suma,t,s)
            para incorporarla luego de la transformada
            inconveniente revisado en version 1.11.1
        '''
        constante = 1
        if termino.is_Mul:
            factor_mul = sym.Mul.make_args(termino)
            for factor_k in factor_mul:
                if not(factor_k.has(t)):
                    constante = constante*factor_k
            termino = termino/constante
        return([termino,constante])
    
    # transformadas de Laplace por términos suma
    ft = sym.expand(ft)  # expresion de sumas
    ft = sym.powsimp(ft) # simplifica exponentes

    term_suma = sym.Add.make_args(ft)
    Fs = 0
    for term_k in term_suma:
        [term_k,constante] = separa_constante(term_k)
        Fsk = sym.laplace_transform(term_k,t,s)
        Fs  = Fs + Fsk[0]*constante
    
    # separa exponenciales constantes
    Fs = sym.expand_power_exp(Fs)
    return(Fs)

Fs = laplace_transform_suma(ft)

# SALIDA
print('\n f(t): ')
sym.pprint(ft)
print('\n F(s): ')
sym.pprint(Fs)

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ]  [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]
..


Ejemplo 5 f(t) con Impulsos unitarios desplazados

Referencia:  Lathi Ej 4.9c p355

Considera la entrada x(t) como una suma de impulsos desplazados en tiempo y de diferente magnitud.

x(t) = \delta (t) - 3 \delta (t-2) + 2 \delta (t-3)

para el algoritmo del ejercicio 4, se modifica la línea de ingreso a:

ft = d - 3*d.subs(t,t-2) + 2*d.subs(t,t-3)

obteniendo como resultado del algoritmo anterior:

 f(t): 
DiracDelta(t) + 2*DiracDelta(t - 3) - 3*DiracDelta(t - 2)

 F(s): 
       -2*s      -3*s
1 - 3*e     + 2*e    
>>> 

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 escalón unitario ] [ ej3 con cos(t) ] [ ej4 impulso unitario y sumas]  [ ej5  impulso unitario desplazado ]


4.2.1 Transformada de Laplace – Integral con Sympy-Python

El algoritmo con Python se desarrolla progresivamente a partir de los ejercicios resueltos de forma analítica, usando las librerias de Sympy-Python. En caso de requerir usar el algoritmo completo,  puede ir directamente al ejercicio 6 que incluye impulsos unitarios y considera el integral unilateral o completo en s.

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]
..


Ejemplo 1. Transformada de Laplace de una exponencial decreciente, un solo termino

Referencia: Lathi ejemplo 4.1. p331, Oppenheim Ejemplo 9.2 p656, Hsu Ejemplo 3.1 p111 Transformada Laplace f(t) Ej01

Para una señal x(t) = e-atμ(t), encuentre la transformada X(s) y su región de convergencia (ROC).

x(t) = e^{-at} μ(t)

Considerando que μ(t)=0 para t<0 y μ(t)=1 para t≥0:

X(s) = \int_{0}^{\infty} e^{-at} e^{-st} \delta t

X(s) = \frac{1}{s+a} Transformada Laplace Fs Ej01

Se puede observar que la región de convergencia, ROC, de X(s) es Re(s)>-a, como se muestra en la gráfica. Para otros valores de s el integral no converge.

Para el ejemplo, siendo a=2, se obtiene con el algoritmo el resultado y las gráficas presentadas.

 dentro de integral f(t)*e(-st):
 -t*(a + s)             
e          *Heaviside(t)

 expresion F(s):
  1  
-----
a + s

 {Q_polos:veces}:  {-a: 1}
 {P_ceros:veces}:  {}

Integral de la Transformada de Laplace con Sympy-Python

Para iniciar se empieza usando una expresión simple como la del ejemplo. Semejante a los algoritmos anteriores realizados, se definen las variables tipo símbolo a usar y el escalón unitario en una expresión ‘u’ para simplificar su escritura posterior.

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

El desarrollo de los integrales se facilita al simplificar la expresión por términos de suma y simplificando los exponentes con la instrucción  sym.powsimp().

>>> fts = ft*sym.exp(-s*t)
>>> fts
exp(-a*t)*exp(-s*t)*Heaviside(t)
>>> fts = sym.powsimp(fts)
>>> fts
exp(-t*(a + s))*Heaviside(t)
>>>

Las instrucciones de las operaciones para el algoritmo serán:

# PROCEDIMIENTO
fts = ft*sym.exp(-s*t) # f(t,s) para integrar
fts = sym.expand(fts)  # expresion de sumas
fts = sym.powsimp(fts) # simplifica exponentes

La integración a realizar es unilateral al considerar principalmente señales causales como se justifica en los textos, por lo que el intervalo de integración es [0,∞].

# integral Laplace unilateral
Fs_L = sym.integrate(fts,(t,0,sym.oo))

El resultado del integral se presenta por por partes o intervalos, de la forma  sym.piecewise() mostrada por facilidad de lectura con sym.pprint(Fs_L).

Fs_L:

/               1                                    pi
|             -----               for |arg(a + s)| < --
|             a + s                                  2 
|                                                      
| oo                                                   
|  /                                                   
< |                                                    
| |   -t*(a + s)                                       
| |  e          *Heaviside(t) dt        otherwise      
| |                                                    
|/                                                     
|0                                                     
\                                

La parte a usar para el desarrollo del ejercicio corresponde la primera expresión, Fs = Fs_L.args[0], que contiene expresión e intervalo. Para quedarse solo con la expresión y omitir el intervalo, se usa Fs[0]

>>> Fs = Fs_L.args[0]
>>> Fs
(1/(a + s), Abs(arg(a + s)) < pi/2)
>>> Fs = Fs[0]
>>> Fs
1/(a + s)
>>>

las instrucciones para el algoritmo serán:

Fs = Fs_L.args[0] # primera ecuacion e intervalo
Fs = Fs[0]        # solo expresion

El comportamiento F(s), puede presentar divisiones para cero para algun valor de s llamados polos, asi como valores de s que resultan en ceros. Se usa la expresión como una fracción de polinomios P(s) y Q(s), tomando el denominador Q(s) para obtener las raíces del polinomio en Q_polos.

# polos y ceros en Fs
[P,Q] = Fs.as_numer_denom()
P = sym.poly(P,s)
Q = sym.poly(Q,s)
P_ceros = sym.roots(P)
Q_polos = sym.roots(Q)

Instrucciones en Python – Ejemplo 1

# Transformada de Laplace
# integral de Laplace unilateral
# blog.espol.edu.ec/telg1001/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

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

ft = sym.exp(-a*t)*u

a_k = 2 # valor de 'a' constante
# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras = 101  # 51 resolucion grafica

# PROCEDIMIENTO
fts = ft*sym.exp(-s*t) # f(t,s) para integrar
fts = sym.expand(fts)  # expresion de sumas
fts = sym.powsimp(fts) # simplifica exponentes

# integral Laplace unilateral
Fs_L = sym.integrate(fts,(t,0,sym.oo))
Fs = Fs_L.args[0] # primera ecuacion e intervalo
Fs = Fs[0]        # solo expresion

# polos y ceros en Fs
[P,Q] = Fs.as_numer_denom()
P = sym.poly(P,s)
Q = sym.poly(Q,s)
P_ceros = sym.roots(P)
Q_polos = sym.roots(Q)

# SALIDA
print(' dentro de integral f(t)*e(-st):')
sym.pprint(fts)
print('\n expresion F(s):')
sym.pprint(Fs)
print('\n {Q_polos:veces}: ',Q_polos)
print(' {P_ceros:veces}: ',P_ceros)

Instrucciones para gráficar la transformada de Laplace

Se sustituye en las ecuaciones resultantes los valores de ‘a‘ con ‘a_k‘  y usando lambdify, aplicarla en en un intervalo [sa , sb] que permita mostrar una gráfica para F(s). Ejemplos de graficas se encuentran en la unidad 1.

El intervalo [sa , sb] considera los extremos de las raices o polos del denominador Q.  Los polos se marcan usando puntos de dispersión con plt.scatter() y una línea vertical plt.axvline() para resaltar en la gráfica el efecto del polo.

# GRAFICA -----------
ft = ft.subs(a,a_k) # a tiene valor a_k
Fs = Fs.subs(a,a_k)

# f(t) # evalua en intervalo
ti  = np.linspace(t_a,t_b,muestras)
f_t = sym.lambdify(t,ft)
fti = f_t(ti)

# Q polinomio denominador con valor de a
# polos en s
[P,Q] = Fs.as_numer_denom()
Q = Q.as_poly(s)
Q_polos = sym.roots(Q)

# estima intervalo para s, raices reales
s_a = 0 ; s_b = 0
polos = list(Q_polos.keys())
s_a = int(min(s_a,min(polos)))-1
s_b = int(max(s_b,max(polos)))+1

# F(s) # evalua en intervalo
F_s = sym.lambdify(s,Fs)
s_i = np.linspace(s_a,s_b,muestras)
Fsi = F_s(s_i) # Revisar cuando s es complejo

# grafica f(t) 
fig_ft, graf_ft = plt.subplots()
plt.plot(ti,fti,label='f(t)')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('f(t) = '+str(ft)+' ; a='+str(a_k))
plt.grid()

# grafica F(s) , corte en plano real
fig_Fs, graf_Fs = plt.subplots()
plt.plot(s_i,Fsi, color='green',label='F(s)')
for raiz in Q_polos.keys():    
    plt.axvline(sym.re(raiz),color='red')
    plt.scatter(sym.re(raiz),sym.im(raiz),
                label='polo:'+str(raiz),
                marker='x',color='red')
plt.axvline(0,color='gray')
plt.legend()
plt.xlabel('s')
plt.ylabel('F(s)')
plt.title('F(s) = '+str(Fs)+' ; a='+str(a_k))
plt.grid()
plt.show()

Se recomienda probar con otros ejemplos desarrollados para la parte analítica y de ser necesario en cada caso mejorar el algoritmo.

Considere también que Sympy dispone de la librerías para transformadas de Laplace con instrucciones simplificadas que se usarán a lo largo de ésta unidad de estudio y que se desarrollan en la siguiente sección.

Fs_L = sym.laplace_transform(ft,t,s)

Transformada de Laplace e Inversas con Sympy-Python

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]

..


Ejemplo 2. Transformada de Laplace con suma de exponenciales

Transformada Laplace f(t) Ej02

Referencia: Oppenheim Ejemplo 9.3 p658

Considere la señal que es la suma de dos exponenciales:

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

Del desarrollo analítico del ejercicio y usando el algoritmo se tiene que:

X(s) = \frac{3}{s+2}-\frac{2}{s+1}

Para observar mejor los polos y ceros de X(s) se simplifica la expresión como factores en el numerador y denominador:

X(s) = \frac{s-1}{(s+2)(s+1)}

La gráfica de F(s) muestra la ubicación de los polos, observando que en estos puntos la función tiende a crecer a infinito positivo o negativo,

Transformada Laplace F(s) Ej02

Para el algoritmo, la función de entrada f(t) se expresa como:

u = sym.Heaviside(t)
ft = 3*sym.exp(-2*t)*u - 2*sym.exp(-t)*u

Usando el algoritmo, añadiendo las partes que consideren el integral de Laplace para cada uno de los términos suma luego agrupando por factores se obtiene:

 dentro de integral f(t)*e(-st):
   -t*(s + 2)                   -t*(s + 1)             
3*e          *Heaviside(t) - 2*e          *Heaviside(t)

 expresion F(s):
     s - 1     
---------------
(s + 1)*(s + 2)

Q_polos{polos:veces}:  {-1: 1, -2: 1}
P_ceros{polos:veces}:  {1: 1}
>>>
X(s) = \frac{3}{s+2}-\frac{2}{s+1}

Considere que en la expresión X(s) no será necesario disponer de la constante ‘a‘, usada para generalizar la respuesta en el ejemplo 1. Por lo que se podría prescindir de las líneas que le hacen referencia en las instrucciones  para la revisión en los polos y los títulos de las gráficas, resumidas como:

# POR COMENTAR O ELIMINAR, si no se usa constante a:
# a = sym.Symbol('a', real=True)

# Para graficar con s
# ft = ft.subs(a,a_k) # a tiene valor a_k
# Fs = Fs.subs(a,a_k)

# Eliminar la 2da evaluacion de Q_polos

# plt.title('f(t) = '+str(ft))
# plt.title('F(s) = '+str(Fs))

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]
..


Ejemplo 3. Transformada de Laplace para suma de términos f(t) con desplazamiento, o función «gate» o compuerta

Referencia: Lathi práctica 4.1.a p337

literal a. Por integración directa, encuentra la transformada X(s) y
Transformada Laplace f(t) Ej03la región de convergencia para la función descrita en la imagen

Analizando la expresión de la función es:

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

y la transformada unilateral de Laplace se expresa como:

X(s) = \frac{1}{s} - \frac{e^{-2s}}{s}

La gráfica de F(s) muestra el resultado de la suma de términos o componentes. La suma de éstos términos hace que se minimice alrededor de cero el efecto creciente de cada componente.

Transformada Laplace F(s) Ej03

Desarrollo con Sympy-Python

Para el caso presentado, la integral se compone de dos términos de una suma, procediendo de forma semejante al ejercicio anterior.

El bloque de ingreso se expresa como:

u = sym.Heaviside(t)
ft = u - u.subs(t,t-2)

El resultado del algoritmo será:

 dentro de integral f(t)*e(-st):
 -s*t                 -s*t                 
e    *Heaviside(t) - e    *Heaviside(t - 2)

 expresion F(s):
     -2*s
1   e    
- - -----
s     s  

Polos y ceros:
término: 1 * 1/s
 Q_polos{polos:veces}:  {0: 1}
 P_ceros{polos:veces}:  {}
término: exp(-2*s) * -1/s
 Q_polos{polos:veces}:  {0: 1}
 P_ceros{polos:veces}:  {}
>>>

Para el caso de polos y ceros el  término con e-s no permite analizar la expresión de forma directa como un polinomio y obtener las raíces. Con la restricción descrita, el análisis se realiza agrupando los términos por cada  e-s llamado componente de suma. Las veces que se repite cada polo será el mayor valor entre cada termino suma.

Instrucciones adicionales para términos con sym.exp(-a*s)

Las operaciones para obtener los polos y ceros en Sympy se realizan como polinomios sym.poly(Q,s), sin embargo ésto no considera o permite usar la expresión de F(s) con términos sym.exp(-a*s), siendo a un valor numérico.

Se obtiene una lista de todos los términos con exponencial de la forma sym.exp(-a*s), al obtener los componentes básicos de toda la expresión con Fs.atoms(sym.exp(-s)).

# separa en lista los sym.exp(-s)
lista_exp = list(Fs.atoms(sym.exp(-s)))

para observar la operación se muestra el resultado de la instrucción sobre F(s)

>>> Fs
1/s - exp(-2*s)/s
>>> Fs.atoms(sym.exp(-s))
{exp(-2*s)}
>>>

Luego se agrupa F(s) por cada elemento de la lista_exp.

# agrupa por cada exp(-s) en lista_exp
if len(lista_exp)>0:
    Fs = sym.expand(Fs,s)
    Fs = sym.collect(Fs,lista_exp)

Los polos y ceros se obtienen para cada grupo de lista_exp(), empezando por aquellos que no tienen sym.exp(-a*s). Los resultado parciales se almacenan en un diccionario polosceros con el objetivo de mostrar lo que sucede con cada componente y su influencia en F(s).

El análisis detallado de los polos y ceros para la expresión se desarrolla en la sección LTI CT Laplace – H(s) Polos reales y complejos con Sympy-Python . Por ahora para mantener el enfoque sobre el desarrollo del integral se usa una función fcnm.busca_polosceros(Fs) obtenida desde telg1001.py que entrega los polos y ceros para cada componente de sym.exp(-a*s).

Instrucciones con Python – Ejemplos 2, 3, 4 y 5

# Transformada de Laplace Ejemplo 2 y 3
# integral de Laplace unilateral con suma de terminos,
#  y desplazamientos en tiempo
# blog.espol.edu.ec/telg1001/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1001 as fcnm

# INGRESO
t = sym.Symbol('t', real=True)
s = sym.Symbol('s')
u = sym.Heaviside(t)

# ft = sym.exp(-2*t)*u
# ft = 3*sym.exp(-2*t)*u - 2*sym.exp(-t)*u
# ft = u
# ft = u.subs(t,t-2)
ft = u - u.subs(t,t-2)
# ft = u.subs(t,t-2) - u.subs(t,t-4)
# ft = sym.exp(-2*t)*u + sym.exp(-t)*sym.cos(3*t)*u

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 10
muestras_t = 101  # 51 resolucion grafica

# PROCEDIMIENTO
fts = ft*sym.exp(-s*t) # f(t,s) para integrar
fts = sym.expand(fts)  # expresion de sumas
fts = sym.powsimp(fts) # simplifica exponentes

# Integral Laplace 
unilateral = True
lim_a = 0 # unilateral
if not(unilateral):
    lim_a = -sym.oo
Fs_L = sym.integrate(fts,(t,lim_a,sym.oo))
Fs = Fs_L.args[0] # primera ecuacion e intervalo
Fs = Fs[0]        # solo expresion

# convierte a Sympy si es solo constante
Fs = sym.sympify(Fs)

# separa en lista los sym.exp(-s)
lista_exp = list(Fs.atoms(sym.exp(-s)))

# agrupa por cada exp(-s) en lista_exp
if len(lista_exp)>0:
    Fs = sym.expand(Fs,s)
    Fs = sym.collect(Fs,lista_exp)
else:
    Fs = sym.factor(Fs,s)

polosceros = fcnm.busca_polosceros(Fs)

# SALIDA
print(' dentro de integral f(t)*e(-st):')
sym.pprint(fts)
print('\n expresion F(s):')
sym.pprint(Fs)
print('\n {Q_polos:veces}: ',polosceros['Q_polos'])
print(' {P_ceros:veces}: ',polosceros['P_ceros'])

# GRAFICA ------------------
# grafica de f(t)
fig_ft = fcnm.graficar_ft(ft,t_a,t_b,
                            muestras = muestras_t,
                            f_nombre='f')
# grafica de polos y ceros
fig_Fs = fcnm.graficar_Fs(Fs,polosceros['Q_polos'],
                          polosceros['P_ceros'],
                          solopolos=False))
plt.show()

Las instrucciones para las gráficas se pueden simplificar para f(t), usando la instrucción de la función fcnm.graficar_ft() desarrollada en la unidad 3 desde 3.2.1 LTI CT y disponible en telg1001.py.

Para analizar en un gráfico el caso de F(s) en detalle, se requiere realizar las líneas de cada componente que tiene un exp() y los polos para cada componente. Estas instrucciones se desarrollan en la sección LTI CT Laplace – H(s) Polos reales y complejos con Sympy-Python por ahora se usará la función gráfica en telg1001.py

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]
..


Ejemplo 4. Transformada de Laplace para suma de términos f(t) con desplazamiento, función «gate» o compuerta causal

Referencia: Lathi práctica 4.1.b p337

literal b. Por integración directa, encuentra la transformada X(s) y la región de convergencia para la función descrita en la imagen

x(t) = \mu (t-2) - \mu (t-4)

Transformada Laplace f(t) Ej04

Con el desarrollo analítico realizado en la pagina anterior se comprueba que el resultado de X(s) es

X(s) = \frac{e^{-2s}}{s}-\frac{e^{-4s}}{s}

Transformada Laplace F(s) Ej04

Desarrollo con Sympy-Python

El bloque de ingreso del ejercicio se actualiza como la suma de escalón unitario μ(t) desplazados,

u = sym.Heaviside(t)
ft = u.subs(t,t-2)-u.subs(t,t-4)

obteniendo con el algoritmo el resultado y gráfica correspondiente.

 dentro de integral f(t)*e(-st):
   -s*t                     -s*t                 
- e    *Heaviside(t - 4) + e    *Heaviside(t - 2)

 expresion F(s):
 -2*s    -4*s
e       e    
----- - -----
  s       s  

Polos y ceros:
término: exp(-2*s) * 1/s
 Q_polos{polos:veces}:  {0: 1}
 P_ceros{polos:veces}:  {}
término: exp(-4*s) * -1/s
 Q_polos{polos:veces}:  {0: 1}
 P_ceros{polos:veces}:  {}  

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]
..


Ejemplo 5. Transformada de Laplace para cos(t)

Referencia: Oppenheim Ejemplo 9.4 p658

Encontrar la transformada de Laplace para:

x(t) = e^{-2t}\mu (t) + e^{-t} \cos (3t) \mu (t)

Transformada Laplace ft Ej05
La expresión para el algoritmo de la señal de entrada es:

ft = sym.exp(-2*t)*u + sym.exp(-t)*sym.cos(3*t)*u

Transformada Laplace F(s) Ej05

el resultado con el algoritmo desarrollado hasta el ejemplo es:

 dentro de integral f(t)*e(-st):
 -t*(s + 2)                 -t*(s + 1)                      
e          *Heaviside(t) + e          *cos(3*t)*Heaviside(t)

 expresion F(s):
       2               
    2*s  + 5*s + 12    
-----------------------
        / 2           \
(s + 2)*\s  + 2*s + 10/

 {Q_polos:veces}:  {-2: 1, -1 - 3*I: 1, -1 + 3*I: 1}
 {P_ceros:veces}:  {-5/4 - sqrt(71)*I/4: 1, -5/4 + sqrt(71)*I/4: 1}

La expresión F(s)  de factores, también puede expresarse como la suma de fracciones parciales con la instrucción sym.apart() que aplica solo a polinomios. En el caso de tener términos con exp(), la instrucción sym.apart() se puede aplicar por cada grupo.

 expresion F(s):
   s + 1         1  
------------ + -----
       2       s + 2
(s + 1)  + 9        

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]
..


Ejemplo 6. Transformada de Laplace con Impulso unitario δ(t) y suma de exponenciales

Referencia: Oppenheim ejemplo 9.5 p661

x(t) = \delta(t) -\frac{4}{3} e^{-t} \mu (t) + \frac{1}{3} e^{2t} \mu (t)

Los resultados de la parte analítica y con algoritmo del ejercicio se muestran en la gráfica.

Transformada Laplace f(t) Ej06
La introducción de un impulso d(t) en el integral de varios componentes ya comienza a generar algun trabajo adicional en la instrucción simbólica sym.integrate(), por lo que es preferible facilitar el desarrollo por cada componente de términos suma. Al entregar a sym.integrate() términos más simples, la función no tiene que analizar varias formas de integración de la expresión, el proceso se acelera y se minimizan errores.

Transformada Laplace F(s) Ej06

Desarrollo y ajustes en el Algoritmo

Inicialmente considere usar como f(t) = δ(t) y revisar el resultado Fs_L, encontrando que es solo una constante,  el resultado de F(s) no es por partes (sym.Piecewise()) como se había realizado hasta ahora en el algoritmo.

Si Fs_L aún es por partes, se puede verificar con la instrucción Fs_L.is_Piecewise que indica True o False para aplicar o no Fs_L[0].

Debe considerar también que el integral unilateral de Laplace presenta el resultado 1/2, en lugar de 1. Revisar las definiciones sobre la función impulso δ(t)  en el enlace DiracDelta en Sympy. Para mantener concordancia con lo desarrollado en los textos, se corrige el resultado a 1.

ft = sym.DiracDelta(t)

dentro de integral f(t)*e(-st):
 -s*t              
e    *DiracDelta(t)

 expresion F(s):
1/2

 {polos:veces}:  {}

Para la convención usada en los libros del curso, se revisará si el término al que se aplica el integral de Laplace tiene un impulso para proceder a multiplicarlo por 2. Si el integral fuese de [-∞,∞] el integral tendria resultado 1, pero en lo realizado es solo unilateral (t≥0) por lo que Sympy se ha aplicado básicamente solo la mitad. Discusiones sobre aquello las puede revisar en el sitio web de Sympy.

unilateral = True

# Integral del impulso en cero es 1
if fts.has(sym.DiracDelta): 
    donde = fcnm.busca_impulso(fts)
    if (0 in donde) and unilateral: # (integral unilateral) x 2
        fts = fts.subs(d,2*d)
    fts = sym.powsimp(fts) # simplifica exponentes

con la variable unilateral también se establece el valor del límite inferior del integral:

    # Integral Laplace de sumas
    lim_a = 0 # unilateral
    if not(unilateral):
        lim_a = -sym.oo

con lo que se actualiza el resultado anterior a:

ft = sym.DiracDelta(t)

dentro de integral f(t)*e(-st):
 -s*t              
e    *DiracDelta(t)

 expresion F(s):
1

 {polos:veces}:  {}

Nota: Sympy hasta la versión 1.11.1, las operaciones en el dominio ‘s’ para la Transformadas Inversas de Laplace se encuentran implementadas para manejar principalmente números enteros y fracciones. Los resultados de expresiones combinadas con coeficientes enteros y coeficientes reales no necesariamente se simplifican entre si, pues se manejan diferentes dominios ‘ZZ’ o ‘QQ’. (Revisión 2022-Nov)

Para optimizar la simplificación de expresiones con coeficientes entre enteros y reales, los números reales se convierten a su aproximación racional con la instrucción sym.Rational(0.333333).limit_denominator(100).

Convirtiendo los coeficientes a racionales como:

k1 = sym.Rational(1/3).limit_denominator(100)
k2 = sym.Rational(4/3).limit_denominator(100)

ft = d - k2*sym.exp(-t)*u + k1*sym.exp(2*t)*u

Uniendo las observaciones, el algoritmo más general se muestra a continuación

Instrucciones en Python Ejercicio 6. Transformada de Laplace con impulso δ(t) y suma de exponenciales

# Transformada de Laplace Ejemplo 6
# integral de Laplace unilateral con suma de terminos,
#  y desplazamientos en tiempo
# blog.espol.edu.ec/telg1001/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1001 as fcnm

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

# ft = u
# ft = u.subs(t,t-2)
# ft = u - u.subs(t,t-2)
# ft = 3*sym.exp(-2*t)*u-2*sym.exp(-t)*u
# ft = u.subs(t,t-2) - u.subs(t,t-4)
# ft = sym.exp(-2*t)*u + sym.exp(-t)*sym.cos(3*t)*u
# ft = d

# coeficientes como racional en dominio 'ZZ' enteros
k1 = sym.Rational(1/3).limit_denominator(100)
k2 = sym.Rational(4/3).limit_denominator(100)

ft = d - k2*sym.exp(-t)*u + k1*sym.exp(2*t)*u
#ft = u.subs(t,t+1)

# Grafica, intervalo tiempo [t_a,t_b]
t_a = -1 ; t_b = 4
muestras_t = 101  # 51 resolucion grafica

# PROCEDIMIENTO
def laplace_integral_sumas(ft,unilateral=True):
    ''' integral de transformada de laplace
        considera impulsos d(t), escalon u(t) desplazados
    '''
    fts = ft*sym.exp(-s*t) # f(t,s) para integrar
    fts = sym.expand(fts)  # expresion de sumas
    fts = sym.powsimp(fts) # simplifica exponentes

    # Integral del impulso en cero es 1
    if fts.has(sym.DiracDelta): 
        donde = fcnm.busca_impulso(fts)
        if (0 in donde) and unilateral: # (integral unilateral) x 2
            fts = fts.subs(d,2*d)
        fts = sym.powsimp(fts) # simplifica exponentes

    # Integral Laplace de sumas
    lim_a = 0 # unilateral
    if not(unilateral):
        lim_a = -sym.oo
    Fs = 0*s
    term_suma = sym.Add.make_args(fts)
    for term_k in term_suma:
        # integral Laplace unilateral
        Fs_L = sym.integrate(term_k,(t,lim_a,sym.oo))
        if Fs_L.is_Piecewise:   # Fs_L es por partes
            Fsk = Fs_L.args[0]  # primera ecuacion e intervalo
            Fsk = Fsk[0]        # solo expresion
        else: # una sola expresión
            Fsk = Fs_L   
        Fs = Fs + Fsk

    # convierte a Sympy si es solo constante
    Fs = sym.sympify(Fs)
    Fs = sym.expand_power_exp(Fs)
    # lista los sym.exp(-s)
    lista_exp = list(Fs.atoms(sym.exp(-s)))

    # agrupa por cada exp(-s) en lista_exp
    if len(lista_exp)>0:
        Fs = sym.expand(Fs,s)
        Fs = sym.collect(Fs,lista_exp)
    else:
        Fs = sym.factor(Fs,s)
    return(Fs)

Fs = laplace_integral_sumas(ft, unilateral=True)
polosceros = fcnm.busca_polosceros(Fs)

# SALIDA
print('\n expresion F(s):')
sym.pprint(Fs)
print('\n {Q_polos:veces}: ',polosceros['Q_polos'])
print(' {P_ceros:veces}: ',polosceros['P_ceros'])

La sección de gráficas ya incorpora desde la unidad 3 el caso que f(t) contiene impulsos unitarios, por lo que las instrucciones se mantienen.

# GRAFICA ------------------
#grafica de f(t)
fig_ft = fcnm.graficar_ft(ft,t_a,t_b,
                            muestras = muestras_t,
                            f_nombre='f')
# grafica de polos y ceros
fig_Fs = fcnm.graficar_Fs(Fs,polosceros['Q_polos'],
                          polosceros['P_ceros'],
                          solopolos=False)
plt.show()

Presentando así, un algoritmo  general para desarrollar integrales de Laplace.

Desde luego existe la forma mas simple de usar con las instrucciones Fs_L = sym.laplace_transform(ft,t,s) mencionadas al inicio de la página. El algoritmo es un ejercicio didáctico de cómo realizar la tabla de pares de transformadas.

En la siguiente sección el algoritmo se simplifica y se muestra el uso de las transformadas de Laplace usando las librerías de Sympy.

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]

..


Ejemplo 7. Transformada de Laplace de Impulsos unitarios desplazados

Referencia:  Lathi Ej 4.9c p355

Considera la entrada x(t) como una suma de impulsos desplazados en tiempo y de diferente magnitud.

x(t) = \delta (t) - 3 \delta (t-2) + 2 \delta (t-3)

para el algoritmo del ejercicio 6, se modifica la línea de ingreso a:

ft = d - 3*d.subs(t,t-2) + 2*d.subs(t,t-3)

obteniendo como resultado del algoritmo anterior:

X(s) = 1 - 3 e^{-2s} + 2 e^{-3s}
 dentro de integral f(t)*e(-st):
 -s*t                    -s*t                        -s*t                  
e    *DiracDelta(t) + 2*e    *DiracDelta(t - 3) - 3*e    *DiracDelta(t - 2)

 expresion F(s):
       -2*s      -3*s
1 - 3*e     + 2*e    

 {Q_polos:veces}:  {}

con las siguientes gráficas:

Transformada Laplace Ej07

y para F(s):

Transformada Laplace Ej07Polos

Transformada de Laplace para: [ ej1 un exponencial ]  [ ej2 suma de términos ]  [ ej3 escalón desplazado ]  [ ej4 desplazados ]  [ ej5 coseno ]  [ ej6 impulso ]  [ ej7 impulsos desplazados ]

4.2 Transformada de Laplace – Integral ejemplos

Referencia: Lathi 4.1. p330. Hsu 3.2.A p110, Oppenheim 9.1 p655

La transformada de Laplace permite simplificar el proceso de solución de ecuaciones integro-diferenciales usando operaciones mas simples al cambiar desde el dominio del tiempo ‘t’ al dominio ‘s’.

Para una señal contínua x(t), la transformada de Laplace esta definida como:

X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} dt

Para el caso de señales contínuas, lineales y causales, se define una señal x(t) que tiene un componente escalón unitario μ(t), por lo que la integral se desarrolla de forma unilateral. Se enfatiza que se entiende como transformada de Laplace unilateral si cada señal x(t) es cero para t<0, y es apropiado indicarlo al multiplicar la señal por el escalón unitario μ(t)

X(s) = \int_{-\infty}^{\infty} x(t) \mu(t) e^{-st} dt = \int_{0}^{\infty} x(t) \mu(t) e^{-st} dt

Para la transformada de Laplace unilateral, existe una transformada inversa de X(s) que es única. En consecuencia, no hay necesidad de especificar la región de convergencia (ROC) de forma explícita. Motivo por el que generalmente no se menciona la ROC para transformadas unilaterales. (Lathi p337).

La señal x(t) es la inversa de la transformada X(s), que se obtiene de la forma:

x(t) = \frac{1}{2πj} \int_{c-j\infty}^{c+j\infty} X(s) e^{st} dt

donde c es una constante seleccionada para asegurar la convergencia de la integral.

Los pares de ecuaciones conocidos como «Pares de la transformada de Laplace» se escriben de forma simbólica:

X(s) \Rightarrow \mathscr{L}[x(t)] x(t) \Rightarrow \mathscr{L}^{-1}[X(s)]

que tienen algunas propiedades de interés para señales y sistemas como la linealidad,

\mathscr{L}[a_1 x_1(t) + a_2 x_2 (t)] = a_1 X_1(s) + a_2 X_2 (s)

Para el desarrollo de ejercicios se usa principalmente la tabla de pares de Transformadas de Laplace, los pares se obtienen al aplicar la definición del integral. El desarrollo el integral se puede también realizar usando Sympy-Python.

Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ]  [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]
..


Ejemplo 1. Integral de la Transformada de Laplace de una exponencial decreciente, un solo termino

Referencia: Lathi ejemplo 4.1. p331, Oppenheim Ejemplo 9.2 p656, Hsu Ejemplo 3.1 p111

Para una señal x(t) = e-atμ(t), encuentre la transformada X(s) y su región de convergencia (ROC)

x(t) = e^{-at} μ(t)

Transformada Laplace f(t) Ej01

X(s) = \int_{-\infty}^{\infty} e^{-at} \mu (t) e^{-st} \delta t

pero tomando en cuenta que μ(t)=0 para t<0 y μ(t)=1 para t≥0:

X(s) = \int_{0}^{\infty} e^{-at} e^{-st} \delta t = \int_{0}^{\infty} e^{-(s+a)t} \delta t = -\frac{1}{s+a} e^{-(s+a)t}\Big|_{0}^{\infty} = \Big[-\frac{1}{s+a} e^{-(s+a)(\infty)}\Big]-\Big[ -\frac{1}{s+a} e^{-(s+a)(0)}\Big] X(s) = \frac{1}{s+a}

con polo s=-a al producir división para cero en la expresión.

Se puede observar que la región de convergencia, ROC, de X(s) es Re(s)>-a, como se muestra en la gráfica. Para otros valores de s el integral no converge.

Para el ejemplo, siendo a=2

Transformada Laplace Fs Ej01


El siguiente video presenta una interpretación gráfica y animada de la transformada de Laplace en el plano s real e imaginario.


Propiedades de la transformada de Laplace

Otra de las propiedades de interés es la diferenciación. Por ejemplo, los sistemas formados por circuitos electrónicos que tienen inductores L, usan la variación de corriente o derivada de la corriente, generando ecuaciones diferenciales ordinarias. La transformada de Laplace presenta una alternativa viable para tratar estos circuitos como se mostrará en los ejemplos de la Unidad.

Propiedad de diferenciación

\frac{\delta x}{\delta t} \Leftrightarrow sX(s) - x(0^{-}) \frac{d^2x}{dt^2} \Leftrightarrow s^2X(s) - sx(0^{-}) - x'(0^{-}) \frac{\delta^3x}{\delta t^3} \Leftrightarrow s^3 X(s) - s^2 x(0^{-}) - sx'(0^{-}) - x''(0^{-})

Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ] [ ej3 escalón Desplazado ] [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]

..


Ejemplo 2. Transformada de Laplace con suma de exponenciales

Referencia: Oppenheim Ejemplo 9.3 p658

Considere la señal que es la suma de dos exponenciales:

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

Transformada Laplace f(t) Ej02

se tiene que:

X(s) = \int_{0}^{\infty} \Big[ 3 e^{-2t} \mu (t) - 2e^{-t} \mu (t) \Big] e^{-st} \delta t X(s) = \int_{0}^{\infty} \Big[ 3 e^{-2t}e^{-st}\mu (t) - 2e^{-t}e^{-st}\mu (t) \Big] \delta t X(s) = \int_{0}^{\infty} \Big[ 3 e^{-(s+2)t}\mu (t) - 2e^{-(s+1)t} \mu (t) \Big] \delta t

En general para tratar este tipo de ejercicios es mejor descomponer la señal o función matemática en varios componentes de suma, siguiendo la propiedad de linealidad de los sistemas.

X(s) = 3 \int_{0}^{\infty} e^{-(s+2)t}\mu (t) \delta t - 2 \int_{0}^{\infty} e^{-(s+1)t} \mu (t) \delta t X(s) = -3 \frac{1}{s+2}e^{-(s+2)t}\Big|_{0}^{\infty} + 2 \frac{1}{s+1}e^{-(s+1)t}\Big|_{0}^{\infty} X(s) = -3 \frac{1}{s+2}\Big( e^{-(s+2)(\infty)} - e^{-(s+2)(0)} \Big) + 2 \frac{1}{s+1} \Big( e^{-(s+1)(\infty)} - e^{-(s+1)(0)} \Big) X(s) = -3 \frac{1}{s+2} \Big( 0-1 \Big) + 2 \frac{1}{s+1} \Big( 0 - 1 \Big) X(s) = 3 \frac{1}{s+2} - 2 \frac{1}{s+1}

con polos en s=-1 y s=-2 al producir división para cero en la expresión

Transformada Laplace F(s) Ej02

Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ] [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]
..


Ejemplo 3. Transformada de Laplace para suma de términos f(t) con desplazamiento, o función «gate» o compuerta

Referencia: Lathi práctica 4.1.a p337

literal a. Por integración directa, encuentra la transformada X(s) y la región de convergencia para la función descrita en la imagen

Transformada Laplace f(t) Ej03
Analizando la expresión de la función es:

x(t) = \mu (t) - \mu (t-2) X(s) = \int_{0}^{\infty} \Big[ \mu (t) - \mu (t-2) \Big] e^{-st} \delta t X(s) = \int_{0}^{\infty} \Big[ \mu (t) e^{-st} - \mu (t-2) e^{-st} \Big] \delta t X(s) = \int_{0}^{\infty} \mu (t) e^{-st} \delta t - \int_{2}^{\infty} \mu (t) e^{-st} \delta t X(s) = -\frac{1}{s} e^{-st} \Big|_{0}^{\infty} +\frac{1}{s} e^{-st} \Big|_{2}^{\infty} X(s) = -\frac{1}{s} \Big( e^{-s(\infty)} -e^{-s(0)}\Big) +\frac{1}{s} \Big( e^{-s(\infty)} -e^{-s(2)} \Big) X(s) = -\frac{1}{s} \Big( 0 - 1 \Big) +\frac{1}{s} \Big( 0 - e^{-2s} \Big) X(s) = \frac{1}{s} -\frac{1}{s} e^{-2s}

con polo en s=0

se observa también que el retraso de tiempo aplicado en μ(t-2) genera en la transformada un termino multiplicado por e(-2s).

Transformada Laplace F(s) Ej03
Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ]  [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]
..


Ejemplo 4. Transformada de Laplace para suma de términos f(t) con desplazamiento, función «gate» o compuerta causal

Referencia: Lathi práctica 4.1.b p337

literal b. Por integración directa, encuentra la transformada X(s) y la región de convergencia para la función descrita en la imagen

Transformada Laplace f(t) Ej04

x(t) = \mu (t-2) - \mu (t-4)

que observando el ejercicio anterior se puede deducir que los retrasos en cada término generan como resultado:

X(s) = \frac{1}{s} e^{-2s}-\frac{1}{s} e^{-4s}

con polo en s=0

Transformada Laplace F(s) Ej04

Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ]  [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]
..


Ejemplo 5. Transformada de Laplace para cos(t)

Referencia: Oppenheim Ejemplo 9.4 p658

Para resolver el ejercicio de transformada de Laplace, observe que la expresión es un polinomio que se desarrolla de forma mas simple aplicando la tabla de transformadas.

x(t) = e^{-2t}\mu (t) + e^{-t} \cos (3t) \mu (t)

Transformada Laplace ft Ej05
Usando el resultado del ejemplo 1 y la tabla de transformadas se tiene que:

X(s) = \frac{1}{s+2} + \frac{s+1}{s^2 +2s+10}

agrupando términos en factores para obsevar mejor los polos

X(s) = \frac{2s^2+5s+12}{(s+2)(s^2 +2s+10)}

Transformada Laplace F(s) Ej05
Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ]  [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]
..


Ejercicio 6. Transformada de Laplace con impulso δ(t) y suma de exponenciales

Referencia: Oppenheim ejemplo 9.5 p661

El ejercicio propuesto contiene un componente de impulso unitario δ(t) aplicado en t=0. Se propone usar la tabla de transformadas de Laplace para resolverlo de forma directa, aunque se propone realizar el integral para comprobar el resultado.

x(t) = \delta(t) -\frac{4}{3} e^{-t} \mu (t) + \frac{1}{3} e^{2t} \mu (t)

Se tiene una función x(t) creciente en el tiempo y no acotada.

Transformada Laplace f(t) Ej06

usando la tabla de transformadas se tiene que:

X(s) = 1 -\frac{4}{3} \frac{1}{s+1} + \frac{1}{3} \frac{1}{s-2}

Para observar los polos y ceros se agrupa X(s) por factores

X(s) =\frac{(s-1)^2}{(s-2)(s+1)}

la gráfica respecto al dominio s, mostrando los polos en s=2 y s=-1 que generan divisiones para cero. Observe que un polo  se encuentra del lado derecho del plano, relacionado con el término creciente en el tiempo y no acotado.

Transformada Laplace F(s) Ej06
Transformada de Laplace para: [ ej1 un Exponencial ]  [ ej2 Suma de términos ]  [ ej3 escalón Desplazado ]  [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]

4.1 LTI Laplace – H(s) Diagramas de bloques, dominio ‘s’ y ecuaciones diferenciales

Referencia: Lathi 4.5 p386, Oppenheim 9.8.2 p708,

Representaciones en diagramas de bloques para los sistemas LTI causales descritos por ecuaciones diferenciales y en el dominio s por funciones racionales H(s). Los diagramas se desarrollan con software abierto como Xcos de SciLab.

[H(s) 1er orden] [H(s) 2do orden combinado de 1er orden] [H(s) 2do Orden componentes Serie o Paralelo]


Ejemplo 1. H(s) Diagrama de bloques y sistema de 1er orden

Referencia: Oppenheim ejemplo 9.28 p708, Lathi ejemplo 4.22a  p392

Considere un sistema LTI causal

\frac{\delta}{\delta t} y(t) + 3y(t) = x(t) sY(s) + 3Y(s) = X(s) sY(s) = X(s) -3Y(s) Y(s) = \frac{1}{s}[X(s) - 3Y(s)]

Se presenta un bloque integrador (1/s) de la señal de entrada X(s) con retroalimentación de salida Y(s).

Agrupando Y(s) para obtener polinomios de s

sY(s) + 3Y(s) = X(s) ( s+3) Y(s) = X(s)

Se puede expresar la relación como una función de transferencia H(s):

H(s)= \frac{Y(s)}{X(s)} =\frac{1}{s+3}

Otra forma de expresar H(s), al mutiplicar el numerador y denominador por 1/s

H(s) = \frac{1}{s+3} = \frac{\frac{1}{s}}{1+\frac{3}{s}}

que confirma que, para dos bloques H1(s) = 1/s y H2(s) = 3, interconectados en el diagrama, se tiene:

\frac{Y(s)}{X(s)} = H(s) = \frac{H_1(s)}{1+H_1(s)H_2(s)}

[H(s) 1er orden] [H(s) 2do orden combinado de 1er orden] [H(s) 2do Orden componentes Serie o Paralelo]


Ejemplo 2. H(s) Diagrama de bloques de sistema de 2do orden como combinación de 1er orden.

Referencia: Oppenheim ejemplo 9.30 p711, Lathi 4.5-3 p393

Considere el sistema causal de segundo orden con la función del sistema:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{1}{s^2 + 3s +2}

Para la ecuación diferencial se usa la forma,

H(s) = \frac{Y(s)}{X(s)} = \frac{1}{s^2 + 3s +2} (s^2 + 3s + 2)Y(s) = X(s) s^2 Y(s) + 3sY(s) + 2Y(s) = X(s)

la entrada x(t) y la salida y(t) para este sistema satisfacen la ecuación diferencial:

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

El diagrama de bloques se obtiene reordenando la ecuación de s para despejar el término de mayor grado.

s^2Y(s) = X(s) - 3sY(s) - 2Y(s) Y(s) = \frac{1}{s^2} \Big[ X(s) - 3sY(s) - 2Y(s)\Big]

reordenando los bloques:

Ejemplo 2.1 Reordenando expresión de H(s) con fracciones parciales

Usando H(s) de la primera expresión del enunciado, se usa fracciones parciales para reordenar la expresión como:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{k_1}{s+1} + \frac{k_2}{s+2}

Para encontrar las constantes k1 y k2, el numerador de la expresión de la derecha debe ser igual a 1 , desarrollando la expresión para el numerador se tiene que:

k_1s + 2k_1 + k_2s + k_2 = (k_1 + k_2)s + (2k_1 +k_2)

igualando al numerador de la parte izquierda que es 1, o expresado como 0s+1, el resultado debería ser:

(k_1 + k_2)s + (2k_1 +k_2) = 0s + 1

Se crean las ecuaciones para cada coeficiente de s en el numerador:

k_1 + k_2 = 0 2k_1 + k_2 = 1

con lo que,

k_1 = 1 , k_2 = -1

al reemplazar,

H(s) = \frac{1}{s+1} -\frac{1}{s+2}

El resultado permite realizar un diagrama equivalente en paralelo de dos sistemas mas simples. La suma de dos componentes de primer orden es un diagrama con bloques de primer orden en paralelo.

En el caso de la expresión racional inicial, se observa que se puede escribir como la multiplicación de dos sistemas de primer orden.

H(s) = \frac{1}{(s+1)(s+2)} = \Big[\frac{1}{s+1}\Big]\Big[\frac{1}{s+2}\Big]

que se representa como bloques en serie o cascada de dos sistemas de primer orden.

que por cierto, también es la convolución de dos sistemas en serie.

Ejemplo 2.2 Resolviendo fracciones parciales usando Sympy-Python

También es posible realizar las fracciones parciales con Sympy ingresando toda la expresión de H(s) de la forma:

Hs = 1/((s+1)*(s+2))

y usando sym.apart() se obtienen las fracciones parciales de la expresión,

H(s):
       1       
---------------
(s + 1)*(s + 2)

 H(s) en  fracciones parciales:
    1       1  
- ----- + -----
  s + 2   s + 1
>>> 

las instrucciones en Sympy-Python a ser usadas son:

# Fracciones parciales con Laplace
# Ps es numerador, Qs es denominador
# Oppenheim 9.30 p711
import sympy as sym

# INGRESO
s = sym.Symbol('s')

Hs = 1/((s+1)*(s+2))

# PROCEDIMIENTO
# fracciones parciales de s
Hs_parcial = sym.apart(Hs,s) 

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en  fracciones parciales:')
sym.pprint(Hs_parcial)

La instruccion sym.apart() se aplica sobre expresiones tipo polinomio, se debe considerar el caso cuando H(s) es solo un componente constante o tiene un desplazamiento de tiempo representado por sym.exp(). El asunto se trata en la siguiente página sobre fracciones parciales.

[H(s) 1er orden] [H(s) 2do orden combinado de 1er orden] [H(s) 2do Orden componentes Serie o Paralelo]


Ejemplo 3. H(s) Diagrama de bloques de sistema 2do orden componentes en serie o paralelo

Referencia: Oppenheim ejemplo 9.31.p712/pdf743, Lathi Ejemplo 4.23 p395

Considere la función del sistema

H(s) = \frac{2s^2 + 4s - 6}{s^2 + 3s + 2}

se puede reescribir como:

H(s) = \Big[\frac{1}{s^2 + 3s + 2} \Big] \Big[ 2s^2 + 4s - 6 \Big]

Para los componentes del numerador P pueden tomar respuestas del bloque denominador Q simplificando el diagrama :

Ejemplo 3.1 H(s) con fracciones parciales

Usando fracciones parciales, se puede convertir H(s) en componentes mas simples en paralelo

H(s) = 2 + \frac{6}{s+2} - \frac{8}{s+1}

Usando las raíces del numerador P se escribe H(s) de la forma en que se obtiene un sistema en serie o cascada,

H(s) = \Big[ 2\Big] \Big[\frac{s-1}{s+2} \Big] \Big[\frac{s+3}{s+1} \Big]

Ejemplo 3.2 Fracciones parciales usando Sympy

H(s):
   2          
2*s  + 4*s - 6
--------------
  2           
 s  + 3*s + 2 

 H(s) en  fracciones parciales:
      6       8  
2 + ----- - -----
    s + 2   s + 1

 H(s) como factores:
2*(s - 1)*(s + 3)
-----------------
 (s + 1)*(s + 2) 

 H(s) simplificada:
  / 2          \
2*\s  + 2*s - 3/
----------------
   2            
  s  + 3*s + 2 
# Fracciones parciales con Laplace, factores
# Ps es numerador, Qs es denominador
# Oppenheim 9.30 p711
import sympy as sym

# INGRESO
s = sym.Symbol('s')

Ps = 2*s**2+4*s-6
Qs = s**2+3*s+2

Hs = Ps/Qs

# PROCEDIMIENTO
# fracciones parciales de s
Hs_parcial = sym.apart(Hs,s)

# expresion con factores de s
Hs_factor = sym.factor(Hs_parcial,s)

# simplificar a la forma inicial
Hs_simple = sym.simplify(Hs_parcial)

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en  fracciones parciales:')
sym.pprint(Hs_parcial)
print('\n H(s) como factores:')
sym.pprint(Hs_factor)
print('\n H(s) simplificada:')
sym.pprint(Hs_simple)

La instrucción sym.factor() aplica sobre expresiones simples. En caso de disponer la expresión como polinomica, puede usar la conversión con Hs_parcial.as_expr(s).

Ejemplo 3.3 Ecuación diferencial de H(s)

H(s) = \frac{Y(s)}{X(s)} = \frac{2s^2 + 4s - 6}{s^2 + 3s + 2}

Se reordena la expresión,

(s^2 + 3s + 2)Y(s) = (2s^2 + 4s - 6) X(s) s^2 Y(s) +3s Y(s) + 2Y(s) = 2 s^2 X(s) + 4s X(s) - 6X(s)

Si considera la forma de la ecuación diferencial, tambien es la de un circuito eléctrico RLC. El primer término Y(t) sería el del inductor L, pero expresado como primera derivada. El segundo término es el resistor y el tercer término es del capacitor. Lo que se puede apreciar dividiendo toda la ecuación para ‘s’.

sY(s) +3 Y(s) + 2\frac{1}{s}Y(s) = 2 s X(s) + 4 X(s) - 6\frac{1}{s}X(s)

que es el circuito de los primeros ejemplos de Sistema – Modelo entrada-salida.

Sustituyendo las s de Laplace por derivadas y Y(s) por y(t) en la expresión de 2do orden,

\frac{\delta ^2}{\delta t^2} y(t) + 3\frac{\delta}{\delta t} y(t) + 2y(t) = 2\frac{\delta ^2}{\delta t^2} x(t) + 4\frac{\delta}{\delta t} x(t) - 6x(t)

En las próximas secciones se analiza la estabilidad del sistema para las señales de salida usando polos y ceros.

Algoritmos y funciones: telg1001.py

Resumen de algoritmos y funciones del curso señales y sistemas. Descargue una copia del archivo como telg1001.py en el directorio de trabajo, importe las funciones y use llamando con fcnm.funcion().

import telg1001 as fcnm

respuesta = fcnm.funcion(parametros)

El contenido de telg1001.py se actualiza cuando al desarrollar un ejercicio se encuentra que se puede mejorar una función manteniendo compatibilidad con lo realizado anteriormente. Se incorpora la referencia al ejercicio como comentario.

# TELG1001 - Señales y Sistemas ver 2017/05/29
# Resumen de problemas resueltos en clases /ESPOL/FCNM-FIEC/
# http://blog.espol.edu.ec/telg1001/
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
s = sym.Symbol('s')
t = sym.Symbol('t',real=True)

# Transformada de Laplace para f(t) con Sympy-Python
# http://blog.espol.edu.ec/telg1001/transformada-de-laplace-para-ft-con-sympy-python/
def separa_constante(termino):
    ''' separa constante antes de usar
        sym.laplace_transform(term_suma,t,s)
        para incorporarla luego de la transformada
        inconveniente revisado en version 1.11.1
    '''
    constante = 1
    if termino.is_Mul:
        for term_mul in termino.args:
            if not(term_mul.has(t)):
                constante = term_mul
        termino = termino/constante
    return([termino,constante])

...

Convolución – Tabla de Propiedades

Referencia: Lathi 2.4-1 p170, Hsu Schaum 2.2.D p57, Oppenheim 2.2 p90 pdf121

1, Conmutativa:

x(t) ⊗ h(t) = h(t) ⊗ x(t)

2. Asociativa

[x(t) ⊗ h1(t)] ⊗ h2(t) = x(t) ⊗ [h1(t) ⊗ h2(t)]

3. Distributiva

x(t) ⊗ [h1(t) + h2(t)] = x(t) ⊗ h1(t) + x(t) ⊗ h2(t)]

4. Desplazamiento
Si: x(t) ⊗ h(t) = c(t)

x1(t) ⊗ x2(t-T) = x1(t-T) ⊗ x2(t) = c(t-T)

x1(t-T1) ⊗ x2(t-T2) = c(t-T1 – T2)

5. Convolución con un impulso

x(t) ⊗ δ(t) = x(t)

6. Duración o ancho de señal

Si la duración (ancho) de x1(t) y x2(t) son finitas dads por T1 y T2, la duración de la convolución de x1(t) ⊗ x2(t) es T1 + T2

Convolución Integrales – Tabla

Referencia: Lathi Tabla 2.1 p176

Respuesta a estado cero y convolución con Sympy-Python

No x1(t) x2(t) x1(t)⊗x2(t) = x2(t)⊗x1(t)
1 x(t) δ(t-T) x(t-T)
2 eλt μ(t) μ(t) \frac{1-e^{\lambda t}}{-\lambda} \mu (t)
3 μ(t) μ(t) t μ(t)
4 eλ1t μ(t) eλ2t μ(t) \frac{e^{\lambda _1 t}-e^{\lambda _2 t}}{\lambda _1 - \lambda _2} \mu(t)
\lambda _1 \neq \lambda _2
5 eλt μ(t) eλt μ(t) t eλt μ(t)
6 t eλt μ(t) eλt μ(t) \frac{1}{2} t^2 e^{\lambda t} \mu(t)
7 t N μ(t) eλt μ(t) \frac{N! e^{\lambda t}}{\lambda^{N+1}} \mu(t) - \sum_{k=0}^{N}\frac{N! t^{N-k}}{\lambda^{N+1} (N-k)!} \mu(t)
8 t M μ(t) t N μ(t) \frac{M! N!}{(M+N+1)!} t^{M+N+1} \mu (t)
9 t eλ1t μ(t) eλ2t μ(t) \frac{e^{\lambda_2 t}- e^{\lambda_1 t} + (\lambda_1-\lambda_2)te^{\lambda_1 t}}{(\lambda_1-\lambda_2)^2} \mu (t)
10 t M eλt μ(t) t N eλt μ(t) \frac{M! N!}{(M+N+1)!} t^{M+N+1} e^{\lambda t}\mu (t)
11 t M eλ1t μ(t) t^N e^{\lambda_2t} \mu (t) \sum_{k=0}^{M}\frac{(-1)^k M! (N+k)! t^{M-k} e^{\lambda_1t}}{k!(M-k)!(\lambda_1 - \lambda_2)^{N+k+1}} \mu (t)
λ1 ≠λ2 + \sum_{k=0}^{N}\frac{(-1)^k N! (M+k)! t^{N-k} e^{\lambda_2t}}{k!(N-k)!(\lambda_2 - \lambda_1)^{M+k+1}} \mu (t)
12 e^{\alpha t} \cos (\beta t + \theta) \mu (t) eλt μ(t) \frac{\cos(\theta - \phi)e^{\lambda t}-e^{-\alpha t} \cos(\beta t + \theta - \phi)}{\sqrt{(\alpha + \lambda)^2 + \beta^2}} \mu (t)
\phi =tan^{-1} \Big[\frac{-\beta}{(\alpha + \lambda})\Big]
13 eλ1t μ(t) eλ2t μ(-t) \frac{e^{\lambda_1 t} \mu (t) + e^{\lambda_2 t} \mu (-t)}{\lambda_2 -\lambda_1}
\text{Re} \lambda_2 > \text{Re} \lambda_1
14 eλ1t μ(-t) eλ2t μ(-t) \frac{e^{\lambda_1 t} -e^{\lambda_2 t}}{\lambda_2 -\lambda_1} \mu (-t)

3.4.5 LTI CT – Respuesta a estado cero – Diagrama Bloques

Se continúa con el ejercicio propuesto en la sección Sistema LTIC-Modelo entrada-salida:

Referencia: 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 continuación del presentado para respuesta a entrada cero, que tiene la ecuación:

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

Al diagrama desarrollado en la sección Respuesta entrada cero – Diagrama Bloques, se añaden los bloques para formar la señal de entrada x(t). La variable tiempo se indica con el Bloque Reloj, el escalón μ(t) se configura para empezar en cero y el valor -3 del exponente se los indica con un bloque de ganancia antes de eu.

Para diferenciación de los bloques por sección, la entrada se usa el color azul, el sistema en color naranja y la salida en color morado.

Para continuar con la simulación, asegúrese que no existan condiciones iniciales en el sistema, es decir los valores de estado inicial para los bloques 1/s se encuentren en cero.

Los osciloscopios de observación, tienen su propia referencia de tiempo con otro reloj. El osciloscopio de la señal de entrada x(t) se muestra primero, y la salida en la gráfica posterior.

La señal de salida se obtiene de forma semejante a la respuesta a impulso, es decir en el punto donde está la primera derivada. Revisar teoría para aplicar la forma: P(D)[y(t)μ(t)], que en éste ejercicio se aplica una D.

El resultado de y(t), solo como referencia, es el que se obtiene en el osciloscopio ubicado a la derecha abajo del diagrama.

Con lo que se muestra que es posible observar resultados usando algunos bloques en un simulador.

3.4.4 LTI CT – Respuesta a estado cero con Scipy-Python y convolución con Numpy

Se continúa con el ejercicio propuesto en la sección Sistema-Modelo entrada-salida:

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

El ejercicio se desarrolla en Python de varias formas para obtener el resultado en una gráfica usando:

1. Un sistema LTI con Scipy.signal.lti() y simulado con Scipy.signal.lsim()
2. La integral de convolución con Numpy.convolve()
3. La integral de convolución con algoritmo numérico

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

..


1. Un sistema LTI con Scipy.signal.lti() y simulado con Scipy.signal.lsim()

Esta forma de describir el problema, simplifica el desarrollar el ejercicio a partir de la descripción de la ecuación del sistema.

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación y condiciones de inicio en cero.

La salida y(t) se calcula con una versión muestreada xi de la función de entrada x(t) a partir del tiempo muestreado en el intervalo [a,b].

respuesta Estado Cero 01 Sympy ZSR

En el algoritmo se incorpora la gráfica de la función de transferencia h(t) como referencia.

# Sistema lineal usando scipy.signal
#   Q(D)y(t)=P(D)x(t)
# ejemplo: (D^2+3D+2)y(t)=(D+0)x(t)
import numpy as np
import scipy.signal as senal

# INGRESO
# tiempo [a,b]
a = 0 ; b = 5 ; muestras=101

# Señal de entrada
u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)

# Sistema LTI
# coeficientes Q  P de la ecuación diferencial
Q = [1., 3., 2.]
P = [1., 0.]
# condiciones de inicio [y'0,y0]
inicio = [0.,0.]

# PROCEDIMIENTO
ti = np.linspace(a, b, muestras)
xi = x(ti)

sistema = senal.lti(P,Q)
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, inicio)
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.subplot(211)
plt.suptitle('Respuesta a estado cero')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')

plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(t_y, yi, color='magenta', label='y(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='y(t)')
# plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

plt.xlabel('t')
plt.legend()
plt.grid()
plt.show()

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

..


2. La integral de convolución con Numpy.convolve()

El algoritmo requiere la expresión de las funciones x(t) y h(t) para determinar y(t). Se realiza el cálculo de la convolución, que para aproximar a la integral requiere multiplicar por dt.

Note que debe mantener el valor dt pequeño,  en el ejercicio se ha definido como 0.01 que corresponden a casi 1000 muestras en el intervalo [-5,5].

Para usar la función np.convolve(), debido a la naturaleza numérica del algoritmo, se requieren muestras en tiempos, simétricos alrededor de cero. Al mostrar la gráfica, se puede prescindir de la sección negativa, sin embargo, para recordarlo se ha dejado la parte negativa del intervalo de tiempo.

# Señales contínuas en convolución
# Compare con el algoritmo en forma discreta
# Lathi ejemplo 2.6 p166
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])

x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)

# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01

# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)

# Integral de Convolucion x[t]*h[t]
# corrección de magnitud por dt para en integral
yi = np.convolve(xi,hi,'same')*dt

# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')

plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()

plt.show()

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

..


3. La integral de convolución con algoritmo numérico

Esta sección desarrolla el algoritmo para la integral de convolución, tomando paso a paso las operaciones para obtener la gráfica del resultado.

El algoritmo se presenta con propósitos didácticos, pues parece ser un poco más lento debido a que se realizan las operaciones de forma básica (interpretada) de Python. La función np.convolve() usa elementos compilados a lenguaje de maquina que se ejecutan de forma más eficiente.

Note que la única sección que cambia es la de convolución, por lo que la gráfica de salida es la misma que la sección anterior:

# Señales contínuas en convolución. Algoritmo numérico.
# compare con el algoritmo en forma discreta
# Lathi ejemplo 2.6 p166
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
u = lambda t: np.heaviside(t,1)
# u = lambda t: np.piecewise(t,t>=0,[1,0])
x = lambda t: 10*np.exp(-3*t)*u(t)
h = lambda t: (2*np.exp(-2*t)-np.exp(-t))*u(t)

# tiempo [a,b] simétrico alrededor de 0
b = 5 ; a = -b
dt = 0.01

# PROCEDIMIENTO
ti = np.arange(a,b+dt,dt)
xi = x(ti)
hi = h(ti)

# Integral de Convolución x[t]*h[t]
muestras = len(xi)
yi = np.zeros(muestras, dtype=float)
for i in range(0,muestras):
    suma=0
    for k in range(0,muestras):
        suma = suma + x(ti[k])*h(ti[i]-ti[k])
    yi[i]= suma*dt

# SALIDA - GRAFICA
plt.suptitle('Integral de Convolución x(t)*h(t)')

plt.subplot(211)
plt.plot(ti,xi,'b', label='x(t)')
plt.plot(ti,hi,'r', label='h(t)')
plt.legend()
plt.grid()

plt.subplot(212)
plt.plot(ti,yi,'m', label='x(t)*h(t)')
plt.xlabel('t')
plt.legend()
plt.grid()

plt.show()

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