4.3.2 H(s) – Polos reales y complejos con Sympy-Python

Referencia: Lathi ejemplo 4 p330, Hsu literal C. p112

La función de transferencia H(s) se escribe como P(s)/Q(s) que es una expresión racional de s. Se considera que en la expresión el grado m del polinomio del numerador P(s)  es mayor que el grado n del polinomio del denominador Q(s). Las raíces del polinomio P(s) del numerador se conocen como ceros, pues convierten la expresión en cero, mientras que para el denominador las raices se denominan polos al convertir en infinito la expresión (Hsu).

\frac{P(s)}{Q(s)} = \frac{a_0(s-ceros(1))\text{...}(s-ceros(m))}{b_0(s - polos(1)) \text{...} (s - polos(n))}

otra forma de escribir los polos y ceros (Lathi)

\frac{P(s)}{Q(s)} = \frac{ceros(1)}{s - polos(1)}+\frac{ceros(2)}{s - polos(2)}+\text{ ... } \text{ ... }+\frac{ceros(n)}{s - polos(n)}+ganancia(s)

Para el desarrollo de un algoritmo con Sympy, el primer paso consiste en determinan los polos de un ejercicio.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo, busca_polosceros(Hs) ]


Ejemplo 1. Polos reales no repetidos y fracciones parciales

Referencia: Hsu problema 3.17a p137, Lathi ejemplo 4.3a p338, Oppenheim ejemplo 9.9 p671

La expresión H(s)  tiene los componentes de P(s) y Q(s) como se indican,

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

Ingreso de H(s) como P(s) y Q(s)

Para la expresión de la función de transferencia, con Sympy se crean los polinomios del numerador Ps y denominador Qs.

# INGRESO
s = sym.Symbol('s')
Ps = 2*s + 4  # 1+0*s cuando es constante
Qs = s**2 + 4*s + 3

Hs = Ps/Qs

Siendo H(s) la forma en que se expresan los ejercicios, se realizan las operaciones para obtener las raíces del numerador y denominador. El resultado buscado con el algoritmo es:

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

 H(s) en fracciones parciales 
  1       1  
----- + -----
s + 3   s + 1

 {Q_polos:veces}: {-1: 1, -3: 1}
 {P_ceros:veces}: {-2: 1}
>>>  

Instrucciones en Python

Hs se convierte en una expresión de Sympy con sym.sympify(Hs)s para el caso en que se de tan solo una constante. En el caso que Hs requiera agrupar términos se usa la instrucción sym.simplify(Hs,inverse=True).

Una vez preparada la expresión Hs, se separa como numerador Ps y denominador Qs en forma de polinomio para obtener las raíces y las veces que ocurren con la instrucción sym.roots().

La expresión de H(s) en fracciones parciales se obtienen usando sym.apart(Hs,s).

# Busca Polos y Ceros de H(s) con Sympy
# Ps es numerador, Qs es denominador
import sympy as sym

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

Ps = 2*s + 4
Qs = s**2 + 4*s + 3
Hs = Ps/Qs

# PROCEDIMIENTO
Hs = sym.sympify(Hs) # convierte a sympy una constante
Hs = sym.simplify(Hs,inverse=True) # agrupa terminos

# polos y ceros de Hs
[P,Q] = Hs.as_numer_denom()
P = P.as_poly(s)  # numerador
Q = Q.as_poly(s)  # denominador
P_ceros = sym.roots(P)
Q_polos = sym.roots(Q)

# en factores
Hs = sym.factor(Hs,s)
# fracciones parciales
Hs_fp = sym.apart(Hs,s)

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en fracciones parciales ')
sym.pprint(Hs_fp)
print('\n {Q_polos:veces}:',Q_polos)
print('\n {P_ceros:veces}:',P_ceros)

Para visualizar el concepto de polos y ceros, se presenta la gráfica de polos y ceros en el dominio s, para un corte en el plano real. Por simplicidad se usará el procedimiento creado para el curso en telg1001.py

Todos los polos y ceros se encuentran en el lado izquierdo del plano, es decir su parte real es negativa.

# GRAFICA -----------
import matplotlib.pyplot as plt
import telg1001 as fcnm
fig_Hs = fcnm.graficar_Fs(Hs,Q_polos,P_ceros,f_nombre='H')
plt.show()

Hs Polos Ceros Ej01

Para revisar los polos de la expresión de fracciones parciales, se considera que los términos de las expresiones podrían tener polos repetidos como en el ejemplo 3. Para polos repetidos, se usa el número de veces mayor cuando aparecen polos de términos suma.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]

..


Ejemplo 2. Grados iguales de P(s) y Q(s), ceros con valores complejos

Referencia: Lathi ejemplo 4.3b p338, Hsu ejercicio 3.20 p140, Oppenheim ejemplo 9.36 p716

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

Ejemplo 2 Desarrollo con Sympy-Python

Al algoritmo del ejercicio anterior se modifica solo el bloque de ingreso con las expresiones para numerador y denominador:

Ps = 2*s**2+5
Qs = s**2 + 3*s + 2
Hs = Ps/Qs

El resultado Hs_fp en fracciones parciales muestra que existe un término constante en las sumas de la expresión. La forma descrita por Lathi al inicio de la página, la constante se la considera como ‘ganancia(s)’.

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

 H(s) en fracciones parciales 
      13      7  
2 - ----- + -----
    s + 2   s + 1


 {Q_polos:veces}: {-1: 1, -2: 1}
 {P_ceros:veces}: {-sqrt(10)*I/2: 1, sqrt(10)*I/2: 1}
>>>  

Por otra parte se tiene que los ceros tienen raíces con componente imaginario representado con el símbolo ‘I‘, que para representar cada raíz en la gráfica se considera el eje vertical también como parte Imag(s). Si es cierto, es diferente a la naturaleza de H(s), perdonen, se usa ésta «libertad» solo para resaltar la influencia de los polos en H(s) como se mostró en el ejercicio anterior, o se realizaría un gráfico con el plano real e imaginario como base y una tercera dimensión para H(s).

Los polos se encuentran en el lado izquierdo del plano, es decir su parte real es negativa.

Hs Polos Ceros Ej02

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]


Ejemplo 3. Polos reales repetidos o raices repetidas en denominador Q(s)

Referencia: Lathi ejemplo 4.3d p342

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

la forma de la expresión para el algoritmo es

Ps = 8*s+10
Qs = (s+1)*((s+2)**3)

Hs = Ps/Qs

En él ejercicio el término del denominador tiene un factor elevado al cubo, lo que muestra que la raiz es repetida y se debe reflejar en la respuesta del algoritmo. En las fracciones parciales se muestra que para el polo repetido existe un término suma para cada grado en el denominador.

H(s):
  2*(4*s + 5)   
----------------
               3
(s + 1)*(s + 2) 

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

 {Q_polos:veces}: {-1: 1, -2: 3}
 {P_ceros:veces}: {-5/4: 1}
>>> 

Todos los polos y ceros se encuentran en el lado izquierdo del plano s.

Se adjunta la gráfica para la interpretación de H(s)

H(s) Polos Ceros Ej03

Si la busqueda de polos es sobre la expresión de fracciones parciales, se tiene que las veces se determinan como el mayor exponente del término del denominador.

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]

..


Ejemplo 4. H(s) con Fracciones parciales y polos de tipo complejo

Referencia: Lathi ejemplo 4.3c p340

H(s) = \frac{6(s+34)}{s(s^2+10s+34)}

la forma de la expresión para el algoritmo es

Ps = 6*(s+34)
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

Incorporando las instrucciones para identificar si el denominador tiene raices complejas, se asignan los valores de los parametros necesarios en un diccionario de variables, obteniendo como resultado,

H(s):
    6*(s + 34)    
------------------
  / 2            \
s*\s  + 10*s + 34/

 H(s) en fracciones parciales 
    6*(s + 9)      6
- -------------- + -
   2               s
  s  + 10*s + 34    

 {Q_polos:veces}: {-5 - 3*I: 1, -5 + 3*I: 1, 0: 1}
 {P_ceros:veces}: {-34: 1}

 h(t): 
                               -5*t                              
- 2*(4*sin(3*t) + 3*cos(3*t))*e    *Heaviside(t) + 6*Heaviside(t)
>>> 

grafica ejemplo 4

H(s) Polos Ceros Ej04

Ejemplo 4 Instrucciones en Python

# Fracciones parciales de H(s) con Sympy
# Ps es numerador, Qs es denominador
# con Transformada Inversa de Laplace
import numpy as np
import sympy as sym

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

Ps = 6*(s+34) 
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

# PROCEDIMIENTO
Hs = sym.simplify(Hs,inverse=True)
# Hs debería ser un solo termino suma
Hs = sym.factor(Hs,s)
# fracciones parciales
Hsp = sym.apart(Hs,s)

# Analiza polos
def polosceros_simple(Hs):
    ''' Busca_polo de un termino sin exp(-a*s)
    '''
    Hs = sym.simplify(Hs,inverse=True)
    # Hs debería ser un solo termino suma
    Hs_factor = sym.factor(Hs,s)
    # polos y ceros de termino Hs
    [P,Q] = Hs_factor.as_numer_denom()
    P = sym.poly(P,s)  # numerador
    Q = sym.poly(Q,s)  # denominador
    P_ceros = sym.roots(P)
    Q_polos = sym.roots(Q)
    respuesta = {'Q_polos' : Q_polos,
                 'P_ceros' : P_ceros}
    return(respuesta)

polosceros = polosceros_simple(Hs)

# SALIDA
print('H(s):')
sym.pprint(Hs)
print('\n H(s) en fracciones parciales ')
sym.pprint(Hsp)
print('\n {Q_polos:veces}:',polosceros['Q_polos'])
print(' {P_ceros:veces}:',polosceros['P_ceros'])
print('\n h(t): ')
sym.pprint(ht)

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]


Ejemplo 5. Con suma de términos y exponenciales de retraso en tiempo

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

la expresión para el algoritmo se escribe como:

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

Para el segundo término que tiene un exponencial, no se puede aplicar directamente la instrucción sym.apart(Hs,s) al no ser reconocida como tipo polinomio. Para continuar se debe separar el término sym.exp() del resto de la expresión, y se podrá aplicar fracciones parciales a esa parte, para luego volver a poner el término.

Para facilitar el desarrollo del algoritmo, se recomienda empezar a analizar una expresión mas simple, como solamente el segundo término de la suma.

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

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

H(s) Polos Ceros Ej05

En el algoritmo se empieza por determinar si la expresión de Hs tiene términos sym.exp(-s), la instrucción que permite listar estos términos es

lista_exp = list(Hs.atoms(sym.exp(-s)))

como la expresión H(s) en los ejercicios puede contener sumas o multiplicaciones de bloques, primero se simplifica la expresión H(s)

    # convierte a Sympy si es solo constante
    Hs = sym.sympify(Hs)
    # simplifica operación Hs1(s)*Hs2(s)
    Hs = sym.expand_power_exp(Hs)
    Hs = sym.expand(Hs, power_exp=False)
    
    term_sum = sym.Add.make_args(Hs)
    Hs_0 = 0*s # reagrupa Fs
    for term_k in term_sum:
        term_k = sym.simplify(term_k)
        term_k = sym.expand_power_exp(term_k)
        Hs_0 = Hs_0 + term_k
    
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_0.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

luego se tabula la lista_exp y se eliminan los términos que pueden ser constantes o no dependen de ‘s’, por ejemplo sym.exp(10).

Luego se separa la expresión por grupos usando sym.collect() creando un diccionario con los elementos de la lista_exp y con las expresiones agrupadas.

    # separados por lista_exp
    separados = sym.collect(Hs_0,lista_exp,evaluate=False)

Con los elementos separados se puede aplicar la búsqueda de polos y ceros para cada elemento de la lista con la función polosceros_simple(Hs_k)

Al final se pueden agrupar los Q_polos y P_ceros de toda la expresión siguiendo criterios como que las veces que ocurre cada polo es el maximo entre los términos separados por lista_exp.

Instrucciones en Python

# Busca Polos y Ceros de H(s) con Sympy
# separa términos con sym.exp(-a*s)
# Ps es numerador, Qs es denominador
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1001 as fcnm

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

# Ps = (2*s + 4)
# Qs = s**2 + 4*s + 3

# Ps = 2*s**2+5
# Qs = s**2 + 3*s + 2

# Ps = 8*s+10
# Qs = (s+1)*((s+2)**3)

# Ps = 6*(s+34)
# Qs = s*(s**2+10*s+34)

# Hs = Ps/Qs
# Hs = sym.exp(-2*s)/(s-2)
# Hs = 1/s - sym.exp(-2*s)/s
# Hs = sym.exp(-2*s)/s -sym.exp(-4*s)/s
# Hs = (2*s**2+5*s+12)/((s+2)*(s**2+2*s+10))
Hs = ((s+2)/(s-1)) + sym.exp(-2*s)/(s-2)

# PROCEDIMIENTO
def busca_polosceros(Hs):
    ''' Busca polos de Hs (divisiones para cero)
        y ceros de Hs,cuenta las veces que aparecen,
        agrupa por exp(-a*s) agrupado en factores.
    '''
    def polosceros_simple(Hs):
        ''' Busca_polo de un termino sin exp(-a*s)
        '''
        Hs = sym.simplify(Hs,inverse=True)
        # Hs debería ser un solo termino suma
        Hs_factor = sym.factor(Hs,s)
        # polos y ceros de termino Hs
        [P,Q] = Hs_factor.as_numer_denom()
        P = sym.poly(P,s)  # numerador
        Q = sym.poly(Q,s)  # denominador
        P_ceros = sym.roots(P)
        Q_polos = sym.roots(Q)
        respuesta = {'Q_polos' : Q_polos,
                     'P_ceros' : P_ceros}
        return(respuesta)
    
    # convierte a Sympy si es solo constante
    Hs = sym.sympify(Hs)
    # simplifica operación Hs1(s)*Hs2(s)
    Hs = sym.expand_power_exp(Hs)
    Hs = sym.expand(Hs, power_exp=False)
    
    term_sum = sym.Add.make_args(Hs)
    Hs_0 = 0*s # reagrupa Fs
    for term_k in term_sum:
        term_k = sym.simplify(term_k)
        term_k = sym.expand_power_exp(term_k)
        Hs_0 = Hs_0 + term_k
    
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_0.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

    # separados por lista_exp
    separados = sym.collect(Hs_0,lista_exp,evaluate=False)

    # polos y ceros por terminos exp(-s) agrupados
    polosceros = {} ; Hs_fp = 0
    for k in separados:
        Hs_k = sym.factor(separados[k],s)
        PZ_k = polosceros_simple(Hs_k)
        PZ_k['Hs_k']  = Hs_k
        polosceros[k] = PZ_k

    # integra polos
    Q_polos = {} ; P_ceros = {}
    for k in polosceros:
        polos = polosceros[k]['Q_polos']
        for unpolo in polos:
            if unpolo in Q_polos:
                veces = max([Q_polos[unpolo],polos[unpolo]])
                Q_polos[unpolo] = veces
            else:
                Q_polos[unpolo] = polos[unpolo]
        ceros = polosceros[k]['P_ceros']
        for uncero in ceros:
            if uncero in P_ceros:
                veces = max([P_ceros[uncero],ceros[uncero]])
                P_ceros[unpolo] = veces
            else:
                P_ceros[uncero] = ceros[uncero]
    # revisa exp(-a*s)
    if len(lista_exp)==0: #sin componentes
        del polosceros['1']
    polosceros['Q_polos'] = Q_polos
    polosceros['P_ceros'] = P_ceros
    return(polosceros)

polosceros = busca_polosceros(Hs)

# SALIDA
print('\n expresion H(s):')
sym.pprint(Hs)
print('\nPolos y ceros:')
lista_PZ = ['Q_polos','P_ceros']
for k in polosceros:
    if not(k in lista_PZ):
        print('término:', k,'*',polosceros[k]['Hs_k'])
        print(' Q_polos{polos:veces}: ',polosceros[k]['Q_polos'])
        print(' P_ceros{polos:veces}: ',polosceros[k]['P_ceros'])
    else:
        print(k,polosceros[k])

La gráfica en éste caso se realiza usando Hs y sus componentes, de ésta forma es posible visualizar el efeto de cada componente con sus polos y ceros.

# GRAFICA ----------------------
def graficar_Fs2(Fs,Q_polos={},P_ceros={},
                s_a=1,s_b=0,muestras=101,f_nombre='F',
                solopolos=False,polosceros={}):
    # grafica Fs
    lista_Hs =[]

    lista_PZ = ['Q_polos','P_ceros']
    lista_term = list(polosceros.keys())
    lista_term.remove('Q_polos')
    lista_term.remove('P_ceros')

    # intervalo para s,y componentes Fs con exp(-a*s)
    s_a = 0 ; s_b = 0 ; Lista_Hs = []
    cond_componente = False
    if len(lista_term)>0:
        cond_componente = True
        if len(polosceros)>1: # tiene componenes con exp(-a*s)
            lista_Hs.append([Hs,{},{}])
        for k in polosceros:
            if not(k in lista_PZ):
                Q_polosk = polosceros[k]['Q_polos']
                P_cerosk = polosceros[k]['P_ceros']
                [s_ak,s_bk] = fcnm.intervalo_s(Q_polosk,P_cerosk)
                s_a = min(s_a,s_ak)
                s_b = max(s_b,s_bk)
                lista_Hs.append([polosceros[k]['Hs_k']*sym.sympify(k),
                                 Q_polosk,P_cerosk])
    else: #sin componenes con exp(-a*s)
        lista_Hs.append([Hs,polosceros['Q_polos'],polosceros['P_ceros']])
        [s_a,s_b] = fcnm.intervalo_s(polosceros['Q_polos'],
                                       polosceros['P_ceros'])
    # graficas por cada componente
    fig_Fs, graf_Fs = plt.subplots()
    # no presenta errores de división para cero en lambdify()
    np.seterr(divide='ignore', invalid='ignore')
    for k in lista_Hs:
        Fs_k = k[0] ;Q_polosk = k[1] ;P_cerosk = k[2]
        
        # convierte a sympy una constante
        Fs_k = sym.sympify(Fs_k,s) 
        if Fs_k.has(s): # no es constante
            F_s = sym.lambdify(s,Fs_k,modules=fcnm.equivalentes)
        else:
            F_s = lambda s: Fs_k + 0*s
        
        s_i = np.linspace(s_a,s_b,muestras)
        Fsi = F_s(s_i) # Revisar cuando s es complejo
        lineaestilo = 'solid'
        if len(Q_polosk)>0 and cond_componente:
            lineaestilo = 'dashed'
            
        graf_Fs.plot(s_i,Fsi,label=Fs_k,linestyle=lineaestilo)
        lineacolor = plt.gca().lines[-1].get_color()

        if len(Q_polosk)>0:
            polo_re = [] ; polo_im = []
            for polos in Q_polosk.keys():    
                graf_Fs.axvline(sym.re(polos),color='red',
                                linestyle='dotted')
                x_polo = sym.re(polos)
                y_polo = sym.im(polos)
                polo_re.append(x_polo)
                polo_im.append(y_polo)
                etiqueta = "("+str(x_polo)+','+str(y_polo)+")"
                graf_Fs.annotate(etiqueta,(x_polo,y_polo))
        
            graf_Fs.scatter(polo_re,polo_im,marker='x',
                            color =lineacolor,
                            label = Q_polosk)
        if len(P_cerosk)>0:
            cero_re = [] ; cero_im = []
            for cero in P_cerosk.keys():    
                x_cero = sym.re(cero)
                y_cero = sym.im(cero)
                cero_re.append(x_cero)
                cero_im.append(y_cero)
                etiqueta = "("+str(x_cero)+','+str(y_cero)+")"
                graf_Fs.annotate(etiqueta,(x_cero,y_cero))
        
            graf_Fs.scatter(cero_re,cero_im,marker='o',
                            color =lineacolor,
                            label = P_cerosk)

    graf_Fs.axvline(0,color='gray')
    graf_Fs.legend()
    graf_Fs.set_xlabel('Real(s)')
    graf_Fs.set_ylabel('F(s) ; Imag(s)')
    graf_Fs.set_title(r'F(s) = $'+str(sym.latex(Hs))+'$')
    graf_Fs.grid()
    return(fig_Fs)

fig_Fs = graficar_Fs2(Hs,polosceros=polosceros)
plt.show()

Parametros cuadráticos

Para el ejercicio 4, se tienen polos con valores complejos y para usar la tabla de transformadas de Laplace se requiere obtener algunos parámetros.

Referencia: Lathi ejemplo 4.3c p340

H(s) = \frac{6(s+34)}{s(s^2+10s+34)}

el desarrollo se puede describir como

H(s) = \frac{6(s+34)}{s(s^2+10s+34)} = \frac{k_1}{s} + \frac{As+B}{s^2+10s+34}

usando el método de Heaviside, se obtiene k1

k_1 = \frac{6(s+34)}{\cancel{s}(s^2+10s+34)}\Big|_{s=0} =\frac{6(0+34)}{(0^2+10(0)+34)} = \frac{6(34)}{34}=6 \frac{6(s+34)}{s(s^2+10s+34)} = \frac{6}{s} + \frac{As+B}{s^2+10s+34}

se simplifican las fracciones al multiplicar ambos lados por s(s2+10s+34)

6(s+34) = 6(s^2+10s+34) + (As+B)s 6s+204 = 6s^2+60s+204 + A s^2+Bs 6s+204 = (6+A)s^2+(60+B)s+204

se puede observar que 6+A=0, por lo que A=-6.
Tambien se tiene que 60+B=6, por lo que B=-54
quedando la expresión como:

H(s) = \frac{6}{s} + \frac{-6s-54}{s^2+10s+34}

Observando la tabla de transformadas de Laplace se compara la expresión con las filas 12c y 12d se usan los parametros de:

\frac{As+B}{s^2+2as+c}

A= -6, B=-64, a=5, c = 34
y se calculan:

b = \sqrt{c-a^2} = \sqrt{34-(5)^2} = 3 r = \sqrt{\frac{A^2 c +B^2 -2ABa}{c-a^2}} r = \sqrt{\frac{(-6)^2 (34) +(-64)^2 -2(-6)(-64)(5)}{(34)-(5)^2}}=10 \theta = \tan ^{-1} \Bigg( \frac{(-6)(5)-(-64)}{(-6)\sqrt{34-(5)^2}} \Bigg) = -0.9272

con lo que se puede escribir la transformada inversa de Laplace para h(t) como

h(t) = [6+10e^{-5t}\cos (3t--0.9272)] \mu (t)

o la otra foma de expresión usando la fila 12d de la tabla

Ejemplo 4 Desarrollo con Sympy-Python

Incorporando las instrucciones para identificar si el denominador tiene raices complejas, se asignan los valores de los parámetros necesarios en un diccionario de variables, obteniendo como resultado,

H(s) en fracciones parciales:
    6*(s + 9)      6
- -------------- + -
   2               s
  s  + 10*s + 34    

 parametros de Q_cuadratico: 
  -6*(s + 9)/(s**2 + 10*s + 34)  :
   {'A': -6.0, 'B': -54.0, 'a': 5.0, 'c': 34.0, 'r': 10.0, 'b': 3.0, 'theta': -0.9272952180016122}
>>> 

mientras que el algoritmo se resume en

# Fracciones parciales de H(s) con Sympy
# Ps es numerador, Qs es denominador
# con Transformada Inversa de Laplace
import numpy as np
import sympy as sym
import telg1001 as fcnm

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

Ps = 6*(s+34) 
Qs = s*(s**2+10*s+34)

Hs = Ps/Qs

# PROCEDIMIENTO
# parametros de termino con Q cuadratico
def Q_cuad_s_parametros(Hs):
    '''Si Q tiene grado=2, obtiene los parametros
       para la tabla de Transformadas Inversas de Laplace
    '''
    def Q_cuad_sin_exp(Hs):
        '''para Hs sin sym.exp(-a*s)
        '''
        # fracciones parciales NO separa Q con factores complejos
        Hs_fp = sym.apart(Hs,s)
        respuesta = {}
        term_suma = sym.Add.make_args(Hs_fp)
        for term_k in term_suma:
            [P,Q]  = term_k.as_numer_denom()
            if sym.degree(Q) == 2 and sym.degree(P) == 1:
                P_coef = P.as_coefficients_dict(s)
                Q_coef = Q.as_coefficients_dict(s)
                Q0 = 1 # Normalizar coefficiente s**2=1
                if Q_coef[s**2]!=1: 
                    Q0 = Q_coef[s**2]
                # Parametros de Q cuadratico
                a = 0 ; c = 0
                if s**1 in Q_coef:
                    a = float(Q_coef[s**1]/Q0)/2
                if s**0 in Q_coef:
                    c = float(Q_coef[s**0]/Q0)
                A = float(P_coef[s**1])
                B = 0
                if s**0 in P_coef:
                    B = float(P_coef[s**0])
                rP = (A**2)*c + B**2 - 2*A*B*a
                rQ = c - a**2
                r  = np.sqrt(rP/rQ)
                b  = np.sqrt(c-a**2)
                thetaP = A*a-B
                thetaQ = A*np.sqrt(c-a**2)
                theta  = np.arctan(thetaP/thetaQ)
                parametro = {'A': A, 'B': B,'a': a, 'c': c,
                             'r': r, 'b': b,'theta':theta}
                respuesta[term_k] = parametro
        return (respuesta)

    Hs = sym.sympify(Hs,s)
    Hs_fp = fcnm.apart_s(Hs) # separa por sym.exp(-a*s)
    # tabula sym.exp(-s) en lista_exp
    lista_exp = list(Hs_fp.atoms(sym.exp(-s)))
    if len(lista_exp)>0: # elimina constantes
        for k in lista_exp:
            if not(k.has(s)): # es constante
                lista_exp.remove(k)

    # separados por lista_exp
    separados = sym.collect(Hs_fp,lista_exp,evaluate=False)

    Qs2 = {} # agrupa parametros Qs2
    for k in separados:
        Qs2_k = Q_cuad_sin_exp(separados[k])
        if len(Qs2_k)>0:
            for term in Qs2_k:
                Qs2[k*term] = Qs2_k[term]
    return(Qs2)

Hs_fp = fcnm.apart_s(Hs)
Qs2 = Q_cuad_s_parametros(Hs)

# SALIDA
print('H(s) en fracciones parciales:')
sym.pprint(Hs_fp)

# parametros de termino con Q cuadratico
if len(Qs2)>0:
    print('\n parametros de Q_cuadratico: ')
    for k in Qs2:
        print(' ',k,' :')
        print('  ',Qs2[k])

[ polos reales no repetidos ]  [ grados iguales P y Q ]  [ polos reales repetidos ]  [polos complejos ]  [ con exponencial s o retraso en tiempo ]

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

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

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

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

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

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


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

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

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

La ecuación diferencial con operadores D es:

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

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

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

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

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

import sympy as sym

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

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

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

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

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

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

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

# ecuacion homogenea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general entrada cero
general = sym.dsolve(homogenea, y(t))
general = general.expand()

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

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

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

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

Se realiza lo mismo para y(0)=0

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

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

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

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

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

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

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

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

resultado del algoritmo:

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

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

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

Instrucciones en Python

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

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

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

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

# PROCEDIMIENTO
edo_tipo = sym.classify_ode(ecuacion, y(t))

# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0

# solucion general de ecuación homogénea
general = sym.dsolve(homogenea, y(t))
general = general.expand()

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

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

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

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

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

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

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

Gráfica de respuesta a entrada cero ZIR

respuesta Entrada Cero 01 Sympy ZIR

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

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

Se añaden las siguientes instrucciones al algoritmo anterior:

# GRAFICA ------------------
import numpy as np
import matplotlib.pyplot as plt

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

# PROCEDIMIENTO
yt = sym.lambdify(t,ZIR, modules=equivalentes) 
ti = np.linspace(t_a,t_b,muestras)
yi = yt(ti)

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

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

..


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

Referencia: Lathi Ejemplo 2.1.b p155

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

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

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

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

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

respuesta Entrada Cero 02 Sympy ZIR

Instrucciones en Python

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

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

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

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

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

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

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

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

# PROCEDIMIENTO
def respuesta_ZIR(ecuacion,cond_inicio=[],t0,
                y = sym.Function('y'),x = sym.Function('x')):
    ''' Sympy: ecuacion: diferencial en forma(LHS,RHS),
        condiciones de inicio t0 y [y'(t0),y(t0)]
        cond_inicio en orden descendente derivada
    '''
    # ecuacion homogenea x(t)=0, entrada cero
    RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
    LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
    homogenea = LHSx0 - RHSx0

    # solucion general entrada cero
    general = sym.dsolve(homogenea,y(t))
    general = general.expand()

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

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

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

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

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

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


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

Referencia: Lathi Ejemplo 2.1.c p155

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

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

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

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

Los resultados con el algoritmo son:

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

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

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

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

respuesta Entrada Cero 03 Sympy ZIR
Instrucciones con Python

# Respuesta a entrada cero ZIR con Sympy-Python
# http://blog.espol.edu.ec/telg1001/lti-ct-respuesta-entrada-cero-con-sympy-python/
# Lathi 2.1.c pdf 155, (D^2+ 4D + 40)y = (D+2)x
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import telg1001 as fcnm

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

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

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

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

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

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

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

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

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

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

iimport telg1001 as fcnm

respuesta = fcnm.funcion(parametros)

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

Mostrar el resultado del diccionario según el tipo de entrada

Para mostrar los resultados de una forma más fácil de leer, se usará sym.pprint() para las ecuaciones, y print() para los demás elementos del resultado

def print_resultado_dict(resultado):
    ''' print de diccionario resultado
        formato de pantalla
    '''
    eq_sistema = ['ZIR','h','ZSR','xh']
    for entrada in resultado:
        tipo = type(resultado[entrada])
        cond = (tipo == sym.core.relational.Equality)
        cond = cond or (entrada in eq_sistema)
        if cond:
            print(entrada,':')
            sym.pprint(resultado[entrada])
        elif tipo==list or tipo==tuple:
            tipoelem = type(resultado[entrada][0])
            if tipoelem==sym.core.relational.Equality:
                print(entrada,':')
                for fila in resultado[entrada]:
                    sym.pprint(fila)
            else:
                print(entrada,':')
                for fila in resultado[entrada]:
                    print(' ',fila)
        else:
            print(entrada,':',resultado[entrada])
    return()