3.3.1 Respuesta a impulso – Python

1. Fórmulas simbólicas con Sympy

Siguiendo las instrucciones de la sección teórica, usando los coeficientes de la ecuación y condiciones de inicio, se modifica el algoritmo para encontrar la respuesta a entrada cero. Las condiciones de inicio se establecen como y(0)=0 y y'(0)=0 para el problema planteado.

La ecuación diferencial: 
                        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
1 = -C1 - 2*C2
{C1: 1, C2: -1}
Solución particular:
        -t    -2*t
y(t) = e   - e    

Solución h(t): 
          -t      -2*t
h(t) = - e   + 2*e    
>>> 

se desarrolla con las instrucciones:

# Lathi 2.1.b pdf 116
# Respuesta a entrada cero con Sympy
# solucion para (D^2+ 3D + 2)y = 0
import sympy as sym

# INGRESO
# Condiciones iniciales
t0 = 0 ; y0 = 0 ; dy0 = 1

t = sym.Symbol('t')
y = sym.Function('y')
ecuacion = sym.Eq(sym.diff(y(t),t,t) + 3*sym.diff(y(t),t)+2*y(t),0)

# PROCEDIMIENTO
general = sym.dsolve(ecuacion, y(t))
general = general.expand()

# Condiciones iniciales
# rhs, righ hand side. Lado derecho de ecuacion
s0 = sym.Eq( y0, general.rhs.subs(t,t0))
s1 = sym.Eq(dy0, general.rhs.diff(t).subs(t,t0))
constante = sym.solve([s0,s1])

# Solucion particular
particular = general
for Ci in constante:
    particular = particular.subs(Ci, constante[Ci])

h = sym.Function('h')
rimpulso = sym.Eq(h(t),particular.rhs.diff(t))

# SALIDA
print('La ecuación diferencial: ')
sym.pprint(ecuacion)
print('Solución general: ')
sym.pprint(general)

print('\nEcuaciones de condiciones iniciales')
sym.pprint(s0)               
sym.pprint(s1)
print(constante)
print('Solución particular:')
sym.pprint(particular)

print('\nSolución h(t): ')
sym.pprint(rimpulso)

y con gráfica:

# GRAFICA
import numpy as np
a = 0 ; b =5 ; muestras = 101
ht = sym.lambdify(t, rimpulso.rhs, 'numpy')

ti = np.linspace(a,b,muestras)
hi = ht(ti)

import matplotlib.pyplot as plt
plt.plot(ti,hi, color='red', label='y(t)')
plt.title('Respuesta impulso unitario: y(t)='+str(rimpulso.rhs))
plt.xlabel('t')
plt.grid()
plt.show()

2. Usando un sistema LTI definido con Scipy.signal

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación y 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].

Del problema del Sistema-Modelo y la ecuacion 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 condicion inicial [1,0] y con la función de scipy, que son las usadas en la gráfica.

# 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: np.piecewise(t,t==0,[1,0])
# 10*np.exp(-3*t)*u(t)

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

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

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

# SALIDA - GRAFICA
import matplotlib.pyplot as plt

plt.subplot(211)
plt.suptitle('Respuesta a Impulso Unitario')
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='yi(t)')
# plt.plot(t_y, yc[:,0], 'g-.', label='yi(t)')
# plt.plot(t_y, yc[:,1], color='orange', label='y0(t)')

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

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:

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

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
import matplotlib.pyplot as plt
plt.plot(ti,zi, color = 'Red', label='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 la teoría antes de aplicar el algoritmo.

Publicado por

Edison Del Rosario

edelros@espol.edu.ec / Profesor del FIEC/FCNM-ESPOL