Examen_2018_2_03_solucion.ipynb

Kernel: Python 3

%pylab inline
from scipy.integrate import odeint
from IPython.display import Math

Populating the interactive namespace from numpy and matplotlib

Elastic pendulum

From Wikipedia

The spring has the rest length $l_0$ and can be stretched to length $x$. The angle of oscillation of the pendulum is $\theta$.

The equations of motion are $$\ddot x =(l_0+x)\dot \theta^2 -\frac{k}{m}x + g\cos \theta\,,$$ and $$\ddot \theta=-\frac{g}{l_0+x}\sin \theta-\frac{2\dot x}{l_0+x}\dot \theta$$

We can define $$\omega=\dot\theta=\frac{d\theta}{dt}\,,\qquad v=\dot x=\frac{dx}{dt}\,,$$

Haciendo $$U=\begin{bmatrix} U_0\\ U_1\\ U_2\\ U_3 \end{bmatrix}=\begin{bmatrix} v\\ \omega\\ x\\ \theta\\ \end{bmatrix}.$$

Podemos obtener $$\begin{align} \frac{\operatorname{d}}{\operatorname{d} t} \begin{bmatrix} x\\ \theta\\ v\\ \omega\\ \end{bmatrix}=& \begin{bmatrix} v\\ \omega\\ (l_0+x)\omega^{2} - \frac{k}{m} x+ g\cos \theta\\ -\frac{1}{l_0+x}\begin{bmatrix} g\sin \theta + 2v\omega \end{bmatrix}&\\ \end{bmatrix}\\ \end{align}$$

1. Coordenadas x, y del movimiento

# definición de la función que da las derivadas de los parámetros por hallar
def du_dt(u0,t,L0=1,k=3.5,m=0.5,g=9.8):
    x,θ,v,ω = u0
    dudt = [v, ω, (L0+x)*ω**2 - k/m*x + g*np.cos(θ), -(g*np.sin(θ) + 2*v*ω)/(L0 + x)]
    return dudt

# definición de constantes
n = 1000   
L0 = 1            # longitud en reposo [m]
k = 4           # constante elástica [N/m]
m = 0.5           # masa [kg]
g = 9.8           # gravedad [m/s²]
u0 = [L0,0.3,0,0]  # condiciones iniciales de L,θ,v,ω

t = np.linspace(0,30,n)
#t = np.linspace(0,40,n)

# solución de las ecuaciones diferenciales: L,θ,v,ω
u = odeint(du_dt,u0,t,args=(L0,k,m,g))

# coordenadas x, y del movimiento
x = L0*np.sin(u[:,1])+u[:,0]*np.sin(u[:,1])
y = -L0*np.cos(u[:,1])-u[:,0]*np.cos(u[:,1])

# gráfico de las coordenadas
plt.plot(x,y,c='mediumslateblue',lw=0.8)
plt.title('Coordenadas del péndulo-resorte')
plt.xlabel('$x(t)$',size=13)
plt.ylabel('$y(t)$',size=13)
#plt.grid()
plt.show()

2. Espacio de configuraciones

# definición de condiciones iniciales aleatorias
n = 100
Lp = np.random.uniform(0,1,n)
θp = np.random.uniform(-np.pi/2,np.pi/2,n)
vp = np.random.uniform(-1,1,n)
ωp = np.random.uniform(-1,1,n)

t = np.linspace(0,30,10000)

# cálculo de L,θ,v,ω en las distintas condiciones iniciales
for L,θ,v,ω in zip(Lp,θp,vp,ωp):
    u0 = [L,θ,v,ω]
    u = odeint(du_dt,u0,t,args=(L0,k,m,g))
    plt.plot(u[:,1],u[:,3],c='b',lw=0.1)
plt.title('Espacio de configuraciones del péndulo-resorte')
plt.xlabel('$θ(t)$',size=13)
plt.ylabel('$ω(t)$',size=13)
#plt.xlim(-0.5,0.5)
#plt.ylim(-0.5,0.5)
#plt.grid()
plt.show()