Ejercicio: 3Eva2021PAOI_T3 EDO 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) =0con valores iniciales: y0(t)=0 , y’0(t) =-5
se puede escribir como:
y"+3 y'+2y = 0 y" = -3y'-2ysustituyendo 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.1t0 = 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()