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

Publicado por

Edison Del Rosario

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