Edison Del Rosario – Página 20 – Señales y Sistemas

2 abril, 201719 abril, 2023

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

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}$

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 [s_a , s_b] que permita mostrar una gráfica para F(s). Ejemplos de graficas se encuentran en la unidad 1.

El intervalo [s_a , s_b] 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

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)

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,

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

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

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.

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

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)

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}

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}:  {}

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)

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

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

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.

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.

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.

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:

y para F(s):

2 abril, 201725 abril, 2023

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)

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

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)

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

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

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

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

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 de Laplace para: [ ej1 un Exponencial ] [ ej2 Suma de términos ] [ ej3 escalón Desplazado ] [ ej4 sumas desplazadas ] [ ej5 coseno ] [ ej6 impulso ]

2 abril, 20174 febrero, 2023

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 H₁(s) = 1/s y H₂(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 k₁ y k₂, 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.

20 marzo, 20174 febrero, 2023

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

...

20 marzo, 20174 febrero, 2023

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) ⊗ h₁(t)] ⊗ h₂(t) = x(t) ⊗ [h₁(t) ⊗ h₂(t)]

3. Distributiva

x(t) ⊗ [h₁(t) + h₂(t)] = x(t) ⊗ h₁(t) + x(t) ⊗ h₂(t)]

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

x₁(t) ⊗ x₂(t-T) = x₁(t-T) ⊗ x₂(t) = c(t-T)

x₁(t-T₁) ⊗ x₂(t-T₂) = c(t-T₁ – T₂)

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 x₁(t) y x₂(t) son finitas dads por T₁ y T₂, la duración de la convolución de x₁(t) ⊗ x₂(t) es T₁ + T₂

20 marzo, 20174 febrero, 2023

Convolución Integrales – Tabla

Referencia: Lathi Tabla 2.1 p176

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

No	x₁(t)	x₂(t)	x₁(t)⊗x₂(t) = x₂(t)⊗x₁(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^λ₁t μ(t)	e^λ₂t μ(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^λ₁t μ(t)	e^λ₂t μ(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^λ₁t μ(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)$
	λ₁ ≠λ₂		$+ \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^λ₁t μ(t)	e^λ₂t μ(-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^λ₁t μ(-t)	e^λ₂t μ(-t)	$\frac{e^{\lambda_1 t} -e^{\lambda_2 t}}{\lambda_2 -\lambda_1} \mu (-t)$

8 marzo, 20172 enero, 2023

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 e^u.

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.

8 marzo, 20172 enero, 2023

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].

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]

8 marzo, 201715 enero, 2023

3.3.3 LTI CT – Respuesta a impulso unitario h(t) – Diagrama Bloques

Para usar el simulador con diagrama de bloques para respuesta a impulso unitario, es necesario usar la expresión:

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

Donde si el orden de P(D) es menor al orden de Q(D) entonces b₀=0.

Para el ejercicio desarrollado en el sistema modelo, P(D)=D, por lo que hay una sola derivada, que se puede obtener de la entrada del primer integrador de y(t).

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

Para obtener h(t), el impulso unitario x(t), es semejante a tener una condicion inicial sin señal de entrada, siguiendo la teoría descrita de sobre condiciones iniciales.

El impulso unitario representa que y'(t) inicia con 1, por ser el término de más alto orden en P(D) y se configura el diagrama de bloques:

Con lo que la respuesta del sistema al impulso unitario que se obtiene es:

El parámetros inicial se configuró en:

Como referencia para observación se muestra la salida y(t):

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

8 marzo, 201715 enero, 2023

3.3.2 LTI CT – Respuesta a impulso unitario h(t) con Scipy-Python o Runge-Kutta

2. Respuesta a impulso con sistema LTI definido en Scipy.signal

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación P, Q y las condiciones de inicio. La respuesta a un impulso unitario δ(t) o función de transferencia h(t) se calcula con una versión muestreada x_i a partir del tiempo muestreado en el intervalo [a,b].

sistema = senal.lti(P,Q)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)

Del problema del «Sistema LTIC – Modelo entrada-salida» y la ecuación de segundo orden, se tiene dos opciones para la respuesta:

– scipy.signal.impulse() donde solo se da el vector ti
– aplicar la teoria de los coeficientes donde solo la derivada de mayor orden de P(D) evaluada en cero, tiene el valor de uno. Todos los demás términos tienen el valor de cero. Revisar la parte conceptual de la sección para respuesta a impulso.

Para el ejercicio y como comprobación se calculan las dos formas, con entrada cero y condición inicial [1,0] y con la función de scipy, que son las usadas en la gráfica.

Instrucciones en Python

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

# INGRESO
# Señal de entrada x(t)
x = lambda t: t*0

# Sistema LTI
# coeficients Q  P de la ecuación diferencial
Q = [1., 3., 2.]
P = [1., 0.]
# condiciones iniciales [y'(0),y(0)]
cond_inicio = [1,0]

# grafica, parametros
a = 0 ; b = 5 # intervalo observacion [a,b] 
muestras = 101

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

sistema = senal.lti(P,Q)

# respuesta entrada cero
[t_y, yi, yc] = senal.lsim(sistema, xi, ti, cond_inicio)
# respuesta a impulso
[t_h, hi] = sistema.impulse(T = ti)

# SALIDA - GRAFICA
plt.subplot(211)
plt.suptitle('Respuesta a Impulso Unitario: h(t) = dy/dt')
plt.plot(ti, xi, color='blue', label='x(t)')
plt.plot(t_h, hi, color='red', label='h(t)')
plt.legend()
plt.grid()

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

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

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

3. Algoritmo de Runge-Kutta d2y/dx2 de Análisis numérico

Para el caso del algoritmo de Runge-Kutta, se aplica las condiciones iniciales de y(0)= 0 y y'(0)=1,

siguiendo lo descrito en la parte teórica y se toman los valores de z=y’ para obtener el resultado.

Las instrucciones aplicadas al algoritmo de entrada cero se muestran a continuación:

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = 1

h = 0.1
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
#plt.plot(ti,yi, label='y(t)')
plt.plot(ti,zi, color = 'Red', label='h(t)=dy/dt')

plt.ylabel('dy/dt')
plt.xlabel('t')
plt.title('Respuesta impulso: Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Recuerde que se aplican las modificaciones de acuerdo al planteamiento del problema, por lo que revise los conceptos antes de aplicar el algoritmo.

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