Ejercicio: 2Eva_2023PAOI_T2 Péndulo vertical amortiguado
literal a
\frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)Se simplifica su forma a:
\frac{d\theta}{dt}= z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' = -\mu z -\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)se usan los valores dados: g = 9.81 m/s2, L = 2 m
f_t(t,\theta,z) = z g_t(t,\theta,z) = - 0.5 z -\frac{9.81}{2}\sin (\theta)y los valores iniciales para la tabla: θ(0) = π/4 rad, θ’ (0) = 0 rad/s, se complementan los valores en la medida que se aplica el desarrollo.
| ti | θ(ti) | θ'(ti)=z |
| 0 | π/4 | 0 |
| 0.2 | 0.7161 | -0.6583 |
| 0.4 | 0.5267 | -0.1156 |
| 0.6 | 0.2579 | -0.1410 |
| … |
literal b
Iteración 1: ti = 0 ; yi = π/4 ; zi = 0
K1y = h * ft(ti,yi,zi)
= 0.2*(0) = 0
K1z = h * gt(ti,yi,zi)
= 0.2*(-0.5(0) -(9.81/2)sin (π/4) = -0.6930
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
= 0.2*(0-0.6930)= -0.1386
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
= 0.2*(-0.5(0-0.6930) -(9.81/2)sin(π/4-0)
= -0.6237
yi = yi + (K1y+K2y)/2
= π/4+ (0+(-0.1386))/2 = 0.7161
zi = zi + (K1z+K2z)/2
= 0+(-0.6930-0.6237)/2 = -0.6583
ti = ti + h = 0 + 0.2 = 0.2
estimado[i] = [0.2,0.7161,-0.6583]
Iteración 2: ti = 0.2 ; yi = 0.7161 ; zi = -0.6583
K1y = h * ft(ti,yi,zi)
= 0.2*(-0.6583) = -0.1317
K1z = h * gt(ti,yi,zi)
= 0.2*(- 0.5 ( -0.6583) -(9.81/2)sin (0.7161)
= -0.5775
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
= 0.2*(-0.6583 -0.5775)= -0.2472
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
= 0.2*(- 0.5 (-0.6583 -0.5775) -(9.81/2)sin(0.7161-0.1317)
= -0.4171
yi = yi + (K1y+K2y)/2
= 0.7161 + (-0.1317-0.2472)/2 = 0.5267
zi = zi + (K1z+K2z)/2
= -0.6583+(-0.5775-0.4171)/2 = -0.1156
ti = ti + h = 0.2 + 0.2 = 0.4
estimado[i] = [0.4,0.5267,-0.1156]
Iteración 3: ti = 0.4 ; yi = 0.5267 ; zi = -1.156
K1y = h * ft(ti,yi,zi)
= 0.2*(-1.156) = -0.2311
K1z = h * gt(ti,yi,zi)
= 0.2*(- 0.5(-1.156) -(9.81/2)sin (0.5267)
= -0.3771
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
= 0.2*(-1.156 -0.3771)= -0.3065
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
= 0.2*(- 0.5 ( -1.156 -0.3771) -(9.81/2)sin(0.5267-0.2311)
= -0.1322
yi = yi + (K1y+K2y)/2
= 0.5267 + (-0.2311-0.3065)/2 = 0.2579
zi = zi + (K1z+K2z)/2
= -1.156+(-0.3771-0.1322)/2 = -1.410
ti = ti + h = 0.4 + 0.2 = 0.6
estimado[i] = [0.6,0.2579,-1.410]
literal c
resultados del algoritmo:
[ t, y, dyi/dti=z, K1y, K1z, K2y, K2z] [[ 0.000e+00 7.854e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00] [ 2.000e-01 7.161e-01 -6.583e-01 0.000e+00 -6.930e-01 -1.386e-01 -6.237e-01] [ 4.000e-01 5.267e-01 -1.156e+00 -1.317e-01 -5.775e-01 -2.472e-01 -4.171e-01] [ 6.000e-01 2.579e-01 -1.410e+00 -2.311e-01 -3.771e-01 -3.065e-01 -1.322e-01] [ 8.000e-01 -3.508e-02 -1.377e+00 -2.820e-01 -1.089e-01 -3.038e-01 1.756e-01] ...
con h=0.2 se tienen 1/0.2 = 5 tramos por segundo, por lo que para 10 segundo serán 50 tramos. La cantidad de muestras = tramos + 1(valor inicial) = 51
con lo que se puede usar el algoritmo en EDO Runge-Kutta d2y/dx2
literal d
Se observa que la respuesta es oscilante y amortiguada en magnitud como se esperaba según el planteamiento. Con el tiempo se estabilizará en cero.
Instrucciones en Python
# 2Eva_2023PAOI_T2 Péndulo vertical amortiguado 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,7),dtype=float) # incluye el punto [x0,y0,z0] estimado[0] = [x0,y0,z0,0,0,0,0] 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,K1y,K1z,K2y,K2z] return(estimado) # INGRESO theta = y g = 9.8 L = 2 ft = lambda t,y,z: z gt = lambda t,y,z: -0.5*z +(-g/L)*np.sin(y) t0 = 0 y0 = np.pi/4 z0 = 0 h = 0.2 muestras = 51 # PROCEDIMIENTO tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras) # SALIDA np.set_printoptions(precision=3) print(' [ t, \t\t y, \t dyi/dti=z, K1y,\t K1z,\t K2y,\t K2z]') print(tabla) # Grafica ti = np.copy(tabla[:,0]) yi = np.copy(tabla[:,1]) zi = np.copy(tabla[:,2]) plt.subplot(121) plt.plot(ti,yi) plt.grid() plt.xlabel('ti') plt.title('yi') plt.subplot(122) plt.plot(ti,zi, color='green') plt.xlabel('ti') plt.title('dyi/dti') plt.grid() plt.show()
