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]