Series de Fourier

Referencia: Lathi 6.1 pdf406, Chapra 5ed Ejemplo 19.2 p548/pdf572

La serie de Fourier aproxima una señal o función contínua mediante una serie infinita de sinusoides.

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

donde los coeficientes de la ecuación se calculan como:

a_k = \frac{2}{T} \int_0^T f(t) \cos(k \omega_0 t) \delta t b_k = \frac{2}{T} \int_0^T f(t) \sin(k \omega_0 t) \delta t

Ejemplo

Utilice la serie de Fourier continua para aproximar la función de onda cuadrada o rectangular. Para el cálculo del ejemplo se usará hasta 4 términos de la serie.

f(t) = \begin{cases} -1 && -T/2<t<-T/4 \\ 1 && -T/4<t<T/4 \\ -1 && T/4<t<T/2 \end{cases}


Coeficientes ak y bk

Para determinar las expresiones de los coeficientes, en Python se usa la libreria Sympy que nos facilita el desarrollo las fórmulas simbólicas ak y bk.

Si requiere revisar el concepto se adjunta el enlace:

Fórmulas simbólicas – Sympy

El desarrollo a papel y lápiz se deja como tarea.

Usando Python se obtiene los siguientes resultados:

 expresión ak:
/1.27323*sin(1.5707*k) - 0.6366*sin(3.1415*k)  for And(k > -oo, k < oo, k != 0)
|--------------------------------------------
<                      k                                
|                                                                  
\          0                                    otherwise 

 expresión bk:
0

Nota: se han truncado los decimales a 4 para facilitar la lectura en la pantalla.

usando formato latex para la expresión, generado por Sympy se obtiene:

a_k = \begin{cases} \frac{1.2732\sin (1.5707 k) - 0.6366 \sin(3.1415 k )}{k} & \text{for}\: k > -\infty \wedge k < \infty \wedge k \neq 0 \\0 & \text{otherwise} \end{cases} b_k = 0

Instrucciones en Python:

# Serie de Fourier, con n coeficientes
# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
T = 2*np.pi

t = sym.Symbol('t')
ft = sym.Piecewise((-1,t <-T/2),
                   (-1,t <-T/4),
                   ( 1,t < T/4),
                   (-1,t < T/2),
                   (-1,True),)
# intervalo
a = -T/2
b = T/2

# número de coeficientes
n = 4

# PROCEDIMIENTO
k = sym.Symbol('k')
w0 = 2*np.pi/T

# Términos ak para coseno
enintegral = ft*sym.cos(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
ak = (2/T)*yaintegrado
ak = sym.simplify(ak)

# Términos bk para seno
enintegral = ft*sym.sin(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
bk = (2/T)*yaintegrado
bk = sym.simplify(bk)

print(' expresión ak:')
sym.pprint(ak)
print('\n ak formato latex')
print(sym.latex(ak))

print('\n expresión bk:')
sym.pprint(bk)
print('\n bk formato latex')
print(sym.latex(bk))

Evaluación de coeficientes

Con las expresiones obtenidas en el bloque anterior, se evaluan los n coeficientes ak y bk.

 coeficientes ak: 
[0, 1.27323954473516, 1.55926873300775e-16, 
  -0.424413181578388]

 coeficientes bk: 
[0, 0, 0, 0]

Instrucciones Python

las instrucciones son adicionales al bloque anterior. La evaluación se mantiene usando las expresiones simbólicas usando la instruccion .subs()

# evalua los coeficientes
a0 = ak.subs(k,0)
b0 = bk.subs(k,0)
aki = [a0]
bki = [b0]
i = 1
while not(i>=n):
    avalor = ak.subs(k,i)
    bvalor = bk.subs(k,i)
    aki.append(avalor)
    bki.append(bvalor)
    i = i+1
print('\n coeficientes ak: ')
print(aki)
print('\n coeficientes bk: ')
print(bki)

Expresión de la Serie de Fourier

Encontrados los coeficientes ak y bk, se los usa en la expresión principal

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

obteniendo la siguiente expresión para la serie de Fourier

 serie de Fourier fs(t): 
1.27323954473516*cos(1.0*t) 
+ 1.55926873300775e-16*cos(2.0*t) 
- 0.424413181578388*cos(3.0*t)

Instrucciones en Python

# serie de Fourier
serieF = a0 + 0*t 
i = 1
while not(i>=n):
    serieF= serie + aki[i]*sym.cos(i*w0*t)
    serie = serie + bki[i]*sym.sin(i*w0*t)
    i = i+1
# serie = sym.simplify(serie)

print('\n serie de Fourier fs(t): ')
sym.pprint(serie)

Gráficas de f(t) y Serie de Fourier

Para comparar la función f(t) con la aproximación en series de Fourier, se usa el intervalo [a,b] con una cantidad de muestras usando la instruccion np.linspace() y guardando el resultado en ti.

Para evaluar los puntos ti en cada expresión, por simplicidad se convierte la expresión de su forma simbólica a numérica lambda, usando la instrucción sym.lambdify()

Intrucciones en Python

# Para evaluación numérica
fsn = sym.lambdify(t,serieF)
ftn = sym.lambdify(t,ft)

# Evaluación para gráfica
muestras = 41
ti = np.linspace(a,b,muestras)
fi = ftn(ti)
fsi = fsn(ti) 

# SALIDA
# Grafica
plt.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
plt.plot(ti,fsi,label = etiqueta)
plt.xlabel('ti')
plt.legend()
plt.title('Serie de Fourier')
plt.show()

3.4.2 Respuesta a estado cero – Diagrama Bloques

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

Referencia: Lathi ejemplo 2.6 pdf129

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 de Bloques desarrollado en la sección 3.2.2 Respuesta a entrada cero, se añaden los bloques para formar la señal de entrada. La variable tiempo se indica con el Bloque Reloj, el escalón μ(t) se configura para empezar en cero y los valores de -3 del exponente se los indica con un bloque de ganancia antes de eu.

Para diferenciación de los bloques, 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 tiene su referencia de tiempo con otro reloj. El osciloscopio de la señal de entrada x(t) se muestra primero, y la salida luego.

la señal de respuesta se obtiene de forma semejante a la respuesta al 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.

3.3.2 Respuesta a impulso – 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):

 

 

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:

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

solución general:

y(t) = C_1 e^{-t} + C_2 e^{-2t}

La respuesta del algoritmo en Python es:

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.

3.4.1 Respuesta a estado cero – Python

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

Referencia:  Lathi ejemplo 2.6 pdf129

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)

Se desarrollan el ejercicio de varias formas para obtener la gráfica de salida usando:

  1. Un sistema 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

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

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

La libreria Scipy.signal permite definir un sistema LTI, usando los coeficientes de la ecuación y condiciones de inicio en cero. La salida y(t) se calcula con una versión muestreada xi de la función de entrada x(t) a partir del tiempo muestreado en el intervalo [a,b].

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

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

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

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

# Sistema
# coeficients 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()

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

Para usar la función convolve se requiere muestras en tiempos, simétricos alrededor de cero debido a la naturaleza del algoritmo. Sección que puede ser eliminada antes de mostrar la gráfica.  Para recordarlo se ha dejado la parte negativa del rango de tiempo.

# Lathi ejemplo 2.6 pdf129
# Señales contínuas en convolución
# Compare con el algoritmo en discretas
import numpy as np

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

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)


# 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
import matplotlib.pyplot as plt

plt.figure(1)
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()

3. La integral de convolución con algoritmo

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.

Parece ser un poco más lento debido a que se realizan las operaciones de forma básica (interpretada) de Python, la función usa elementos compilados más eficientes.

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

# Lathi ejemplo 2.6 pdf129
# Señales contínuas en convolución
# Compare con el algoritmo en discretas
import numpy as np

# INGRESO
# Rango [a,b], simétrico a 0
b = 5 ; a = -b; dt =0.01

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)


# 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
import matplotlib.pyplot as plt

plt.figure(1)
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()

3.2.1 Respuesta entrada cero – Python

1. Fórmulas simbólicas con Sympy

Para éste desarrollo se toma la función original , y se reordena de la forma:

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

Se define la variable t como símbolo, a y como función, el programa presenta la fórmula a la que se pueden aplicar las condiciones iniciales y una solución particular.

La ecuación diferencial:

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

solución general:

y(t) = C_1 e^{-t} + C_2 e^{-2t}
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
-5 = -C1 - 2*C2
{C1: -5, C2: 5}
Solución particular:
            -t      -2*t
y(t) = - 5*e   + 5*e    
>>>

Las intrucciones en Python son:

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

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

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

Para obtener la solución particular con las condiciones iniciales, se generan dos ecuaciones sutituyendo del lado izquierdo (lhs) los valores iniciales, y en el dado derecho (rhs) los valores de t=0 en la ecuación general y en la derivada del lado derecho.

Adicionalmente se puede realizar la gráfica, convirtiendo la ecuación de forma simbólica a forma numérica numpy con sym.lambdify(t,particular,’numpy’). Se establece el rango de observacion para evaluar la función y se envia a graficar. Se añaden las siguientes instrucciones:

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

ti = np.linspace(a,b,muestras)
yi = yt(ti)

import matplotlib.pyplot as plt
plt.plot(ti,yi, color='blue', label='y(t)')
plt.title('Respuesta entrada cero: y(t)='+str(particular.rhs))
plt.xlabel('t')
plt.grid()
plt.show()

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

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

La función Scipy.signal.lti() permite definir un sistema LTI, usando los coeficientes de la ecuación como numerador y denominador P(D)/Q(D).

La función scipy.signal.lsim(), usando:
– un intervalo [a,b] de tiempo ti muestreado,
– una señal de entrada x(t) como vector xi  y
– las condiciones de inicio  como vector [y'(0),y(0)]
determina la señal de salida en dos formas:
– la respuesta del sistema yi
– los componentes en forma de columnas yc=[yi(t),y0(t)]

usando los datos del ejercicio se grafica la respuesta a entrada cero:

# 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: t*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 = [-5,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 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()

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

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)

derivadas sutituyendo:

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

con los valores iniciales de t0=0, y0=0, z0=-5 evaluado en el rango de tiempo [0,tn=5] y como tamaño de paso h=dt=0.1

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

# 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 = -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('t, y, z')
print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(ti,yi, label='y(t)')

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

Más detalles del algoritmo en :

8.2.2 Runge-Kutta d2y/dx2

3.2.2 Respuesta entrada cero – Diagrama Bloques

Desarrollo con diagrama de bloques y un simulador

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). Con problema planteado ésta sección se cambia la forma de la ecuación:

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)

se despeja la respuesta y(t):

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

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.


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