s3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

Ejercicio: 3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

la ecuación a resolver es:

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

con valores iniciales: y0(t)=0 , y’0(t) =-5

se puede escribir como:

y"+3 y'+2y = 0 y" = -3y'-2y

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

y’ = z = f(x)

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=6] y definir el tamaño de paso o resolución del resultado

h = \delta t = \frac{b-a}{n} = \frac{6-0}{60} = 0.1

t0 = 0, y0 = 0,  y’0 = z0 = -5

iteración 1

K1y = h * f(ti,yi,zi) = 0.1 (-5) = -0.5

K1z = h * g(ti,yi,zi) ) = 0.1*(-3(-5)-2(0)) = 1.5

K2y = h * f(ti+h, yi + K1y, zi + K1z) = 0.1(-5+1.5) = -0.35

K2z = h * g(ti+h, yi + K1y, zi + K1z) = 0.1 ( -3(-5+1.5) – 2(0-0.5)) = 1.15

yi = yi + (K1y+K2y)/2 =0+(-0.5-0.35)/2=-0.425

zi = zi + (K1z+K2z)/2 = -5+(1.5+1.15)/2 = -3.675

ti = ti + h = 0+0.1 = 0.1

iteración 2

K1y = 0.1 (-3.675) = -0.3675

K1z = 0.1*(-3(-3.675)-2(-0.425)) = 1.1875

K2y = 0.1(-3.675+ 1.1875) = -0.24875

K2z = 0.1 ( -3(-3.675+ 1.1875) – 2(-0.425-0.3675)) = 0.90475

yi = -0.425+(-0.3675-0.24875)/2=-0.7331

zi = -3.675+( 1.1875+0.90475)/2 = -2.6288

ti =0.1+0.1 = 0.2

iteración 3

continuar como ejercicio

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

los valores de las iteraciones son:

t, y, z
[[ 0.        0.       -5.      ]
 [ 0.1      -0.425    -3.675   ]
 [ 0.2      -0.733125 -2.628875]
 [ 0.3      -0.949248 -1.807592]
 [ 0.4      -1.093401 -1.167208]
 [ 0.5      -1.18168  -0.67202 ]
 [ 0.6      -1.226984 -0.293049]
 [ 0.7      -1.239624 -0.006804]
 [ 0.8      -1.227806  0.205735]
 [ 0.9      -1.19804   0.359943]
 [ 1.       -1.155465  0.468225]
 [ 1.1      -1.104111  0.540574]
 [ 1.2      -1.047121  0.585021]
 [ 1.3      -0.986923  0.608001]
 [ 1.4      -0.925374  0.614658]
 [ 1.5      -0.863874  0.609087]
 [ 1.6      -0.803463  0.594537]
 [ 1.7      -0.744893  0.573574]
 [ 1.8      -0.68869   0.548208]
 [ 1.9      -0.635205  0.520011]
 [ 2.       -0.584652  0.490193]
 [ 2.1      -0.53714   0.459683]
 [ 2.2      -0.492695  0.42918 ]
 [ 2.3      -0.451288  0.399206]
 [ 2.4      -0.412843  0.370135]
 [ 2.5      -0.377253  0.342233]
 [ 2.6      -0.34439   0.315674]
 [ 2.7      -0.314114  0.290567]
 [ 2.8      -0.286275  0.266966]
 [ 2.9      -0.26072   0.244887]
 [ 3.       -0.237297  0.224314]
 [ 3.1      -0.215858  0.205211]
 [ 3.2      -0.196256  0.187526]
 [ 3.3      -0.178354  0.171195]
 [ 3.4      -0.162019  0.156149]
 [ 3.5      -0.147126  0.142312]
 [ 3.6      -0.133558  0.129611]
 [ 3.7      -0.121206  0.117969]
 [ 3.8      -0.109966  0.107312]
 [ 3.9      -0.099745  0.097569]
 [ 4.       -0.090454  0.08867 ]
 [ 4.1      -0.082013  0.080549]
 [ 4.2      -0.074346  0.073146]
 [ 4.3      -0.067385  0.066401]
 [ 4.4      -0.061067  0.06026 ]
 [ 4.5      -0.055334  0.054673]
 [ 4.6      -0.050134  0.049591]
 [ 4.7      -0.045417  0.044972]
 [ 4.8      -0.04114   0.040776]
 [ 4.9      -0.037263  0.036964]
 [ 5.       -0.033748  0.033503]
 [ 5.1      -0.030563  0.030362]
 [ 5.2      -0.027677  0.027512]
 [ 5.3      -0.025062  0.024926]
 [ 5.4      -0.022692  0.022581]
 [ 5.5      -0.020546  0.020455]
 [ 5.6      -0.018602  0.018527]
 [ 5.7      -0.016841  0.01678 ]
 [ 5.8      -0.015246  0.015196]
 [ 5.9      -0.013802  0.013761]
 [ 6.       -0.012494  0.012461]]

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 = 6
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.ylabel('y(t)')
plt.xlabel('t')
plt.title('Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()