Convolución – Tabla de Propiedades

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

1, Conmutativa:

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

2. Asociativa

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

3. Distributiva

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

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

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

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

5. Convolución con un impulso

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

6. Duración o ancho de señal

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

Convolución Integrales – Tabla

Referencia: Lathi Tabla 2.1 p176

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

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

3.4.5 LTI CT – Respuesta a estado cero – Diagrama Bloques

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

Referencia: Lathi ejemplo 2.6 p166

En la entrada de sistema se aplica:

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

el sistema tiene respuesta a impulso:

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

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

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

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

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

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

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

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

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

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

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

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

Referencia Lathi ejemplo 2.6 p166

En la entrada de sistema se aplica:

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

El sistema tiene respuesta a impulso:

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

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

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

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

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

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

..


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

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

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

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

respuesta Estado Cero 01 Sympy ZSR

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

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

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

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

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

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

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

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

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

plt.legend()
plt.grid()

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

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

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

..


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

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

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

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

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

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

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

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

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

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

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

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

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

plt.show()

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

..


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

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

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

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

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

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

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

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

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

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

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

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

plt.show()

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

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 b0=0.

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

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

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

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

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

El parámetros inicial se configuró en:

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

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

 

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 xi a partir del tiempo muestreado en el intervalo [a,b].

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

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]

3.2.3 LTI CT – Respuesta a entrada cero – Diagrama Bloques

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

Un diferenciador D representa un sistema inestable, pues para una entrada acotada como el escalón μ(t) resulta en una salida no acotada como el impulso unitario δ(t). Un sistema con diferenciadores expuesto a ruido amplifica el ruido.

Para los diagramas de bloques y simuladores se prefiere usar integradores (1/D) ≈ (1/s).

Con el problema planteado, en ésta sección se cambia la expresión de la ecuación, ubicando el operador D de mayor grado a la izquierda:

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

La expresión es semejante a lo realizado para el método de Runge-Kutta presentado en la sección con Python.

Se despeja la respuesta para y(t) del lado izquierdo:

\frac{D^2 y(t)}{D^2} = \frac{-3Dy(t) - 2y(t)}{D^2} y(t) = -3 \frac{1}{D} y(t) - 2 \frac{1}{D^2} y(t)

La forma de la ecuación permite crear el diagrama del sistema, usando 1/D como un integrador y 1/D2 representa una segunda integración en serie.

Los coeficientes de los términos se ubican como ganancias retroalimentadas.

Los valores iniciales se configuran en los integradores observando cual salida se usará en el sumador.

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

..


3.1 Usando el Simulador XCOS-Scilab

Para XCOS, se añade un visor de señal de salida (Sinks/Scope) y un reloj de muestreo. Para y'(0)=5 siguiendo el diagrama se ubica y'(t)=dy/dt como la salida del primer integrador . El valor se cambia dando doble click al elemento marcado en rojo y escribiendo el valor de -5.

La gráfica se obtiene usando el osciloscopio (scope)

El tiempo de observación se establece con el valor de «Simulatión/Final Integration time» a 5 segundos.

De forma semejante se puede realizar el ejercicio en otros simuladores, como por ejemplo Simulink-Matlab.


3.2 Desarrollo como circuito en simulador

La simulación en forma de circuito usando las herramientas de Simulink/Matlab y SimPowerSystems, para cada componente (Elements) se crea con «Series RLC Branch».

Para extraer las lecturas de las señal se usa un medidor (Measurements) que se ubica en serie con los demás componentes y se visualiza con un osciloscopio (Simulink/Commonly Used Blocks/Scope).

Las conexiones son similares a las que se usarían en el laboratorio. Se puede observar que las respuestas serán las mismas.

Las condiciones iniciales se pueden establecer en los componentes que pueden tener un valor inicial como el capacitor en 5 voltios.

El signo de la corriente de salida ya está incluido en la inversión de polaridad respecto a la polaridad de la fuente.

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

3.2.2 LTI CT – Respuesta a entrada cero con Scipy-Python o Runge-Kutta

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

Siguiendo con el ejemplo de la referencia planteado en «Sistema LTIC – Modelo entrada-salida» y las instrucciones del desarrollo analítico para «LTIC – Respuesta entrada cero«, se tiene la expresión:

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

Se prefiere escribir el termino de mayor orden primero y para entrada cero.Se tiene que x(t)=0, con condiciones iniciales y(0)=0, y'(0)=-5.

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

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

..


2. Respuesta entrada cero – desarrollo numérico con Scipy-Python

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

Scipy define un sistema LTI usando los coeficientes de la ecuación ordenados como la expresión P(D)/Q(D con la función Scipy.signal.lti().

Se describe una señal de entrada x(t), muestreada en un intervalo [a,b] junto a las condiciones iniciales en un vector [y'(0),y(0)].

El resultado será la simulación del paso de la señal x(t) por el sistema LTI, usando la función scipy.signal.lsim(), con lo que se obtiene la respuesta del sistema y(t) y los componentes yc=[yi(t),y0(t)], donde y0(t) representa la respuesta a entrada cero.

El algoritmo tiene el siguiente resultado:

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
# tiempo [a,b] intervalo observación
a = 0 ; b = 5 # intervalo tiempo [a,b]
muestras = 41

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

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

# 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 Entrada 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 entrada cero: [Desarrollo Analítico]  [Sympy-Python]  [Scipy-Python ]  [Runge-Kutta d2y/dx2]  [Simulador]
..


3. Respuesta entrada cero – desarrollo numérico con Runge-Kutta de Métodos numéricos en Python

Se reordena el sistema de ecuaciones para plantear una solución con el algoritmo Runge-Kutta de 2do orden para ecuaciones diferenciales con segundas derivadas

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

Se ubica el término de mayor grado a la izquierda de la ecuación:

D^2 y(t) = - 3Dy(t) - 2y(t)

sutituyendo las expresiones de las derivadas como las funciones f(x) por z y g(x) por z’:

Dy = y’ = z = f(x)
D2y = y" = z’= -3z – 2y = g(x)

Los valores iniciales de t0=0, y0=0, z0=-5 se usan en el algoritmo.

En este caso también se requiere conocer un intervalo de tiempo a observar [0,tn=5] y definir el tamaño de paso o resolución del resultado h=dt=0.1

El algoritmo permite obtener la gráfica y la tabla de datos.

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

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('ti, yi, zi')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y0(t)')

plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Entrada cero con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Más detalles del algoritmo de Runge-Kutta para segunda derivada en :

6.2.2 EDO Runge-Kutta d2y/dx2

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

3.8 Gráfica animada para interpretar el Integral de convolución con matplotlib-Python

Para interpretar mejor la operación se presenta una animación gráfica de h(t-τ) y y(t), para diferentes valores de t en el intervalo de observación [t_a,t_b] se incorpora la función graf_animada_xh_y().

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

Integral de Convolucion 01 animado

Otros ejemplos de gif animado de gráficas con matplotlib, se pueden revisar en: Una partícula en movimiento circular del curso CCPG1001 Fundamentos de Programación.

# GRAFICA CON ANIMACION ------------
def graf_animada_xh_y(xt,ht,yt,t_a,t_b,
                      muestras=101,y_nombre='y',
                      reprod_x = 4,retardo  = 200,
                      archivo_nombre = ''):
    '''grafica animada convolucionx(t) y h(t)
       en dos subgráficas con Parametros de animación trama/foto
        y_nombre = 'ZSR' # o y nombre de resultado convolución
        reprod_x = 4     # velocidad de reproducción
        retardo  = 200   # milisegundos entre tramas
        archivo_nombre = '' # crea gif animado si hay nombre
    '''
    # grafica evaluación numerica
    x_t = sym.lambdify(t,xt,modules=equivalentes)
    h_t = sym.lambdify(t,ht,modules=equivalentes)
    y_t = sym.lambdify(t,yt,modules=equivalentes)

    ti = np.linspace(t_a,t_b,muestras)
    xi = x_t(ti)
    hi = h_t(ti)
    yi = y_t(ti)

    import matplotlib.animation as animation
    # h(t-tau) para cada t
    ht_tau   = []
    for tau in range(0,muestras,reprod_x):
        ht_tau.append(h_t(ti[tau]-ti))
    tramas = len(ht_tau) # tramas creadas

    # figura con dos sub-graficas
    fig_anim = plt.figure()
    graf_a1  = fig_anim.add_subplot(211)
    graf_a2  = fig_anim.add_subplot(212)

    # grafico superior
    x_linea, = graf_a1.plot(ti,xi, color='blue',
                            label=r'$x(\tau)$')
    h_linea, = graf_a1.plot(ti,hi,color='magenta',
                             linestyle='dashed',
                             label=r'$h(\tau)$')
    htau_linea, = graf_a1.plot(ti,ht_tau[0],
                               color='magenta',
                               label=r'$h(t-\tau)$')
    punto1, = graf_a1.plot(0,0, color='magenta',marker=6)

    # grafico inferior
    color_y = 'green'
    if y_nombre=='ZSR':
        color_y ='dodgerblue'
    y_linea, = graf_a2.plot(ti,yi, color=color_y,
                            label=y_nombre+'(t)')
    punto2,  = graf_a2.plot(0,0, color=color_y,marker=6)
    y_sombra, = graf_a2.plot(ti,yi, color=color_y)
    y_sombra.set_visible(False) # Para fijar leyend()


    # Configura gráfica
    titulo = r''+y_nombre+'(t)= x(t)$\circledast$h(t)'
    graf_a1.set_title(titulo)
    ymax1 = np.max([np.max(xi),np.max(hi)])*1.11
    ymin1 = np.min([np.min(xi),np.min(hi)])-0.1*ymax1
    graf_a1.set_xlim([t_a,t_b])
    graf_a1.set_ylim([ymin1,ymax1])
    graf_a1.set_xlabel(r'$\tau$')
    graf_a1.legend()
    graf_a1.grid()

    ymax2 = np.max(yi)*1.1
    ymin2 = np.min(yi)-0.1*ymax2
    graf_a2.set_xlim([t_a,t_b])
    graf_a2.set_ylim([ymin2,ymax2])
    graf_a2.set_xlabel('t')
    graf_a2.legend()
    graf_a2.grid()

    # cuadros de texto en gráfico
    txt_x = (t_b+t_a)/2
    txt_y = ymax1*(1-0.09)
    txt_tau = graf_a1.text(txt_x,txt_y,'t='+str(t_a),
                   horizontalalignment='center')

    def trama_actualizar(i,ti,ht_tau):
        # actualiza cada linea
        htau_linea.set_xdata(ti)
        htau_linea.set_ydata(ht_tau[i])

        hasta   = i*reprod_x
        porusar = (muestras-reprod_x*(i+1))
        if porusar>=reprod_x: # en intervalo
            y_linea.set_xdata(ti[0:hasta])
            y_linea.set_ydata(yi[0:hasta])
            punto1.set_xdata(ti[hasta])
            punto1.set_ydata(0)
            punto2.set_xdata(ti[hasta])
            punto2.set_ydata(0)
        else: # insuficientes en intervalo
            y_linea.set_xdata(ti)
            y_linea.set_ydata(yi)
            punto1.set_xdata(ti[-1])
            punto1.set_ydata(0)
            punto2.set_xdata(ti[-1])
            punto2.set_ydata(0)

        # actualiza texto
        t_trama = np.around(ti[i*reprod_x],4)
        txt_tau.set_text('t= '+str(t_trama))
        
        return(htau_linea,y_linea,punto1,punto2,txt_tau)

    def trama_limpiar(): # Limpia Trama anterior
        htau_linea.set_ydata(np.ma.array(ti, mask=True))
        y_linea.set_ydata(np.ma.array(ti, mask=True))
        punto1.set_ydata(np.ma.array(ti, mask=True))
        punto2.set_ydata(np.ma.array(ti, mask=True))
        txt_tau.set_text('')
        return(htau_linea,y_linea,punto1,punto2,txt_tau)

    i   = np.arange(0,tramas,1) # Trama contador
    ani = animation.FuncAnimation(fig_anim,trama_actualizar,i ,
                                  fargs = (ti,ht_tau),
                                  init_func = trama_limpiar,
                                  interval = retardo,
                                  blit=True)
    # Guarda archivo GIF animado o video
    if len(archivo_nombre)>0:
        ani.save(archivo_nombre+'_animado.gif',
                 writer='imagemagick')
        #ani.save(archivo_nombre+'_video.mp4')
    plt.draw()
    #plt.show()
    return(ani)

# grafica animada de convolución
n_archivo = '' # sin crear archivo gif animado 
n_archivo = 'convolucion01' # requiere 'imagemagick'
figura_animada = graf_animada_xh_y(x,h,y,t_a,t_b,
                      muestras, reprod_x = 4,
                      archivo_nombre = n_archivo)
plt.show()

Para obtener un gif animado se debe asignar un 'archivo_nombre' para identificar el ejercicio; si se mantiene el nombre vacio '', solamente se crea la gráfica. El resultado se almacena en el archivo_nombre_animado.gif del directorio de trabajo,

La función también se incorpora  a telg1001.py para su uso posterior en los ejercicios

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

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

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

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

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

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

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

>>> sym.__version__
'1.11.1'

Simplificar multiplicación de impulso unitario δ(t)

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

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

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

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

xh = x*h

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

simplifica impulso Sympy 01

otro caso a considerar es:

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

xh = x*h

simplifica impulso Sympy 02

Instrucciones con Python

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

# Simplificar multiplicacion de impulso unitario d(t)*d(t-1)
# http://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import telg1001 as fcnm

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

# entrada x(t), respuesta impulso h(t)
# x = d
# h = d.subs(t,t-1)
# h = d # h= d**2

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

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

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

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

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

xh = x*h
ft = simplifica_impulso(xh)

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

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

 


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

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

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

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

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

xh = x*h

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

simplifica escalon Sympy 01

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

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

xh = x*h

simplifica escalon Sympy 02

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

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

Instrucciones con Python

# Simplificar multiplicacion de escalon unitario u(t)*u(t-1)
# http://blog.espol.edu.ec/telg1001/simplificar-multiplicacion-impulso-o-escalon-unitario-sympy/
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
import telg1001 as fcnm

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

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

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

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

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

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


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

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

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

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

    return(respuesta)

xh = x*h
ft = simplifica_escalon(xh)

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

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

..